E2SFCA GEE-OSM Model-Func 2.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "442b3921",
   "metadata": {},
   "outputs": [],
   "source": [
    "# system packages\n",
    "import sys\n",
    "import time\n",
    "import warnings\n",
    "import os\n",
    "\n",
    "# non-geo numeric packages\n",
    "import numpy as np\n",
    "import math\n",
    "from itertools import product, combinations\n",
    "import pandas as pd\n",
    "\n",
    "# network and OSM packages\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "city_geo = ox.geocoder.geocode_to_gdf\n",
    "\n",
    "# Earth engine packages\n",
    "import ee\n",
    "import geemap\n",
    "\n",
    "# General geo-packages\n",
    "import libpysal\n",
    "import rasterio\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "07277cdb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=ysO-w_I_ZrnW6bb7tIWC7iPqTXsjagMdEamSp1IkN-c&tc=pP62VB4mx4S9VBgdyrwvDQV-1zhTmgVntc1KDp-6GFI&cc=9gIXLbrGLsHogtsg27yNcBKnBdtzqG1GunD9OvEMukM>https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=ysO-w_I_ZrnW6bb7tIWC7iPqTXsjagMdEamSp1IkN-c&tc=pP62VB4mx4S9VBgdyrwvDQV-1zhTmgVntc1KDp-6GFI&cc=9gIXLbrGLsHogtsg27yNcBKnBdtzqG1GunD9OvEMukM</a></p>\n",
       "        <p>The authorization workflow will generate a code, which you should paste in the box below.</p>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter verification code: 4/1AVHEtk54M9lXpJZmq6XFNBSHtB9a5SsmpGNtVe-ZuHFqFzT-Wq5hRHF5p70\n",
      "\n",
      "Successfully saved authorization token.\n"
     ]
    }
   ],
   "source": [
    "# Authenticate and Initialize Google Earth Engine\n",
    "ee.Authenticate()\n",
    "ee.Initialize()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "47f4b61b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 0 cities and thresholds\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m\n",
    "\n",
    "# Extract iso-3166 country codes\n",
    "iso = pd.read_excel('iso_countries.xlsx')\n",
    "\n",
    "# Extract cities list\n",
    "cities = pd.read_excel('cities.xlsx') # all cities\n",
    "\n",
    "# 'cities_adj' serves by default as city-input for functions\n",
    "# cities_adj = cities\n",
    "# cities_adj = cities[cities['Included (Y/N)'] == 'Y']\n",
    "cities_adj = cities[cities['City'].isin(['Addis Ababa','Dhaka','Shijiazhuang','Damascus'])]\n",
    "cities_adj = cities_adj.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fa5e18f7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/f0de31b02bb07953d66a210523e93d01-a18cb8bb1b94fd93363a46bc57c91e8e:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\ETH_Addis Ababa_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/ff61bd25735c10bedb131185a12a5a39-9bb3905e6784cef45ce03ed7960665c8:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\SYR_Damascus_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/e374c5e555e31457cffe53390382762f-cdb1a3867aae40ff788480c6db729411:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\BGD_Dhaka_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/c00d2dec9576c7941a5faf8cf3e6eb04-657492e08299f81a7f4fd40fe2362577:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\CHN_Shijiazhuang_2020.tif\n",
      "Wall time: 13.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 1. Required preprocess for information extraction\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# In essence, we use Google Earth Engine to extract a country's grid raster and clip it with the city's preferred OSM area\n",
    "# Predifine in Excel: the (1) city name as \"City\" and (2) the OSM area that needs to be extracted as \"OSM_area\"\n",
    "# i.e. City = \"Los Angeles\" and OSM_area = \"Los Angeles county, Orange county CA\"\n",
    "files = gee_worldpop_extract(cities_adj, iso, 'C:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "df72d9e2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Addis Ababa 0.4 mns\n",
      "Damascus 0.53 mns\n",
      "Dhaka 0.79 mns\n",
      "Shijiazhuang 1.14 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Addis Ababa done 1.65 mns\n",
      "Damascus done 2.22 mns\n",
      "Dhaka done 3.3 mns\n",
      "Shijiazhuang done 3.67 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Addis Ababa done\n",
      "Damascus done\n",
      "Dhaka done\n",
      "Shijiazhuang done\n",
      "Wall time: 4min 50s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 2. Information extraction\n",
    "\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print(' ')\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_networks(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\n",
    "\n",
    "# Extract urban greenspace (UGS)\n",
    "UGS = urban_greenspace(cities_adj, \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": 38,
   "id": "6a673017",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get fake UGS entry points\n",
      "Addis Ababa 0.0 % done 0.01  mns\n",
      "Addis Ababa 73.5 % done 0.39  mns\n",
      "Addis Ababa 100 % done 0.53  mns\n",
      "Damascus 0.0 % done 0.53  mns\n",
      "Damascus 41.2 % done 0.67  mns\n",
      "Damascus 82.3 % done 0.81  mns\n",
      "Damascus 100 % done 0.88  mns\n",
      "Dhaka 0.0 % done 0.89  mns\n",
      "Dhaka 19.4 % done 1.12  mns\n",
      "Dhaka 38.8 % done 1.33  mns\n",
      "Dhaka 58.1 % done 1.56  mns\n",
      "Dhaka 77.5 % done 1.77  mns\n",
      "Dhaka 96.9 % done 2.01  mns\n",
      "Dhaka 100 % done 2.05  mns\n",
      "Shijiazhuang 0.0 % done 2.05  mns\n",
      "Shijiazhuang 100 % done 2.14  mns\n",
      " \n",
      "get potential (Euclidean) suitible combinations\n",
      "Addis Ababa\n",
      "in chunk 1 / 7 87293 suitible comb.\n",
      "in chunk 2 / 7 3449 suitible comb.\n",
      "in chunk 3 / 7 9022 suitible comb.\n",
      "in chunk 4 / 7 1081 suitible comb.\n",
      "in chunk 5 / 7 6206 suitible comb.\n",
      "in chunk 6 / 7 36604 suitible comb.\n",
      "in chunk 7 / 7 53040 suitible comb.\n",
      "total combinations within distance 196695\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 2.48  mns\n",
      "Damascus\n",
      "in chunk 1 / 4 41477 suitible comb.\n",
      "in chunk 2 / 4 83841 suitible comb.\n",
      "in chunk 3 / 4 160101 suitible comb.\n",
      "in chunk 4 / 4 124914 suitible comb.\n",
      "total combinations within distance 410333\n",
      "0.0 % gridentry done 0.0  mns\n",
      "60.9 % gridentry done 0.4  mns\n",
      "100 % gridentry done 3.93  mns\n",
      "Dhaka\n",
      "in chunk 1 / 4 13479 suitible comb.\n",
      "in chunk 2 / 4 14607 suitible comb.\n",
      "in chunk 3 / 4 8244 suitible comb.\n",
      "in chunk 4 / 4 6374 suitible comb.\n",
      "total combinations within distance 42704\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 4.83  mns\n",
      "Shijiazhuang\n",
      "in chunk 1 / 3 6195 suitible comb.\n",
      "in chunk 2 / 3 24362 suitible comb.\n",
      "in chunk 3 / 3 7373 suitible comb.\n",
      "total combinations within distance 37930\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 5.38  mns\n",
      " \n",
      "Check grids within UGS\n",
      "0 0.02  mns\n",
      "100 0.67  mns\n",
      "Check grids within UGS\n",
      "0 0.89  mns\n",
      "100 1.05  mns\n",
      "200 1.2  mns\n",
      "Check grids within UGS\n",
      "0 1.26  mns\n",
      "100 1.55  mns\n",
      "200 1.79  mns\n",
      "300 2.03  mns\n",
      "400 2.31  mns\n",
      "500 2.61  mns\n",
      "Check grids within UGS\n",
      "0 2.71  mns\n",
      "Wall time: 10min 49s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 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_networks['nodes'], \n",
    "                           cities_adj['City'],\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_networks['nodes'], # For finding nearest grid entry points\n",
    "                                 thresholds,\n",
    "                                 cities_adj['City'],\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": 39,
   "id": "b98ff692",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addis Ababa 1 / 1 range 0 - 169540\n",
      "0.0 % done 0.0 mns\n",
      "5.9 % done 0.1 mns\n",
      "11.8 % done 0.23 mns\n",
      "17.69 % done 0.39 mns\n",
      "23.59 % done 0.54 mns\n",
      "29.49 % done 0.66 mns\n",
      "35.39 % done 0.81 mns\n",
      "41.29 % done 1.03 mns\n",
      "47.19 % done 1.2 mns\n",
      "53.08 % done 1.35 mns\n",
      "index 96883 No route\n",
      "index 96884 No route\n",
      "index 96885 No route\n",
      "index 96886 No route\n",
      "index 96982 No route\n",
      "index 96983 No route\n",
      "index 96985 No route\n",
      "index 96986 No route\n",
      "index 96987 No route\n",
      "58.98 % done 1.95 mns\n",
      "64.88 % done 2.03 mns\n",
      "70.78 % done 2.12 mns\n",
      "76.68 % done 2.21 mns\n",
      "82.58 % done 2.31 mns\n",
      "88.47 % done 2.39 mns\n",
      "94.37 % done 2.52 mns\n",
      "for 32 routes nearest nodes found\n",
      "100.0 % pathfinding done 2.61 mns\n",
      "formatting done 3.53 mns\n",
      "dissolving done 4.51 mns\n",
      "Damascus 1 / 2 range 0 - 250000\n",
      "0.0 % done 4.52 mns\n",
      "2.48 % done 9.2 mns\n",
      "4.97 % done 9.48 mns\n",
      "7.45 % done 10.13 mns\n",
      "9.93 % done 10.25 mns\n",
      "12.42 % done 10.35 mns\n",
      "14.9 % done 10.51 mns\n",
      "17.38 % done 10.69 mns\n",
      "19.87 % done 10.87 mns\n",
      "22.35 % done 11.1 mns\n",
      "24.83 % done 11.23 mns\n",
      "27.32 % done 11.41 mns\n",
      "29.8 % done 11.53 mns\n",
      "32.28 % done 11.66 mns\n",
      "34.76 % done 11.76 mns\n",
      "37.25 % done 11.85 mns\n",
      "39.73 % done 11.96 mns\n",
      "42.21 % done 12.08 mns\n",
      "44.7 % done 12.19 mns\n",
      "47.18 % done 12.31 mns\n",
      "49.66 % done 12.46 mns\n",
      "52.15 % done 12.57 mns\n",
      "54.63 % done 12.68 mns\n",
      "57.11 % done 12.79 mns\n",
      "59.6 % done 13.18 mns\n",
      "for 1790 routes nearest nodes found\n",
      "62.08 % pathfinding done 13.26 mns\n",
      "formatting done 14.66 mns\n",
      "dissolving done 16.14 mns\n",
      "Damascus 2 / 2 range 250000 - 402704\n",
      "62.08 % done 16.14 mns\n",
      "64.56 % done 16.24 mns\n",
      "67.05 % done 16.33 mns\n",
      "69.53 % done 16.44 mns\n",
      "72.01 % done 16.52 mns\n",
      "74.5 % done 16.82 mns\n",
      "76.98 % done 16.95 mns\n",
      "79.46 % done 17.28 mns\n",
      "index 329255 No route\n",
      "index 329256 No route\n",
      "index 329257 No route\n",
      "index 329258 No route\n",
      "index 329259 No route\n",
      "index 329260 No route\n",
      "index 329261 No route\n",
      "index 329262 No route\n",
      "index 329263 No route\n",
      "index 329264 No route\n",
      "index 329305 No route\n",
      "index 329306 No route\n",
      "index 329307 No route\n",
      "index 329308 No route\n",
      "index 329309 No route\n",
      "index 329310 No route\n",
      "index 329311 No route\n",
      "index 329312 No route\n",
      "index 329313 No route\n",
      "index 329314 No route\n",
      "index 329315 No route\n",
      "index 329316 No route\n",
      "index 329317 No route\n",
      "index 329318 No route\n",
      "index 329319 No route\n",
      "index 329320 No route\n",
      "index 329321 No route\n",
      "index 329322 No route\n",
      "index 329323 No route\n",
      "index 329324 No route\n",
      "index 329325 No route\n",
      "index 329326 No route\n",
      "index 329327 No route\n",
      "index 329328 No route\n",
      "index 329329 No route\n",
      "index 329330 No route\n",
      "index 329331 No route\n",
      "index 329332 No route\n",
      "index 329333 No route\n",
      "index 329334 No route\n",
      "index 329335 No route\n",
      "index 329336 No route\n",
      "index 329337 No route\n",
      "index 329338 No route\n",
      "index 329339 No route\n",
      "index 329340 No route\n",
      "index 329341 No route\n",
      "index 329342 No route\n",
      "index 329343 No route\n",
      "index 329344 No route\n",
      "index 329345 No route\n",
      "index 329346 No route\n",
      "index 329347 No route\n",
      "index 329348 No route\n",
      "index 329349 No route\n",
      "index 329350 No route\n",
      "index 329351 No route\n",
      "index 329352 No route\n",
      "index 329353 No route\n",
      "index 329354 No route\n",
      "index 329355 No route\n",
      "index 329356 No route\n",
      "index 329357 No route\n",
      "index 329358 No route\n",
      "index 329359 No route\n",
      "index 329360 No route\n",
      "index 329361 No route\n",
      "index 329362 No route\n",
      "index 329363 No route\n",
      "index 329364 No route\n",
      "index 329365 No route\n",
      "index 329366 No route\n",
      "index 329367 No route\n",
      "index 329368 No route\n",
      "index 329369 No route\n",
      "index 329370 No route\n",
      "index 329371 No route\n",
      "index 329372 No route\n",
      "index 329373 No route\n",
      "index 329374 No route\n",
      "index 329375 No route\n",
      "index 329376 No route\n",
      "index 329377 No route\n",
      "index 329378 No route\n",
      "index 329379 No route\n",
      "index 329380 No route\n",
      "index 329381 No route\n",
      "index 329382 No route\n",
      "index 329383 No route\n",
      "index 329384 No route\n",
      "index 329385 No route\n",
      "index 329386 No route\n",
      "index 329387 No route\n",
      "index 329388 No route\n",
      "index 329389 No route\n",
      "index 329400 No route\n",
      "index 329401 No route\n",
      "index 329402 No route\n",
      "index 329403 No route\n",
      "index 329404 No route\n",
      "index 329405 No route\n",
      "index 329406 No route\n",
      "index 329407 No route\n",
      "index 329408 No route\n",
      "index 329409 No route\n",
      "index 329410 No route\n",
      "index 329411 No route\n",
      "index 329412 No route\n",
      "index 329413 No route\n",
      "index 329414 No route\n",
      "index 329415 No route\n",
      "index 329416 No route\n",
      "index 329417 No route\n",
      "index 329418 No route\n",
      "index 329419 No route\n",
      "index 329420 No route\n",
      "index 329421 No route\n",
      "index 329422 No route\n",
      "index 329423 No route\n",
      "index 329424 No route\n",
      "index 329425 No route\n",
      "index 329446 No route\n",
      "index 329447 No route\n",
      "index 329448 No route\n",
      "index 329449 No route\n",
      "index 329450 No route\n",
      "index 329451 No route\n",
      "index 329452 No route\n",
      "index 329453 No route\n",
      "index 329454 No route\n",
      "index 329455 No route\n",
      "index 329456 No route\n",
      "index 329457 No route\n",
      "index 329458 No route\n",
      "index 329459 No route\n",
      "index 329460 No route\n",
      "index 329461 No route\n",
      "index 329462 No route\n",
      "index 329463 No route\n",
      "index 329464 No route\n",
      "index 329465 No route\n",
      "index 329466 No route\n",
      "index 329467 No route\n",
      "index 329468 No route\n",
      "index 329469 No route\n",
      "index 329470 No route\n",
      "index 329471 No route\n",
      "index 329472 No route\n",
      "index 329473 No route\n",
      "index 329474 No route\n",
      "index 329475 No route\n",
      "index 329476 No route\n",
      "index 329477 No route\n",
      "index 329478 No route\n",
      "index 329479 No route\n",
      "index 329480 No route\n",
      "index 329481 No route\n",
      "index 329482 No route\n",
      "index 329483 No route\n",
      "index 329484 No route\n",
      "index 329485 No route\n",
      "index 329486 No route\n",
      "index 329487 No route\n",
      "index 329488 No route\n",
      "index 329489 No route\n",
      "81.95 % done 22.08 mns\n",
      "index 339174 No route\n",
      "index 339175 No route\n",
      "index 339176 No route\n",
      "index 339177 No route\n",
      "index 339178 No route\n",
      "index 339179 No route\n",
      "index 339180 No route\n",
      "index 339181 No route\n",
      "index 339182 No route\n",
      "index 339183 No route\n",
      "index 339184 No route\n",
      "index 339185 No route\n",
      "index 339186 No route\n",
      "index 339187 No route\n",
      "index 339188 No route\n",
      "index 339189 No route\n",
      "index 339190 No route\n",
      "index 339191 No route\n",
      "index 339192 No route\n",
      "index 339193 No route\n",
      "index 339194 No route\n",
      "index 339195 No route\n",
      "index 339308 No route\n",
      "index 339309 No route\n",
      "index 339316 No route\n",
      "index 339317 No route\n",
      "index 339318 No route\n",
      "index 339319 No route\n",
      "index 339320 No route\n",
      "index 339321 No route\n",
      "index 339322 No route\n",
      "index 339323 No route\n",
      "index 339324 No route\n",
      "index 339325 No route\n",
      "index 339326 No route\n",
      "index 339327 No route\n",
      "index 339328 No route\n",
      "index 339329 No route\n",
      "index 339330 No route\n",
      "index 339331 No route\n",
      "index 339332 No route\n",
      "index 339333 No route\n",
      "index 339334 No route\n",
      "index 339335 No route\n",
      "index 339336 No route\n",
      "index 339337 No route\n",
      "index 339338 No route\n",
      "index 339339 No route\n",
      "index 339340 No route\n",
      "index 339341 No route\n",
      "index 339342 No route\n",
      "index 339343 No route\n",
      "index 339344 No route\n",
      "index 339345 No route\n",
      "index 339346 No route\n",
      "index 339347 No route\n",
      "index 339348 No route\n",
      "index 339349 No route\n",
      "index 339350 No route\n",
      "index 339351 No route\n",
      "index 339352 No route\n",
      "index 339353 No route\n",
      "index 339354 No route\n",
      "index 339355 No route\n",
      "index 339356 No route\n",
      "index 339357 No route\n",
      "index 339432 No route\n",
      "index 339433 No route\n",
      "84.43 % done 28.32 mns\n",
      "86.91 % done 28.38 mns\n",
      "89.4 % done 28.45 mns\n",
      "91.88 % done 29.03 mns\n",
      "94.36 % done 29.1 mns\n",
      "96.85 % done 29.2 mns\n",
      "99.33 % done 29.47 mns\n",
      "index 402558 No route\n",
      "index 402559 No route\n",
      "index 402560 No route\n",
      "index 402561 No route\n",
      "index 402562 No route\n",
      "index 402563 No route\n",
      "index 402564 No route\n",
      "index 402565 No route\n",
      "index 402566 No route\n",
      "index 402567 No route\n",
      "index 402568 No route\n",
      "index 402569 No route\n",
      "index 402570 No route\n",
      "index 402571 No route\n",
      "index 402572 No route\n",
      "index 402573 No route\n",
      "index 402574 No route\n",
      "index 402575 No route\n",
      "index 402576 No route\n",
      "index 402577 No route\n",
      "index 402578 No route\n",
      "index 402579 No route\n",
      "index 402580 No route\n",
      "index 402581 No route\n",
      "index 402582 No route\n",
      "index 402583 No route\n",
      "index 402584 No route\n",
      "index 402585 No route\n",
      "index 402586 No route\n",
      "index 402587 No route\n",
      "index 402588 No route\n",
      "index 402589 No route\n",
      "index 402590 No route\n",
      "index 402591 No route\n",
      "index 402592 No route\n",
      "index 402593 No route\n",
      "index 402594 No route\n",
      "index 402595 No route\n",
      "index 402596 No route\n",
      "index 402597 No route\n",
      "index 402598 No route\n",
      "index 402599 No route\n",
      "index 402600 No route\n",
      "index 402601 No route\n",
      "index 402602 No route\n",
      "index 402603 No route\n",
      "index 402604 No route\n",
      "index 402605 No route\n",
      "index 402606 No route\n",
      "index 402607 No route\n",
      "index 402608 No route\n",
      "index 402609 No route\n",
      "index 402610 No route\n",
      "index 402611 No route\n",
      "index 402612 No route\n",
      "index 402613 No route\n",
      "index 402614 No route\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "index 402615 No route\n",
      "index 402616 No route\n",
      "index 402617 No route\n",
      "index 402618 No route\n",
      "index 402619 No route\n",
      "index 402620 No route\n",
      "index 402621 No route\n",
      "index 402636 No route\n",
      "index 402637 No route\n",
      "index 402638 No route\n",
      "index 402639 No route\n",
      "index 402640 No route\n",
      "index 402641 No route\n",
      "index 402642 No route\n",
      "index 402643 No route\n",
      "index 402644 No route\n",
      "index 402645 No route\n",
      "index 402646 No route\n",
      "index 402647 No route\n",
      "index 402648 No route\n",
      "index 402649 No route\n",
      "index 402650 No route\n",
      "index 402651 No route\n",
      "index 402652 No route\n",
      "index 402653 No route\n",
      "index 402654 No route\n",
      "index 402655 No route\n",
      "index 402656 No route\n",
      "index 402657 No route\n",
      "index 402658 No route\n",
      "index 402659 No route\n",
      "index 402660 No route\n",
      "index 402661 No route\n",
      "index 402662 No route\n",
      "index 402663 No route\n",
      "index 402664 No route\n",
      "index 402665 No route\n",
      "index 402666 No route\n",
      "index 402667 No route\n",
      "index 402668 No route\n",
      "index 402669 No route\n",
      "index 402670 No route\n",
      "index 402671 No route\n",
      "index 402672 No route\n",
      "index 402673 No route\n",
      "index 402674 No route\n",
      "index 402675 No route\n",
      "index 402676 No route\n",
      "index 402677 No route\n",
      "index 402678 No route\n",
      "index 402679 No route\n",
      "index 402680 No route\n",
      "index 402681 No route\n",
      "index 402682 No route\n",
      "index 402683 No route\n",
      "index 402684 No route\n",
      "index 402685 No route\n",
      "index 402686 No route\n",
      "index 402687 No route\n",
      "index 402688 No route\n",
      "index 402689 No route\n",
      "index 402690 No route\n",
      "index 402691 No route\n",
      "index 402692 No route\n",
      "index 402693 No route\n",
      "for 3431 routes nearest nodes found\n",
      "100.0 % pathfinding done 30.66 mns\n",
      "formatting done 31.38 mns\n",
      "dissolving done 32.21 mns\n",
      "Dhaka 1 / 1 range 0 - 41282\n",
      "0.0 % done 32.23 mns\n",
      "24.22 % done 32.27 mns\n",
      "48.45 % done 32.3 mns\n",
      "72.67 % done 32.34 mns\n",
      "96.89 % done 32.39 mns\n",
      "for 0 routes nearest nodes found\n",
      "100.0 % pathfinding done 32.4 mns\n",
      "formatting done 32.55 mns\n",
      "dissolving done 32.74 mns\n",
      "Shijiazhuang 1 / 1 range 0 - 35050\n",
      "0.0 % done 32.75 mns\n",
      "28.53 % done 32.78 mns\n",
      "57.06 % done 32.82 mns\n",
      "85.59 % done 32.86 mns\n",
      "for 234 routes nearest nodes found\n",
      "100.0 % pathfinding done 33.28 mns\n",
      "formatting done 33.36 mns\n",
      "dissolving done 33.51 mns\n",
      "Wall time: 33min 30s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 4. Finding shortest routes.\n",
    "Routes = route_finding (road_networks['graphs'], # graphs of the road networks\n",
    "               suitible_InOut_UGS, # potential suitible routes with grid-UGS comb. separated in or out UGS.\n",
    "               road_networks['nodes'], \n",
    "               road_networks['edges'], \n",
    "               cities_adj['City'], \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": 48,
   "id": "a49890ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "300 Addis Ababa\n",
      "600 Addis Ababa\n",
      "1000 Addis Ababa\n",
      "300 Damascus\n",
      "600 Damascus\n",
      "1000 Damascus\n",
      "300 Dhaka\n",
      "600 Dhaka\n",
      "1000 Dhaka\n",
      "300 Shijiazhuang\n",
      "600 Shijiazhuang\n",
      "1000 Shijiazhuang\n",
      "Wall time: 56.3 s\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>City</th>\n",
       "      <th>Addis Ababa</th>\n",
       "      <th>Damascus</th>\n",
       "      <th>Dhaka</th>\n",
       "      <th>Shijiazhuang</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>population</th>\n",
       "      <td>3,552,822.00</td>\n",
       "      <td>2,075,713.00</td>\n",
       "      <td>12,219,271.00</td>\n",
       "      <td>2,979,832.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 300</th>\n",
       "      <td>72.47</td>\n",
       "      <td>1.30</td>\n",
       "      <td>3.79</td>\n",
       "      <td>38.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 300</th>\n",
       "      <td>6.83</td>\n",
       "      <td>40.38</td>\n",
       "      <td>15.92</td>\n",
       "      <td>5.09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 300</th>\n",
       "      <td>8,329,597.37</td>\n",
       "      <td>14,559.79</td>\n",
       "      <td>35,965.75</td>\n",
       "      <td>29,968.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 300</th>\n",
       "      <td>74.40</td>\n",
       "      <td>1.05</td>\n",
       "      <td>4.16</td>\n",
       "      <td>41.60</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 300</th>\n",
       "      <td>3.74</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.23</td>\n",
       "      <td>9.15</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 600</th>\n",
       "      <td>64.06</td>\n",
       "      <td>1.30</td>\n",
       "      <td>3.13</td>\n",
       "      <td>36.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 600</th>\n",
       "      <td>30.58</td>\n",
       "      <td>128.85</td>\n",
       "      <td>72.75</td>\n",
       "      <td>19.87</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 600</th>\n",
       "      <td>9,501,572.99</td>\n",
       "      <td>17,957.38</td>\n",
       "      <td>55,841.92</td>\n",
       "      <td>46,858.40</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 600</th>\n",
       "      <td>66.01</td>\n",
       "      <td>0.87</td>\n",
       "      <td>3.28</td>\n",
       "      <td>41.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 600</th>\n",
       "      <td>3.13</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.14</td>\n",
       "      <td>8.52</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 1000</th>\n",
       "      <td>58.99</td>\n",
       "      <td>1.31</td>\n",
       "      <td>3.00</td>\n",
       "      <td>35.52</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 1000</th>\n",
       "      <td>84.73</td>\n",
       "      <td>285.91</td>\n",
       "      <td>132.60</td>\n",
       "      <td>50.30</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 1000</th>\n",
       "      <td>10,952,377.93</td>\n",
       "      <td>20,125.89</td>\n",
       "      <td>69,327.12</td>\n",
       "      <td>67,570.82</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 1000</th>\n",
       "      <td>62.26</td>\n",
       "      <td>0.73</td>\n",
       "      <td>2.82</td>\n",
       "      <td>39.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 1000</th>\n",
       "      <td>2.86</td>\n",
       "      <td>0.01</td>\n",
       "      <td>0.15</td>\n",
       "      <td>7.90</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                    Addis Ababa             Damascus                Dhaka  \\\n",
       "population             3,552,822.00         2,075,713.00        12,219,271.00   \n",
       "Sc-access 300                 72.47                 1.30                 3.79   \n",
       "M-dist 300                     6.83                40.38                15.92   \n",
       "M-area 300             8,329,597.37            14,559.79            35,965.75   \n",
       "M-supply 300                  74.40                 1.05                 4.16   \n",
       "Sc-norm 300                    3.74                 0.01                 0.23   \n",
       "Sc-access 600                 64.06                 1.30                 3.13   \n",
       "M-dist 600                    30.58               128.85                72.75   \n",
       "M-area 600             9,501,572.99            17,957.38            55,841.92   \n",
       "M-supply 600                  66.01                 0.87                 3.28   \n",
       "Sc-norm 600                    3.13                 0.01                 0.14   \n",
       "Sc-access 1000                58.99                 1.31                 3.00   \n",
       "M-dist 1000                   84.73               285.91               132.60   \n",
       "M-area 1000           10,952,377.93            20,125.89            69,327.12   \n",
       "M-supply 1000                 62.26                 0.73                 2.82   \n",
       "Sc-norm 1000                   2.86                 0.01                 0.15   \n",
       "\n",
       "City                   Shijiazhuang  \n",
       "population             2,979,832.00  \n",
       "Sc-access 300                 38.17  \n",
       "M-dist 300                     5.09  \n",
       "M-area 300                29,968.11  \n",
       "M-supply 300                  41.60  \n",
       "Sc-norm 300                    9.15  \n",
       "Sc-access 600                 36.90  \n",
       "M-dist 600                    19.87  \n",
       "M-area 600                46,858.40  \n",
       "M-supply 600                  41.00  \n",
       "Sc-norm 600                    8.52  \n",
       "Sc-access 1000                35.52  \n",
       "M-dist 1000                   50.30  \n",
       "M-area 1000               67,570.82  \n",
       "M-supply 1000                 39.37  \n",
       "Sc-norm 1000                   7.90  "
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, \n",
    "                             population_grids, \n",
    "                             thresholds, \n",
    "                             cities_adj['City'], \n",
    "                             save_path = 'C:/Dumps/GEE-WP Scores/E2SFCA/', \n",
    "                             grid_size = 100)\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "313b0bd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gee_worldpop_extract (city_file, iso, save_path = None):\n",
    "    \n",
    "    cities = city_file\n",
    "    \n",
    "    # Get included city areas\n",
    "    OSM_incl = [cities[cities['City'] == city]['OSM_area'].tolist()[0].rsplit(', ') for city in cities['City'].tolist()]\n",
    "\n",
    "    # Get the city geoms\n",
    "    obj = [city_geo(city).dissolve()['geometry'].tolist()[0] for city in OSM_incl]\n",
    "\n",
    "    # Get the city countries\n",
    "    obj_displ = [city_geo(city).dissolve()['display_name'].tolist()[0].rsplit(', ')[-1]for city in OSM_incl]\n",
    "    obj_displ = np.where(pd.Series(obj_displ).str.contains(\"Ivoire\"),\"CIte dIvoire\",obj_displ)\n",
    "\n",
    "    # Get the country's iso-code\n",
    "    iso_list = [iso[iso['name'] == ob]['alpha3'].tolist()[0] for ob in obj_displ]\n",
    "\n",
    "    # Based on the iso-code return the worldpop 2020\n",
    "    ee_worldpop = [ee.ImageCollection(\"WorldPop/GP/100m/pop\")\\\n",
    "        .filter(ee.Filter.date('2020'))\\\n",
    "        .filter(ee.Filter.inList('country', [io])).first() for io in iso_list]\n",
    "\n",
    "    # Clip the countries with the city geoms.\n",
    "    clipped = [ee_worldpop[i].clip(shapely.geometry.mapping(obj[i])) for i in range(0,len(obj))]\n",
    "\n",
    "    # Create path if non-existent\n",
    "    if save_path == None:\n",
    "        path = ''\n",
    "    else:\n",
    "        path = save_path\n",
    "        if not os.path.exists(path):\n",
    "                    os.makedirs(path)\n",
    "\n",
    "    # Export as TIFF file.\n",
    "    # Stored in form path + USA_Los Angeles_2020.tif\n",
    "    filenames = [path+iso_list[i]+'_'+cities['City'][i]+'_2020.tif' for i in range(len(obj))]\n",
    "    [geemap.ee_export_image(clipped[i], filename = filenames[i]) for i in range(0,len(obj))]\n",
    "    return(filenames)\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "e5a4c548",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 2 population grids extraction\n",
    "def city_grids_format(city_grids, cities_area, 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(city_grids)):\n",
    "        \n",
    "        # Open the raster file\n",
    "        with rasterio.open(city_grids[i]) as src:\n",
    "            band= src.read() # the population values\n",
    "            aff = src.transform # the raster bounds and size (affine)\n",
    "        \n",
    "        # Get the rowwise arrays, get a 2D dataframe\n",
    "        grid = pd.DataFrame()\n",
    "        for b in enumerate(band[0]):\n",
    "            grid = pd.concat([grid, pd.Series(b[1],name=b[0])],axis=1)\n",
    "        grid= grid.unstack().reset_index()\n",
    "        \n",
    "        # Unstack df to columns\n",
    "        grid.columns = ['row','col','value']\n",
    "        grid['minx'] = aff[2]+aff[0]*grid['col']\n",
    "        grid['miny'] = aff[5]+aff[4]*grid['row']\n",
    "        grid['maxx'] = aff[2]+aff[0]*grid['col']+aff[0]\n",
    "        grid['maxy'] = aff[5]+aff[4]*grid['row']+aff[4]\n",
    "        \n",
    "        # Create polygon from affine bounds and row/col indices\n",
    "        grid['geometry'] = [Polygon([(grid.minx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.maxy[i]),\n",
    "                                   (grid.minx[i],grid.maxy[i])])\\\n",
    "                          for i in range(len(grid))]\n",
    "        \n",
    "        # Set the df as geo-df\n",
    "        grid = gpd.GeoDataFrame(grid, crs = 4326) \n",
    "\n",
    "        # Get dissolvement_key for dissolvement. \n",
    "        grid['row3'] = np.floor(grid['row']/(grid_size/100)).astype(int)\n",
    "        grid['col3'] = np.floor(grid['col']/(grid_size/100)).astype(int)\n",
    "        grid['dissolve_key'] = grid['row3'].astype(str) +'-'+ grid['col3'].astype(str)\n",
    "        \n",
    "        # Define a city's OSM area as Polygon.\n",
    "        geo_ls = gpd.GeoSeries(city_geo(cities_area[i].split(', ')).dissolve().geometry)\n",
    "        \n",
    "        # Intersect grids with the city boundary Polygon.\n",
    "        insec = grid.intersection(geo_ls.tolist()[0])\n",
    "        \n",
    "        # Exclude grids outside the specified city boundaries\n",
    "        insec = insec[insec.area > 0]\n",
    "        \n",
    "        # Join in other information.\n",
    "        insec = gpd.GeoDataFrame(geometry = insec, crs = 4326).join(grid.loc[:, grid.columns != 'geometry'])\n",
    "        \n",
    "        # Dissolve into block by block grids\n",
    "        popgrid = insec[['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(insec.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), coords error.\n",
    "        # Therefore recreate the 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(city_grids[i].rsplit('_')[3], round((time.time() - start_time)/60,2),'mns')\n",
    "    return(grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5a8007ac",
   "metadata": {},
   "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 enumerate(cities['OSM_area']):\n",
    "        # Get graph, road nodes and edges\n",
    "        road_node = pd.DataFrame()\n",
    "        roads = pd.DataFrame()\n",
    "        \n",
    "        # For each included OSM_area get the roads\n",
    "        for district in i[1].rsplit(', '):\n",
    "            graph = ox.graph_from_place(district, network_type = \"all\", buffer_dist = (np.max(thresholds)+1000))\n",
    "            node, edge = ox.graph_to_gdfs(graph)\n",
    "            road_node = pd.concat([road_node, node], axis = 0)\n",
    "            roads = pd.concat([roads, edge], axis = 0)\n",
    "        \n",
    "        # Eliminate lists in the df which prevents drop of duplicate columns\n",
    "        road_edge = pd.DataFrame([[c[0] if isinstance(c,list) else c for c in roads[col]]\\\n",
    "                              for col in roads]).transpose()\n",
    "        road_edge.columns = roads.columns\n",
    "        road_edge.index = roads.index\n",
    "        road_edge = gpd.GeoDataFrame(road_edge, crs = 4326)\n",
    "        \n",
    "        # Return the unique nodes and edges of the (often) adjacent OSM_areas.\n",
    "        road_node = road_node.drop_duplicates()\n",
    "        road_edge = road_edge.drop_duplicates()\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['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(cities['City'][i[0]].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": 23,
   "id": "de012d0d",
   "metadata": {},
   "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 enumerate(cities['OSM_area']):\n",
    "        # Tags seen as Urban Greenspace (UGS) require the following:\n",
    "        # 1. Tag represent an area\n",
    "        # 2. The area is outdoor\n",
    "        # 3. The area is (semi-)publically available\n",
    "        # 4. The area is likely to contain trees, grass and/or greenery\n",
    "        # 5. The area can reasonable be used for walking or recreational activities\n",
    "        tags = {'landuse':['allotments','forest','greenfield','village_green'],\\\n",
    "                'leisure':['garden','fitness_station','nature_reserve','park','playground'],\\\n",
    "                'natural':'grassland'}\n",
    "        gdf = ox.geometries_from_place(i[1].rsplit(', '),tags = tags,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(cities['City'][i[0]], '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": 24,
   "id": "06ac19c9",
   "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'].to_crs(4326))\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": 25,
   "id": "0c537f03",
   "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": 26,
   "id": "e354607f",
   "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": 27,
   "id": "6058fb80",
   "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": 28,
   "id": "5f6944b4",
   "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'].to_crs(4326)\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": 29,
   "id": "6d5b61e2",
   "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": 30,
   "id": "854c5949",
   "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('index',i+block*k,'No route')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "3ecd5703",
   "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": 32,
   "id": "50b89e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Improve: 2-to-2 instead of 1-to-all.\n",
    "\n",
    "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",
    "                    try:\n",
    "                        found_route.append(nx.shortest_path(graph, near_grid, near_UGS, \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",
    "                        try:\n",
    "                            found_route.append(nx.shortest_path(graph, near_UGS, near_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(1)\n",
    "                        except:\n",
    "                            pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "1cf68b13",
   "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": 34,
   "id": "4e347555",
   "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": 35,
   "id": "1ce9ee7b",
   "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": 36,
   "id": "3467e686",
   "metadata": {},
   "outputs": [],
   "source": [
    "def E2SCFA_scores(min_gridUGS_comb, grids, thresholds, cities, save_path = 'C:/Dumps/GEE-WP Scores/E2SFCA/', grid_size = 100):\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",
    "        E2SFCA_score = E2SFCA_score.fillna(0)\n",
    "        \n",
    "        if not os.path.exists(save_path+str(grid_size)+'m grids'+'/grid_geoms/'):\n",
    "            os.makedirs(save_path+str(grid_size)+'m grids'+'/grid_geoms/')\n",
    "        \n",
    "        E2SFCA_score.to_file(save_path+str(grid_size)+'m grids'+'/grid_geoms/'+cities[i]+'.gpkg') # 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(save_path):\n",
    "            os.makedirs(save_path)\n",
    "        \n",
    "        E2SFCA_score.loc[:, E2SFCA_score.columns != 'geometry'].to_csv(save_path+cities[i]+'.csv')\n",
    "    E2SFCA_summary.columns = cities\n",
    "    \n",
    "    if not os.path.exists(save_path):\n",
    "        os.makedirs(save_path)\n",
    "    \n",
    "    E2SFCA_summary.to_csv(save_path+str(grid_size)+'m grids'+'all_cities.csv')\n",
    "    E2SFCA_summary\n",
    "    return({'score summary':E2SFCA_summary,'score detail':E2SFCA_cities})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc2146e7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0710956a",
   "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
}