E2SFCA WorldPoP-OSM Model 2.1-Func.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "675eb1c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import libpysal\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "import time\n",
    "import os\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon\n",
    "from shapely.ops import nearest_points, polygonize\n",
    "import shapely\n",
    "from itertools import product, combinations\n",
    "import math\n",
    "import warnings\n",
    "import socket\n",
    "from wpgpDownload.utils.dl import wpFtp\n",
    "from wpgpDownload.utils.isos import Countries\n",
    "from wpgpDownload.utils.convenience_functions import download_country_covariates as dl\n",
    "from wpgpDownload.utils.wpcsv import Product\n",
    "import georasters as gr\n",
    "from wpgpDownload.utils.convenience_functions import refresh_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "50e92f97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 0 cities and assumptions\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "cities = ['Tel Aviv','Philadelphia','Washington DC','Ghent','Dhaka Metropolitan']\n",
    "\n",
    "# idea to convert to dask-pandas and dask-geopandas\n",
    "# https://towardsdatascience.com/pandas-with-dask-for-an-ultra-fast-notebook-e2621c3769f\n",
    "# Or with Koalas (Spark-like pandas)\n",
    "\n",
    "# Assumptions\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "d841a73e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "if prefer dwnl from terminal: \n",
      "wpgpDownload download -i ISR --id 5089\n",
      "wpgpDownload download -i USA --id 4983\n",
      "wpgpDownload download -i BEL --id 5007\n",
      "wpgpDownload download -i BGD --id 5004\n",
      " \n",
      "downloaded:\n",
      "ISR downloaded 0.01 mns\n",
      "USA downloaded 2.05 mns\n",
      "BEL downloaded 2.14 mns\n",
      "BGD downloaded 2.17 mns\n"
     ]
    }
   ],
   "source": [
    "# 1. Required preprocess for information extraction\n",
    "\n",
    "# Let's ignore depreciation warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Get the city boundaries\n",
    "bound_df = ox.geocoder.geocode_to_gdf(cities) # gets city boundaries from OSM\n",
    "\n",
    "# Get unique iso-codes of selected cities (only load country raster once)\n",
    "unique_iso = iso_countries(bound_df, # Finding the country of the bounded city\n",
    "                           cities)\n",
    "print(' ')\n",
    "\n",
    "print('downloaded:')\n",
    "# Get raster of countries (if automatic download is preferred (standard))\n",
    "raster = countries_grids(unique_iso,\n",
    "                         r'D:\\Dumps\\WorldPoP_Grids') # custom path, where grid files can be stored without downloading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "9e396ef0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Tel Aviv 0.04 mns\n",
      "Philadelphia 0.27 mns\n",
      "Washington DC 0.36 mns\n",
      "Ghent 0.47 mns\n",
      "Dhaka Metropolitan 0.6 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Tel Aviv done 0.65 mns\n",
      "Philadelphia done 3.15 mns\n",
      "Washington DC done 5.25 mns\n",
      "Ghent done 6.03 mns\n",
      "Dhaka Metropolitan done 6.93 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Tel Aviv done\n",
      "Philadelphia done\n",
      "Washington DC done\n",
      "Ghent done\n",
      "Dhaka Metropolitan done\n"
     ]
    }
   ],
   "source": [
    "# 2. Information extraction\n",
    "\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(bound_df, # city boundaries\n",
    "                                     unique_iso,\n",
    "                                     raster, # country raster\n",
    "                                     cities, \n",
    "                                     grid_size = 100)\n",
    "print(' ')\n",
    "\n",
    "# Get road networks\n",
    "road_network = road_networks(cities, # Get 'all' (drive,walk,bike) network\n",
    "                                 thresholds,\n",
    "                                 undirected = True)\n",
    "\n",
    "\n",
    "print(' ')\n",
    "# Extracting UGS\n",
    "UGS = urban_greenspace(cities, \n",
    "                       thresholds,\n",
    "                       one_UGS_buf = 25, # buffer at which UGS is seen as one\n",
    "                       min_UGS_size = 400) # WHO sees this as minimum UGS size (400m2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "eb027bb5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get fake UGS entry points\n",
      "Tel Aviv 0.0 % done 0.0  mns\n",
      "Tel Aviv 12.7 % done 0.1  mns\n",
      "Tel Aviv 25.3 % done 0.18  mns\n",
      "Tel Aviv 38.0 % done 0.27  mns\n",
      "Tel Aviv 50.6 % done 0.35  mns\n",
      "Tel Aviv 63.3 % done 0.43  mns\n",
      "Tel Aviv 75.9 % done 0.51  mns\n",
      "Tel Aviv 88.6 % done 0.58  mns\n",
      "Tel Aviv 100 % done 0.7  mns\n",
      "Philadelphia 0.0 % done 0.72  mns\n",
      "Philadelphia 11.1 % done 1.09  mns\n",
      "Philadelphia 22.3 % done 1.36  mns\n",
      "Philadelphia 33.4 % done 1.65  mns\n",
      "Philadelphia 44.5 % done 1.95  mns\n",
      "Philadelphia 55.7 % done 2.21  mns\n",
      "Philadelphia 66.8 % done 2.47  mns\n",
      "Philadelphia 78.0 % done 2.72  mns\n",
      "Philadelphia 89.1 % done 2.97  mns\n",
      "Philadelphia 100 % done 3.3  mns\n",
      "Washington DC 0.0 % done 3.3  mns\n",
      "Washington DC 9.9 % done 3.77  mns\n",
      "Washington DC 19.8 % done 4.07  mns\n",
      "Washington DC 29.6 % done 4.31  mns\n",
      "Washington DC 39.5 % done 4.54  mns\n",
      "Washington DC 49.4 % done 4.79  mns\n",
      "Washington DC 59.3 % done 5.11  mns\n",
      "Washington DC 69.2 % done 5.39  mns\n",
      "Washington DC 79.1 % done 5.63  mns\n",
      "Washington DC 88.9 % done 5.88  mns\n",
      "Washington DC 98.8 % done 6.16  mns\n",
      "Washington DC 100 % done 6.27  mns\n",
      "Ghent 0.0 % done 6.27  mns\n",
      "Ghent 21.3 % done 6.39  mns\n",
      "Ghent 42.6 % done 6.49  mns\n",
      "Ghent 63.8 % done 6.59  mns\n",
      "Ghent 85.1 % done 6.69  mns\n",
      "Ghent 100 % done 6.78  mns\n",
      "Dhaka Metropolitan 0.0 % done 6.78  mns\n",
      "Dhaka Metropolitan 29.8 % done 6.89  mns\n",
      "Dhaka Metropolitan 59.5 % done 7.0  mns\n",
      "Dhaka Metropolitan 89.3 % done 7.1  mns\n",
      "Dhaka Metropolitan 100 % done 7.15  mns\n",
      " \n",
      "get potential (Euclidean) suitible combinations\n",
      "Tel Aviv\n",
      "in chunk 1 / 4 221868 suitible comb.\n",
      "in chunk 2 / 4 258024 suitible comb.\n",
      "in chunk 3 / 4 263214 suitible comb.\n",
      "in chunk 4 / 4 284885 suitible comb.\n",
      "total combinations within distance 1027991\n",
      "0.0 % gridentry done 0.0  mns\n",
      "24.3 % gridentry done 0.37  mns\n",
      "48.6 % gridentry done 0.77  mns\n",
      "73.0 % gridentry done 1.17  mns\n",
      "97.3 % gridentry done 1.57  mns\n",
      "100 % gridentry done 2.78  mns\n",
      "Philadelphia\n",
      "in chunk 1 / 21 22427 suitible comb.\n",
      "in chunk 2 / 21 27735 suitible comb.\n",
      "in chunk 3 / 21 39002 suitible comb.\n",
      "in chunk 4 / 21 41583 suitible comb.\n",
      "in chunk 5 / 21 41232 suitible comb.\n",
      "in chunk 6 / 21 42745 suitible comb.\n",
      "in chunk 7 / 21 40047 suitible comb.\n",
      "in chunk 8 / 21 21229 suitible comb.\n",
      "in chunk 9 / 21 24934 suitible comb.\n",
      "in chunk 10 / 21 38874 suitible comb.\n",
      "in chunk 11 / 21 46802 suitible comb.\n",
      "in chunk 12 / 21 73938 suitible comb.\n",
      "in chunk 13 / 21 115472 suitible comb.\n",
      "in chunk 14 / 21 54155 suitible comb.\n",
      "in chunk 15 / 21 42204 suitible comb.\n",
      "in chunk 16 / 21 3549 suitible comb.\n",
      "in chunk 17 / 21 7429 suitible comb.\n",
      "in chunk 18 / 21 13983 suitible comb.\n",
      "in chunk 19 / 21 16260 suitible comb.\n",
      "in chunk 20 / 21 16033 suitible comb.\n",
      "in chunk 21 / 21 21199 suitible comb.\n",
      "total combinations within distance 750832\n",
      "0.0 % gridentry done 0.0  mns\n",
      "33.3 % gridentry done 0.38  mns\n",
      "66.6 % gridentry done 0.78  mns\n",
      "99.9 % gridentry done 1.15  mns\n",
      "100 % gridentry done 11.57  mns\n",
      "Washington DC\n",
      "in chunk 1 / 14 79816 suitible comb.\n",
      "in chunk 2 / 14 85313 suitible comb.\n",
      "in chunk 3 / 14 83642 suitible comb.\n",
      "in chunk 4 / 14 75099 suitible comb.\n",
      "in chunk 5 / 14 36093 suitible comb.\n",
      "in chunk 6 / 14 25261 suitible comb.\n",
      "in chunk 7 / 14 35141 suitible comb.\n",
      "in chunk 8 / 14 49627 suitible comb.\n",
      "in chunk 9 / 14 45442 suitible comb.\n",
      "in chunk 10 / 14 49858 suitible comb.\n",
      "in chunk 11 / 14 52554 suitible comb.\n",
      "in chunk 12 / 14 42279 suitible comb.\n",
      "in chunk 13 / 14 47738 suitible comb.\n",
      "in chunk 14 / 14 58537 suitible comb.\n",
      "total combinations within distance 766400\n",
      "0.0 % gridentry done 0.0  mns\n",
      "32.6 % gridentry done 0.37  mns\n",
      "65.2 % gridentry done 0.75  mns\n",
      "97.9 % gridentry done 1.14  mns\n",
      "100 % gridentry done 17.34  mns\n",
      "Ghent\n",
      "in chunk 1 / 9 66203 suitible comb.\n",
      "in chunk 2 / 9 154803 suitible comb.\n",
      "in chunk 3 / 9 227987 suitible comb.\n",
      "in chunk 4 / 9 257571 suitible comb.\n",
      "in chunk 5 / 9 278548 suitible comb.\n",
      "in chunk 6 / 9 241102 suitible comb.\n",
      "in chunk 7 / 9 98915 suitible comb.\n",
      "in chunk 8 / 9 5059 suitible comb.\n",
      "in chunk 9 / 9 9822 suitible comb.\n",
      "total combinations within distance 1340010\n",
      "0.0 % gridentry done 0.0  mns\n",
      "18.7 % gridentry done 0.39  mns\n",
      "37.3 % gridentry done 0.76  mns\n",
      "56.0 % gridentry done 1.13  mns\n",
      "74.6 % gridentry done 1.52  mns\n",
      "93.3 % gridentry done 1.91  mns\n",
      "100 % gridentry done 22.32  mns\n",
      "Dhaka Metropolitan\n",
      "in chunk 1 / 3 9399 suitible comb.\n",
      "in chunk 2 / 3 9452 suitible comb.\n",
      "in chunk 3 / 3 3863 suitible comb.\n",
      "total combinations within distance 22714\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 22.85  mns\n",
      " \n",
      "Check grids within UGS\n",
      "0 0.0  mns\n",
      "100 0.07  mns\n",
      "200 0.14  mns\n",
      "300 0.2  mns\n",
      "Check grids within UGS\n",
      "0 0.27  mns\n",
      "100 0.71  mns\n",
      "200 0.99  mns\n",
      "300 1.24  mns\n",
      "400 1.49  mns\n",
      "Check grids within UGS\n",
      "0 1.65  mns\n",
      "100 2.05  mns\n",
      "200 2.2  mns\n",
      "300 2.35  mns\n",
      "400 2.5  mns\n",
      "500 2.64  mns\n",
      "Check grids within UGS\n",
      "0 2.67  mns\n",
      "100 2.95  mns\n",
      "200 3.21  mns\n",
      "Check grids within UGS\n",
      "0 3.31  mns\n",
      "100 3.56  mns\n"
     ]
    }
   ],
   "source": [
    "# 3. Preprocess information for route finding\n",
    "\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_network['nodes'], \n",
    "                           cities, \n",
    "                           UGS_entry_buf = 25, # road nodes within 25 meters are seen as fake entry points\n",
    "                           walk_radius = 500, # assume that the average person only views a UGS up to 500m in radius\n",
    "                                                # more attractive\n",
    "                           entry_point_merge = 0) # merges closeby fake UGS entry points within X meters \n",
    "                                                    # what may be done for performance\n",
    "print(' ')\n",
    "# Checks all potential suitible combinations (points that fall within max threshold Euclidean distance from the ego)\n",
    "suitible = suitible_combinations(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_network['nodes'], # For finding nearest grid entry points\n",
    "                                 thresholds,\n",
    "                                 cities,\n",
    "                                 chunk_size = 10000000) # calculating per chunk of num UGS entry points * num pop_grids\n",
    "                                                        # Preventing normal PC meltdown, set lower if PC gets stuck\n",
    "print(' ')\n",
    "# Checks if grids are already in a UGS\n",
    "suitible_InOut_UGS = grids_in_UGS (suitible, UGS, population_grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "40c30628",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tel Aviv 1 / 4 range 0 - 250000\n",
      "0.0 % done 0.01 mns\n",
      "1.01 % done 0.28 mns\n",
      "2.02 % done 0.57 mns\n",
      "3.04 % done 0.78 mns\n",
      "4.05 % done 1.09 mns\n",
      "5.06 % done 1.29 mns\n",
      "6.07 % done 1.44 mns\n",
      "7.09 % done 1.59 mns\n",
      "8.1 % done 1.73 mns\n",
      "9.11 % done 1.88 mns\n",
      "10.12 % done 2.01 mns\n",
      "11.14 % done 2.19 mns\n",
      "12.15 % done 2.41 mns\n",
      "13.16 % done 2.63 mns\n",
      "14.17 % done 2.81 mns\n",
      "15.19 % done 3.02 mns\n",
      "16.2 % done 3.2 mns\n",
      "17.21 % done 3.37 mns\n",
      "18.22 % done 3.56 mns\n",
      "19.24 % done 3.7 mns\n",
      "20.25 % done 3.95 mns\n",
      "21.26 % done 4.42 mns\n",
      "22.27 % done 4.67 mns\n",
      "23.29 % done 4.81 mns\n",
      "24.3 % done 4.95 mns\n",
      "for 228 routes nearest nodes found\n",
      "25.31 % pathfinding done 5.08 mns\n",
      "formatting done 6.67 mns\n",
      "dissolving done 8.15 mns\n",
      "Tel Aviv 2 / 4 range 250000 - 500000\n",
      "25.31 % done 8.15 mns\n",
      "26.32 % done 8.3 mns\n",
      "27.33 % done 8.46 mns\n",
      "28.35 % done 8.64 mns\n",
      "29.36 % done 8.8 mns\n",
      "30.37 % done 8.96 mns\n",
      "31.38 % done 9.11 mns\n",
      "32.4 % done 9.25 mns\n",
      "33.41 % done 9.39 mns\n",
      "34.42 % done 9.55 mns\n",
      "35.43 % done 9.69 mns\n",
      "36.45 % done 9.84 mns\n",
      "37.46 % done 9.98 mns\n",
      "38.47 % done 10.14 mns\n",
      "39.48 % done 10.33 mns\n",
      "40.5 % done 10.51 mns\n",
      "41.51 % done 10.67 mns\n",
      "42.52 % done 10.81 mns\n",
      "43.53 % done 10.98 mns\n",
      "44.55 % done 11.11 mns\n",
      "45.56 % done 11.77 mns\n",
      "46.57 % done 11.97 mns\n",
      "47.58 % done 12.12 mns\n",
      "48.6 % done 12.28 mns\n",
      "49.61 % done 12.44 mns\n",
      "for 173 routes nearest nodes found\n",
      "50.62 % pathfinding done 12.6 mns\n",
      "formatting done 14.2 mns\n",
      "dissolving done 15.61 mns\n",
      "Tel Aviv 3 / 4 range 500000 - 750000\n",
      "50.62 % done 15.61 mns\n",
      "51.63 % done 15.77 mns\n",
      "52.64 % done 15.93 mns\n",
      "53.66 % done 16.07 mns\n",
      "54.67 % done 16.19 mns\n",
      "55.68 % done 16.32 mns\n",
      "56.69 % done 16.45 mns\n",
      "57.71 % done 16.61 mns\n",
      "58.72 % done 21.75 mns\n",
      "59.73 % done 21.94 mns\n",
      "60.74 % done 22.22 mns\n",
      "61.76 % done 22.43 mns\n",
      "62.77 % done 22.57 mns\n",
      "63.78 % done 22.76 mns\n",
      "64.79 % done 22.95 mns\n",
      "65.81 % done 23.17 mns\n",
      "66.82 % done 23.38 mns\n",
      "67.83 % done 23.63 mns\n",
      "68.84 % done 23.79 mns\n",
      "69.86 % done 23.97 mns\n",
      "70.87 % done 24.15 mns\n",
      "71.88 % done 24.33 mns\n",
      "72.89 % done 24.49 mns\n",
      "73.91 % done 24.66 mns\n",
      "74.92 % done 24.85 mns\n",
      "for 513 routes nearest nodes found\n",
      "75.93 % pathfinding done 25.05 mns\n",
      "formatting done 26.7 mns\n",
      "dissolving done 28.13 mns\n",
      "Tel Aviv 4 / 4 range 750000 - 987751\n",
      "75.93 % done 28.13 mns\n",
      "76.94 % done 28.34 mns\n",
      "77.95 % done 28.55 mns\n",
      "78.97 % done 28.73 mns\n",
      "79.98 % done 28.92 mns\n",
      "80.99 % done 29.1 mns\n",
      "82.0 % done 29.36 mns\n",
      "83.02 % done 29.84 mns\n",
      "84.03 % done 30.01 mns\n",
      "85.04 % done 30.19 mns\n",
      "86.05 % done 30.33 mns\n",
      "87.07 % done 30.55 mns\n",
      "88.08 % done 30.72 mns\n",
      "89.09 % done 30.87 mns\n",
      "90.1 % done 31.02 mns\n",
      "91.12 % done 31.2 mns\n",
      "92.13 % done 31.34 mns\n",
      "93.14 % done 31.55 mns\n",
      "94.15 % done 31.74 mns\n",
      "95.17 % done 31.92 mns\n",
      "96.18 % done 32.13 mns\n",
      "97.19 % done 32.32 mns\n",
      "98.2 % done 32.54 mns\n",
      "99.22 % done 32.74 mns\n",
      "for 198 routes nearest nodes found\n",
      "100.0 % pathfinding done 32.93 mns\n",
      "formatting done 34.35 mns\n",
      "dissolving done 35.66 mns\n",
      "Philadelphia 1 / 3 range 0 - 250000\n",
      "0.0 % done 35.68 mns\n",
      "1.37 % done 35.75 mns\n",
      "2.75 % done 35.87 mns\n",
      "4.12 % done 35.97 mns\n",
      "5.49 % done 36.08 mns\n",
      "6.87 % done 36.18 mns\n",
      "8.24 % done 36.31 mns\n",
      "9.61 % done 36.43 mns\n",
      "10.99 % done 36.56 mns\n",
      "12.36 % done 36.65 mns\n",
      "13.73 % done 36.78 mns\n",
      "15.11 % done 36.88 mns\n",
      "16.48 % done 37.02 mns\n",
      "17.85 % done 37.18 mns\n",
      "19.23 % done 37.3 mns\n",
      "20.6 % done 37.47 mns\n",
      "21.97 % done 37.75 mns\n",
      "23.35 % done 38.38 mns\n",
      "24.72 % done 38.45 mns\n",
      "26.09 % done 38.6 mns\n",
      "27.47 % done 38.68 mns\n",
      "28.84 % done 38.82 mns\n",
      "30.22 % done 38.92 mns\n",
      "31.59 % done 39.0 mns\n",
      "32.96 % done 39.11 mns\n",
      "for 166 routes nearest nodes found\n",
      "34.34 % pathfinding done 39.19 mns\n",
      "formatting done 40.41 mns\n",
      "dissolving done 41.67 mns\n",
      "Philadelphia 2 / 3 range 250000 - 500000\n",
      "34.34 % done 41.67 mns\n",
      "35.71 % done 41.74 mns\n",
      "37.08 % done 41.81 mns\n",
      "38.46 % done 41.9 mns\n",
      "39.83 % done 42.03 mns\n",
      "41.2 % done 42.1 mns\n",
      "42.58 % done 42.19 mns\n",
      "43.95 % done 42.32 mns\n",
      "45.32 % done 42.4 mns\n",
      "46.7 % done 42.54 mns\n",
      "48.07 % done 42.64 mns\n",
      "49.44 % done 42.74 mns\n",
      "50.82 % done 42.86 mns\n",
      "52.19 % done 42.96 mns\n",
      "53.56 % done 43.12 mns\n",
      "54.94 % done 43.31 mns\n",
      "56.31 % done 43.47 mns\n",
      "57.68 % done 43.64 mns\n",
      "59.06 % done 43.78 mns\n",
      "60.43 % done 43.92 mns\n",
      "61.8 % done 44.06 mns\n",
      "63.18 % done 44.22 mns\n",
      "64.55 % done 44.41 mns\n",
      "65.92 % done 44.57 mns\n",
      "67.3 % done 44.73 mns\n",
      "for 0 routes nearest nodes found\n",
      "68.67 % pathfinding done 44.91 mns\n",
      "formatting done 46.12 mns\n",
      "dissolving done 47.44 mns\n",
      "Philadelphia 3 / 3 range 500000 - 728113\n",
      "68.67 % done 47.45 mns\n",
      "70.04 % done 47.63 mns\n",
      "71.42 % done 47.8 mns\n",
      "72.79 % done 48.02 mns\n",
      "74.16 % done 48.18 mns\n",
      "75.54 % done 48.33 mns\n",
      "76.91 % done 48.67 mns\n",
      "78.28 % done 48.8 mns\n",
      "79.66 % done 48.91 mns\n",
      "81.03 % done 49.02 mns\n",
      "82.4 % done 49.12 mns\n",
      "83.78 % done 49.21 mns\n",
      "85.15 % done 49.49 mns\n",
      "86.53 % done 49.62 mns\n",
      "87.9 % done 49.72 mns\n",
      "89.27 % done 50.41 mns\n",
      "657313 No route 10\n",
      "657314 No route 10\n",
      "657317 No route 10\n",
      "657321 No route 10\n",
      "657322 No route 10\n",
      "658764 No route 10\n",
      "658765 No route 10\n",
      "658766 No route 10\n",
      "658767 No route 10\n",
      "658768 No route 10\n",
      "658769 No route 10\n",
      "658770 No route 10\n",
      "658771 No route 10\n",
      "658772 No route 10\n",
      "658773 No route 10\n",
      "658774 No route 10\n",
      "658775 No route 10\n",
      "658776 No route 10\n",
      "658777 No route 10\n",
      "658778 No route 10\n",
      "658779 No route 10\n",
      "658780 No route 10\n",
      "658781 No route 10\n",
      "658782 No route 10\n",
      "658783 No route 10\n",
      "658784 No route 10\n",
      "658785 No route 10\n",
      "658786 No route 10\n",
      "658787 No route 10\n",
      "658788 No route 10\n",
      "658789 No route 10\n",
      "658790 No route 10\n",
      "658791 No route 10\n",
      "658792 No route 10\n",
      "658793 No route 10\n",
      "658794 No route 10\n",
      "658795 No route 10\n",
      "658796 No route 10\n",
      "658797 No route 10\n",
      "658798 No route 10\n",
      "658799 No route 10\n",
      "658800 No route 10\n",
      "658801 No route 10\n",
      "658802 No route 10\n",
      "658803 No route 10\n",
      "658832 No route 10\n",
      "90.65 % done 53.7 mns\n",
      "92.02 % done 53.98 mns\n",
      "93.39 % done 54.17 mns\n",
      "94.77 % done 54.88 mns\n",
      "96.14 % done 54.97 mns\n",
      "97.51 % done 55.15 mns\n",
      "98.89 % done 55.23 mns\n",
      "for 572 routes nearest nodes found\n",
      "100.0 % pathfinding done 55.31 mns\n",
      "formatting done 56.41 mns\n",
      "dissolving done 57.68 mns\n",
      "Washington DC 1 / 3 range 0 - 250000\n",
      "0.0 % done 57.73 mns\n",
      "1.37 % done 57.92 mns\n",
      "2.75 % done 58.07 mns\n",
      "4.12 % done 58.31 mns\n",
      "5.5 % done 58.58 mns\n",
      "6.87 % done 58.86 mns\n",
      "8.25 % done 58.96 mns\n",
      "9.62 % done 59.13 mns\n",
      "11.0 % done 59.26 mns\n",
      "12.37 % done 59.49 mns\n",
      "13.74 % done 59.77 mns\n",
      "15.12 % done 60.05 mns\n",
      "16.49 % done 60.32 mns\n",
      "17.87 % done 60.43 mns\n",
      "19.24 % done 60.58 mns\n",
      "20.62 % done 60.72 mns\n",
      "21.99 % done 60.9 mns\n",
      "23.36 % done 61.05 mns\n",
      "24.74 % done 61.27 mns\n",
      "26.11 % done 61.46 mns\n",
      "27.49 % done 61.56 mns\n",
      "28.86 % done 61.66 mns\n",
      "30.24 % done 61.79 mns\n",
      "31.61 % done 61.92 mns\n",
      "32.99 % done 62.03 mns\n",
      "for 28 routes nearest nodes found\n",
      "34.36 % pathfinding done 62.15 mns\n",
      "formatting done 63.69 mns\n",
      "dissolving done 65.23 mns\n",
      "Washington DC 2 / 3 range 250000 - 500000\n",
      "34.36 % done 65.23 mns\n",
      "35.73 % done 65.38 mns\n",
      "37.11 % done 65.51 mns\n",
      "38.48 % done 65.64 mns\n",
      "39.86 % done 65.72 mns\n",
      "41.23 % done 65.8 mns\n",
      "42.61 % done 65.94 mns\n",
      "43.98 % done 66.04 mns\n",
      "45.35 % done 66.1 mns\n",
      "46.73 % done 66.19 mns\n",
      "48.1 % done 66.3 mns\n",
      "49.48 % done 66.38 mns\n",
      "50.85 % done 66.49 mns\n",
      "52.23 % done 66.55 mns\n",
      "53.6 % done 66.7 mns\n",
      "54.98 % done 66.82 mns\n",
      "56.35 % done 66.91 mns\n",
      "57.72 % done 67.05 mns\n",
      "59.1 % done 67.15 mns\n",
      "60.47 % done 67.26 mns\n",
      "61.85 % done 67.35 mns\n",
      "63.22 % done 67.45 mns\n",
      "64.6 % done 67.54 mns\n",
      "65.97 % done 67.6 mns\n",
      "67.35 % done 67.71 mns\n",
      "for 0 routes nearest nodes found\n",
      "68.72 % pathfinding done 67.8 mns\n",
      "formatting done 68.98 mns\n",
      "dissolving done 70.36 mns\n",
      "Washington DC 3 / 3 range 500000 - 727595\n",
      "68.72 % done 70.37 mns\n",
      "70.09 % done 70.47 mns\n",
      "71.47 % done 70.56 mns\n",
      "72.84 % done 70.71 mns\n",
      "74.22 % done 70.84 mns\n",
      "75.59 % done 70.94 mns\n",
      "76.97 % done 71.05 mns\n",
      "78.34 % done 71.17 mns\n",
      "79.71 % done 71.3 mns\n",
      "81.09 % done 71.45 mns\n",
      "82.46 % done 71.55 mns\n",
      "83.84 % done 71.67 mns\n",
      "85.21 % done 71.79 mns\n",
      "86.59 % done 71.9 mns\n",
      "87.96 % done 72.01 mns\n",
      "89.34 % done 72.08 mns\n",
      "90.71 % done 72.2 mns\n",
      "92.08 % done 72.33 mns\n",
      "93.46 % done 72.43 mns\n",
      "94.83 % done 72.53 mns\n",
      "96.21 % done 72.72 mns\n",
      "97.58 % done 72.88 mns\n",
      "98.96 % done 73.07 mns\n",
      "for 0 routes nearest nodes found\n",
      "100.0 % pathfinding done 73.15 mns\n",
      "formatting done 74.28 mns\n",
      "dissolving done 75.55 mns\n",
      "Ghent 1 / 6 range 0 - 250000\n",
      "0.0 % done 76.2 mns\n",
      "0.76 % done 77.02 mns\n",
      "1.52 % done 77.11 mns\n",
      "2.28 % done 77.21 mns\n",
      "3.04 % done 77.36 mns\n",
      "3.81 % done 77.48 mns\n",
      "4.57 % done 77.59 mns\n",
      "5.33 % done 77.74 mns\n",
      "6.09 % done 77.92 mns\n",
      "6.85 % done 78.08 mns\n",
      "7.61 % done 78.26 mns\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.37 % done 78.48 mns\n",
      "9.13 % done 78.72 mns\n",
      "9.89 % done 78.88 mns\n",
      "10.66 % done 79.03 mns\n",
      "11.42 % done 79.16 mns\n",
      "12.18 % done 79.3 mns\n",
      "12.94 % done 79.43 mns\n",
      "13.7 % done 79.58 mns\n",
      "14.46 % done 79.72 mns\n",
      "15.22 % done 79.84 mns\n",
      "15.98 % done 79.96 mns\n",
      "16.74 % done 80.12 mns\n",
      "17.51 % done 80.26 mns\n",
      "18.27 % done 80.49 mns\n",
      "for 656 routes nearest nodes found\n",
      "19.03 % pathfinding done 80.71 mns\n",
      "formatting done 82.46 mns\n",
      "dissolving done 84.04 mns\n",
      "Ghent 2 / 6 range 250000 - 500000\n",
      "19.03 % done 84.05 mns\n",
      "19.79 % done 84.28 mns\n",
      "20.55 % done 84.52 mns\n",
      "21.31 % done 84.79 mns\n",
      "22.07 % done 85.06 mns\n",
      "22.83 % done 85.17 mns\n",
      "23.6 % done 85.26 mns\n",
      "24.36 % done 85.42 mns\n",
      "25.12 % done 85.59 mns\n",
      "25.88 % done 85.85 mns\n",
      "26.64 % done 86.06 mns\n",
      "27.4 % done 86.32 mns\n",
      "28.16 % done 86.62 mns\n",
      "28.92 % done 86.88 mns\n",
      "29.68 % done 87.05 mns\n",
      "30.45 % done 87.24 mns\n",
      "31.21 % done 87.44 mns\n",
      "31.97 % done 87.63 mns\n",
      "32.73 % done 87.9 mns\n",
      "33.49 % done 88.07 mns\n",
      "34.25 % done 88.23 mns\n",
      "35.01 % done 88.51 mns\n",
      "35.77 % done 88.79 mns\n",
      "36.53 % done 89.07 mns\n",
      "37.3 % done 89.34 mns\n",
      "for 0 routes nearest nodes found\n",
      "38.06 % pathfinding done 89.59 mns\n",
      "formatting done 91.3 mns\n",
      "dissolving done 92.95 mns\n",
      "Ghent 3 / 6 range 500000 - 750000\n",
      "38.06 % done 92.96 mns\n",
      "38.82 % done 93.21 mns\n",
      "39.58 % done 93.47 mns\n",
      "40.34 % done 93.73 mns\n",
      "41.1 % done 93.94 mns\n",
      "41.86 % done 94.18 mns\n",
      "42.62 % done 94.35 mns\n",
      "43.38 % done 94.5 mns\n",
      "44.15 % done 94.64 mns\n",
      "44.91 % done 94.97 mns\n",
      "45.67 % done 95.27 mns\n",
      "46.43 % done 95.57 mns\n",
      "47.19 % done 95.86 mns\n",
      "47.95 % done 96.21 mns\n",
      "48.71 % done 96.43 mns\n",
      "49.47 % done 96.65 mns\n",
      "50.23 % done 96.86 mns\n",
      "51.0 % done 97.14 mns\n",
      "51.76 % done 97.41 mns\n",
      "52.52 % done 97.62 mns\n",
      "53.28 % done 97.77 mns\n",
      "54.04 % done 98.04 mns\n",
      "54.8 % done 98.31 mns\n",
      "55.56 % done 98.58 mns\n",
      "56.32 % done 98.83 mns\n",
      "for 0 routes nearest nodes found\n",
      "57.08 % pathfinding done 99.17 mns\n",
      "formatting done 100.87 mns\n",
      "dissolving done 102.49 mns\n",
      "Ghent 4 / 6 range 750000 - 1000000\n",
      "57.08 % done 102.5 mns\n",
      "57.85 % done 102.85 mns\n",
      "58.61 % done 103.15 mns\n",
      "59.37 % done 103.5 mns\n",
      "60.13 % done 103.78 mns\n",
      "60.89 % done 104.07 mns\n",
      "61.65 % done 104.36 mns\n",
      "62.41 % done 104.64 mns\n",
      "63.17 % done 104.92 mns\n",
      "63.93 % done 105.16 mns\n",
      "64.7 % done 105.43 mns\n",
      "65.46 % done 105.71 mns\n",
      "66.22 % done 105.99 mns\n",
      "66.98 % done 106.2 mns\n",
      "67.74 % done 106.4 mns\n",
      "68.5 % done 106.6 mns\n",
      "69.26 % done 106.92 mns\n",
      "70.02 % done 107.27 mns\n",
      "70.79 % done 107.59 mns\n",
      "71.55 % done 107.83 mns\n",
      "72.31 % done 107.95 mns\n",
      "73.07 % done 108.09 mns\n",
      "965561 No route 10\n",
      "965562 No route 10\n",
      "965563 No route 10\n",
      "965564 No route 10\n",
      "965565 No route 10\n",
      "965566 No route 10\n",
      "965567 No route 10\n",
      "965568 No route 10\n",
      "965569 No route 10\n",
      "965570 No route 10\n",
      "965571 No route 10\n",
      "965572 No route 10\n",
      "965573 No route 10\n",
      "965574 No route 10\n",
      "965575 No route 10\n",
      "965576 No route 10\n",
      "965577 No route 10\n",
      "965578 No route 10\n",
      "965579 No route 10\n",
      "965580 No route 10\n",
      "965581 No route 10\n",
      "965582 No route 10\n",
      "965583 No route 10\n",
      "965584 No route 10\n",
      "965585 No route 10\n",
      "965586 No route 10\n",
      "965587 No route 10\n",
      "965588 No route 10\n",
      "965589 No route 10\n",
      "965590 No route 10\n",
      "965591 No route 10\n",
      "965592 No route 10\n",
      "965593 No route 10\n",
      "965594 No route 10\n",
      "965595 No route 10\n",
      "965596 No route 10\n",
      "965597 No route 10\n",
      "965598 No route 10\n",
      "965599 No route 10\n",
      "965600 No route 10\n",
      "965601 No route 10\n",
      "965602 No route 10\n",
      "965603 No route 10\n",
      "965604 No route 10\n",
      "965605 No route 10\n",
      "965606 No route 10\n",
      "965607 No route 10\n",
      "965608 No route 10\n",
      "965609 No route 10\n",
      "965610 No route 10\n",
      "965611 No route 10\n",
      "965612 No route 10\n",
      "965613 No route 10\n",
      "965614 No route 10\n",
      "965615 No route 10\n",
      "965616 No route 10\n",
      "965617 No route 10\n",
      "965618 No route 10\n",
      "965619 No route 10\n",
      "965620 No route 10\n",
      "965621 No route 10\n",
      "965622 No route 10\n",
      "965623 No route 10\n",
      "965624 No route 10\n",
      "965625 No route 10\n",
      "965626 No route 10\n",
      "965627 No route 10\n",
      "965628 No route 10\n",
      "965629 No route 10\n",
      "965630 No route 10\n",
      "965631 No route 10\n",
      "965632 No route 10\n",
      "965633 No route 10\n",
      "965634 No route 10\n",
      "965635 No route 10\n",
      "965636 No route 10\n",
      "965637 No route 10\n",
      "965638 No route 10\n",
      "965639 No route 10\n",
      "965640 No route 10\n",
      "965641 No route 10\n",
      "965642 No route 10\n",
      "965643 No route 10\n",
      "965644 No route 10\n",
      "965645 No route 10\n",
      "965646 No route 10\n",
      "965647 No route 10\n",
      "965648 No route 10\n",
      "965649 No route 10\n",
      "965650 No route 10\n",
      "965651 No route 10\n",
      "965652 No route 10\n",
      "965653 No route 10\n",
      "965654 No route 10\n",
      "965655 No route 10\n",
      "965656 No route 10\n",
      "965657 No route 10\n",
      "965946 No route 10\n",
      "965947 No route 10\n",
      "965948 No route 10\n",
      "73.83 % done 110.2 mns\n",
      "74.59 % done 110.42 mns\n",
      "75.35 % done 115.27 mns\n",
      "for 1984 routes nearest nodes found\n",
      "76.11 % pathfinding done 115.47 mns\n",
      "formatting done 117.16 mns\n",
      "dissolving done 118.77 mns\n",
      "Ghent 5 / 6 range 1000000 - 1250000\n",
      "76.11 % done 118.77 mns\n",
      "76.87 % done 118.99 mns\n",
      "77.64 % done 119.25 mns\n",
      "78.4 % done 119.52 mns\n",
      "79.16 % done 119.86 mns\n",
      "79.92 % done 120.18 mns\n",
      "80.68 % done 120.36 mns\n",
      "81.44 % done 120.56 mns\n",
      "82.2 % done 120.65 mns\n",
      "82.96 % done 120.74 mns\n",
      "83.72 % done 121.02 mns\n",
      "84.49 % done 121.3 mns\n",
      "85.25 % done 121.57 mns\n",
      "86.01 % done 121.83 mns\n",
      "86.77 % done 122.09 mns\n",
      "87.53 % done 122.4 mns\n",
      "88.29 % done 122.73 mns\n",
      "89.05 % done 123.0 mns\n",
      "89.81 % done 123.22 mns\n",
      "90.57 % done 123.5 mns\n",
      "91.34 % done 123.73 mns\n",
      "1206550 No route 10\n",
      "1206551 No route 10\n",
      "1206552 No route 10\n",
      "1206553 No route 10\n",
      "92.1 % done 124.66 mns\n",
      "92.86 % done 124.73 mns\n",
      "93.62 % done 124.8 mns\n",
      "94.38 % done 124.91 mns\n",
      "for 425 routes nearest nodes found\n",
      "95.14 % pathfinding done 125.11 mns\n",
      "formatting done 126.78 mns\n",
      "dissolving done 128.35 mns\n",
      "Ghent 6 / 6 range 1250000 - 1313836\n",
      "95.14 % done 128.36 mns\n",
      "95.9 % done 128.6 mns\n",
      "96.66 % done 128.79 mns\n",
      "97.42 % done 128.99 mns\n",
      "98.19 % done 129.97 mns\n",
      "98.95 % done 130.74 mns\n",
      "99.71 % done 131.77 mns\n",
      "for 1000 routes nearest nodes found\n",
      "100.0 % pathfinding done 132.64 mns\n",
      "formatting done 133.11 mns\n",
      "dissolving done 133.57 mns\n",
      "Dhaka Metropolitan 1 / 1 range 0 - 21615\n",
      "0.0 % done 133.73 mns\n",
      "46.26 % done 133.77 mns\n",
      "92.53 % done 133.81 mns\n",
      "for 0 routes nearest nodes found\n",
      "100.0 % pathfinding done 133.81 mns\n",
      "formatting done 134.02 mns\n",
      "dissolving done 134.25 mns\n"
     ]
    }
   ],
   "source": [
    "# 4. Finding shortest routes.\n",
    "Routes = route_finding (road_network['graphs'], # graphs of the road networks\n",
    "               suitible_InOut_UGS, # potential suitible routes with grid-UGS comb. separated in or out UGS.\n",
    "               road_network['nodes'], \n",
    "               road_network['edges'], \n",
    "               cities, \n",
    "               block_size = 250000, # Chunk to spread dataload.\n",
    "               nn_iter = 10) # max amount of nearest nodes to be found (both for UGS entry and grid-centroid road entries)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "43e3b723",
   "metadata": {},
   "outputs": [],
   "source": [
    "gpd.GeoDataFrame(UGS_entry[4][['geometry_x']], geometry = 'geometry_x').to_file('D:/Dumps/Scores output WP2-OSM/Park entry points/Dhaka.shp')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08decce9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300 Tel Aviv\n",
      "600 Tel Aviv\n",
      "1000 Tel Aviv\n",
      "300 Philadelphia\n",
      "600 Philadelphia\n",
      "1000 Philadelphia\n",
      "300 Washington DC\n",
      "600 Washington DC\n",
      "1000 Washington DC\n",
      "300 Ghent\n",
      "600 Ghent\n",
      "1000 Ghent\n",
      "300 Dhaka Metropolitan\n",
      "600 Dhaka Metropolitan\n",
      "1000 Dhaka Metropolitan\n"
     ]
    }
   ],
   "source": [
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, population_grids, thresholds, cities)\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e7c8c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def iso_countries(bounds, cities):\n",
    "    # bound_df = ox.geocoder.geocode_to_gdf(cities)\n",
    "    # The 'Countries' is a list of iso-countries and descriptions from the package wpgpDownload.utils.isos\n",
    "    C = pd.DataFrame(Countries)\n",
    "    start_time = time.time()\n",
    "    iso_countries = []\n",
    "    print('if prefer dwnl from terminal: ')\n",
    "    \n",
    "    # Check the display name in the city boundaries to get the country name (enabling only specifying city in front)\n",
    "    for i in bounds['display_name']:\n",
    "        country = i.rsplit(',')[-1][1:]\n",
    "        iso = C[C['name'] == country].iloc[0,1]\n",
    "        # Get unique ISO countries, so all country-grids are only loaded once\n",
    "        if iso not in iso_countries:\n",
    "            iso_countries.append(iso)\n",
    "            \n",
    "            # List data and extract raster file download string with 2020 population (if download manually is preferred)\n",
    "            products = Product(iso)\n",
    "            Results = products.description_contains('people per grid-cell 2020')\n",
    "            list1 = []\n",
    "            for p in Results:\n",
    "                prints = '%s/%s\\t%s\\t%s' % (p.idx, p.country_name,p.dataset_name,p.path)\n",
    "                list1.append(prints)\n",
    "            print('wpgpDownload download -i',iso,'--id',list1[0].split(\"\\t\")[0].split('/')[0])\n",
    "    \n",
    "    return(iso_countries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c68fbf7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def countries_grids(iso_countries, download_dir = ' '):\n",
    "    start_time = time.time()\n",
    "    blocks = []\n",
    "    for iso in iso_countries:\n",
    "        # Check if raster files already exist on the system path or a manually specified path\n",
    "        path1 = os.getcwd() +'\\\\'+ iso.lower() + '_ppp_2020.tif'\n",
    "        path2 = download_dir +'\\\\'+ iso.lower() + '_ppp_2020.tif'\n",
    "        # First check the manual path\n",
    "        if os.path.exists(path2): \n",
    "            block = gr.from_file(path2)\n",
    "            blocks.append(block)\n",
    "        else:\n",
    "            # Then the system path\n",
    "            if os.path.exists(path1): \n",
    "                block = gr.from_file(path1)\n",
    "                blocks.append(block)\n",
    "            else:\n",
    "                # Otherwise run a suprocess (spr.run) command to download via the terminal in notebook.\n",
    "                runstr = 'wpgpDownload download -i '+ iso+ ' -f people --datasets'\n",
    "                p1 = spr.run('wpgpDownload download -i '+ iso+ ' -f people --datasets', \n",
    "                                    shell = True, \n",
    "                                    capture_output = True)\n",
    "                # decode the output to a list of available datasets from WorldPoP\n",
    "                datasets = p1.stdout.decode().rsplit('\\n')\n",
    "\n",
    "                # The first population raster grid (id-sorted) is the general one, without specifying to demographic groups\n",
    "                for i in enumerate(datasets):\n",
    "                    if '2020' in i[1]:\n",
    "                        ds = datasets[i[0]].rsplit('\\t')[0]\n",
    "                        print(ds)\n",
    "                        # if we found the file, we can stop the loop (we don't need the demograhically specified files)\n",
    "                        break\n",
    "                # Construct the download string\n",
    "                dwnl = 'wpgpDownload download -i '+iso+' --id '+str(ds)\n",
    "                # Get the specified file (terminal)\n",
    "                spr.run(dwnl, shell = True)\n",
    "                # Extract the file\n",
    "                block = gr.from_file(path1)\n",
    "                blocks.append(block)\n",
    "        print(iso,'downloaded', round((time.time() - start_time)/60,2),'mns')\n",
    "    return(blocks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a44b205a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 2 population grids extraction\n",
    "def city_grids_format(bounds, iso_countries, country_grids, cities, grid_size = 100):\n",
    "    start_time = time.time()\n",
    "    grids = []\n",
    "    print(str(grid_size) + 'm resolution grids extraction')\n",
    "    for i in range(len(cities)):\n",
    "        C = pd.DataFrame(Countries)\n",
    "        iso = C[bounds['display_name'][i].rsplit(',')[-1][1:] == C['name']].iloc[0,1]\n",
    "        contains = [j for j, x in enumerate(iso_countries) if x == iso][0]\n",
    "\n",
    "        # Clip the city from the country\n",
    "        clipped = country_grids[contains].clip(bounds['geometry'][i])\n",
    "        clipped = clipped[0].to_geopandas()\n",
    "\n",
    "        # Get dissolvement_key for dissolvement. \n",
    "        clipped['row3'] = np.floor(clipped['row']/(grid_size/100)).astype(int)\n",
    "        clipped['col3'] = np.floor(clipped['col']/(grid_size/100)).astype(int)\n",
    "        clipped['dissolve_key'] = clipped['row3'].astype(str) +'-'+ clipped['col3'].astype(str)\n",
    "\n",
    "        # Dissolve into block by block grids\n",
    "        popgrid = clipped[['dissolve_key','geometry','row3','col3']].dissolve('dissolve_key')\n",
    "\n",
    "        # Get those grids populations and area. Only blocks with population and full blocks\n",
    "        popgrid['population'] = round(clipped.groupby('dissolve_key')['value'].sum()).astype(int)\n",
    "        popgrid['area_m'] = round(gpd.GeoSeries(popgrid['geometry'], crs = 4326).to_crs(3043).area).astype(int)\n",
    "        popgrid = popgrid[popgrid['population'] > 0]\n",
    "        popgrid = popgrid[popgrid['area_m'] / popgrid['area_m'].max() > 0.95]\n",
    "\n",
    "        # Get centroids and coords\n",
    "        popgrid['centroid'] = popgrid['geometry'].centroid\n",
    "        popgrid['centroid_m'] = gpd.GeoSeries(popgrid['centroid'], crs = 4326).to_crs(3043)\n",
    "        popgrid['grid_lon'] = popgrid['centroid_m'].x\n",
    "        popgrid['grid_lat'] = popgrid['centroid_m'].y\n",
    "        popgrid = popgrid.reset_index()\n",
    "\n",
    "        minx = popgrid.bounds['minx']\n",
    "        maxx = popgrid.bounds['maxx']\n",
    "        miny = popgrid.bounds['miny']\n",
    "        maxy = popgrid.bounds['maxy']\n",
    "\n",
    "        # Some geometries result in a multipolygon when dissolving (like i.e. 0.05 meters) which is in my mind an coords error\n",
    "        # I therefore create one polygon\n",
    "        Poly = []\n",
    "        for k in range(len(popgrid)):\n",
    "            Poly.append(Polygon([(minx[k],maxy[k]),(maxx[k],maxy[k]),(maxx[k],miny[k]),(minx[k],miny[k])]))\n",
    "        popgrid['geometry'] = Poly\n",
    "\n",
    "        grids.append(popgrid)\n",
    "\n",
    "        print(cities[i].rsplit(',')[0], round((time.time() - start_time)/60,2),'mns')\n",
    "    return(grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5bc1aa68",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 3 Road networks\n",
    "def road_networks (cities, thresholds, undirected = False):\n",
    "    print('get road networks from OSM')\n",
    "    start_time = time.time()\n",
    "    graphs = list()\n",
    "    road_nodes = list()\n",
    "    road_edges = list()\n",
    "    road_conn = list()\n",
    "\n",
    "    for i in cities:\n",
    "        # Get graph, road nodes and edges\n",
    "        graph = ox.graph_from_place(i, network_type = \"all\", buffer_dist = (np.max(thresholds)+1000))\n",
    "        #graphs.append(graph)\n",
    "\n",
    "        road_node, road_edge = ox.graph_to_gdfs(graph)\n",
    "\n",
    "        # Road nodes format\n",
    "        road_node = road_node.to_crs(4326)\n",
    "        road_node['geometry_m'] = gpd.GeoSeries(road_node['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_node['osmid_var'] = road_node.index\n",
    "        road_node = gpd.GeoDataFrame(road_node, geometry = 'geometry', crs = 4326)\n",
    "\n",
    "        # format road edges\n",
    "        road_edge = road_edge.to_crs(4326)\n",
    "        road_edge['geometry_m'] = gpd.GeoSeries(road_edge['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_edge = road_edge.reset_index()\n",
    "        road_edge.rename(columns={'u':'from', 'v':'to', 'key':'keys'}, inplace=True)\n",
    "        road_edge['key'] = road_edge['from'].astype(str) + '-' + road_edge['to'].astype(str)\n",
    "        \n",
    "        if undirected == True:\n",
    "            # Apply one-directional to both for walking\n",
    "            both = road_edge[road_edge['oneway'] == False]\n",
    "            one = road_edge[road_edge['oneway'] == True]\n",
    "            rev = pd.DataFrame()\n",
    "            rev[['from','to']] = one[['to','from']]\n",
    "            rev = pd.concat([rev,one.iloc[:,2:]],axis = 1)\n",
    "            edge_bidir = pd.concat([both, one, rev])\n",
    "            edge_bidir = edge_bidir.reset_index()\n",
    "            edge_bidir['oneway'] = False\n",
    "        else:\n",
    "            edge_bidir = road_edge\n",
    "\n",
    "        # Exclude highways and ramps on edges    \n",
    "        edge_filter = edge_bidir[(edge_bidir['highway'].str.contains('motorway') | \n",
    "              (edge_bidir['highway'].str.contains('trunk') & \n",
    "               edge_bidir['maxspeed'].astype(str).str.contains(\n",
    "                   '40 mph|45 mph|50 mph|55 mph|60 mph|65|70|75|80|85|90|95|100|110|120|130|140'))) == False]\n",
    "        road_edges.append(edge_filter)\n",
    "\n",
    "        # Exclude isolated nodes\n",
    "        fltrnodes = pd.Series(list(edge_filter['from']) + list(edge_filter['to'])).unique()\n",
    "        newnodes = road_node[road_node['osmid_var'].isin(fltrnodes)]\n",
    "        road_nodes.append(newnodes)\n",
    "\n",
    "        # Get only necessary road connections columns for network performance\n",
    "        road_con = edge_filter[['osmid','key','length','geometry']]\n",
    "        road_con = road_con.set_index('key')\n",
    "\n",
    "        road_conn.append(road_con)\n",
    "\n",
    "        # formatting to graph again.\n",
    "        newnodes = newnodes.loc[:, ~newnodes.columns.isin(['geometry_m', 'osmid_var'])]\n",
    "        edge_filter = edge_filter.set_index(['from','to','keys'])\n",
    "        edge_filter = edge_filter.loc[:, ~edge_filter.columns.isin(['geometry_m', 'key'])]\n",
    "\n",
    "        graph2 = ox.graph_from_gdfs(newnodes, edge_filter)\n",
    "\n",
    "        graphs.append(graph2)\n",
    "        print(i.rsplit(',')[0], 'done', round((time.time() - start_time) / 60,2),'mns')\n",
    "    return({'graphs':graphs,'nodes':road_nodes,'edges':road_conn,'edges long':road_edges})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3d3ceef5",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 4 city greenspace\n",
    "def urban_greenspace (cities, thresholds, one_UGS_buf = 25, min_UGS_size = 400):\n",
    "    print('get urban greenspaces from OSM')\n",
    "    parks_in_range = list()\n",
    "    for i in cities:\n",
    "        gdf = ox.geometries_from_place(i, tags={'leisure':'park'}, buffer_dist = np.max(thresholds))\n",
    "        gdf = gdf[(gdf.geom_type == 'Polygon') | (gdf.geom_type == 'MultiPolygon')]\n",
    "        greenspace = gdf.reset_index()    \n",
    "        warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "        green_buffer = gpd.GeoDataFrame(geometry = greenspace.to_crs(3043).buffer(one_UGS_buf).to_crs(4326))\n",
    "        greenspace['geometry_w_buffer'] = green_buffer\n",
    "        greenspace['geometry_w_buffer'] = gpd.GeoSeries(greenspace['geometry_w_buffer'], crs = 4326)\n",
    "        greenspace['geom buffer diff'] = greenspace['geometry_w_buffer'].difference(greenspace['geometry'])\n",
    "\n",
    "        # This function group components in itself that overlap (with the buffer set of 25 metres)\n",
    "        # https://stackoverflow.com/questions/68036051/geopandas-self-intersection-grouping\n",
    "        W = libpysal.weights.fuzzy_contiguity(greenspace['geometry_w_buffer'])\n",
    "        greenspace['components'] = W.component_labels\n",
    "        parks = greenspace.dissolve('components')\n",
    "\n",
    "        # Exclude parks below 0.04 ha.\n",
    "        parks = parks[parks.to_crs(3043).area > min_UGS_size]\n",
    "        print(i, 'done')\n",
    "        parks = parks.reset_index()\n",
    "        parks['geometry_m'] = parks['geometry'].to_crs(3043)\n",
    "        parks['park_area'] = parks['geometry_m'].area\n",
    "        parks_in_range.append(parks)\n",
    "    return(parks_in_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cdc51e59",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5 park entry points\n",
    "def UGS_fake_entry(UGS, road_nodes, cities, UGS_entry_buf = 25, walk_radius = 500, entry_point_merge = 0):\n",
    "    print('get fake UGS entry points')\n",
    "    start_time = time.time()\n",
    "    ParkRoads = list()\n",
    "    for j in range(len(cities)):\n",
    "        ParkRoad = pd.DataFrame()\n",
    "        mat = list()\n",
    "        # For all\n",
    "        for i in range(len(UGS[j])):\n",
    "            dist = road_nodes[j]['geometry'].to_crs(3043).distance(UGS[j]['geometry'].to_crs(\n",
    "                3043)[i])\n",
    "            buf_nodes = road_nodes[j][(dist < UGS_entry_buf) & (dist > 0)]\n",
    "            mat.append(list(np.repeat(i, len(buf_nodes))))\n",
    "            ParkRoad = pd.concat([ParkRoad, buf_nodes])\n",
    "            if i % 100 == 0: print(cities[j].rsplit(',')[0], round(i/len(UGS[j])*100,1),'% done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "        # Park no list conversion\n",
    "        mat_u = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat) for i in b]\n",
    "\n",
    "        # Format\n",
    "        ParkRoad['Park_No'] = mat_u\n",
    "        ParkRoad = ParkRoad.reset_index()\n",
    "        ParkRoad['park_lon'] = ParkRoad['geometry_m'].x\n",
    "        ParkRoad['park_lat'] = ParkRoad['geometry_m'].y\n",
    "        \n",
    "        # Get the road nodes intersecting with the parks' buffer\n",
    "        ParkRoad = pd.merge(ParkRoad, UGS[j][['geometry','park_area']], left_on = 'Park_No', right_index = True)\n",
    "\n",
    "        # Get the walkable park size\n",
    "        ParkRoad['park_size_walkable'] = ParkRoad['geometry_m'].buffer(walk_radius).to_crs(4326).intersection(ParkRoad['geometry_y'])\n",
    "        ParkRoad['walk_area'] = ParkRoad['park_size_walkable'].to_crs(3043).area\n",
    "        #ParkRoad['park_area'] = ParkRoad['geometry_y'].to_crs(3043).area\n",
    "        ParkRoad['share_walked'] = ParkRoad['walk_area'] / ParkRoad['park_area']\n",
    "                \n",
    "        # Merge fake UGS entry points if within X meters of each other for better system performance\n",
    "        # Standard no merging\n",
    "        ParkRoad = simplify_UGS_entry(ParkRoad, entry_point_merge = 0)\n",
    "                \n",
    "        ParkRoads.append(ParkRoad)\n",
    "\n",
    "        print(cities[j].rsplit(',')[0],'100 % done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "    return(ParkRoads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "af76feed",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5.5 (not in use, buffer is 0, thus retains all the park entry points as is)\n",
    "def simplify_UGS_entry(fake_UGS_entry, entry_point_merge = 0):\n",
    "    # Get buffer of nodes close to each other.\n",
    "    # Get the buffer\n",
    "    ParkComb = fake_UGS_entry\n",
    "    ParkComb['geometry_m_buffer'] = ParkComb['geometry_m'].buffer(entry_point_merge)\n",
    "\n",
    "    # Get and merge components\n",
    "    M = libpysal.weights.fuzzy_contiguity(ParkComb['geometry_m_buffer'])\n",
    "    ParkComb['components'] = M.component_labels\n",
    "\n",
    "    # Take centroid of merged components\n",
    "    centr = gpd.GeoDataFrame(ParkComb, geometry = 'geometry_x', crs = 4326).dissolve('components')['geometry_x'].centroid\n",
    "    centr = gpd.GeoDataFrame(centr)\n",
    "    centr.columns = ['comp_centroid']\n",
    "\n",
    "    # Get node closest to the centroid of all merged nodes, which accesses the road network.\n",
    "    ParkComb = pd.merge(ParkComb, centr, left_on = 'components', right_index = True)\n",
    "    ParkComb['centr_dist'] = ParkComb['geometry_x'].distance(ParkComb['comp_centroid'])\n",
    "    ParkComb = ParkComb.iloc[ParkComb.groupby('components')['centr_dist'].idxmin()]\n",
    "    return(ParkComb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "19711d72",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 6 grid-parkentry combinations within euclidean threshold distance\n",
    "def suitible_combinations(UGS_entry, pop_grids, road_nodes, thresholds, cities, chunk_size = 10000000):\n",
    "    print('get potential (Euclidean) suitible combinations')\n",
    "    start_time = time.time()\n",
    "    RoadComb = list()\n",
    "    for l in range(len(cities)):\n",
    "        #blockA = block_combinations\n",
    "        print(cities[l])\n",
    "        len1 = len(pop_grids[l])\n",
    "        len2 = len(UGS_entry[l])\n",
    "\n",
    "        # Reduce the size of combinations per iteration\n",
    "        len4 = 1\n",
    "        len5 = len1 * len2\n",
    "        blockC = len5\n",
    "        while blockC > chunk_size:\n",
    "            blockC = len5 / len4\n",
    "            #print(blockC, len4)\n",
    "            len4 = len4+1\n",
    "\n",
    "        # Amount of grids taken per iteration block\n",
    "        block = round(len1 / len4)\n",
    "\n",
    "        output = pd.DataFrame()\n",
    "        # Checking all the combinations at once is too performance intensive, it is broken down per 1000 (or what you want)\n",
    "        for i in range(len4):\n",
    "            # Check all grid-park combinations per block\n",
    "            l1, l2 = range(i*block,(i+1)*block), range(0,len2)\n",
    "            listed = pd.DataFrame(list(product(l1, l2)))\n",
    "\n",
    "            # Merge grid and park information\n",
    "            grid_merged = pd.merge(listed, \n",
    "                                   pop_grids[l][['grid_lon','grid_lat','centroid','centroid_m']],\n",
    "                                   left_on = 0, right_index = True)\n",
    "            node_merged = pd.merge(grid_merged, \n",
    "                                   UGS_entry[l][['Park_No','osmid','geometry_x','geometry_y','geometry_m','park_lon','park_lat',\n",
    "                                       'share_walked','park_area','walk_area']], \n",
    "                                   left_on = 1, right_index = True)\n",
    "\n",
    "            # Preset index for merging\n",
    "            node_merged['key'] = range(0,len(node_merged))\n",
    "            node_merged = node_merged.set_index('key')\n",
    "            node_merged = node_merged.loc[:, ~node_merged.columns.isin(['index'])]\n",
    "\n",
    "            # Create lists for better computational performance\n",
    "            glon = list(node_merged['grid_lon'])\n",
    "            glat = list(node_merged['grid_lat'])\n",
    "            plon = list(node_merged['park_lon'])\n",
    "            plat = list(node_merged['park_lat'])\n",
    "\n",
    "            # Get the euclidean distances\n",
    "            mat = list()\n",
    "            for j in range(len(node_merged)):\n",
    "                mat.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2))\n",
    "\n",
    "            # Check if distances are within 1000m and join remaining info and concat in master df per 1000.\n",
    "            mat_df = pd.DataFrame(mat)[(np.array(mat) <= np.max(thresholds))]\n",
    "\n",
    "            # join the other gravity euclidean scores and other information\n",
    "            mat_df.columns = ['Euclidean']    \n",
    "            mat_df = mat_df.join(node_merged)\n",
    "\n",
    "            output = pd.concat([output, mat_df])\n",
    "\n",
    "            print('in chunk',(i+1),'/',len4,len(mat_df),'suitible comb.')\n",
    "        # Renaming columns\n",
    "        print('total combinations within distance',len(output))\n",
    "\n",
    "        output.columns = ['Euclidean','Grid_No','Park_entry_No','grid_lon','grid_lat','Grid_coords_centroid','Grid_m_centroid',\n",
    "                      'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid','park_lon',\n",
    "                      'park_lat','parkshare_walked','park_area','walk_area_m2']\n",
    "\n",
    "        output = output[['Euclidean','Grid_No','Park_entry_No','Grid_coords_centroid','Grid_m_centroid','walk_area_m2',\n",
    "                     'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid','park_area']]\n",
    "\n",
    "        # Reinstate geographic elements\n",
    "        output = gpd.GeoDataFrame(output, geometry = 'Grid_coords_centroid', crs = 4326)\n",
    "        output['Grid_m_centroid'] = gpd.GeoSeries(output['Grid_m_centroid'], crs = 3043)\n",
    "        output['Parkroad_coords_centroid'] = gpd.GeoSeries(output['Parkroad_coords_centroid'], crs = 4326)\n",
    "        output['Parkroad_m_centroid'] = gpd.GeoSeries(output['Parkroad_m_centroid'], crs = 3043)\n",
    "\n",
    "        # Get the nearest entrance point for the grid centroids\n",
    "        output = gridroad_entry(output, road_nodes[l])\n",
    "\n",
    "        print('100 % gridentry done', round((time.time() - start_time) / 60,2),' mns')\n",
    "        RoadComb.append(output)\n",
    "    return (RoadComb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e9d2c955",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gridroad_entry (suitible_comb, road_nodes):    \n",
    "    start_time = time.time()\n",
    "    mat5 = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        try:\n",
    "            nearest = int(road_nodes['geometry'].sindex.nearest(suitible_comb['Grid_coords_centroid'].iloc[i])[1])\n",
    "            mat5.append(road_nodes['osmid_var'].iloc[nearest])\n",
    "        except: \n",
    "            # sometimes two nodes are the exact same distance, then the first in the list is taken.\n",
    "            nearest = int(road_nodes['geometry'].sindex.nearest(suitible_comb['Grid_coords_centroid'].iloc[i])[1][0])\n",
    "            mat5.append(road_nodes['osmid_var'].iloc[nearest])\n",
    "        if i % 250000 == 0: print(round(i/len(suitible_comb)*100,1),'% gridentry done', round((time.time() - start_time) / 60,2),' mns')\n",
    "    # format resulting dataframe\n",
    "    suitible_comb['grid_osm'] = mat5\n",
    "    suitible_comb = pd.merge(suitible_comb, road_nodes['geometry'], left_on = 'grid_osm', right_index = True)\n",
    "    suitible_comb['geometry_m'] = gpd.GeoSeries(suitible_comb['geometry'], crs = 4326).to_crs(3043)\n",
    "    suitible_comb = suitible_comb.reset_index()\n",
    "    return(suitible_comb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "a8f10468",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check grids in or out of UGS\n",
    "def grids_in_UGS (suitible_comb, UGS, pop_grid): \n",
    "    start_time = time.time()\n",
    "    RoadInOut = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        UGS_geoms = UGS[i]['geometry']\n",
    "        grid = pop_grid[i]['centroid']\n",
    "        lst = list()\n",
    "        print('Check grids within UGS')\n",
    "        for l in enumerate(UGS_geoms):\n",
    "            lst.append(grid.intersection(l[1]).is_empty == False)\n",
    "            if l[0] % 100 == 0: print(l[0], round((time.time() - start_time) / 60,2),' mns')\n",
    "\n",
    "        dfGrUGS = pd.DataFrame(pd.DataFrame(np.array(lst)).unstack())\n",
    "        dfGrUGS.columns = ['in_out_UGS']\n",
    "        merged = pd.merge(suitible_comb[i], dfGrUGS, left_on = ['Grid_No','Park_No'], right_index = True, how = 'left')\n",
    "        RoadInOut.append(merged)\n",
    "    return(RoadInOut)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "b44555aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 7 calculate route networks of all grid-parkentry combinations within euclidean threshold distance\n",
    "def route_finding (graphs, combinations, road_nodes, road_edges, cities, block_size = 250000, nn_iter = 10):\n",
    "\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    start_time = time.time()\n",
    "\n",
    "    Routes = list()\n",
    "    Routes_detail = list()\n",
    "    for j in range(len(cities)):\n",
    "        Graph = graphs[j]\n",
    "        suit_raw = combinations[j] # iloc to test the iteration speed.\n",
    "        nodes = road_nodes[j]\n",
    "\n",
    "        In_UGS = suit_raw[suit_raw['in_out_UGS'] == True] # Check if a grid centroid is in an UGS\n",
    "        suitible = suit_raw[suit_raw['in_out_UGS'] == False].reset_index(drop = True) # recreate a subsequential index\n",
    "                                                                                      # for the other grids outside UGS\n",
    "        block = block_size # Execute with chunks for performance improvement.\n",
    "\n",
    "        Route_parts = pd.DataFrame()\n",
    "        Route_dparts = pd.DataFrame()\n",
    "        len2 = int(np.ceil(len(suitible)/block))\n",
    "        # Divide in chunks of block for computational load\n",
    "        for k in range(len2):    \n",
    "            suitible_chunk = suitible.iloc[k*block:k*block+block] # Select chunk\n",
    "\n",
    "            parknode = list(suitible_chunk['Parkroad_osmid'])\n",
    "            gridnode = list(suitible_chunk['grid_osm'])\n",
    "\n",
    "            s_mat = list([]) # origin (normally grid) osmid\n",
    "            s_mat1 = list([]) # destination (normally UGS) osmid\n",
    "            s_mat2 = list([]) # route id\n",
    "            s_mat3 = list([]) # step id\n",
    "            s_mat4 = list([]) # way calculated\n",
    "            s_mat5 = list([]) # way calculated id\n",
    "            mat_nn = [] # found nearest nodes by block\n",
    "            len1 = len(suitible_chunk)\n",
    "\n",
    "            print(cities[j].rsplit(',')[0], k+1,'/',len2,'range',k*block,'-',k*block+np.where(k*block+block >= len1,len1,block))\n",
    "            for i in range(len(suitible_chunk)):\n",
    "                try: \n",
    "                    # from grid to UGS.\n",
    "                    shortest = nx.shortest_path(Graph, gridnode[i], parknode[i], 'travel_dist', method = 'dijkstra')\n",
    "                    s_mat.append(shortest)\n",
    "                    shortest_to = list(shortest[1:len(shortest)])\n",
    "                    shortest_to.append(-1)\n",
    "                    s_mat1.append(shortest_to)\n",
    "                    s_mat2.append(list(np.repeat(i+block*k, len(shortest))))\n",
    "                    s_mat3.append(list(np.arange(0, len(shortest))))\n",
    "                    s_mat4.append('normal way')\n",
    "                    s_mat5.append(1)\n",
    "                except:\n",
    "                    try:\n",
    "                        # Check the reverse\n",
    "                        shortest = nx.shortest_path(Graph, parknode[i], gridnode[i], 'travel_dist', method = 'dijkstra')\n",
    "                        s_mat.append(shortest)\n",
    "                        shortest_to = list(shortest[1:len(shortest)])\n",
    "                        shortest_to.append(-1)\n",
    "                        s_mat1.append(shortest_to)\n",
    "                        s_mat2.append(list(np.repeat(i+block*k, len(shortest))))\n",
    "                        s_mat3.append(list(np.arange(0, len(shortest))))\n",
    "                        s_mat4.append('reverse way')\n",
    "                        s_mat5.append(0)\n",
    "                    except:\n",
    "                        # Otherwise find nearest nodes (grid and UGS) and try to find routes between them\n",
    "                        nn_route_finding(Graph, suitible_chunk, nodes, s_mat, s_mat1, s_mat2, s_mat3,\n",
    "                                             s_mat4, s_mat5, mat_nn, i, block, k, nn_iter)\n",
    "                        \n",
    "                if i % 10000 == 0: print(round((i+block*k)/len(suitible)*100,2),'% done',\n",
    "                                         round((time.time() - start_time) / 60,2),'mns')\n",
    "            print('for', len(mat_nn),'routes nearest nodes found')\n",
    "\n",
    "            print(round((i+block*k)/len(suitible)*100,2),'% pathfinding done', round((time.time() - start_time) / 60,2),'mns')\n",
    "\n",
    "            # Formats route information by route and step (detailed)\n",
    "            routes = route_formatting(s_mat, s_mat1, s_mat2, s_mat3, road_edges[j]) # Formats lists to routes detail.\n",
    "            print('formatting done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            # Summarizes information by route\n",
    "            routes2 = route_summarization(routes, suitible_chunk, road_nodes[j], s_mat4, s_mat5) # formats routes to summary\n",
    "            print('dissolving done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            Route_parts = pd.concat([Route_parts, routes2])\n",
    "            Route_dparts = pd.concat([Route_dparts, routes])\n",
    "\n",
    "        # Format grids in UGS to enable smooth df concat\n",
    "        In_UGS = In_UGS.set_geometry(In_UGS['Grid_coords_centroid'])\n",
    "        In_UGS = In_UGS[['geometry','Grid_No','grid_osm','Park_No','Park_entry_No','Parkroad_osmid',\n",
    "                                   'Grid_m_centroid','walk_area_m2',\n",
    "                                   'Euclidean','geometry_m']]\n",
    "\n",
    "        In_UGS['realG_osmid'] = suit_raw['Parkroad_osmid']\n",
    "        In_UGS['realP_osmid'] = suit_raw['grid_osm']\n",
    "        In_UGS['way_calc'] = 'grid in UGS'\n",
    "\n",
    "        Route_parts = pd.concat([Route_parts,In_UGS])\n",
    "        Route_parts = Route_parts.reset_index(drop = True)\n",
    "\n",
    "        Route_parts['gridpark_no'] = Route_parts['Grid_No'].astype(str) +'-'+ Route_parts['Park_No'].astype(str)\n",
    "\n",
    "        # All fill value 0 because no routes are calculated for grid centroids in UGSs\n",
    "        to_fill = ['way-id','route_cost','steps','real_G-entry','Tcost']                                   \n",
    "        Route_parts[to_fill] = Route_parts[to_fill].fillna(0)  \n",
    "            \n",
    "        Routes.append(Route_parts)\n",
    "        Routes_detail.append(Route_dparts)\n",
    "    return(Routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "a680b152",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_route_finding(graph, suitible_chunk, nodes, mat_from, mat_to, mat_route, mat_step,\n",
    "                                             mat_way, mat_wbin, mat_nn, i, block, k, nn_iter):\n",
    "                        \n",
    "    # Order in route for nearest node:\n",
    "    # 1. gridnode to nearest to the original failed parknode\n",
    "    # 2. The reverse of 1.\n",
    "    # 3. nearest gridnode to the failed one and route to park\n",
    "    # 4. The reverse of 3.\n",
    "                        \n",
    "    gridosm = suitible_chunk['grid_osm'] # grid osmid\n",
    "    UGSosm = suitible_chunk['Parkroad_osmid'] # UGS osmid\n",
    "    nodeosm = nodes['osmid_var'] # road node osmid\n",
    "    nodegeom = nodes['geometry'] # road node geometry\n",
    "                        \n",
    "    len3 = 0\n",
    "    alt_route = list([])\n",
    "    while len3 < nn_iter and len(alt_route) < 1: # If a route is found (alt_route == 1) or until max iterations\n",
    "\n",
    "        len3 = len3 +1\n",
    "                            \n",
    "        nn = nn_finding(gridosm, UGSosm, nodeosm, nodegeom, nodes, i, len3) # finds nearest node.\n",
    "\n",
    "        nn_routing (graph, nn['currUGS'], nn['nearUGS'], nn['currgrid'], nn['neargrid'], \n",
    "                                        mat_way, mat_wbin, len3, alt_route) # executes route finding in try order.\n",
    "    if len(alt_route) == 0: \n",
    "        alt = alt_route \n",
    "    else: \n",
    "        alt = alt_route[0]\n",
    "    len4 = len(alt)\n",
    "    if len4 > 0: # If a route is found\n",
    "        mat_nn.append(i+block*k)\n",
    "        mat_from.append(alt)\n",
    "        shortest_to = list(alt[1:len(alt)])\n",
    "        shortest_to.append(-1)\n",
    "        mat_to.append(shortest_to)\n",
    "        mat_route.append(list(np.repeat(i+block*k,len4)))\n",
    "        mat_step.append(list(np.arange(0, len4)))\n",
    "    else: # If a route is not found\n",
    "        mat_from.append(-1)\n",
    "        mat_to.append(-1)\n",
    "        mat_route.append(i+block*k)\n",
    "        mat_step.append(-1)\n",
    "        mat_way.append('no way')\n",
    "        mat_wbin.append(2)\n",
    "        print(i+block*k,'No route',nn_iter)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "695e9681",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_finding (gridosm, UGSosm, nodeosm, nodegeom, nodes, i, nn_i): \n",
    "    # Grid nearest\n",
    "    g_geom = nodegeom[nodeosm == int(gridosm[i:i+1])] # Get geom of current node UGS\n",
    "    g_nearest = pd.DataFrame((abs(float(g_geom.x) - nodegeom.x)**2 # Check distance UGS\n",
    "    +abs(float(g_geom.y) - nodegeom.y)**2)**(1/2)\n",
    "                            ).join(nodeosm).sort_values(0) # sort by distance ascending UGS\n",
    "\n",
    "    g_grid = g_nearest.iloc[nn_i,1] # get the nearest node according to the nn_iter UGS entry\n",
    "    g_park = list(UGSosm)[i] # current node\n",
    "        \n",
    "    p_geom = nodegeom[nodeosm == int(UGSosm[i:i+1])] # get the geom of the current node grid\n",
    "    p_nearest = pd.DataFrame((abs(float(p_geom.x) - nodegeom.x)**2 # Check distance grid\n",
    "    +abs(float(p_geom.y) - nodegeom.y)**2)**(1/2)\n",
    "                            ).join(nodeosm).sort_values(0) # sort by distance ascending grid\n",
    "\n",
    "    p_grid = list(gridosm)[i] # current node\n",
    "    p_park = p_nearest.iloc[nn_i,1] # get the nearest node to the nn_iter grid\n",
    "    return({'currUGS':p_grid, 'nearUGS':p_park,'currgrid':g_park, 'neargrid':g_grid})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "65f7c352",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_routing (graph, curr_UGS, near_UGS, curr_grid, near_grid, mat_way, mat_wbin, nn_i, found_route):\n",
    "    try:\n",
    "        found_route.append(nx.shortest_path(graph, curr_UGS, near_UGS, \n",
    "                                          'travel_dist', method = 'dijkstra'))\n",
    "        mat_way.append(str(nn_i)+'grid > n-park') # grid to nearest unseen UGS node\n",
    "        mat_wbin.append(1)\n",
    "    except:\n",
    "        try:\n",
    "            found_route.append(nx.shortest_path(graph, near_UGS, curr_UGS, \n",
    "                                              'travel_dist', method = 'dijkstra'))\n",
    "            mat_way.append(str(nn_i)+'n-park > grid') # nearest unseen UGS node to grid\n",
    "            mat_wbin.append(0)\n",
    "        except:\n",
    "            try:\n",
    "                found_route.append(nx.shortest_path(graph, curr_grid, near_grid, \n",
    "                                                  'travel_dist', method = 'dijkstra'))\n",
    "                mat_way.append(str(nn_i)+'n-grid > park') # nearest grid node to UGS\n",
    "                mat_wbin.append(1)\n",
    "            except:\n",
    "                try:\n",
    "                    found_route.append(nx.shortest_path(graph, near_grid, curr_grid, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                    mat_way.append(str(nn_i)+'park > n-grid') # UGS to nearest grid node\n",
    "                    mat_wbin.append(0)\n",
    "                except:\n",
    "                    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "9dd0bea7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def route_formatting(mat_from, mat_to, mat_route, mat_step, road_edges):\n",
    "    # Unpack lists\n",
    "    s_mat_u = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_from) for i in b]\n",
    "    s_mat_u1 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_to) for i in b]\n",
    "    s_mat_u2 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_route) for i in b]\n",
    "    s_mat_u3 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_step) for i in b]\n",
    "\n",
    "    # Format df\n",
    "    routes = pd.DataFrame([s_mat_u,s_mat_u1,s_mat_u2,s_mat_u3]).transpose()\n",
    "    routes.columns = ['from','to','route','step']\n",
    "    mat_key = list([])\n",
    "    for n in range(len(routes)): # get key of origin and destination\n",
    "        mat_key.append(str(int(s_mat_u[n])) + '-' + str(int(s_mat_u1[n])))\n",
    "    routes['key'] = mat_key\n",
    "    routes = routes.set_index('key')\n",
    "\n",
    "    # Add route information\n",
    "    routes = routes.join(road_edges, how = 'left') # to add road node information\n",
    "    routes = gpd.GeoDataFrame(routes, geometry = 'geometry', crs = 4326)\n",
    "    routes = routes.sort_values(by = ['route','step'])\n",
    "    return(routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "84241ed6",
   "metadata": {},
   "outputs": [],
   "source": [
    "def route_summarization(routes, suitible_comb, road_nodes, mat_way, mat_wbin):\n",
    "    # dissolve route\n",
    "    routes2 = routes[['route','geometry']].dissolve('route')\n",
    "\n",
    "    # get used grid- and parkosm. Differs at NN-route.\n",
    "    route_reset = routes.reset_index()\n",
    "    origin = route_reset['from'].iloc[list(route_reset.groupby('route')['step'].idxmin()),]\n",
    "    origin = origin.reset_index().iloc[:,-1]\n",
    "    dest = route_reset['from'].iloc[list(route_reset.groupby('route')['step'].idxmax()),]\n",
    "    dest = dest.reset_index().iloc[:,-1]\n",
    "\n",
    "    # grid > park = 1, park > grid = 0, no way = 2, detailed way in way_calc.\n",
    "    routes2['way-id'] = mat_wbin\n",
    "    routes2['realG_osmid'] = np.where(routes2['way-id'] == 1, origin, dest)\n",
    "    routes2['realP_osmid'] = np.where(routes2['way-id'] == 1, dest, origin)\n",
    "    routes2['way_calc'] = mat_way\n",
    "\n",
    "    # get route cost, steps, additional information.\n",
    "    routes2['route_cost'] = routes.groupby('route')['length'].sum()\n",
    "    routes2['steps'] = routes.groupby('route')['step'].max()\n",
    "    routes2['index'] = suitible_comb.index\n",
    "    routes2 = routes2.set_index(['index'])\n",
    "    routes2.index = routes2.index.astype(int)\n",
    "    routes2 = pd.merge(routes2, suitible_comb[['Grid_No','grid_osm','Park_No','Park_entry_No','Parkroad_osmid',\n",
    "                                          'Grid_m_centroid','walk_area_m2','Euclidean']],\n",
    "                                            left_index = True, right_index = True)\n",
    "    routes2 = pd.merge(routes2, road_nodes['geometry_m'], how = 'left', left_on = 'realG_osmid', right_index = True)\n",
    "    # calculate distance of used road-entry for grid-centroid.\n",
    "    routes2['real_G-entry'] = round(gpd.GeoSeries(routes2['Grid_m_centroid'], crs = 3043).distance(routes2['geometry_m']),3)\n",
    "                                    \n",
    "    # Calculcate total route cost for the four gravity variants\n",
    "    routes2['Tcost'] = routes2['route_cost'] + routes2['real_G-entry']\n",
    "    return(routes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "e355b43e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_gridUGS_comb (routes, grids, UGS):\n",
    "    gp_nearest = []\n",
    "    for i in range(len(routes)):\n",
    "        gp_nn = routes[i][routes[i]['Tcost'] <= max(thresholds)]\n",
    "        gp_nn = pd.merge(gp_nn, grids[i]['population'], left_on='Grid_No', right_index = True)\n",
    "        gp_nn = pd.merge(gp_nn, UGS[i]['park_area'], left_on = 'Park_No', right_index = True)\n",
    "        gp_nn = gp_nn.reset_index()\n",
    "\n",
    "        gp_nn = gp_nn.iloc[gp_nn.groupby('gridpark_no')['Tcost'].idxmin()]\n",
    "        gp_nn.index.name = 'idx'\n",
    "        gp_nn = gp_nn.sort_values('idx')\n",
    "        gp_nn = gp_nn.reset_index()\n",
    "        gp_nearest.append(gp_nn)\n",
    "    gp_nearest[0].sort_values('Grid_No')\n",
    "    return(gp_nearest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "ebb4c053",
   "metadata": {},
   "outputs": [],
   "source": [
    "def E2SCFA_scores(min_gridUGS_comb, grids, thresholds, cities):\n",
    "    pd.options.display.float_format = '{:20,.2f}'.format\n",
    "    E2SFCA_cities = []\n",
    "    E2SFCA_summary = pd.DataFrame()\n",
    "    for i in range(len(cities)):\n",
    "        E2SFCA_score = grids[i][['population','geometry']]\n",
    "        for j in range(len(thresholds)):\n",
    "            subset = min_gridUGS_comb[i][min_gridUGS_comb[i]['Tcost'] <= thresholds[j]]\n",
    "\n",
    "            # use gussian distribution: let v= 923325, then the weight for 800m is 0.5\n",
    "            v = -thresholds[j]**2/np.log(0.5)\n",
    "\n",
    "            # add a column of weight: apply the decay function on distance\n",
    "            subset['weight'] = np.exp(-(subset['Tcost']**2/v)).astype(float)\n",
    "            subset['pop_weight'] = subset['weight'] * subset['population']\n",
    "\n",
    "            # get the sum of weighted population each green space has to serve.\n",
    "            s_w_p = pd.DataFrame(subset.groupby('Park_No').sum('pop_weight')['pop_weight'])\n",
    "\n",
    "            # delete other columns, because they are useless after groupby\n",
    "            s_w_p = s_w_p.rename({'pop_weight':'pop_weight_sum'},axis = 1)\n",
    "            middle = pd.merge(subset,s_w_p, how = 'left', on = 'Park_No' )\n",
    "\n",
    "            # calculate the supply-demand ratio for each green space\n",
    "            middle['green_supply'] = middle['park_area']/middle['pop_weight_sum']\n",
    "\n",
    "            # caculate the accessbility score for each green space that each population grid cell could reach\n",
    "            middle['Sc-access'] = middle['weight'] * middle['green_supply']\n",
    "            # add the scores for each population grid cell\n",
    "            pop_score_df = pd.DataFrame(middle.groupby('Grid_No').sum('Sc-access')['Sc-access'])\n",
    "\n",
    "            # calculate the mean distance of all the green space each population grid cell could reach\n",
    "            mean_dist = middle.groupby('Grid_No').mean('Tcost')['Tcost']\n",
    "            pop_score_df['M-dist'] = mean_dist\n",
    "\n",
    "            # calculate the mean area of all the green space each population grid cell could reach\n",
    "            mean_area = middle.groupby('Grid_No').mean('park_area')['park_area']\n",
    "            pop_score_df['M-area'] = mean_area\n",
    "\n",
    "            # calculate the mean supply_demand ratio of all the green space each population grid cell could reach\n",
    "            mean_supply = middle.groupby('Grid_No').mean('green_supply')['green_supply']\n",
    "            pop_score_df['M-supply'] = mean_supply\n",
    "\n",
    "            pop_score = pop_score_df\n",
    "\n",
    "            pop_score_df = pop_score_df.join(grids[i]['population'], how = 'right')\n",
    "            pop_score_df['Sc-norm'] = pop_score_df['Sc-access'] / pop_score_df['population']\n",
    "\n",
    "            pop_score_df = pop_score_df.loc[:, pop_score_df.columns != 'population']\n",
    "            pop_score_df = pop_score_df.add_suffix(' '+str(thresholds[j]))\n",
    "            E2SFCA_score = E2SFCA_score.join(pop_score_df, how = 'left')\n",
    "\n",
    "            print(thresholds[j], cities[i])\n",
    "            \n",
    "        if not os.path.exists('D:Dumps/E2SFCA-OD/Scores/grid_geoms/'):\n",
    "            os.makedirs('D:Dumps/E2SFCA-OD/Scores/grid_geoms/')\n",
    "\n",
    "        E2SFCA_score = E2SFCA_score.fillna(0)\n",
    "        E2SFCA_score.to_file('D:Dumps/E2SFCA-OD/Scores/grid_geoms/'+cities[i]+'.shp') # Detailed scores\n",
    "        pop_sum = pd.Series(E2SFCA_score['population'].sum()).astype(int)\n",
    "        pop_sum.index = ['population']\n",
    "        mean_metrics = E2SFCA_score.loc[:, E2SFCA_score.columns != 'population'].mean()\n",
    "        E2SFCA_sum = pd.concat([pop_sum, mean_metrics])\n",
    "        E2SFCA_summary = pd.concat([E2SFCA_summary, E2SFCA_sum], axis = 1) # summarized results\n",
    "        E2SFCA_cities.append(E2SFCA_score)\n",
    "        \n",
    "        if not os.path.exists('D:/Dumps/E2SFCA-OD/Scores/'):\n",
    "            os.makedirs('D:/Dumps/E2SFCA-OD/Scores/')\n",
    "        \n",
    "        E2SFCA_score.loc[:, E2SFCA_score.columns != 'geometry'].to_csv('D:/Dumps/E2SFCA-OD/Scores/'+cities[i]+'.csv')\n",
    "    E2SFCA_summary.columns = cities\n",
    "    \n",
    "    if not os.path.exists('D:/Dumps/E2SFCA-OD/Scores/'):\n",
    "        os.makedirs('D:/Dumps/E2SFCA-OD/Scores/')\n",
    "    \n",
    "    E2SFCA_summary.to_csv('D:/Dumps/E2SFCA-OD/Scores/all_cities.csv')\n",
    "    E2SFCA_summary\n",
    "    return({'score summary':E2SFCA_summary,'score detail':E2SFCA_cities})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "d90c64ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.2 mns\n"
     ]
    }
   ],
   "source": [
    "print(round((time.time() - start) / 60,2),'mns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "665a9824",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82686f3c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}