Gravity GEE-OSM Model-Func 2.1.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "675eb1c0",
   "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": 3,
   "id": "d8720c72",
   "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=6Nyq081WQ2fxQk4vf_G50hV-El-gqb9XQ7gCdDos9Zs&tc=zRf1u4bSDil4SW7BAcLzjPIQLU3eRyU-IKx1CYwyy9U&cc=w8N3s9CfLVmyb5vRLFu-42NMmY-mOXNhNyKWWW4c8P4>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=6Nyq081WQ2fxQk4vf_G50hV-El-gqb9XQ7gCdDos9Zs&tc=zRf1u4bSDil4SW7BAcLzjPIQLU3eRyU-IKx1CYwyy9U&cc=w8N3s9CfLVmyb5vRLFu-42NMmY-mOXNhNyKWWW4c8P4</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/1AVHEtk575uHD8-D1tg494NBnIX6eXe9QnW-xLoeWVMhBNniS2D26aLNfXSg\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": 4,
   "id": "674831cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Thresholds and cities\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m\n",
    "\n",
    "# Extract cities list\n",
    "iso = pd.read_excel('iso_countries.xlsx')\n",
    "cities = pd.read_excel('cities.xlsx')\n",
    "cities_adj = cities[cities['City'].isin(['Addis Ababa','Dhaka','Shijiazhuang','Santo Domingo'])]\n",
    "cities_adj = cities_adj.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7de9fb99",
   "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-76a258af2bfa11770a4695c13a1f5b2c: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/e374c5e555e31457cffe53390382762f-45c7cc7aa94fed3f50d21c2f19e15a2e: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/aa1d5bb09020fa6804d00e3e57c83081-aba785176e22ca563efe5ea4e8d70269:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\DOM_Santo Domingo_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/c00d2dec9576c7941a5faf8cf3e6eb04-879651c71c4de2a1d0bce11577d82895:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\CHN_Shijiazhuang_2020.tif\n",
      "Wall time: 15.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 1. Required preprocess for information extraction\n",
    "warnings.filterwarnings('ignore')\n",
    "\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": 35,
   "id": "9e396ef0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Addis Ababa 0.42 mns\n",
      "Dhaka 0.68 mns\n",
      "Santo Domingo 1.78 mns\n",
      "Shijiazhuang 2.14 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Addis Ababa done 1.66 mns\n",
      "Dhaka done 2.68 mns\n",
      "Santo Domingo done 4.66 mns\n",
      "Shijiazhuang done 5.04 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Addis Ababa done\n",
      "Dhaka done\n",
      "Santo Domingo done\n",
      "Shijiazhuang done\n",
      "Wall time: 7min 14s\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": 36,
   "id": "eb027bb5",
   "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.42  mns\n",
      "Addis Ababa 100 % done 0.56  mns\n",
      "Dhaka 0.0 % done 0.57  mns\n",
      "Dhaka 19.4 % done 0.81  mns\n",
      "Dhaka 38.8 % done 1.04  mns\n",
      "Dhaka 58.1 % done 1.3  mns\n",
      "Dhaka 77.5 % done 1.55  mns\n",
      "Dhaka 96.9 % done 1.81  mns\n",
      "Dhaka 100 % done 1.85  mns\n",
      "Santo Domingo 0.0 % done 1.85  mns\n",
      "Santo Domingo 25.6 % done 2.22  mns\n",
      "Santo Domingo 51.2 % done 2.52  mns\n",
      "Santo Domingo 76.7 % done 2.84  mns\n",
      "Santo Domingo 100 % done 3.16  mns\n",
      "Shijiazhuang 0.0 % done 3.16  mns\n",
      "Shijiazhuang 100 % done 3.25  mns\n",
      " \n",
      "get potential (Euclidean) suitible combinations\n",
      "Addis Ababa\n",
      "chunk 1 / 7 131913 suitible comb.\n",
      "chunk 2 / 7 18619 suitible comb.\n",
      "chunk 3 / 7 42384 suitible comb.\n",
      "chunk 4 / 7 17050 suitible comb.\n",
      "chunk 5 / 7 40319 suitible comb.\n",
      "chunk 6 / 7 251030 suitible comb.\n",
      "chunk 7 / 7 202336 suitible comb.\n",
      "total combinations within distance 703651\n",
      "0.0 % gridentry done 0.01  mns\n",
      "35.5 % gridentry done 0.4  mns\n",
      "71.1 % gridentry done 0.77  mns\n",
      "100 % gridentry done 6.39  mns\n",
      "Dhaka\n",
      "chunk 1 / 4 32312 suitible comb.\n",
      "chunk 2 / 4 43932 suitible comb.\n",
      "chunk 3 / 4 18061 suitible comb.\n",
      "chunk 4 / 4 11801 suitible comb.\n",
      "total combinations within distance 106106\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 8.53  mns\n",
      "Santo Domingo\n",
      "chunk 1 / 9 39743 suitible comb.\n",
      "chunk 2 / 9 40758 suitible comb.\n",
      "chunk 3 / 9 35444 suitible comb.\n",
      "chunk 4 / 9 27934 suitible comb.\n",
      "chunk 5 / 9 27941 suitible comb.\n",
      "chunk 6 / 9 88294 suitible comb.\n",
      "chunk 7 / 9 11727 suitible comb.\n",
      "chunk 8 / 9 21231 suitible comb.\n",
      "chunk 9 / 9 26519 suitible comb.\n",
      "total combinations within distance 319591\n",
      "0.0 % gridentry done 0.02  mns\n",
      "78.2 % gridentry done 0.43  mns\n",
      "100 % gridentry done 16.53  mns\n",
      "Shijiazhuang\n",
      "chunk 1 / 3 9055 suitible comb.\n",
      "chunk 2 / 3 28934 suitible comb.\n",
      "chunk 3 / 3 11558 suitible comb.\n",
      "total combinations within distance 49547\n",
      "0.0 % gridentry done 0.01  mns\n",
      "100 % gridentry done 17.85  mns\n",
      " \n",
      "Check grids within UGS\n",
      "0 0.03  mns\n",
      "100 0.69  mns\n",
      "Check grids within UGS\n",
      "0 0.93  mns\n",
      "100 1.2  mns\n",
      "200 1.46  mns\n",
      "300 1.71  mns\n",
      "400 1.95  mns\n",
      "500 2.21  mns\n",
      "Check grids within UGS\n",
      "0 2.28  mns\n",
      "100 2.71  mns\n",
      "200 3.03  mns\n",
      "300 3.37  mns\n",
      "Check grids within UGS\n",
      "0 3.74  mns\n",
      "Wall time: 25min 21s\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'], \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": 37,
   "id": "40c30628",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "comb. by city\n",
      "Addis Ababa 703651\n",
      "Dhaka 106106\n",
      "Santo Domingo 319591\n",
      "Shijiazhuang 49555\n",
      " \n",
      "Addis Ababa 1 / 3 range 0 - 250000\n",
      "0.0 % done 0.08 mns\n",
      "1.77 % done 0.54 mns\n",
      "3.55 % done 1.33 mns\n",
      "5.32 % done 1.84 mns\n",
      "7.09 % done 2.19 mns\n",
      "8.87 % done 2.68 mns\n",
      "10.64 % done 2.83 mns\n",
      "12.41 % done 3.04 mns\n",
      "14.19 % done 3.2 mns\n",
      "15.96 % done 3.32 mns\n",
      "17.73 % done 3.46 mns\n",
      "19.5 % done 3.67 mns\n",
      "21.28 % done 3.81 mns\n",
      "23.05 % done 3.92 mns\n",
      "24.82 % done 4.39 mns\n",
      "26.6 % done 4.68 mns\n",
      "28.37 % done 4.99 mns\n",
      "index 166318 No route\n",
      "index 166319 No route\n",
      "index 166320 No route\n",
      "index 166321 No route\n",
      "index 166322 No route\n",
      "index 166323 No route\n",
      "index 166324 No route\n",
      "index 166325 No route\n",
      "index 166326 No route\n",
      "index 166327 No route\n",
      "index 166328 No route\n",
      "index 166329 No route\n",
      "index 166330 No route\n",
      "index 166331 No route\n",
      "index 166332 No route\n",
      "index 166333 No route\n",
      "index 166334 No route\n",
      "index 166335 No route\n",
      "index 166336 No route\n",
      "index 166337 No route\n",
      "index 166338 No route\n",
      "index 166339 No route\n",
      "index 166340 No route\n",
      "index 166341 No route\n",
      "index 166342 No route\n",
      "index 166343 No route\n",
      "index 166344 No route\n",
      "index 166345 No route\n",
      "index 166346 No route\n",
      "index 166347 No route\n",
      "index 166348 No route\n",
      "index 166349 No route\n",
      "index 166350 No route\n",
      "index 166351 No route\n",
      "index 166352 No route\n",
      "index 166353 No route\n",
      "index 166490 No route\n",
      "index 166491 No route\n",
      "index 166492 No route\n",
      "index 166493 No route\n",
      "index 166494 No route\n",
      "index 166495 No route\n",
      "index 166496 No route\n",
      "index 166497 No route\n",
      "index 166498 No route\n",
      "index 166499 No route\n",
      "index 166500 No route\n",
      "index 166501 No route\n",
      "index 166502 No route\n",
      "index 166503 No route\n",
      "index 166504 No route\n",
      "index 166505 No route\n",
      "index 166506 No route\n",
      "index 166507 No route\n",
      "index 166508 No route\n",
      "index 166509 No route\n",
      "index 166510 No route\n",
      "index 166511 No route\n",
      "index 166512 No route\n",
      "index 166513 No route\n",
      "index 166514 No route\n",
      "index 166515 No route\n",
      "index 166516 No route\n",
      "index 166517 No route\n",
      "index 166518 No route\n",
      "index 166519 No route\n",
      "index 166520 No route\n",
      "index 166521 No route\n",
      "index 166522 No route\n",
      "index 166523 No route\n",
      "index 166524 No route\n",
      "30.14 % done 8.98 mns\n",
      "31.92 % done 9.21 mns\n",
      "33.69 % done 9.39 mns\n",
      "35.46 % done 9.61 mns\n",
      "37.24 % done 9.86 mns\n",
      "39.01 % done 10.07 mns\n",
      "40.78 % done 10.32 mns\n",
      "42.56 % done 10.53 mns\n",
      "240 nearest nodes found\n",
      "44.33 % pathfinding done 10.76 mns\n",
      "formatting done 12.71 mns\n",
      "dissolving done 14.41 mns\n",
      "Addis Ababa 2 / 3 range 250000 - 500000\n",
      "44.33 % done 14.42 mns\n",
      "46.1 % done 14.73 mns\n",
      "47.88 % done 15.09 mns\n",
      "49.65 % done 15.52 mns\n",
      "51.42 % done 15.99 mns\n",
      "53.19 % done 16.56 mns\n",
      "54.97 % done 17.16 mns\n",
      "56.74 % done 17.57 mns\n",
      "58.51 % done 17.92 mns\n",
      "60.29 % done 18.41 mns\n",
      "62.06 % done 19.01 mns\n",
      "63.83 % done 19.35 mns\n",
      "65.61 % done 19.92 mns\n",
      "67.38 % done 20.35 mns\n",
      "69.15 % done 20.69 mns\n",
      "70.93 % done 21.0 mns\n",
      "72.7 % done 21.35 mns\n",
      "74.47 % done 21.65 mns\n",
      "76.25 % done 21.97 mns\n",
      "78.02 % done 22.33 mns\n",
      "79.79 % done 22.72 mns\n",
      "81.57 % done 23.0 mns\n",
      "83.34 % done 23.36 mns\n",
      "85.11 % done 23.71 mns\n",
      "86.88 % done 23.96 mns\n",
      "0 nearest nodes found\n",
      "88.66 % pathfinding done 24.31 mns\n",
      "formatting done 26.58 mns\n",
      "dissolving done 28.46 mns\n",
      "Addis Ababa 3 / 3 range 500000 - 563966\n",
      "88.66 % done 28.47 mns\n",
      "90.43 % done 28.87 mns\n",
      "92.2 % done 29.26 mns\n",
      "93.98 % done 29.82 mns\n",
      "95.75 % done 30.34 mns\n",
      "97.52 % done 30.78 mns\n",
      "99.3 % done 31.06 mns\n",
      "0 nearest nodes found\n",
      "100.0 % pathfinding done 31.27 mns\n",
      "formatting done 31.83 mns\n",
      "dissolving done 32.3 mns\n",
      "Dhaka 1 / 1 range 0 - 103381\n",
      "0.0 % done 32.34 mns\n",
      "9.67 % done 32.43 mns\n",
      "19.35 % done 32.5 mns\n",
      "29.02 % done 32.61 mns\n",
      "38.69 % done 32.73 mns\n",
      "48.36 % done 32.81 mns\n",
      "58.04 % done 32.96 mns\n",
      "67.71 % done 33.13 mns\n",
      "77.38 % done 33.23 mns\n",
      "87.06 % done 33.37 mns\n",
      "96.73 % done 33.43 mns\n",
      "0 nearest nodes found\n",
      "100.0 % pathfinding done 33.49 mns\n",
      "formatting done 34.0 mns\n",
      "dissolving done 34.63 mns\n",
      "Santo Domingo 1 / 2 range 0 - 250000\n",
      "0.0 % done 34.64 mns\n",
      "3.39 % done 34.85 mns\n",
      "6.79 % done 34.96 mns\n",
      "10.18 % done 35.03 mns\n",
      "13.58 % done 35.27 mns\n",
      "16.97 % done 35.36 mns\n",
      "20.37 % done 35.43 mns\n",
      "23.76 % done 35.64 mns\n",
      "27.16 % done 36.49 mns\n",
      "30.55 % done 37.49 mns\n",
      "index 91443 No route\n",
      "index 91444 No route\n",
      "index 91445 No route\n",
      "index 91446 No route\n",
      "index 91447 No route\n",
      "33.95 % done 37.99 mns\n",
      "37.34 % done 38.11 mns\n",
      "40.74 % done 41.64 mns\n",
      "44.13 % done 41.79 mns\n",
      "47.53 % done 42.05 mns\n",
      "50.92 % done 42.21 mns\n",
      "54.32 % done 42.31 mns\n",
      "57.71 % done 45.45 mns\n",
      "61.11 % done 45.63 mns\n",
      "64.5 % done 46.06 mns\n",
      "67.9 % done 46.58 mns\n",
      "71.29 % done 46.78 mns\n",
      "74.69 % done 47.25 mns\n",
      "78.08 % done 48.34 mns\n",
      "81.48 % done 48.98 mns\n",
      "1446 nearest nodes found\n",
      "84.87 % pathfinding done 49.19 mns\n",
      "formatting done 50.73 mns\n",
      "dissolving done 52.31 mns\n",
      "Santo Domingo 2 / 2 range 250000 - 294553\n",
      "84.87 % done 52.32 mns\n",
      "88.27 % done 52.49 mns\n",
      "91.66 % done 52.58 mns\n",
      "95.06 % done 53.01 mns\n",
      "98.45 % done 53.09 mns\n",
      "10 nearest nodes found\n",
      "100.0 % pathfinding done 53.19 mns\n",
      "formatting done 53.5 mns\n",
      "dissolving done 53.76 mns\n",
      "Shijiazhuang 1 / 1 range 0 - 46393\n",
      "0.0 % done 53.8 mns\n",
      "21.55 % done 53.82 mns\n",
      "43.11 % done 53.88 mns\n",
      "64.66 % done 53.93 mns\n",
      "86.22 % done 53.96 mns\n",
      "752 nearest nodes found\n",
      "100.0 % pathfinding done 57.81 mns\n",
      "formatting done 57.95 mns\n",
      "dissolving done 58.18 mns\n",
      "Wall time: 58min 10s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 4. Finding shortest routes.\n",
    "\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "08decce9",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addis Ababa\n",
      "entrance 0.09 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 2.08 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 4.41 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 6.33 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Addis Ababa done 8.04 mns\n",
      "Dhaka\n",
      "entrance 8.09 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 9.07 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 10.19 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 11.23 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Dhaka done 12.15 mns\n",
      "Santo Domingo\n",
      "entrance 12.21 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 13.59 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 15.12 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 16.56 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Santo Domingo done 17.87 mns\n",
      "Shijiazhuang\n",
      "entrance 17.9 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 19.31 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 20.72 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 22.11 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Shijiazhuang done 23.68 mns\n",
      "Wall time: 23min 41s\n",
      "Parser   : 114 ms\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></th>\n",
       "      <th>City</th>\n",
       "      <th>Addis Ababa</th>\n",
       "      <th>Dhaka</th>\n",
       "      <th>Santo Domingo</th>\n",
       "      <th>Shijiazhuang</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.054715</td>\n",
       "      <td>0.018963</td>\n",
       "      <td>0.048196</td>\n",
       "      <td>0.017725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.028417</td>\n",
       "      <td>0.038374</td>\n",
       "      <td>0.061172</td>\n",
       "      <td>0.020078</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.037736</td>\n",
       "      <td>0.080649</td>\n",
       "      <td>0.092689</td>\n",
       "      <td>0.034942</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.879132</td>\n",
       "      <td>0.862014</td>\n",
       "      <td>0.797944</td>\n",
       "      <td>0.927256</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.072605</td>\n",
       "      <td>0.063743</td>\n",
       "      <td>0.111227</td>\n",
       "      <td>0.019429</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.057362</td>\n",
       "      <td>0.105380</td>\n",
       "      <td>0.120699</td>\n",
       "      <td>0.053866</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.088570</td>\n",
       "      <td>0.152223</td>\n",
       "      <td>0.147924</td>\n",
       "      <td>0.061785</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.781464</td>\n",
       "      <td>0.678655</td>\n",
       "      <td>0.620150</td>\n",
       "      <td>0.864920</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.107542</td>\n",
       "      <td>0.165167</td>\n",
       "      <td>0.207606</td>\n",
       "      <td>0.022579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.098166</td>\n",
       "      <td>0.139371</td>\n",
       "      <td>0.156812</td>\n",
       "      <td>0.092906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.136279</td>\n",
       "      <td>0.149666</td>\n",
       "      <td>0.131158</td>\n",
       "      <td>0.110465</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.658013</td>\n",
       "      <td>0.545795</td>\n",
       "      <td>0.504424</td>\n",
       "      <td>0.774049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.052455</td>\n",
       "      <td>0.030043</td>\n",
       "      <td>0.051165</td>\n",
       "      <td>0.017331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.045486</td>\n",
       "      <td>0.068396</td>\n",
       "      <td>0.095850</td>\n",
       "      <td>0.019573</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.069503</td>\n",
       "      <td>0.130404</td>\n",
       "      <td>0.133106</td>\n",
       "      <td>0.029872</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.832556</td>\n",
       "      <td>0.771156</td>\n",
       "      <td>0.719879</td>\n",
       "      <td>0.933223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.075727</td>\n",
       "      <td>0.119634</td>\n",
       "      <td>0.134388</td>\n",
       "      <td>0.018031</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.106110</td>\n",
       "      <td>0.149924</td>\n",
       "      <td>0.178782</td>\n",
       "      <td>0.049080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.134679</td>\n",
       "      <td>0.161483</td>\n",
       "      <td>0.168948</td>\n",
       "      <td>0.059303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.683484</td>\n",
       "      <td>0.568958</td>\n",
       "      <td>0.517882</td>\n",
       "      <td>0.873586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.143390</td>\n",
       "      <td>0.243144</td>\n",
       "      <td>0.248310</td>\n",
       "      <td>0.019352</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.165758</td>\n",
       "      <td>0.173180</td>\n",
       "      <td>0.203543</td>\n",
       "      <td>0.085560</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.153518</td>\n",
       "      <td>0.099939</td>\n",
       "      <td>0.123444</td>\n",
       "      <td>0.101709</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.537334</td>\n",
       "      <td>0.483737</td>\n",
       "      <td>0.424702</td>\n",
       "      <td>0.793379</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.050487</td>\n",
       "      <td>0.021061</td>\n",
       "      <td>0.046970</td>\n",
       "      <td>0.017384</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.035933</td>\n",
       "      <td>0.051803</td>\n",
       "      <td>0.074203</td>\n",
       "      <td>0.018999</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.050630</td>\n",
       "      <td>0.102159</td>\n",
       "      <td>0.112594</td>\n",
       "      <td>0.031190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.862950</td>\n",
       "      <td>0.824976</td>\n",
       "      <td>0.766233</td>\n",
       "      <td>0.932427</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.068493</td>\n",
       "      <td>0.086046</td>\n",
       "      <td>0.110608</td>\n",
       "      <td>0.018218</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.077895</td>\n",
       "      <td>0.124249</td>\n",
       "      <td>0.152799</td>\n",
       "      <td>0.049919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.111782</td>\n",
       "      <td>0.172700</td>\n",
       "      <td>0.169600</td>\n",
       "      <td>0.057317</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.741830</td>\n",
       "      <td>0.617005</td>\n",
       "      <td>0.566993</td>\n",
       "      <td>0.874546</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.109936</td>\n",
       "      <td>0.203525</td>\n",
       "      <td>0.225546</td>\n",
       "      <td>0.020040</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.136287</td>\n",
       "      <td>0.168548</td>\n",
       "      <td>0.183420</td>\n",
       "      <td>0.086654</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.153250</td>\n",
       "      <td>0.127683</td>\n",
       "      <td>0.148020</td>\n",
       "      <td>0.104655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.600527</td>\n",
       "      <td>0.500243</td>\n",
       "      <td>0.443014</td>\n",
       "      <td>0.788651</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.050894</td>\n",
       "      <td>0.019009</td>\n",
       "      <td>0.046231</td>\n",
       "      <td>0.017464</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.031417</td>\n",
       "      <td>0.043458</td>\n",
       "      <td>0.064434</td>\n",
       "      <td>0.018655</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.042942</td>\n",
       "      <td>0.089820</td>\n",
       "      <td>0.100090</td>\n",
       "      <td>0.032568</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.874747</td>\n",
       "      <td>0.847713</td>\n",
       "      <td>0.789245</td>\n",
       "      <td>0.931313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.068050</td>\n",
       "      <td>0.068707</td>\n",
       "      <td>0.102754</td>\n",
       "      <td>0.018544</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.064916</td>\n",
       "      <td>0.114810</td>\n",
       "      <td>0.136369</td>\n",
       "      <td>0.050761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.095388</td>\n",
       "      <td>0.167107</td>\n",
       "      <td>0.160670</td>\n",
       "      <td>0.057671</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.771645</td>\n",
       "      <td>0.649377</td>\n",
       "      <td>0.600207</td>\n",
       "      <td>0.873023</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.104168</td>\n",
       "      <td>0.177191</td>\n",
       "      <td>0.212447</td>\n",
       "      <td>0.020690</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.110630</td>\n",
       "      <td>0.159517</td>\n",
       "      <td>0.170063</td>\n",
       "      <td>0.087984</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.147923</td>\n",
       "      <td>0.146339</td>\n",
       "      <td>0.149748</td>\n",
       "      <td>0.110847</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.637279</td>\n",
       "      <td>0.516953</td>\n",
       "      <td>0.467743</td>\n",
       "      <td>0.780479</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                          Addis Ababa     Dhaka  Santo Domingo  \\\n",
       "entrance_300        1 high       0.054715  0.018963       0.048196   \n",
       "                    2 medium     0.028417  0.038374       0.061172   \n",
       "                    3 low        0.037736  0.080649       0.092689   \n",
       "                    4 no         0.879132  0.862014       0.797944   \n",
       "entrance_600        1 high       0.072605  0.063743       0.111227   \n",
       "                    2 medium     0.057362  0.105380       0.120699   \n",
       "                    3 low        0.088570  0.152223       0.147924   \n",
       "                    4 no         0.781464  0.678655       0.620150   \n",
       "entrance_1000       1 high       0.107542  0.165167       0.207606   \n",
       "                    2 medium     0.098166  0.139371       0.156812   \n",
       "                    3 low        0.136279  0.149666       0.131158   \n",
       "                    4 no         0.658013  0.545795       0.504424   \n",
       "gravity**(1/2)_300  1 high       0.052455  0.030043       0.051165   \n",
       "                    2 medium     0.045486  0.068396       0.095850   \n",
       "                    3 low        0.069503  0.130404       0.133106   \n",
       "                    4 no         0.832556  0.771156       0.719879   \n",
       "gravity**(1/2)_600  1 high       0.075727  0.119634       0.134388   \n",
       "                    2 medium     0.106110  0.149924       0.178782   \n",
       "                    3 low        0.134679  0.161483       0.168948   \n",
       "                    4 no         0.683484  0.568958       0.517882   \n",
       "gravity**(1/2)_1000 1 high       0.143390  0.243144       0.248310   \n",
       "                    2 medium     0.165758  0.173180       0.203543   \n",
       "                    3 low        0.153518  0.099939       0.123444   \n",
       "                    4 no         0.537334  0.483737       0.424702   \n",
       "gravity**(1/3)_300  1 high       0.050487  0.021061       0.046970   \n",
       "                    2 medium     0.035933  0.051803       0.074203   \n",
       "                    3 low        0.050630  0.102159       0.112594   \n",
       "                    4 no         0.862950  0.824976       0.766233   \n",
       "gravity**(1/3)_600  1 high       0.068493  0.086046       0.110608   \n",
       "                    2 medium     0.077895  0.124249       0.152799   \n",
       "                    3 low        0.111782  0.172700       0.169600   \n",
       "                    4 no         0.741830  0.617005       0.566993   \n",
       "gravity**(1/3)_1000 1 high       0.109936  0.203525       0.225546   \n",
       "                    2 medium     0.136287  0.168548       0.183420   \n",
       "                    3 low        0.153250  0.127683       0.148020   \n",
       "                    4 no         0.600527  0.500243       0.443014   \n",
       "gravity**(1/5)_300  1 high       0.050894  0.019009       0.046231   \n",
       "                    2 medium     0.031417  0.043458       0.064434   \n",
       "                    3 low        0.042942  0.089820       0.100090   \n",
       "                    4 no         0.874747  0.847713       0.789245   \n",
       "gravity**(1/5)_600  1 high       0.068050  0.068707       0.102754   \n",
       "                    2 medium     0.064916  0.114810       0.136369   \n",
       "                    3 low        0.095388  0.167107       0.160670   \n",
       "                    4 no         0.771645  0.649377       0.600207   \n",
       "gravity**(1/5)_1000 1 high       0.104168  0.177191       0.212447   \n",
       "                    2 medium     0.110630  0.159517       0.170063   \n",
       "                    3 low        0.147923  0.146339       0.149748   \n",
       "                    4 no         0.637279  0.516953       0.467743   \n",
       "\n",
       "City                          Shijiazhuang  \n",
       "entrance_300        1 high        0.017725  \n",
       "                    2 medium      0.020078  \n",
       "                    3 low         0.034942  \n",
       "                    4 no          0.927256  \n",
       "entrance_600        1 high        0.019429  \n",
       "                    2 medium      0.053866  \n",
       "                    3 low         0.061785  \n",
       "                    4 no          0.864920  \n",
       "entrance_1000       1 high        0.022579  \n",
       "                    2 medium      0.092906  \n",
       "                    3 low         0.110465  \n",
       "                    4 no          0.774049  \n",
       "gravity**(1/2)_300  1 high        0.017331  \n",
       "                    2 medium      0.019573  \n",
       "                    3 low         0.029872  \n",
       "                    4 no          0.933223  \n",
       "gravity**(1/2)_600  1 high        0.018031  \n",
       "                    2 medium      0.049080  \n",
       "                    3 low         0.059303  \n",
       "                    4 no          0.873586  \n",
       "gravity**(1/2)_1000 1 high        0.019352  \n",
       "                    2 medium      0.085560  \n",
       "                    3 low         0.101709  \n",
       "                    4 no          0.793379  \n",
       "gravity**(1/3)_300  1 high        0.017384  \n",
       "                    2 medium      0.018999  \n",
       "                    3 low         0.031190  \n",
       "                    4 no          0.932427  \n",
       "gravity**(1/3)_600  1 high        0.018218  \n",
       "                    2 medium      0.049919  \n",
       "                    3 low         0.057317  \n",
       "                    4 no          0.874546  \n",
       "gravity**(1/3)_1000 1 high        0.020040  \n",
       "                    2 medium      0.086654  \n",
       "                    3 low         0.104655  \n",
       "                    4 no          0.788651  \n",
       "gravity**(1/5)_300  1 high        0.017464  \n",
       "                    2 medium      0.018655  \n",
       "                    3 low         0.032568  \n",
       "                    4 no          0.931313  \n",
       "gravity**(1/5)_600  1 high        0.018544  \n",
       "                    2 medium      0.050761  \n",
       "                    3 low         0.057671  \n",
       "                    4 no          0.873023  \n",
       "gravity**(1/5)_1000 1 high        0.020690  \n",
       "                    2 medium      0.087984  \n",
       "                    3 low         0.110847  \n",
       "                    4 no          0.780479  "
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# 5. summarize scores\n",
    "grid_scores = grid_score_summary (Routes['route summary'], # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'C:/Dumps/GEE-WP Scores/Gravity/\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "60b25137",
   "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": 11,
   "id": "a44b205a",
   "metadata": {
    "scrolled": false
   },
   "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": 12,
   "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 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": 13,
   "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 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": 14,
   "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']], 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",
    "        # Get size inflation factors for the gravity model\n",
    "        ParkRoad['size_infl_factor'] = ParkRoad['walk_area'] / ParkRoad['walk_area'].median()\n",
    "        ParkRoad['size_infl_sqr2'] = ParkRoad['size_infl_factor']**(1/2)\n",
    "        ParkRoad['size_infl_sqr3'] = ParkRoad['size_infl_factor']**(1/3)\n",
    "        ParkRoad['size_infl_sqr5'] = ParkRoad['size_infl_factor']**(1/5)\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",
    "        \n",
    "    return(ParkRoads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "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": 16,
   "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",
    "                                       'size_infl_sqr2','size_infl_sqr3','size_infl_sqr5','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",
    "            infl2 = list(node_merged['size_infl_sqr2'])\n",
    "            infl3 = list(node_merged['size_infl_sqr3'])\n",
    "            infl5 = list(node_merged['size_infl_sqr5'])\n",
    "\n",
    "            # Get the euclidean distances\n",
    "            mat = list()\n",
    "            mat2 = list()\n",
    "            mat3 = list()\n",
    "            mat4 = 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",
    "                mat2.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl2[j])\n",
    "                mat3.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl3[j])\n",
    "                mat4.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl5[j])\n",
    "\n",
    "            # Check if distances are within 1000m and join remaining info and concat in master df per 1000.\n",
    "            mat_df = pd.DataFrame(mat3)[(np.array(mat) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat2) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat3) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat4) <= np.max(thresholds))]\n",
    "\n",
    "            # join the other gravity euclidean scores and other information\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat), lsuffix='_infl', rsuffix='_entr', how = 'left')\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat2), lsuffix='_entry', rsuffix='_pwr', how = 'left')\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat4), lsuffix='_pwr', rsuffix='_root', how = 'left')\n",
    "            mat_df.columns = ['size_infl_eucl2','raw euclidean','size_infl_eucl3','size_infl_eucl5']    \n",
    "            mat_df = mat_df.join(node_merged)\n",
    "\n",
    "            output = pd.concat([output, mat_df])\n",
    "\n",
    "            print('chunk',(i+1),'/',len4,len(mat_df),'suitible comb.')\n",
    "        # Renaming columns\n",
    "        print('total combinations within distance',len(output))\n",
    "\n",
    "        output.columns = ['size_infl_eucl3','raw euclidean','size_infl_eucl2','size_infl_eucl5',\n",
    "                          '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',\n",
    "                          'park_lon','park_lat','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5',\n",
    "                          'parkshare_walked','park_area','walk_area_m2']\n",
    "\n",
    "        output = output[['raw euclidean','size_infl_eucl2','size_infl_eucl3','size_infl_eucl5',\n",
    "                         'Grid_No','Park_entry_No','Grid_coords_centroid','Grid_m_centroid',\n",
    "                          'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid',\n",
    "                         'walk_area_m2','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5']]\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": 17,
   "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": 1,
   "id": "a8f10468",
   "metadata": {},
   "outputs": [],
   "source": [
    "def grids_in_UGS (suitible_comb, UGS, pop_grid): \n",
    "    print('grids in UGS')\n",
    "    start_time = time.time()\n",
    "    RoadInOut = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        gridUGS = pop_grid[i]['centroid'].intersection(UGS[i].dissolve().geometry[0]).is_empty == False\n",
    "        gridUGS.name = 'in_out_UGS'\n",
    "        merged = pd.merge(suitible_comb[i], gridUGS, left_on = 'Grid_No', right_index = True)\n",
    "        RoadInOut.append(merged)\n",
    "        print(i)\n",
    "    return(RoadInOut) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "07b96b77",
   "metadata": {
    "scrolled": false
   },
   "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",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    start_time = time.time()\n",
    "    \n",
    "    print('comb. by city')\n",
    "    for n in enumerate(cities): # Know how much comb. need to be calculcated.\n",
    "        print(n[1], len(combinations[n[0]]))\n",
    "    print(' ')\n",
    "    \n",
    "    Routes = list()\n",
    "    Routes_detail = list()\n",
    "    for j in range(len(cities)):\n",
    "        suit_raw = combinations[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",
    "        \n",
    "        len2 = int(np.ceil(len(suitible)/block_size)) # get number of blocks (chunks)\n",
    "        Route_parts = pd.DataFrame()\n",
    "        Route_dparts = pd.DataFrame()\n",
    "\n",
    "        # Divide in chunks of block for computational load\n",
    "        for k in range(len2):    \n",
    "            suitible_chunk = suitible.iloc[k*block_size:k*block_size+block_size] # Get block ids\n",
    "\n",
    "            parknode = list(suitible_chunk['Parkroad_osmid']) # UGS road entry ids\n",
    "            gridnode = list(suitible_chunk['grid_osm']) # grid centroid road entry ids\n",
    "\n",
    "            s_mat = list([]) # osmid from\n",
    "            s_mat1 = list([]) # osmid to\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 = [] # sums number of routes containing nearest nodes.\n",
    "            len1 = len(suitible_chunk)\n",
    "\n",
    "            print(cities[j].rsplit(',')[0], k+1,'/',len2, \n",
    "                  'range',k*block_size,'-',k*block_size+np.where(k*block_size+block_size >= len1,len1,block_size))\n",
    "            \n",
    "            for i in range(len(suitible_chunk)):\n",
    "                try:\n",
    "                    shortest = nx.shortest_path(graphs[j], 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_size*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(graphs[j], 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_size*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 the nearest node is taken, which is iterated X times at max, check assumptions, block #0 \n",
    "                        nn_route_finding(graphs[j], suitible_chunk, road_nodes[j],\n",
    "                                   s_mat, s_mat1, s_mat2, s_mat3, s_mat4, s_mat5, mat_nn, # matrice info see above\n",
    "                                   it = i, block = k, block_size = block_size, \n",
    "                                         nn_iter = 10) # max nearest nodes to be found\n",
    "                        \n",
    "                if i % 10000 == 0: print(round((i+block_size*k)/len(suitible)*100,2),'% done',\n",
    "                                         round((time.time() - start_time) / 60,2),'mns')\n",
    "            print(len(mat_nn),'nearest nodes found')\n",
    "\n",
    "            print(round((i+block_size*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])\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)\n",
    "            print('dissolving done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            # Concats chunk with others already calculated\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','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5',\n",
    "                                   'raw 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','raw_Tcost','grav2_Tcost','grav3_Tcost','grav5_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({'route summary':Routes,'route detail':Routes_detail})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "460e1c8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_route_finding (Graph, comb, nodes, \n",
    "                mat_from, mat_to, mat_route, mat_step, mat_way, mat_wbin, mat_nn, \n",
    "                it, block, block_size = 250000, nn_iter = 10):\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",
    "    len3 = 0\n",
    "    alt_route = list([])\n",
    "    \n",
    "    gosm = comb['grid_osm'] # grid osmids (origin)\n",
    "    posm = comb['Parkroad_osmid'] # UGS osmids (destination)\n",
    "    node = nodes['geometry'] # road node geoms\n",
    "    node_osm = nodes['osmid_var'] # road node osmids\n",
    "    \n",
    "    while len3 < nn_iter and len(alt_route) < 1: # continue if no more than 10 nearest nodes or if a route is found\n",
    "        \n",
    "        len3 = len3 +1\n",
    "        # Finds nearest node per iteration.\n",
    "        nn = nn_finding(gosm, posm, node, node_osm, it, len3)\n",
    "        \n",
    "         # routing within graph and current and found nearest nodes of grids and UGS\n",
    "        nn_routing(Graph, nn['curr_park'], nn['near_park'], nn['curr_grid'], nn['near_grid'],\n",
    "                        mat_way, mat_wbin, alt_route, len3)\n",
    "        \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 append\n",
    "        mat_nn.append(it+block_size*block)\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(it+block_size*block,len4)))\n",
    "        mat_step.append(list(np.arange(0, len4)))\n",
    "    else: # if no route is found fill values.\n",
    "        mat_from.append(-1)\n",
    "        mat_to.append(-1)\n",
    "        mat_route.append(it+block_size*block)\n",
    "        mat_step.append(-1)\n",
    "        mat_way.append('no way')\n",
    "        mat_wbin.append(2)\n",
    "        print('index',it+block_size*block,'No route')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b7b34fa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_finding (grid_osmid, UGS_osmid, node_geom, node_osmid, it, nn_i):\n",
    "    # Grid nearest\n",
    "    g_geom = node_geom[node_osmid == int(grid_osmid[it:it+1])] # Get current grid road entry geometry\n",
    "    g_nearest = pd.DataFrame((abs(float(g_geom.x) - node_geom.x)**2 # Find nearest.\n",
    "                              +abs(float(g_geom.y) - node_geom.y)**2)**(1/2)\n",
    "                            ).join(node_osmid).sort_values(0)\n",
    "\n",
    "    g_grid = g_nearest.iloc[nn_i,1] # Take '1' because 0 will get the current node with distance 0.\n",
    "    g_park = list(UGS_osmid)[it]\n",
    "\n",
    "    p_geom = node_geom[node_osmid == int(UGS_osmid[it:it+1])] # Get current UGS raod entry geometry\n",
    "    p_nearest = pd.DataFrame((abs(float(p_geom.x) - node_geom.x)**2 # Find nearest\n",
    "                              +abs(float(p_geom.y) - node_geom.y)**2)**(1/2)\n",
    "                            ).join(node_osmid).sort_values(0)\n",
    "\n",
    "    p_grid = list(grid_osmid)[it]\n",
    "    p_park = p_nearest.iloc[nn_i,1] # Take '1' because 0 will get the current node with distance 0.\n",
    "    \n",
    "    return({'curr_park':p_grid, 'near_park':p_park, 'curr_grid':g_park, 'near_grid':g_grid}) # return as dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "79bfe490",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_routing (Graph, curr_park, near_park, curr_grid, near_grid, mat_way, mat_wbin, found_route, nn_i):\n",
    "    try: # First try from current grid to nearest UGS id.\n",
    "        found_route.append(nx.shortest_path(Graph, curr_park, near_park, \n",
    "                                          'travel_dist', method = 'dijkstra'))\n",
    "        mat_way.append(str(nn_i)+'grid > n-park')\n",
    "        mat_wbin.append(1)\n",
    "    except:\n",
    "        try: # Else try the reverse.\n",
    "            found_route.append(nx.shortest_path(Graph, near_park, curr_park, \n",
    "                                              'travel_dist', method = 'dijkstra'))\n",
    "            mat_way.append(str(nn_i)+'n-park > grid')\n",
    "            mat_wbin.append(0)\n",
    "        except:\n",
    "            try: # If no success try from current UGS id to nearest grid id\n",
    "                found_route.append(nx.shortest_path(Graph, near_grid, curr_grid, \n",
    "                                                  'travel_dist', method = 'dijkstra'))\n",
    "                mat_way.append(str(nn_i)+'n-grid > park')\n",
    "                mat_wbin.append(1)\n",
    "            except:\n",
    "                try: # Else try the reverse\n",
    "                    found_route.append(nx.shortest_path(Graph, curr_grid, near_grid, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                    mat_way.append(str(nn_i)+'park > n-grid')\n",
    "                    mat_wbin.append(0)\n",
    "                except: # if no routes are found pass.\n",
    "                    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f4984b15",
   "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)):\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')\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": 24,
   "id": "90524c2b",
   "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','size_infl_sqr2','size_infl_sqr3',\n",
    "                                          'size_infl_sqr5','raw 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['raw_Tcost'] = routes2['route_cost'] + routes2['real_G-entry']\n",
    "    routes2['grav2_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr2']\n",
    "    routes2['grav3_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr3']\n",
    "    routes2['grav5_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr5']\n",
    "    return(routes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7ffd4567",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 8 determine best parkentry points from each grid, then calculate grid scores\n",
    "# and finally aggregate city access in categories (high, medium, low and no access)\n",
    "def grid_score_summary (routes, cities, pop_grids, ext = '', grid_size = 100, save_path = 'C:/Dumps/GEE-WP Scores/Gravity/'):\n",
    "    start_time = time.time()\n",
    "    popg_acc = pd.DataFrame()\n",
    "    grid_scores = list([])\n",
    "    gridpark = list([])\n",
    "    for n in range(len(cities)):    \n",
    "        print(cities[n])\n",
    "\n",
    "        # For the four distance decay variants regarding park size.\n",
    "        l1 = list(['raw','grav2','grav3','grav5'])\n",
    "        m1 = list(['entrance','gravity**(1/2)','gravity**(1/3)','gravity**(1/5)'])\n",
    "        grid_score = list([])\n",
    "        gridparks = list([])\n",
    "        gridpark.append(gridparks)\n",
    "        popgrid_access = pd.DataFrame()\n",
    "        for i in range(len(l1)):\n",
    "            # Get the lowest indices grouped by a key consisting of grid no and park no (best entry point from a grid to a park)\n",
    "            var_best_routes = best_gridpark_comb (routes[n], l1[i], pop_grids[n])\n",
    "\n",
    "            grdsc = pd.DataFrame()\n",
    "            gridsc = pd.DataFrame()\n",
    "            print(m1[i], round((time.time() - start_time) / 60,2), 'mns')\n",
    "\n",
    "            # For each threshold given, calculate a score\n",
    "            for k in range(len(thresholds)):\n",
    "                \n",
    "                t = thresholds[k]\n",
    "                score = 'tr_'+ str(t)\n",
    "                scores = determine_scores(var_best_routes, pop_grids[n], thresholds[k], l1[i], cities[n], \n",
    "                                          save_path, grid_size = 100)\n",
    "                \n",
    "                grdsc = pd.concat([grdsc, scores['score_w_route']], axis = 1)\n",
    "                gridsc = pd.concat([gridsc, scores['grid_score']])\n",
    "                                \n",
    "                # Group according to the categories just created and sum the populations living in those grids\n",
    "                popgacc = pd.DataFrame()\n",
    "                popgacc[m1[i]+'_'+str(t)] = scores['score_w_route'].groupby(score+'_access')['population'].sum()\n",
    "                popgrid_access = pd.concat([popgrid_access, popgacc],axis=1)   \n",
    "\n",
    "                print('grid ',t)\n",
    "\n",
    "            grid_score.append(grdsc)\n",
    "\n",
    "            gridsc = gridsc.join(pop_grids[n]['geometry'])\n",
    "            gridsc = gpd.GeoDataFrame(gridsc, geometry = 'geometry', crs = 4326)\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",
    "            gridsc.to_file(save_path+str(grid_size)+'m grids/Grid_geoms/gridscore_'+ l1[i] + '_' + cities[n] + '.gpkg')\n",
    "\n",
    "            # Detailed scores to files number of cities * ways to measure = number of files.\n",
    "            # Different threshold-scores are in the same dataframe\n",
    "            gridsc = gridsc.loc[:, gridsc.columns!='geometry']\n",
    "\n",
    "            if not os.path.exists(save_path+str(grid_size)+'m grids/Grid_csv/'):\n",
    "                os.makedirs(save_path+str(grid_size)+'m grids/Grid_csv/')\n",
    "\n",
    "            gridsc.to_csv(save_path+str(grid_size)+'m grids/Grid_csv/gridscore_'+ l1[i] + '_' + cities[n] + '.csv')\n",
    "            gridparks.append(var_best_routes)\n",
    "\n",
    "        grid_scores.append(grid_score)\n",
    "\n",
    "        # For each city, divide the population access by group by the total to get its share.\n",
    "        popgrid_access = popgrid_access / popgrid_access.sum()\n",
    "        popgrid_access = pd.DataFrame(popgrid_access.unstack())\n",
    "        popg_acc = pd.concat([popg_acc, popgrid_access], axis = 1)\n",
    "\n",
    "        print(cities[n],'done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "    popg_acc.columns = cities\n",
    "    popg_acc.to_csv(save_path+str(grid_size)+'m grids/popgrid_access.csv')\n",
    "    return(popg_acc)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "8c00432c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_gridpark_comb (routes, var_abbr, pop_grid):\n",
    "    Rclean = routes[routes['way_calc'] != 'no way'].reset_index()\n",
    "    str1 = 'gridpark_' + var_abbr\n",
    "    locals()[str1] = Rclean.iloc[Rclean.groupby('gridpark_no')[(str(var_abbr) +'_Tcost')].idxmin()]  \n",
    "\n",
    "    # Get grid information\n",
    "    locals()[str1] = pd.merge(locals()[str1], pop_grid[['population','geometry']],\n",
    "                            left_on = 'Grid_No', right_index = True, how = 'outer')\n",
    "    locals()[str1] = locals()[str1].reset_index()\n",
    "\n",
    "    # formatting\n",
    "    locals()[str1]['Park_No'] = locals()[str1]['Park_No'].fillna(-1)\n",
    "    locals()[str1]['Park_No'] = locals()[str1]['Park_No'].astype(int)\n",
    "    locals()[str1]['Park_entry_No'] = locals()[str1]['Park_entry_No'].fillna(-1)\n",
    "    locals()[str1]['Park_entry_No'] = locals()[str1]['Park_entry_No'].astype(int)\n",
    "    return(locals()[str1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d3f1ff5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def determine_scores(var_df, pop_grid, thresholds, var_abbr, city, save_path, grid_size = 100):\n",
    "    t = thresholds\n",
    "    str2 = str(t)\n",
    "    score = 'tr_'+ str2\n",
    "\n",
    "    #Only get routes within the threshold given (it loops over every threshold) and calculate the scores\n",
    "    thold = var_df[var_df[var_abbr + '_Tcost'] <= t]\n",
    "    thold[score] = t - thold[var_abbr + '_Tcost']\n",
    "    thold['pop' + score] = thold[score] * thold['population']\n",
    "    thold['walk_area_ha' + str2] = var_df['walk_area_m2'] /10000\n",
    "    thold['walkha_person' + str2] = thold['population'] / thold['walk_area_ha' + str2]\n",
    "\n",
    "    # Join the gridpark information from before.\n",
    "    var_df = var_df.join(thold[[score,'pop' + score,'walk_area_ha' + str2, 'walkha_person' + str2]])\n",
    "    # get the grid_scores\n",
    "    gs = pd.DataFrame()\n",
    "    gs[[score,'pop_' + score,'walkha_' + str2]] = var_df.groupby(\n",
    "            'Grid_No')[score,'pop' + score, 'walk_area_ha' + str2].sum()\n",
    "\n",
    "    gs['walkha_person_' + score] = var_df.groupby('Grid_No')['walkha_person' + str2].mean()\n",
    "\n",
    "    trstr = var_df[var_df[score] > 0]\n",
    "    gs[score + '_parks'] = trstr.groupby('Grid_No')['gridpark_no'].count()\n",
    "\n",
    "    # Add the routes as a dissolved line_geom\n",
    "    gs[score + '_routes'] = gpd.GeoDataFrame(trstr[['Grid_No','geometry_x']],\n",
    "                                                  geometry = 'geometry_x', crs = 4326).dissolve('Grid_No')\n",
    "\n",
    "    # Add parks which grids have access to with its closest access point\n",
    "    gs[score+'Park:entry'] = trstr[trstr['Park_No'] >=0].groupby('Grid_No')['Park_No'].apply(list).astype(str\n",
    "    ) + ':' + trstr[trstr['Park_entry_No'] >=0].groupby('Grid_No')['Park_entry_No'].apply(list).astype(str)\n",
    "                \n",
    "    # determine the thresholds category-score. \n",
    "    # High >= threshold (perfect score to one park), medium is above half perfect, \n",
    "    # low is below this and no is no access to a park for a certain grid within the threshold given\n",
    "    gs[score+'_access'] = np.select([gs[score] >= t, (gs[score] < t) & (\n",
    "    gs[score]>= t/2), (gs[score] < t/2) & (gs[score]> 0), gs[score] <= 0],\n",
    "          ['1 high','2 medium','3 low','4 no'])\n",
    "    gs = gs.join(pop_grid['population'], how = 'outer')\n",
    "            \n",
    "    gs = gpd.GeoDataFrame(gs, geometry = score + '_routes', crs = 4326)\n",
    "            \n",
    "    if not os.path.exists(save_path+str(grid_size)+'m grids/Grid_lines/'):\n",
    "        os.makedirs(save_path+str(grid_size)+'m grids/Grid_lines/')\n",
    "                \n",
    "    gs.to_file(save_path+str(grid_size)+'m grids/Grid_lines/gridscore_'+ var_abbr + '_' + str2 + '_' + city + '.gpkg')\n",
    "            \n",
    "    gsc = gs.loc[:,~gs.columns.isin([score + '_routes'])]\n",
    "\n",
    "    return({'grid_score':gsc,'score_w_route':gs})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a28fb0a1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e0d37c8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81453915",
   "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.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}