Spatial_data_analysis_hw2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Анализ пространственных данных. Домашнее задание №2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Мягкий дедлайн: __4 ноября 2020 г. 23:59__\n",
    "\n",
    "Жесткий дедлайн (со штрафом в _50%_ от количества набранных вами за ДЗ баллов): __5 ноября 2020 г. 08:59__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Визуализация \"чего-либо\" __без__ выполненного основного задания оценивается в __0 баллов__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ФИО: `Шевчук Анастасия Олеговна`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Группа: `MADE-DS-12`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Задание №1. Горячая точка (алгоритм - 10 баллов, визуализация - 10 баллов)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Генерируйте рандомные точки на планете Земля до тех пор, пока не попадете на территорию ``Афганистана``"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Вы можете использовать функции принадлжености точки полигону и расстояния от точки до полигона (в метрах)\n",
    "2. Предложите не наивный алгоритм поиска (генерировать __напрямую__ точку из полигона границ Афганистана __запрещено__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Визуализируйте пошагово предложенный алгоритм при помощи ``Folium``"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import folium\n",
    "import requests\n",
    "import json \n",
    "\n",
    "from random import uniform\n",
    "import geopandas as gpd\n",
    "from datetime import datetime, timedelta\n",
    "from shapely.geometry import Point, Polygon\n",
    "from haversine import haversine, Unit\n",
    "from folium import plugins"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get data\n",
    "url_prefix = 'https://raw.githubusercontent.com/python-visualization/folium/master/examples/data'\n",
    "source_url = f'{url_prefix}/world-countries.json'\n",
    "\n",
    "response = requests.get(source_url)\n",
    "all_countries_data = json.loads(response.content.decode())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract Afghanistan data\n",
    "for item in all_countries_data['features']:\n",
    "    if item['properties']['name'] == 'Afghanistan':\n",
    "        afghanistan_data = item\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for generating a random point from the area bounded by the given lat and lon\n",
    "def random_point(lat_1, lat_2, lon_1, lon_2):\n",
    "    return Point(uniform(lat_1, lat_2), uniform(lon_1, lon_2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_dist_to_afghan_center_in_meters(x, y):\n",
    "    return haversine((x, y), afghanistan_center_coordinates, unit = Unit.METERS)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_distances(point):\n",
    "    result = []\n",
    "    result.append(calculate_dist_to_afghan_center_in_meters(point.x + abs(lat_2 - lat_1) // 4, point.y))\n",
    "    result.append(calculate_dist_to_afghan_center_in_meters(point.x,  point.y - abs(lon_1 - lon_2) // 4))\n",
    "    result.append(calculate_dist_to_afghan_center_in_meters(point.x, point.y + abs(lon_1 - lon_2) // 4))\n",
    "    result.append(calculate_dist_to_afghan_center_in_meters(point.x - abs(lat_2 - lat_1) // 4, point.y))\n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "afghanistan_polygon = Polygon(afghanistan_data['geometry']['coordinates'][0])\n",
    "afghanistan_center_coordinates = [afghanistan_polygon.centroid.y, afghanistan_polygon.centroid.x]\n",
    "lat_1, lat_2, lon_1, lon_2 = -90, 90, -180, 180\n",
    "point = random_point(lat_1, lat_2, lon_1, lon_2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = []\n",
    "while (not Point(point.y,point.x).within(afghanistan_polygon)):\n",
    "    points.append(point)\n",
    "    north, west, east, south = calculate_distances(point)\n",
    "    if (west < east):\n",
    "        lon_2 = point.y\n",
    "    else:\n",
    "        lon_1 = point.y\n",
    "    if (north < south):\n",
    "        lat_1 = point.x\n",
    "    else:\n",
    "        lat_2 = point.x\n",
    "    point = random_point(lat_1, lat_2, lon_1, lon_2)\n",
    "points.append(point)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_afghan_point = point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create map\n",
    "map_ = folium.Map(location = afghanistan_center_coordinates,\n",
    "               zoom_start = 2)\n",
    "folium.GeoJson(afghanistan_data).add_to(map_)\n",
    "folium.LayerControl().add_to(map_);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Preparing for visualization with Timestamped Geo Json\n",
    "date = datetime.now()\n",
    "dates = [(date + timedelta(minutes = delta * 5)).strftime(\"%Y-%m-%dT%H:%M:%S\") for delta in range(len(points))]\n",
    "features = [\n",
    "    {\n",
    "        'type': 'Feature',\n",
    "        'geometry': {\n",
    "            'type': 'LineString',\n",
    "            'coordinates':  [[point.y, point.x]],\n",
    "        },\n",
    "        'properties': {\n",
    "            'times': [date]\n",
    "        }\n",
    "    }\n",
    "    for point, date in zip (points, dates)\n",
    "]\n",
    "\n",
    "plugins.TimestampedGeoJson({\n",
    "    'type': 'FeatureCollection',\n",
    "    'features': features,\n",
    "}, period='PT1M', add_last_point=True).add_to(map_);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style=\"width:100%;\"><div style=\"position:relative;width:100%;height:0;padding-bottom:60%;\"><span style=\"color:#565656\">Make this Notebook Trusted to load map: File -> Trust Notebook</span><iframe src=\"about:blank\" style=\"position:absolute;width:100%;height:100%;left:0;top:0;border:none !important;\" data-html=PCFET0NUWVBFIGh0bWw+CjxoZWFkPiAgICAKICAgIDxtZXRhIGh0dHAtZXF1aXY9ImNvbnRlbnQtdHlwZSIgY29udGVudD0idGV4dC9odG1sOyBjaGFyc2V0PVVURi04IiAvPgogICAgCiAgICAgICAgPHNjcmlwdD4KICAgICAgICAgICAgTF9OT19UT1VDSCA9IGZhbHNlOwogICAgICAgICAgICBMX0RJU0FCTEVfM0QgPSBmYWxzZTsKICAgICAgICA8L3NjcmlwdD4KICAgIAogICAgPHNjcmlwdCBzcmM9Imh0dHBzOi8vY2RuLmpzZGVsaXZyLm5ldC9ucG0vbGVhZmxldEAxLjYuMC9kaXN0L2xlYWZsZXQuanMiPjwvc2NyaXB0PgogICAgPHNjcmlwdCBzcmM9Imh0dHBzOi8vY29kZS5qcXVlcnkuY29tL2pxdWVyeS0xLjEyLjQubWluLmpzIj48L3NjcmlwdD4KICAgIDxzY3JpcHQgc3JjPSJodHRwczovL21heGNkbi5ib290c3RyYXBjZG4uY29tL2Jvb3RzdHJhcC8zLjIuMC9qcy9ib290c3RyYXAubWluLmpzIj48L3NjcmlwdD4KICAgIDxzY3JpcHQgc3JjPSJodHRwczovL2NkbmpzLmNsb3VkZmxhcmUuY29tL2FqYXgvbGlicy9MZWFmbGV0LmF3ZXNvbWUtbWFya2Vycy8yLjAuMi9sZWFmbGV0LmF3ZXNvbWUtbWFya2Vycy5qcyI+PC9zY3JpcHQ+CiAgICA8bGluayByZWw9InN0eWxlc2hlZXQiIGhyZWY9Imh0dHBzOi8vY2RuLmpzZGVsaXZyLm5ldC9ucG0vbGVhZmxldEAxLjYuMC9kaXN0L2xlYWZsZXQuY3NzIi8+CiAgICA8bGluayByZWw9InN0eWxlc2hlZXQiIGhyZWY9Imh0dHBzOi8vbWF4Y2RuLmJvb3RzdHJhcGNkbi5jb20vYm9vdHN0cmFwLzMuMi4wL2Nzcy9ib290c3RyYXAubWluLmNzcyIvPgogICAgPGxpbmsgcmVsPSJzdHlsZXNoZWV0IiBocmVmPSJodHRwczovL21heGNkbi5ib290c3RyYXBjZG4uY29tL2Jvb3RzdHJhcC8zLjIuMC9jc3MvYm9vdHN0cmFwLXRoZW1lLm1pbi5jc3MiLz4KICAgIDxsaW5rIHJlbD0ic3R5bGVzaGVldCIgaHJlZj0iaHR0cHM6Ly9tYXhjZG4uYm9vdHN0cmFwY2RuLmNvbS9mb250LWF3ZXNvbWUvNC42LjMvY3NzL2ZvbnQtYXdlc29tZS5taW4uY3NzIi8+CiAgICA8bGluayByZWw9InN0eWxlc2hlZXQiIGhyZWY9Imh0dHBzOi8vY2RuanMuY2xvdWRmbGFyZS5jb20vYWpheC9saWJzL0xlYWZsZXQuYXdlc29tZS1tYXJrZXJzLzIuMC4yL2xlYWZsZXQuYXdlc29tZS1tYXJrZXJzLmNzcyIvPgogICAgPGxpbmsgcmVsPSJzdHlsZXNoZWV0IiBocmVmPSJodHRwczovL3Jhd2Nkbi5naXRoYWNrLmNvbS9weXRob24tdmlzdWFsaXphdGlvbi9mb2xpdW0vbWFzdGVyL2ZvbGl1bS90ZW1wbGF0ZXMvbGVhZmxldC5hd2Vzb21lLnJvdGF0ZS5jc3MiLz4KICAgIDxzdHlsZT5odG1sLCBib2R5IHt3aWR0aDogMTAwJTtoZWlnaHQ6IDEwMCU7bWFyZ2luOiAwO3BhZGRpbmc6IDA7fTwvc3R5bGU+CiAgICA8c3R5bGU+I21hcCB7cG9zaXRpb246YWJzb2x1dGU7dG9wOjA7Ym90dG9tOjA7cmlnaHQ6MDtsZWZ0OjA7fTwvc3R5bGU+CiAgICAKICAgICAgICAgICAgPG1ldGEgbmFtZT0idmlld3BvcnQiIGNvbnRlbnQ9IndpZHRoPWRldmljZS13aWR0aCwKICAgICAgICAgICAgICAgIGluaXRpYWwtc2NhbGU9MS4wLCBtYXhpbXVtLXNjYWxlPTEuMCwgdXNlci1zY2FsYWJsZT1ubyIgLz4KICAgICAgICAgICAgPHN0eWxlPgogICAgICAgICAgICAgICAgI21hcF8xMzM3NWJiNDUxNjU0MDY3OTUzYjMyMjZmZjdiMzM5NyB7CiAgICAgICAgICAgICAgICAgICAgcG9zaXRpb246IHJlbGF0aXZlOwogICAgICAgICAgICAgICAgICAgIHdpZHRoOiAxMDAuMCU7CiAgICAgICAgICAgICAgICAgICAgaGVpZ2h0OiAxMDAuMCU7CiAgICAgICAgICAgICAgICAgICAgbGVmdDogMC4wJTsKICAgICAgICAgICAgICAgICAgICB0b3A6IDAuMCU7CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIDwvc3R5bGU+CiAgICAgICAgCiAgICA8c2NyaXB0IHNyYz0iaHR0cHM6Ly9jZG5qcy5jbG91ZGZsYXJlLmNvbS9hamF4L2xpYnMvanF1ZXJ5LzIuMC4wL2pxdWVyeS5taW4uanMiPjwvc2NyaXB0PgogICAgPHNjcmlwdCBzcmM9Imh0dHBzOi8vY2RuanMuY2xvdWRmbGFyZS5jb20vYWpheC9saWJzL2pxdWVyeXVpLzEuMTAuMi9qcXVlcnktdWkubWluLmpzIj48L3NjcmlwdD4KICAgIDxzY3JpcHQgc3JjPSJodHRwczovL3Jhd2Nkbi5naXRoYWNrLmNvbS9uZXphc2EvaXNvODYwMS1qcy1wZXJpb2QvbWFzdGVyL2lzbzg2MDEubWluLmpzIj48L3NjcmlwdD4KICAgIDxzY3JpcHQgc3JjPSJodHRwczovL3Jhd2Nkbi5naXRoYWNrLmNvbS9zb2NpYi9MZWFmbGV0LlRpbWVEaW1lbnNpb24vbWFzdGVyL2Rpc3QvbGVhZmxldC50aW1lZGltZW5zaW9uLm1pbi5qcyI+PC9zY3JpcHQ+CiAgICA8bGluayByZWw9InN0eWxlc2hlZXQiIGhyZWY9Imh0dHBzOi8vY2RuanMuY2xvdWRmbGFyZS5jb20vYWpheC9saWJzL2hpZ2hsaWdodC5qcy84LjQvc3R5bGVzL2RlZmF1bHQubWluLmNzcyIvPgogICAgPGxpbmsgcmVsPSJzdHlsZXNoZWV0IiBocmVmPSJodHRwczovL3Jhd2Nkbi5naXRoYWNrLmNvbS9zb2NpYi9MZWFmbGV0LlRpbWVEaW1lbnNpb24vbWFzdGVyL2Rpc3QvbGVhZmxldC50aW1lZGltZW5zaW9uLmNvbnRyb2wubWluLmNzcyIvPgogICAgPHNjcmlwdCBzcmM9Imh0dHBzOi8vY2RuanMuY2xvdWRmbGFyZS5jb20vYWpheC9saWJzL21vbWVudC5qcy8yLjE4LjEvbW9tZW50Lm1pbi5qcyI+PC9zY3JpcHQ+CjwvaGVhZD4KPGJvZHk+ICAgIAogICAgCiAgICAgICAgICAgIDxkaXYgY2xhc3M9ImZvbGl1bS1tYXAiIGlkPSJtYXBfMTMzNzViYjQ1MTY1NDA2Nzk1M2IzMjI2ZmY3YjMzOTciID48L2Rpdj4KICAgICAgICAKPC9ib2R5Pgo8c2NyaXB0PiAgICAKICAgIAogICAgICAgICAgICB2YXIgbWFwXzEzMzc1YmI0NTE2NTQwNjc5NTNiMzIyNmZmN2IzMzk3ID0gTC5tYXAoCiAgICAgICAgICAgICAgICAibWFwXzEzMzc1YmI0NTE2NTQwNjc5NTNiMzIyNmZmN2IzMzk3IiwKICAgICAgICAgICAgICAgIHsKICAgICAgICAgICAgICAgICAgICBjZW50ZXI6IFszMy44NTYzOTkyNDMzMTg1MSwgNjYuMDg2NjkwMTczNDQ2MzldLAogICAgICAgICAgICAgICAgICAgIGNyczogTC5DUlMuRVBTRzM4NTcsCiAgICAgICAgICAgICAgICAgICAgem9vbTogMiwKICAgICAgICAgICAgICAgICAgICB6b29tQ29udHJvbDogdHJ1ZSwKICAgICAgICAgICAgICAgICAgICBwcmVmZXJDYW52YXM6IGZhbHNlLAogICAgICAgICAgICAgICAgfQogICAgICAgICAgICApOwoKICAgICAgICAgICAgCgogICAgICAgIAogICAgCiAgICAgICAgICAgIHZhciB0aWxlX2xheWVyXzk0NjY0ZjZlYjdmMTQ1NTk5ZGFjZDY1YWJiZDE2OTVlID0gTC50aWxlTGF5ZXIoCiAgICAgICAgICAgICAgICAiaHR0cHM6Ly97c30udGlsZS5vcGVuc3RyZWV0bWFwLm9yZy97en0ve3h9L3t5fS5wbmciLAogICAgICAgICAgICAgICAgeyJhdHRyaWJ1dGlvbiI6ICJEYXRhIGJ5IFx1MDAyNmNvcHk7IFx1MDAzY2EgaHJlZj1cImh0dHA6Ly9vcGVuc3RyZWV0bWFwLm9yZ1wiXHUwMDNlT3BlblN0cmVldE1hcFx1MDAzYy9hXHUwMDNlLCB1bmRlciBcdTAwM2NhIGhyZWY9XCJodHRwOi8vd3d3Lm9wZW5zdHJlZXRtYXAub3JnL2NvcHlyaWdodFwiXHUwMDNlT0RiTFx1MDAzYy9hXHUwMDNlLiIsICJkZXRlY3RSZXRpbmEiOiBmYWxzZSwgIm1heE5hdGl2ZVpvb20iOiAxOCwgIm1heFpvb20iOiAxOCwgIm1pblpvb20iOiAwLCAibm9XcmFwIjogZmFsc2UsICJvcGFjaXR5IjogMSwgInN1YmRvbWFpbnMiOiAiYWJjIiwgInRtcyI6IGZhbHNlfQogICAgICAgICAgICApLmFkZFRvKG1hcF8xMzM3NWJiNDUxNjU0MDY3OTUzYjMyMjZmZjdiMzM5Nyk7CiAgICAgICAgCiAgICAKICAgICAgICBmdW5jdGlvbiBnZW9fanNvbl9kZGMwMDQ3ZTIyNTk0NDM2YWU4N2JhZWVkZTAwODM5Nl9vbkVhY2hGZWF0dXJlKGZlYXR1cmUsIGxheWVyKSB7CiAgICAgICAgICAgIGxheWVyLm9uKHsKICAgICAgICAgICAgICAgIGNsaWNrOiBmdW5jdGlvbihlKSB7CiAgICAgICAgICAgICAgICAgICAgbWFwXzEzMzc1YmI0NTE2NTQwNjc5NTNiMzIyNmZmN2IzMzk3LmZpdEJvdW5kcyhlLnRhcmdldC5nZXRCb3VuZHMoKSk7CiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgIH0pOwogICAgICAgIH07CiAgICAgICAgdmFyIGdlb19qc29uX2RkYzAwNDdlMjI1OTQ0MzZhZTg3YmFlZWRlMDA4Mzk2ID0gTC5nZW9Kc29uKG51bGwsIHsKICAgICAgICAgICAgICAgIG9uRWFjaEZlYXR1cmU6IGdlb19qc29uX2RkYzAwNDdlMjI1OTQ0MzZhZTg3YmFlZWRlMDA4Mzk2X29uRWFjaEZlYXR1cmUsCiAgICAgICAgICAgIAogICAgICAgIH0pLmFkZFRvKG1hcF8xMzM3NWJiNDUxNjU0MDY3OTUzYjMyMjZmZjdiMzM5Nyk7CgogICAgICAgIGZ1bmN0aW9uIGdlb19qc29uX2RkYzAwNDdlMjI1OTQ0MzZhZTg3YmFlZWRlMDA4Mzk2X2FkZCAoZGF0YSkgewogICAgICAgICAgICBnZW9fanNvbl9kZGMwMDQ3ZTIyNTk0NDM2YWU4N2JhZWVkZTAwODM5Ni5hZGREYXRhKGRhdGEpOwogICAgICAgIH0KICAgICAgICAgICAgZ2VvX2pzb25fZGRjMDA0N2UyMjU5NDQzNmFlODdiYWVlZGUwMDgzOTZfYWRkKHsiZ2VvbWV0cnkiOiB7ImNvb3JkaW5hdGVzIjogW1tbNjEuMjEwODE3LCAzNS42NTAwNzJdLCBbNjIuMjMwNjUxLCAzNS4yNzA2NjRdLCBbNjIuOTg0NjYyLCAzNS40MDQwNDFdLCBbNjMuMTkzNTM4LCAzNS44NTcxNjZdLCBbNjMuOTgyODk2LCAzNi4wMDc5NTddLCBbNjQuNTQ2NDc5LCAzNi4zMTIwNzNdLCBbNjQuNzQ2MTA1LCAzNy4xMTE4MThdLCBbNjUuNTg4OTQ4LCAzNy4zMDUyMTddLCBbNjUuNzQ1NjMxLCAzNy42NjExNjRdLCBbNjYuMjE3Mzg1LCAzNy4zOTM3OV0sIFs2Ni41MTg2MDcsIDM3LjM2Mjc4NF0sIFs2Ny4wNzU3ODIsIDM3LjM1NjE0NF0sIFs2Ny44MywgMzcuMTQ0OTk0XSwgWzY4LjEzNTU2MiwgMzcuMDIzMTE1XSwgWzY4Ljg1OTQ0NiwgMzcuMzQ0MzM2XSwgWzY5LjE5NjI3MywgMzcuMTUxMTQ0XSwgWzY5LjUxODc4NSwgMzcuNjA4OTk3XSwgWzcwLjExNjU3OCwgMzcuNTg4MjIzXSwgWzcwLjI3MDU3NCwgMzcuNzM1MTY1XSwgWzcwLjM3NjMwNCwgMzguMTM4Mzk2XSwgWzcwLjgwNjgyMSwgMzguNDg2MjgyXSwgWzcxLjM0ODEzMSwgMzguMjU4OTA1XSwgWzcxLjIzOTQwNCwgMzcuOTUzMjY1XSwgWzcxLjU0MTkxOCwgMzcuOTA1Nzc0XSwgWzcxLjQ0ODY5MywgMzcuMDY1NjQ1XSwgWzcxLjg0NDYzOCwgMzYuNzM4MTcxXSwgWzcyLjE5MzA0MSwgMzYuOTQ4Mjg4XSwgWzcyLjYzNjg5LCAzNy4wNDc1NThdLCBbNzMuMjYwMDU2LCAzNy40OTUyNTddLCBbNzMuOTQ4Njk2LCAzNy40MjE1NjZdLCBbNzQuOTgwMDAyLCAzNy40MTk5OV0sIFs3NS4xNTgwMjgsIDM3LjEzMzAzMV0sIFs3NC41NzU4OTMsIDM3LjAyMDg0MV0sIFs3NC4wNjc1NTIsIDM2LjgzNjE3Nl0sIFs3Mi45MjAwMjUsIDM2LjcyMDAwN10sIFs3MS44NDYyOTIsIDM2LjUwOTk0Ml0sIFs3MS4yNjIzNDgsIDM2LjA3NDM4OF0sIFs3MS40OTg3NjgsIDM1LjY1MDU2M10sIFs3MS42MTMwNzYsIDM1LjE1MzIwM10sIFs3MS4xMTUwMTksIDM0LjczMzEyNl0sIFs3MS4xNTY3NzMsIDM0LjM0ODkxMV0sIFs3MC44ODE4MDMsIDMzLjk4ODg1Nl0sIFs2OS45MzA1NDMsIDM0LjAyMDEyXSwgWzcwLjMyMzU5NCwgMzMuMzU4NTMzXSwgWzY5LjY4NzE0NywgMzMuMTA1NDk5XSwgWzY5LjI2MjUyMiwgMzIuNTAxOTQ0XSwgWzY5LjMxNzc2NCwgMzEuOTAxNDEyXSwgWzY4LjkyNjY3NywgMzEuNjIwMTg5XSwgWzY4LjU1NjkzMiwgMzEuNzEzMzFdLCBbNjcuNzkyNjg5LCAzMS41ODI5M10sIFs2Ny42ODMzOTQsIDMxLjMwMzE1NF0sIFs2Ni45Mzg4OTEsIDMxLjMwNDkxMV0sIFs2Ni4zODE0NTgsIDMwLjczODg5OV0sIFs2Ni4zNDY0NzMsIDI5Ljg4Nzk0M10sIFs2NS4wNDY4NjIsIDI5LjQ3MjE4MV0sIFs2NC4zNTA0MTksIDI5LjU2MDAzMV0sIFs2NC4xNDgwMDIsIDI5LjM0MDgxOV0sIFs2My41NTAyNjEsIDI5LjQ2ODMzMV0sIFs2Mi41NDk4NTcsIDI5LjMxODU3Ml0sIFs2MC44NzQyNDgsIDI5LjgyOTIzOV0sIFs2MS43ODEyMjIsIDMwLjczNTg1XSwgWzYxLjY5OTMxNCwgMzEuMzc5NTA2XSwgWzYwLjk0MTk0NSwgMzEuNTQ4MDc1XSwgWzYwLjg2MzY1NSwgMzIuMTgyOTJdLCBbNjAuNTM2MDc4LCAzMi45ODEyNjldLCBbNjAuOTYzNywgMzMuNTI4ODMyXSwgWzYwLjUyODQzLCAzMy42NzY0NDZdLCBbNjAuODAzMTkzLCAzNC40MDQxMDJdLCBbNjEuMjEwODE3LCAzNS42NTAwNzJdXV0sICJ0eXBlIjogIlBvbHlnb24ifSwgImlkIjogIkFGRyIsICJwcm9wZXJ0aWVzIjogeyJuYW1lIjogIkFmZ2hhbmlzdGFuIn0sICJ0eXBlIjogIkZlYXR1cmUifSk7CiAgICAgICAgCiAgICAKICAgICAgICAgICAgdmFyIGxheWVyX2NvbnRyb2xfYmU1Nzg4YmQ5MDdhNGE2YWI2OTlkNDVmNjAzMDNiNDAgPSB7CiAgICAgICAgICAgICAgICBiYXNlX2xheWVycyA6IHsKICAgICAgICAgICAgICAgICAgICAib3BlbnN0cmVldG1hcCIgOiB0aWxlX2xheWVyXzk0NjY0ZjZlYjdmMTQ1NTk5ZGFjZDY1YWJiZDE2OTVlLAogICAgICAgICAgICAgICAgfSwKICAgICAgICAgICAgICAgIG92ZXJsYXlzIDogIHsKICAgICAgICAgICAgICAgICAgICAibWFjcm9fZWxlbWVudF9kZGMwMDQ3ZTIyNTk0NDM2YWU4N2JhZWVkZTAwODM5NiIgOiBnZW9fanNvbl9kZGMwMDQ3ZTIyNTk0NDM2YWU4N2JhZWVkZTAwODM5NiwKICAgICAgICAgICAgICAgIH0sCiAgICAgICAgICAgIH07CiAgICAgICAgICAgIEwuY29udHJvbC5sYXllcnMoCiAgICAgICAgICAgICAgICBsYXllcl9jb250cm9sX2JlNTc4OGJkOTA3YTRhNmFiNjk5ZDQ1ZjYwMzAzYjQwLmJhc2VfbGF5ZXJzLAogICAgICAgICAgICAgICAgbGF5ZXJfY29udHJvbF9iZTU3ODhiZDkwN2E0YTZhYjY5OWQ0NWY2MDMwM2I0MC5vdmVybGF5cywKICAgICAgICAgICAgICAgIHsiYXV0b1pJbmRleCI6IHRydWUsICJjb2xsYXBzZWQiOiB0cnVlLCAicG9zaXRpb24iOiAidG9wcmlnaHQifQogICAgICAgICAgICApLmFkZFRvKG1hcF8xMzM3NWJiNDUxNjU0MDY3OTUzYjMyMjZmZjdiMzM5Nyk7CiAgICAgICAgCiAgICAKICAgICAgICAgICAgTC5Db250cm9sLlRpbWVEaW1lbnNpb25DdXN0b20gPSBMLkNvbnRyb2wuVGltZURpbWVuc2lvbi5leHRlbmQoewogICAgICAgICAgICAgICAgX2dldERpc3BsYXlEYXRlRm9ybWF0OiBmdW5jdGlvbihkYXRlKXsKICAgICAgICAgICAgICAgICAgICB2YXIgbmV3ZGF0ZSA9IG5ldyBtb21lbnQoZGF0ZSk7CiAgICAgICAgICAgICAgICAgICAgY29uc29sZS5sb2cobmV3ZGF0ZSkKICAgICAgICAgICAgICAgICAgICByZXR1cm4gbmV3ZGF0ZS5mb3JtYXQoIllZWVktTU0tREQgSEg6bW06c3MiKTsKICAgICAgICAgICAgICAgIH0KICAgICAgICAgICAgfSk7CiAgICAgICAgICAgIG1hcF8xMzM3NWJiNDUxNjU0MDY3OTUzYjMyMjZmZjdiMzM5Ny50aW1lRGltZW5zaW9uID0gTC50aW1lRGltZW5zaW9uKAogICAgICAgICAgICAgICAgewogICAgICAgICAgICAgICAgICAgIHBlcmlvZDogIlBUMU0iLAogICAgICAgICAgICAgICAgfQogICAgICAgICAgICApOwogICAgICAgICAgICB2YXIgdGltZURpbWVuc2lvbkNvbnRyb2wgPSBuZXcgTC5Db250cm9sLlRpbWVEaW1lbnNpb25DdXN0b20oCiAgICAgICAgICAgICAgICB7ImF1dG9QbGF5IjogdHJ1ZSwgImxvb3BCdXR0b24iOiBmYWxzZSwgIm1heFNwZWVkIjogMTAsICJtaW5TcGVlZCI6IDAuMSwgInBsYXllck9wdGlvbnMiOiB7Imxvb3AiOiB0cnVlLCAic3RhcnRPdmVyIjogdHJ1ZSwgInRyYW5zaXRpb25UaW1lIjogMjAwfSwgInBvc2l0aW9uIjogImJvdHRvbWxlZnQiLCAidGltZVNsaWRlckRyYWdVcGRhdGUiOiBmYWxzZX0KICAgICAgICAgICAgKTsKICAgICAgICAgICAgbWFwXzEzMzc1YmI0NTE2NTQwNjc5NTNiMzIyNmZmN2IzMzk3LmFkZENvbnRyb2wodGhpcy50aW1lRGltZW5zaW9uQ29udHJvbCk7CgogICAgICAgICAgICB2YXIgZ2VvSnNvbkxheWVyID0gTC5nZW9Kc29uKHsidHlwZSI6ICJGZWF0dXJlQ29sbGVjdGlvbiIsICJmZWF0dXJlcyI6IFt7InR5cGUiOiAiRmVhdHVyZSIsICJnZW9tZXRyeSI6IHsidHlwZSI6ICJMaW5lU3RyaW5nIiwgImNvb3JkaW5hdGVzIjogW1s0OS43OTM2NjE3MzkxMTA4LCA3LjI3NjQ2NTkyMTQyNDEwNV1dfSwgInByb3BlcnRpZXMiOiB7InRpbWVzIjogWyIyMDIwLTExLTA0VDA1OjA5OjQwIl19fSwgeyJ0eXBlIjogIkZlYXR1cmUiLCAiZ2VvbWV0cnkiOiB7InR5cGUiOiAiTGluZVN0cmluZyIsICJjb29yZGluYXRlcyI6IFtbNjUuMzI2NDMwNTAxMDkwODQsIDY2LjQxNzAzMDkyMDU1NTFdXX0sICJwcm9wZXJ0aWVzIjogeyJ0aW1lcyI6IFsiMjAyMC0xMS0wNFQwNToxNDo0MCJdfX0sIHsidHlwZSI6ICJGZWF0dXJlIiwgImdlb21ldHJ5IjogeyJ0eXBlIjogIkxpbmVTdHJpbmciLCAiY29vcmRpbmF0ZXMiOiBbWzEyOS42NzAzODk1OTMxNDgyMywgNjUuODE1MzE0MDA4MjIxMDFdXX0sICJwcm9wZXJ0aWVzIjogeyJ0aW1lcyI6IFsiMjAyMC0xMS0wNFQwNToxOTo0MCJdfX0sIHsidHlwZSI6ICJGZWF0dXJlIiwgImdlb21ldHJ5IjogeyJ0eXBlIjogIkxpbmVTdHJpbmciLCAiY29vcmRpbmF0ZXMiOiBbWzEwNC43MTk3MTUwNDU0NDkxNiwgMzMuNjEwMjMzMjgzNDcwNjM1XV19LCAicHJvcGVydGllcyI6IHsidGltZXMiOiBbIjIwMjAtMTEtMDRUMDU6MjQ6NDAiXX19LCB7InR5cGUiOiAiRmVhdHVyZSIsICJnZW9tZXRyeSI6IHsidHlwZSI6ICJMaW5lU3RyaW5nIiwgImNvb3JkaW5hdGVzIjogW1s2OS41NzMyODkxOTQxMDkxMSwgNTcuMzc2MTY2MDA2NzQwMzldXX0sICJwcm9wZXJ0aWVzIjogeyJ0aW1lcyI6IFsiMjAyMC0xMS0wNFQwNToyOTo0MCJdfX0sIHsidHlwZSI6ICJGZWF0dXJlIiwgImdlb21ldHJ5IjogeyJ0eXBlIjogIkxpbmVTdHJpbmciLCAiY29vcmRpbmF0ZXMiOiBbWzY4LjAyMDcxMjEwODA5ODQ2LCA1NS4xMzY3NzQ4NjQyODMzXV19LCAicHJvcGVydGllcyI6IHsidGltZXMiOiBbIjIwMjAtMTEtMDRUMDU6MzQ6NDAiXX19LCB7InR5cGUiOiAiRmVhdHVyZSIsICJnZW9tZXRyeSI6IHsidHlwZSI6ICJMaW5lU3RyaW5nIiwgImNvb3JkaW5hdGVzIjogW1s2Ny4xMjc3Nzk0MTgyMzA1NywgNDMuNzgzNjQ2ODI0OTIyNDg1XV19LCAicHJvcGVydGllcyI6IHsidGltZXMiOiBbIjIwMjAtMTEtMDRUMDU6Mzk6NDAiXX19LCB7InR5cGUiOiAiRmVhdHVyZSIsICJnZW9tZXRyeSI6IHsidHlwZSI6ICJMaW5lU3RyaW5nIiwgImNvb3JkaW5hdGVzIjogW1s2Ny4zMTgxNjYzMTExNDYyNCwgMzkuMTQ0MjI5ODU5NTg0MzVdXX0sICJwcm9wZXJ0aWVzIjogeyJ0aW1lcyI6IFsiMjAyMC0xMS0wNFQwNTo0NDo0MCJdfX0sIHsidHlwZSI6ICJGZWF0dXJlIiwgImdlb21ldHJ5IjogeyJ0eXBlIjogIkxpbmVTdHJpbmciLCAiY29vcmRpbmF0ZXMiOiBbWzY3LjYwNzgyODMyMzA4MjYzLCAzNC44NDcxMTk4MjgxNTc3Nl1dfSwgInByb3BlcnRpZXMiOiB7InRpbWVzIjogWyIyMDIwLTExLTA0VDA1OjQ5OjQwIl19fV19LCB7CiAgICAgICAgICAgICAgICAgICAgcG9pbnRUb0xheWVyOiBmdW5jdGlvbiAoZmVhdHVyZSwgbGF0TG5nKSB7CiAgICAgICAgICAgICAgICAgICAgICAgIGlmIChmZWF0dXJlLnByb3BlcnRpZXMuaWNvbiA9PSAnbWFya2VyJykgewogICAgICAgICAgICAgICAgICAgICAgICAgICAgaWYoZmVhdHVyZS5wcm9wZXJ0aWVzLmljb25zdHlsZSl7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmV0dXJuIG5ldyBMLk1hcmtlcihsYXRMbmcsIHsKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgaWNvbjogTC5pY29uKGZlYXR1cmUucHJvcGVydGllcy5pY29uc3R5bGUpfSk7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAvL2Vsc2UKICAgICAgICAgICAgICAgICAgICAgICAgICAgIHJldHVybiBuZXcgTC5NYXJrZXIobGF0TG5nKTsKICAgICAgICAgICAgICAgICAgICAgICAgfQogICAgICAgICAgICAgICAgICAgICAgICBpZiAoZmVhdHVyZS5wcm9wZXJ0aWVzLmljb24gPT0gJ2NpcmNsZScpIHsKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGlmIChmZWF0dXJlLnByb3BlcnRpZXMuaWNvbnN0eWxlKSB7CiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgcmV0dXJuIG5ldyBMLmNpcmNsZU1hcmtlcihsYXRMbmcsIGZlYXR1cmUucHJvcGVydGllcy5pY29uc3R5bGUpCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgfTsKICAgICAgICAgICAgICAgICAgICAgICAgICAgIC8vZWxzZQogICAgICAgICAgICAgICAgICAgICAgICAgICAgcmV0dXJuIG5ldyBMLmNpcmNsZU1hcmtlcihsYXRMbmcpOwogICAgICAgICAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICAgICAgICAgICAgIC8vZWxzZQoKICAgICAgICAgICAgICAgICAgICAgICAgcmV0dXJuIG5ldyBMLk1hcmtlcihsYXRMbmcpOwogICAgICAgICAgICAgICAgICAgIH0sCiAgICAgICAgICAgICAgICAgICAgc3R5bGU6IGZ1bmN0aW9uIChmZWF0dXJlKSB7CiAgICAgICAgICAgICAgICAgICAgICAgIHJldHVybiBmZWF0dXJlLnByb3BlcnRpZXMuc3R5bGU7CiAgICAgICAgICAgICAgICAgICAgfSwKICAgICAgICAgICAgICAgICAgICBvbkVhY2hGZWF0dXJlOiBmdW5jdGlvbihmZWF0dXJlLCBsYXllcikgewogICAgICAgICAgICAgICAgICAgICAgICBpZiAoZmVhdHVyZS5wcm9wZXJ0aWVzLnBvcHVwKSB7CiAgICAgICAgICAgICAgICAgICAgICAgIGxheWVyLmJpbmRQb3B1cChmZWF0dXJlLnByb3BlcnRpZXMucG9wdXApOwogICAgICAgICAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICAgICAgICAgfQogICAgICAgICAgICAgICAgfSkKCiAgICAgICAgICAgIHZhciB0aW1lc3RhbXBlZF9nZW9fanNvbl8zMzE0YzQzM2VmNzg0MjBmYTAxNGNmY2NkYmY2YmExOSA9IEwudGltZURpbWVuc2lvbi5sYXllci5nZW9Kc29uKAogICAgICAgICAgICAgICAgZ2VvSnNvbkxheWVyLAogICAgICAgICAgICAgICAgewogICAgICAgICAgICAgICAgICAgIHVwZGF0ZVRpbWVEaW1lbnNpb246IHRydWUsCiAgICAgICAgICAgICAgICAgICAgYWRkbGFzdFBvaW50OiB0cnVlLAogICAgICAgICAgICAgICAgICAgIGR1cmF0aW9uOiB1bmRlZmluZWQsCiAgICAgICAgICAgICAgICB9CiAgICAgICAgICAgICkuYWRkVG8obWFwXzEzMzc1YmI0NTE2NTQwNjc5NTNiMzIyNmZmN2IzMzk3KTsKICAgICAgICAKPC9zY3JpcHQ+ onload=\"this.contentDocument.open();this.contentDocument.write(atob(this.getAttribute('data-html')));this.contentDocument.close();\" allowfullscreen webkitallowfullscreen mozallowfullscreen></iframe></div></div>"
      ],
      "text/plain": [
       "<folium.folium.Map at 0x60b29f0>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "map_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Задание №2. Качество жизни (20 баллов)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Расстояние от точки до 5 ближайших __*__ банкоматов, находящихся в стране с наибольшим количеством объектов жилой недвижимости\n",
    "2. Расстояние от точки до 5 ближайших школ, находящихся в стране с наибольшим количеством аптек в столице\n",
    "3. Расстояние от точки до 5 ближайших кинотеатров, наодящихся в стране с самым большим отношением числа железнодорожных станций к автобусным остановкам в южной части __**__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для измерения показателя качества жизни в точке, найденной в предыдущем задании, вам необходимо рассчитать следующую сумму расстояний (в метрах):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__*__ При поиске _N_ ближайших объектов обязательно использовать ``R-tree``\n",
    "\n",
    "__**__ Южной частью страны является территория, находящаяся к югу от множества точек, равноудаленных от самой северной и самой южной точек страны"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from OSMPythonTools.overpass import overpassQueryBuilder, Overpass\n",
    "from OSMPythonTools.nominatim import Nominatim\n",
    "from itertools import chain \n",
    "from rtree import index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "overpass = Overpass()\n",
    "nominatim = Nominatim()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Расстояние от точки до 5 ближайших банкоматов, находящихся в стране с наибольшим количеством объектов жилой недвижимости"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all countries\n",
    "response = overpass.query('relation[admin_level=2][type=boundary][boundary=administrative];out;', timeout = 120)\n",
    "# Filter data\n",
    "countries_data = {}\n",
    "for country in response.elements():\n",
    "    if country.id() % 10 == 0:\n",
    "        countries_data[country.id() + 3600000000] = country.tags()['name:en'] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "residential_buildings_tags = \"apartments|bungalow|cabin|detached|dormitory|farm|ger|hotel|house|houseboat|residential|semidetached_house|static_caravan|terrace\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Count residential buildings in target countries\n",
    "buildings_count = dict.fromkeys(countries_data.keys())\n",
    "for id_ in countries_data.keys():\n",
    "    query = overpassQueryBuilder(area=id_, elementType=['node'], selector=f'\"building\"~\"{residential_buildings_tags}\"')\n",
    "    counter = overpass.query(query, timeout=120).countElements()\n",
    "    buildings_count[id_] = counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Malawi - country with the biggest number of resudential buildings\n"
     ]
    }
   ],
   "source": [
    "# Get country with the biggest number of resudential buildings\n",
    "country_id = max(buildings_count, key = lambda k: buildings_count[k])\n",
    "print(countries_data[country_id], '- country with the biggest number of resudential buildings')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all atms in target country\n",
    "response = overpass.query(f'area({country_id});(nwr[amenity=atm](area););out ids geom;', timeout=120)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find 5 closest atms to afghan random point\n",
    "idx = index.Index()\n",
    "for point in response.elements():\n",
    "    idx.insert(point.id(), point.geometry().coordinates)\n",
    "nearest_atms = idx.nearest((random_afghan_point.x, random_afghan_point.y), 5, objects = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summ distances\n",
    "result = sum([haversine((p.bbox[0], p.bbox[1]), (random_afghan_point.x, random_afghan_point.y), unit = Unit.METERS) for p in nearest_atms])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result sum 34920891.24153749\n"
     ]
    }
   ],
   "source": [
    "print('Result sum', result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Расстояние от точки до 5 ближайших школ, находящихся в стране с наибольшим количеством аптек в столице"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Unable to find area id for Canberra\n",
      "Unable to find area id for Mbabane\n",
      "Unable to find capital city for São Tomé and Príncipe\n",
      "Unable to find capital city for Tokelau\n",
      "Taiwan - country with the biggest number of pharmacies in capital\n"
     ]
    }
   ],
   "source": [
    "max_pharmacy = 0\n",
    "for id_ in countries_data:\n",
    "    # Get ciuntries capital\n",
    "    query = f'area({id_});(nwr[capital=yes](area); nwr[capital=country](area); nwr[capital=2](area); nwr[admin_level=2](area);); out body;'\n",
    "    res = overpass.query(query)\n",
    "    if (len(res.elements()) == 0):\n",
    "        print(f'Unable to find capital city for {countries_data[id_]}')\n",
    "        continue\n",
    "    capital = res.elements()[0].tags()['name:en']\n",
    "    # For each capital city get it's area id\n",
    "    areaId = nominatim.query(capital).areaId()\n",
    "    if (areaId == None):\n",
    "        print(f'Unable to find area id for {capital}')\n",
    "        continue\n",
    "    # Get all pharmacies in capital\n",
    "    query = f'area({id_});(nwr[shop=chemist](area);nwr[amenity=pharmacy](area);); out body;'\n",
    "    result = overpass.query(query).countElements()\n",
    "    if (result > max_pharmacy):\n",
    "        max_pharmacy = result\n",
    "        country_id = id_ \n",
    "print(countries_data[country_id], '- country with the biggest number of pharmacies in capital')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all schools in target country\n",
    "query = overpassQueryBuilder(area=country_id, elementType=['node', 'way'], selector='amenity=school',  includeGeometry=True)\n",
    "response = overpass.query(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find 5 closest schools to afghan random point\n",
    "idx = index.Index()\n",
    "for el in response.elements():\n",
    "    if (el.geometry().type == 'Pont'):\n",
    "        idx.insert(el.id(), el.geometry().coordinates)\n",
    "    if (el.geometry().type == 'Polygon'):\n",
    "        polygon_center = Polygon(list(chain.from_iterable(el.geometry().coordinates))).centroid\n",
    "        idx.insert(el.id(), [polygon_center.x, polygon_center.y])\n",
    "nearest_schools = idx.nearest((random_afghan_point.x, random_afghan_point.y), 5, objects = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summ distances\n",
    "result += sum([haversine((p.bbox[0], p.bbox[1]), (random_afghan_point.x, random_afghan_point.y), unit = Unit.METERS) for p in nearest_schools])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result sum 42972593.44071838\n"
     ]
    }
   ],
   "source": [
    "print('Result sum', result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 3. Расстояние от точки до 5 ближайших кинотеатров, наодящихся в стране с самым большим отношением числа железнодорожных станций к автобусным остановкам в южной части "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to get array of points of country bounds\n",
    "def get_county_bounds(county_name):\n",
    "    resp = overpass.query(f'relation[admin_level=2][type=boundary][boundary=administrative][\"name:en\"=\"{county_name}\"];out geom;')\n",
    "    if resp.elements()[0].geometry().type == 'Polygon':\n",
    "        bound_points = np.array(resp.elements()[0].geometry().coordinates[0])\n",
    "        bound_points.reshape(bound_points.shape[-2], 2)\n",
    "    elif resp.elements()[0].geometry().type == 'MultiPolygon':\n",
    "        bound_points = np.array([point for polygon in resp.elements()[0].geometry().coordinates for point in polygon[0]])\n",
    "    return bound_points"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to build triangle using extreme points ща сщгтекн ищгтвы\n",
    "def get_south_part_polygon(county_name):\n",
    "    bound_points = get_county_bounds(county_name)\n",
    "    max_lon = max(bound_points, key=lambda point:point[0])\n",
    "    max_lat = max(bound_points, key=lambda point:point[1])\n",
    "    min_lon = min(bound_points, key=lambda point:point[0])\n",
    "    min_lat = min(bound_points, key=lambda point:point[1])\n",
    "    return Polygon([[max_lon[0], max_lat[1]], [max_lon[0], min_lat[1]], [min_lon[0], min_lat[1]]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to count all points within polygon\n",
    "def count_points_within_polygon(response, polygon):\n",
    "    counter = 0\n",
    "    for el in response.elements():\n",
    "        if el.geometry().type == 'Point':\n",
    "            point = Point(el.geometry().coordinates[0], el.geometry().coordinates[1])\n",
    "        elif el.geometry().type == 'Polygon':\n",
    "            centroid = Polygon(el.geometry().coordinates[0]).centroid\n",
    "            point = Point(centroid.x, centroid.y)\n",
    "        if (point.within(polygon)):\n",
    "            counter += 1\n",
    "    return counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Australia - is the result country\n"
     ]
    }
   ],
   "source": [
    "max_value = -1\n",
    "for id_ in countries_data:\n",
    "    south_part = get_south_part_polygon(countries_data[id_])\n",
    "    # Get all bus stations\n",
    "    bus_query = overpassQueryBuilder(area=id_, elementType=['node', 'way'], selector='amenity=bus_station', out='geom')\n",
    "    # Get all train stations\n",
    "    train_query = overpassQueryBuilder(area=id_, elementType=['node', 'way'], selector='building=train_station', out='geom')\n",
    "    bus_stations_amount = count_points_within_polygon(overpass.query(bus_query), south_part)\n",
    "    train_stations_amount = count_points_within_polygon(overpass.query(train_query), south_part)\n",
    "    current_value = train_stations_amount / (bus_stations_amount + 0.001)\n",
    "    if (current_value > max_value):\n",
    "        max_value = current_value\n",
    "        country_id = id_\n",
    "print(countries_data[country_id], '- is the result country')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all cinemas in target country\n",
    "query = overpassQueryBuilder(area=country_id, elementType=['node', 'way'], selector='amenity=cinema',  includeGeometry=True)\n",
    "response = overpass.query(query)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find 5 closest cinemas to afghan random point\n",
    "idx = index.Index()\n",
    "for el in response.elements():\n",
    "    if (el.geometry().type == 'Pont'):\n",
    "        idx.insert(el.id(), el.geometry().coordinates)\n",
    "    if (el.geometry().type == 'Polygon'):\n",
    "        polygon_center = Polygon(list(chain.from_iterable(el.geometry().coordinates))).centroid\n",
    "        idx.insert(el.id(), [polygon_center.x, polygon_center.y])\n",
    "nearest_cinemas = idx.nearest((random_afghan_point.x, random_afghan_point.y), 5, objects = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Summ distances\n",
    "result += sum([haversine((p.bbox[0], p.bbox[1]), (random_afghan_point.x, random_afghan_point.y), unit = Unit.METERS) for p in nearest_cinemas])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result sum 78577549.44103174\n"
     ]
    }
   ],
   "source": [
    "print('Result sum', result)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}