Gravity GEE-OSM Model 3.1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "71639573",
   "metadata": {},
   "outputs": [],
   "source": [
    "# system packages\n",
    "import sys\n",
    "import time\n",
    "import warnings\n",
    "import os\n",
    "\n",
    "# non-geo numeric packages\n",
    "import numpy as np\n",
    "import math\n",
    "from itertools import product, combinations\n",
    "import pandas as pd\n",
    "\n",
    "# network and OSM packages\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "city_geo = ox.geocoder.geocode_to_gdf\n",
    "\n",
    "# Earth engine packages\n",
    "import ee\n",
    "import geemap\n",
    "\n",
    "# General geo-packages\n",
    "import libpysal\n",
    "import rasterio\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c23309aa",
   "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=tsCuP-qjW4JeUmhgwpYpwO9vWjmZUZy-MIMkci8oBgE&tc=2Yc3pLjw9noTLpdqKBAlseUCcuUo4viITq7QnrOSIBk&cc=-zV7pBn4zGMQxv0BZy4quUZlRJPCvRKagcWBSV05TUA>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=tsCuP-qjW4JeUmhgwpYpwO9vWjmZUZy-MIMkci8oBgE&tc=2Yc3pLjw9noTLpdqKBAlseUCcuUo4viITq7QnrOSIBk&cc=-zV7pBn4zGMQxv0BZy4quUZlRJPCvRKagcWBSV05TUA</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/1AbUR2VNHdBMb5WnD-iCCuJknJxVQknFuYm5GFjzJ5UyAWaU3grycF-9Pakg\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": 12,
   "id": "6cc4fa99",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['China']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/747b1673c81bd682790926900b1b8d50-cf06eb95fab2260ad60ca4e9e3ce799e:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\CHN_Zibo_2020.tif\n",
      "get road networks from OSM\n",
      "Zibo done 0.11 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Zibo done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Zibo 0.65 mns\n",
      "\n",
      "Zibo\n",
      "0.0 % fake entry points done 0.0  mns\n",
      "Zibo 100 % fake entry points done 0.06  mns\n",
      "Zibo 100 % done 0.07  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "100 % finding combinations done\n",
      "Zibo 185966 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Zibo\n",
      "0.0 % done 0.03 mns\n",
      "100 % done 0.05 mns\n",
      "\n",
      "Zibo\n",
      "0.0 % 0.12 mns\n",
      "48.46 % 0.69 mns\n",
      "96.93 % 1.2 mns\n",
      "100 % done 1.7 mns\n",
      "\n",
      "Zibo\n",
      "entrance 0.04 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 1.04 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 2.17 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 3.18 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Zibo done 4.12 mns\n",
      "CPU times: total: 5min 27s\n",
      "Wall time: 7min 1s\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>Zibo</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.061045</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.010304</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.028739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.899912</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.070557</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.035164</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.067924</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.826355</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.099716</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.061056</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.154426</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.684802</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.060585</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.017913</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.043578</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.877925</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.067843</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.060586</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.165472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.706099</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.114434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.171624</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.250295</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.463647</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.060516</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.013619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.032726</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.893139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.067402</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.044538</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.120300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.767761</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.097625</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.114762</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.220840</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.566774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.060577</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.012341</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.029410</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.897672</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.067514</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.040294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.090859</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.801333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.095977</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.080162</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.192615</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.631246</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                              Zibo\n",
       "entrance_300        1 high    0.061045\n",
       "                    2 medium  0.010304\n",
       "                    3 low     0.028739\n",
       "                    4 no      0.899912\n",
       "entrance_600        1 high    0.070557\n",
       "                    2 medium  0.035164\n",
       "                    3 low     0.067924\n",
       "                    4 no      0.826355\n",
       "entrance_1000       1 high    0.099716\n",
       "                    2 medium  0.061056\n",
       "                    3 low     0.154426\n",
       "                    4 no      0.684802\n",
       "gravity**(1/2)_300  1 high    0.060585\n",
       "                    2 medium  0.017913\n",
       "                    3 low     0.043578\n",
       "                    4 no      0.877925\n",
       "gravity**(1/2)_600  1 high    0.067843\n",
       "                    2 medium  0.060586\n",
       "                    3 low     0.165472\n",
       "                    4 no      0.706099\n",
       "gravity**(1/2)_1000 1 high    0.114434\n",
       "                    2 medium  0.171624\n",
       "                    3 low     0.250295\n",
       "                    4 no      0.463647\n",
       "gravity**(1/3)_300  1 high    0.060516\n",
       "                    2 medium  0.013619\n",
       "                    3 low     0.032726\n",
       "                    4 no      0.893139\n",
       "gravity**(1/3)_600  1 high    0.067402\n",
       "                    2 medium  0.044538\n",
       "                    3 low     0.120300\n",
       "                    4 no      0.767761\n",
       "gravity**(1/3)_1000 1 high    0.097625\n",
       "                    2 medium  0.114762\n",
       "                    3 low     0.220840\n",
       "                    4 no      0.566774\n",
       "gravity**(1/5)_300  1 high    0.060577\n",
       "                    2 medium  0.012341\n",
       "                    3 low     0.029410\n",
       "                    4 no      0.897672\n",
       "gravity**(1/5)_600  1 high    0.067514\n",
       "                    2 medium  0.040294\n",
       "                    3 low     0.090859\n",
       "                    4 no      0.801333\n",
       "gravity**(1/5)_1000 1 high    0.095977\n",
       "                    2 medium  0.080162\n",
       "                    3 low     0.192615\n",
       "                    4 no      0.631246"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Zibo'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Zibo', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bd476950",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['India']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/6c681418fb83e1a2bfd116d8664a4eaf-d39ade2fd44e70c79976ee87c8eeb3ae:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\IND_Indore_2020.tif\n",
      "get road networks from OSM\n",
      "Indore done 1.32 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Indore done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Indore 1.57 mns\n",
      "\n",
      "Indore\n",
      "0.0 % fake entry points done 0.0  mns\n",
      "19.3 % fake entry points done 0.39  mns\n",
      "38.5 % fake entry points done 0.76  mns\n",
      "57.8 % fake entry points done 1.12  mns\n",
      "77.1 % fake entry points done 1.49  mns\n",
      "96.3 % fake entry points done 1.87  mns\n",
      "Indore 100 % fake entry points done 1.94  mns\n",
      "Indore 100 % done 1.95  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "33.88 % 0.6 mns\n",
      "67.75 % 1.24 mns\n",
      "100 % finding combinations done\n",
      "Indore 280376 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Indore\n",
      "0.0 % done 0.9 mns\n",
      "33.88 % done 0.95 mns\n",
      "67.75 % done 1.51 mns\n",
      "100 % done 2.09 mns\n",
      "\n",
      "Indore\n",
      "0.0 % 0.09 mns\n",
      "10.16 % 0.72 mns\n",
      "20.33 % 1.38 mns\n",
      "30.49 % 2.01 mns\n",
      "40.65 % 2.59 mns\n",
      "50.81 % 3.15 mns\n",
      "60.98 % 3.71 mns\n",
      "71.14 % 4.3 mns\n",
      "81.3 % 4.96 mns\n",
      "91.46 % 5.62 mns\n",
      "100 % done 6.12 mns\n",
      "\n",
      "Indore\n",
      "entrance 0.1 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 2.75 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 6.08 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 8.83 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Indore done 11.28 mns\n",
      "CPU times: total: 19min 48s\n",
      "Wall time: 26min 17s\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>Indore</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.083574</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.034204</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.089766</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.792456</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.175276</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.081642</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.144365</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.598717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.305662</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.099871</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.117090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.477377</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.099466</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.061854</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.157258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.681422</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.267546</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.155744</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.144601</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.432109</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.476678</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.080923</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.062692</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.379707</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.088046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.041824</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.109602</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.760528</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.197517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.105679</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.176369</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.520436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.388061</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.118423</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.100937</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.392579</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.084744</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.036039</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.095337</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.783880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.181455</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.088498</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.152329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.577717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.332228</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.111418</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.132371</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.423983</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                            Indore\n",
       "entrance_300        1 high    0.083574\n",
       "                    2 medium  0.034204\n",
       "                    3 low     0.089766\n",
       "                    4 no      0.792456\n",
       "entrance_600        1 high    0.175276\n",
       "                    2 medium  0.081642\n",
       "                    3 low     0.144365\n",
       "                    4 no      0.598717\n",
       "entrance_1000       1 high    0.305662\n",
       "                    2 medium  0.099871\n",
       "                    3 low     0.117090\n",
       "                    4 no      0.477377\n",
       "gravity**(1/2)_300  1 high    0.099466\n",
       "                    2 medium  0.061854\n",
       "                    3 low     0.157258\n",
       "                    4 no      0.681422\n",
       "gravity**(1/2)_600  1 high    0.267546\n",
       "                    2 medium  0.155744\n",
       "                    3 low     0.144601\n",
       "                    4 no      0.432109\n",
       "gravity**(1/2)_1000 1 high    0.476678\n",
       "                    2 medium  0.080923\n",
       "                    3 low     0.062692\n",
       "                    4 no      0.379707\n",
       "gravity**(1/3)_300  1 high    0.088046\n",
       "                    2 medium  0.041824\n",
       "                    3 low     0.109602\n",
       "                    4 no      0.760528\n",
       "gravity**(1/3)_600  1 high    0.197517\n",
       "                    2 medium  0.105679\n",
       "                    3 low     0.176369\n",
       "                    4 no      0.520436\n",
       "gravity**(1/3)_1000 1 high    0.388061\n",
       "                    2 medium  0.118423\n",
       "                    3 low     0.100937\n",
       "                    4 no      0.392579\n",
       "gravity**(1/5)_300  1 high    0.084744\n",
       "                    2 medium  0.036039\n",
       "                    3 low     0.095337\n",
       "                    4 no      0.783880\n",
       "gravity**(1/5)_600  1 high    0.181455\n",
       "                    2 medium  0.088498\n",
       "                    3 low     0.152329\n",
       "                    4 no      0.577717\n",
       "gravity**(1/5)_1000 1 high    0.332228\n",
       "                    2 medium  0.111418\n",
       "                    3 low     0.132371\n",
       "                    4 no      0.423983"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Indore'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Indore', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "c7a11a85",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Guatemala']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/0e4df133f539bd506383073a7f23eee0-ab4e233e7091b25445264bad507663a0:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\GTM_Guatemala City_2020.tif\n",
      "get road networks from OSM\n",
      "Guatemala City done 1.25 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Guatemala City done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Guatemala City 1.88 mns\n",
      "\n",
      "Guatemala City\n",
      "0.0 % fake entry points done 0.0  mns\n",
      "11.5 % fake entry points done 0.43  mns\n",
      "23.0 % fake entry points done 0.84  mns\n",
      "34.5 % fake entry points done 1.27  mns\n",
      "46.0 % fake entry points done 1.72  mns\n",
      "57.5 % fake entry points done 2.18  mns\n",
      "69.0 % fake entry points done 2.66  mns\n",
      "80.6 % fake entry points done 3.07  mns\n",
      "92.1 % fake entry points done 3.5  mns\n",
      "Guatemala City 100 % fake entry points done 3.8  mns\n",
      "Guatemala City 100 % done 3.82  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "31.85 % 0.35 mns\n",
      "63.69 % 0.66 mns\n",
      "95.54 % 0.98 mns\n",
      "100 % finding combinations done\n",
      "Guatemala City 86042 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Guatemala City\n",
      "0.0 % done 0.82 mns\n",
      "31.85 % done 0.85 mns\n",
      "63.69 % done 1.4 mns\n",
      "95.54 % done 1.95 mns\n",
      "100 % done 2.46 mns\n",
      "\n",
      "Guatemala City\n",
      "0.0 % 0.09 mns\n",
      "53.94 % 0.67 mns\n",
      "100 % done 1.17 mns\n",
      "\n",
      "Guatemala City\n",
      "entrance 0.03 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 1.17 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 2.37 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 3.55 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Guatemala City done 4.7 mns\n",
      "CPU times: total: 13min 19s\n",
      "Wall time: 16min 24s\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>Guatemala City</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.115183</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.003307</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.026636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.854874</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.126035</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.031104</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.070580</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.772281</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.149469</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.063332</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.068504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.718695</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.114965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.008013</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.044122</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.832900</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.127667</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.051648</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.080059</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.740626</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.152504</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.094569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.046102</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.706825</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.114848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.004513</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.029720</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.850919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.123450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.036877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.087718</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.751955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.149661</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.080530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.062399</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.707410</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.114994</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.003096</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.025052</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.856858</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.122476</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.031353</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.081530</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.764640</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.148229</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.070518</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.073420</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.707834</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                          Guatemala City\n",
       "entrance_300        1 high          0.115183\n",
       "                    2 medium        0.003307\n",
       "                    3 low           0.026636\n",
       "                    4 no            0.854874\n",
       "entrance_600        1 high          0.126035\n",
       "                    2 medium        0.031104\n",
       "                    3 low           0.070580\n",
       "                    4 no            0.772281\n",
       "entrance_1000       1 high          0.149469\n",
       "                    2 medium        0.063332\n",
       "                    3 low           0.068504\n",
       "                    4 no            0.718695\n",
       "gravity**(1/2)_300  1 high          0.114965\n",
       "                    2 medium        0.008013\n",
       "                    3 low           0.044122\n",
       "                    4 no            0.832900\n",
       "gravity**(1/2)_600  1 high          0.127667\n",
       "                    2 medium        0.051648\n",
       "                    3 low           0.080059\n",
       "                    4 no            0.740626\n",
       "gravity**(1/2)_1000 1 high          0.152504\n",
       "                    2 medium        0.094569\n",
       "                    3 low           0.046102\n",
       "                    4 no            0.706825\n",
       "gravity**(1/3)_300  1 high          0.114848\n",
       "                    2 medium        0.004513\n",
       "                    3 low           0.029720\n",
       "                    4 no            0.850919\n",
       "gravity**(1/3)_600  1 high          0.123450\n",
       "                    2 medium        0.036877\n",
       "                    3 low           0.087718\n",
       "                    4 no            0.751955\n",
       "gravity**(1/3)_1000 1 high          0.149661\n",
       "                    2 medium        0.080530\n",
       "                    3 low           0.062399\n",
       "                    4 no            0.707410\n",
       "gravity**(1/5)_300  1 high          0.114994\n",
       "                    2 medium        0.003096\n",
       "                    3 low           0.025052\n",
       "                    4 no            0.856858\n",
       "gravity**(1/5)_600  1 high          0.122476\n",
       "                    2 medium        0.031353\n",
       "                    3 low           0.081530\n",
       "                    4 no            0.764640\n",
       "gravity**(1/5)_1000 1 high          0.148229\n",
       "                    2 medium        0.070518\n",
       "                    3 low           0.073420\n",
       "                    4 no            0.707834"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Guatemala City'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Guatemala City', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "053c71e1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Nigeria']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/9dd5fa1c7e53874dedb5b233adbd80ae-a2c6dceb08a2300a65a81a799020bb10:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\NGA_Lagos_2020.tif\n",
      "get road networks from OSM\n",
      "Lagos done 3.41 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Lagos done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Lagos 3.94 mns\n",
      "\n",
      "Lagos\n",
      "0.0 % fake entry points done 0.0  mns\n",
      "90.9 % fake entry points done 0.93  mns\n",
      "Lagos 100 % fake entry points done 1.07  mns\n",
      "Lagos 100 % done 1.07  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "91.24 % 0.98 mns\n",
      "100 % finding combinations done\n",
      "Lagos 281939 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Lagos\n",
      "0.0 % done 0.55 mns\n",
      "91.24 % done 0.6 mns\n",
      "100 % done 1.11 mns\n",
      "\n",
      "Lagos\n",
      "0.0 % 0.07 mns\n",
      "14.55 % 0.72 mns\n",
      "29.1 % 1.35 mns\n",
      "43.66 % 1.99 mns\n",
      "58.21 % 2.57 mns\n",
      "72.76 % 3.16 mns\n",
      "87.31 % 3.75 mns\n",
      "100 % done 4.25 mns\n",
      "\n",
      "Lagos\n",
      "entrance 0.06 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 5.58 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 11.07 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 16.43 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Lagos done 22.35 mns\n",
      "CPU times: total: 33min 4s\n",
      "Wall time: 37min 20s\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>Lagos</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.019116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.001731</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.005822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.973331</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019403</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.007535</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.022363</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.950699</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.021123</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.020488</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.055212</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.903177</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.001577</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.004368</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.974944</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019151</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.005975</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.015787</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.959087</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019395</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.015779</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.033207</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.931619</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.001370</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.004075</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.975454</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019145</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.005436</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.014338</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.961080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.013951</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.033367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.933296</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019097</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.001273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.004153</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.975477</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019165</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.005403</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.015215</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.960217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.019473</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.014473</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.036325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.929729</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                             Lagos\n",
       "entrance_300        1 high    0.019116\n",
       "                    2 medium  0.001731\n",
       "                    3 low     0.005822\n",
       "                    4 no      0.973331\n",
       "entrance_600        1 high    0.019403\n",
       "                    2 medium  0.007535\n",
       "                    3 low     0.022363\n",
       "                    4 no      0.950699\n",
       "entrance_1000       1 high    0.021123\n",
       "                    2 medium  0.020488\n",
       "                    3 low     0.055212\n",
       "                    4 no      0.903177\n",
       "gravity**(1/2)_300  1 high    0.019110\n",
       "                    2 medium  0.001577\n",
       "                    3 low     0.004368\n",
       "                    4 no      0.974944\n",
       "gravity**(1/2)_600  1 high    0.019151\n",
       "                    2 medium  0.005975\n",
       "                    3 low     0.015787\n",
       "                    4 no      0.959087\n",
       "gravity**(1/2)_1000 1 high    0.019395\n",
       "                    2 medium  0.015779\n",
       "                    3 low     0.033207\n",
       "                    4 no      0.931619\n",
       "gravity**(1/3)_300  1 high    0.019101\n",
       "                    2 medium  0.001370\n",
       "                    3 low     0.004075\n",
       "                    4 no      0.975454\n",
       "gravity**(1/3)_600  1 high    0.019145\n",
       "                    2 medium  0.005436\n",
       "                    3 low     0.014338\n",
       "                    4 no      0.961080\n",
       "gravity**(1/3)_1000 1 high    0.019386\n",
       "                    2 medium  0.013951\n",
       "                    3 low     0.033367\n",
       "                    4 no      0.933296\n",
       "gravity**(1/5)_300  1 high    0.019097\n",
       "                    2 medium  0.001273\n",
       "                    3 low     0.004153\n",
       "                    4 no      0.975477\n",
       "gravity**(1/5)_600  1 high    0.019165\n",
       "                    2 medium  0.005403\n",
       "                    3 low     0.015215\n",
       "                    4 no      0.960217\n",
       "gravity**(1/5)_1000 1 high    0.019473\n",
       "                    2 medium  0.014473\n",
       "                    3 low     0.036325\n",
       "                    4 no      0.929729"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Lagos'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Lagos', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "c2a18ebf",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[\"Côte d'Ivoire\"]\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/133236acd52df493cbb311880356fb52-9d9f1f1d949ea30476fc728970d3891a:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\CIV_Abidjan_2020.tif\n",
      "get road networks from OSM\n",
      "Abidjan done 2.32 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Abidjan done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Abidjan 2.97 mns\n",
      "\n",
      "Abidjan\n",
      "0.0 % fake entry points done 0.01  mns\n",
      "32.4 % fake entry points done 0.58  mns\n",
      "64.7 % fake entry points done 1.2  mns\n",
      "97.1 % fake entry points done 1.92  mns\n",
      "Abidjan 100 % fake entry points done 2.0  mns\n",
      "Abidjan 100 % done 2.03  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "36.66 % 0.58 mns\n",
      "73.31 % 1.22 mns\n",
      "100 % finding combinations done\n",
      "Abidjan 1004661 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Abidjan\n",
      "0.0 % done 0.9 mns\n",
      "36.66 % done 1.0 mns\n",
      "73.31 % done 1.67 mns\n",
      "100 % done 2.3 mns\n",
      "\n",
      "Abidjan\n",
      "0.0 % 0.15 mns\n",
      "4.7 % 0.78 mns\n",
      "9.39 % 1.4 mns\n",
      "14.09 % 2.02 mns\n",
      "18.78 % 2.76 mns\n",
      "23.48 % 3.46 mns\n",
      "28.17 % 4.28 mns\n",
      "32.87 % 4.9 mns\n",
      "37.56 % 5.67 mns\n",
      "42.26 % 6.28 mns\n",
      "46.96 % 6.97 mns\n",
      "51.65 % 7.58 mns\n",
      "56.35 % 8.23 mns\n",
      "61.04 % 9.23 mns\n",
      "65.74 % 10.29 mns\n",
      "70.43 % 11.42 mns\n",
      "75.13 % 14.15 mns\n",
      "79.82 % 14.88 mns\n",
      "84.52 % 15.56 mns\n",
      "89.21 % 16.32 mns\n",
      "93.91 % 17.16 mns\n",
      "98.61 % 17.7 mns\n",
      "100 % done 18.24 mns\n",
      "\n",
      "Abidjan\n",
      "entrance 0.28 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 7.64 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 11.18 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 13.25 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Abidjan done 15.1 mns\n",
      "CPU times: total: 23min 41s\n",
      "Wall time: 53min 21s\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>Abidjan</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.072288</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.014692</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.043270</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.869749</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.082141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.056726</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.132901</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.728232</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.126597</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.132357</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.237717</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.503329</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.071434</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.014058</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.042598</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.871910</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.075725</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.056474</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.144057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.723744</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.098300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.152791</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.281451</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.467458</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.071418</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.011184</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.033816</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.883582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.074609</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.045042</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.113934</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.766414</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.092339</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.119038</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.225830</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.562792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.071435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.010711</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.033386</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.884468</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.075771</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.043119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.108906</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.772204</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.097217</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.111303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.216417</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.575063</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                           Abidjan\n",
       "entrance_300        1 high    0.072288\n",
       "                    2 medium  0.014692\n",
       "                    3 low     0.043270\n",
       "                    4 no      0.869749\n",
       "entrance_600        1 high    0.082141\n",
       "                    2 medium  0.056726\n",
       "                    3 low     0.132901\n",
       "                    4 no      0.728232\n",
       "entrance_1000       1 high    0.126597\n",
       "                    2 medium  0.132357\n",
       "                    3 low     0.237717\n",
       "                    4 no      0.503329\n",
       "gravity**(1/2)_300  1 high    0.071434\n",
       "                    2 medium  0.014058\n",
       "                    3 low     0.042598\n",
       "                    4 no      0.871910\n",
       "gravity**(1/2)_600  1 high    0.075725\n",
       "                    2 medium  0.056474\n",
       "                    3 low     0.144057\n",
       "                    4 no      0.723744\n",
       "gravity**(1/2)_1000 1 high    0.098300\n",
       "                    2 medium  0.152791\n",
       "                    3 low     0.281451\n",
       "                    4 no      0.467458\n",
       "gravity**(1/3)_300  1 high    0.071418\n",
       "                    2 medium  0.011184\n",
       "                    3 low     0.033816\n",
       "                    4 no      0.883582\n",
       "gravity**(1/3)_600  1 high    0.074609\n",
       "                    2 medium  0.045042\n",
       "                    3 low     0.113934\n",
       "                    4 no      0.766414\n",
       "gravity**(1/3)_1000 1 high    0.092339\n",
       "                    2 medium  0.119038\n",
       "                    3 low     0.225830\n",
       "                    4 no      0.562792\n",
       "gravity**(1/5)_300  1 high    0.071435\n",
       "                    2 medium  0.010711\n",
       "                    3 low     0.033386\n",
       "                    4 no      0.884468\n",
       "gravity**(1/5)_600  1 high    0.075771\n",
       "                    2 medium  0.043119\n",
       "                    3 low     0.108906\n",
       "                    4 no      0.772204\n",
       "gravity**(1/5)_1000 1 high    0.097217\n",
       "                    2 medium  0.111303\n",
       "                    3 low     0.216417\n",
       "                    4 no      0.575063"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Abidjan'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Abidjan', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bfb90317",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['South Africa']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/8b604efbbcc47906bd0ec1f825dd5a07-94db0194128512b4409e41d829b10f5c:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\ZAF_Johannesburg_2020.tif\n",
      "get road networks from OSM\n",
      "Johannesburg done 4.98 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Johannesburg done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Johannesburg 9.64 mns\n",
      "\n",
      "Johannesburg\n",
      "0.0 % fake entry points done 0.02  mns\n",
      "15.1 % fake entry points done 1.64  mns\n",
      "30.3 % fake entry points done 3.12  mns\n",
      "45.4 % fake entry points done 4.59  mns\n",
      "60.5 % fake entry points done 6.05  mns\n",
      "75.6 % fake entry points done 7.61  mns\n",
      "90.8 % fake entry points done 9.16  mns\n",
      "Johannesburg 100 % fake entry points done 10.14  mns\n",
      "Johannesburg 100 % done 10.18  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "11.62 % 1.37 mns\n",
      "23.25 % 2.78 mns\n",
      "34.87 % 4.27 mns\n",
      "46.49 % 5.79 mns\n",
      "58.11 % 7.35 mns\n",
      "69.74 % 8.93 mns\n",
      "81.36 % 10.57 mns\n",
      "92.98 % 12.21 mns\n",
      "100 % finding combinations done\n",
      "Johannesburg 2744633 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Johannesburg\n",
      "0.0 % done 8.99 mns\n",
      "11.62 % done 9.07 mns\n",
      "23.25 % done 9.69 mns\n",
      "34.87 % done 10.35 mns\n",
      "46.49 % done 11.03 mns\n",
      "58.11 % done 11.75 mns\n",
      "69.74 % done 12.52 mns\n",
      "81.36 % done 13.31 mns\n",
      "92.98 % done 14.11 mns\n",
      "100 % done 14.81 mns\n",
      "\n",
      "Johannesburg\n",
      "0.0 % 0.26 mns\n",
      "1.54 % 0.86 mns\n",
      "3.07 % 1.44 mns\n",
      "4.61 % 2.07 mns\n",
      "6.15 % 2.68 mns\n",
      "7.68 % 3.36 mns\n",
      "9.22 % 4.07 mns\n",
      "10.76 % 4.77 mns\n",
      "12.29 % 5.36 mns\n",
      "13.83 % 6.07 mns\n",
      "15.37 % 6.76 mns\n",
      "16.9 % 7.46 mns\n",
      "18.44 % 8.06 mns\n",
      "19.98 % 8.64 mns\n",
      "21.51 % 9.22 mns\n",
      "23.05 % 9.77 mns\n",
      "24.59 % 10.41 mns\n",
      "26.12 % 10.99 mns\n",
      "27.66 % 11.55 mns\n",
      "29.2 % 12.13 mns\n",
      "30.73 % 12.69 mns\n",
      "32.27 % 13.31 mns\n",
      "33.81 % 13.91 mns\n",
      "35.34 % 14.52 mns\n",
      "36.88 % 15.15 mns\n",
      "38.42 % 15.8 mns\n",
      "39.95 % 16.39 mns\n",
      "41.49 % 16.98 mns\n",
      "43.03 % 17.57 mns\n",
      "44.56 % 18.17 mns\n",
      "46.1 % 18.79 mns\n",
      "47.64 % 19.43 mns\n",
      "49.17 % 20.06 mns\n",
      "50.71 % 20.77 mns\n",
      "52.25 % 21.35 mns\n",
      "53.78 % 21.98 mns\n",
      "55.32 % 22.57 mns\n",
      "56.86 % 23.2 mns\n",
      "58.39 % 23.83 mns\n",
      "59.93 % 24.45 mns\n",
      "61.47 % 25.1 mns\n",
      "63.0 % 25.76 mns\n",
      "64.54 % 26.39 mns\n",
      "66.08 % 27.06 mns\n",
      "67.61 % 27.71 mns\n",
      "69.15 % 28.4 mns\n",
      "70.69 % 29.05 mns\n",
      "72.22 % 29.66 mns\n",
      "73.76 % 30.27 mns\n",
      "75.3 % 30.92 mns\n",
      "76.83 % 31.54 mns\n",
      "78.37 % 32.13 mns\n",
      "79.91 % 32.72 mns\n",
      "81.44 % 33.35 mns\n",
      "82.98 % 33.96 mns\n",
      "84.52 % 34.58 mns\n",
      "86.05 % 35.2 mns\n",
      "87.59 % 35.83 mns\n",
      "89.13 % 36.44 mns\n",
      "90.66 % 37.08 mns\n",
      "92.2 % 37.73 mns\n",
      "93.74 % 38.32 mns\n",
      "95.27 % 38.96 mns\n",
      "96.81 % 39.6 mns\n",
      "98.35 % 40.28 mns\n",
      "99.88 % 40.8 mns\n",
      "100 % done 41.35 mns\n",
      "\n",
      "Johannesburg\n",
      "entrance 0.33 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 6.18 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 12.88 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 22.86 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Johannesburg done 29.13 mns\n",
      "CPU times: total: 1h 22min 43s\n",
      "Wall time: 2h 3min 34s\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>Johannesburg</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.023651</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.011473</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.029957</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.934919</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.033872</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.039216</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.076822</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.850090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.070556</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.071597</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.118249</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.739598</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.023450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.010536</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.028201</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.937813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.026192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.039883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.075343</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.858582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.052881</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.077931</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.136987</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.732201</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.023435</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.009185</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.025046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.942334</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.026190</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.035334</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.068071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.870405</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.052249</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.067883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.113313</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.766554</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.023475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.009082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.024925</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.942517</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.027660</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.034707</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.068490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.869142</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.055687</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.065740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.108727</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.769845</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                          Johannesburg\n",
       "entrance_300        1 high        0.023651\n",
       "                    2 medium      0.011473\n",
       "                    3 low         0.029957\n",
       "                    4 no          0.934919\n",
       "entrance_600        1 high        0.033872\n",
       "                    2 medium      0.039216\n",
       "                    3 low         0.076822\n",
       "                    4 no          0.850090\n",
       "entrance_1000       1 high        0.070556\n",
       "                    2 medium      0.071597\n",
       "                    3 low         0.118249\n",
       "                    4 no          0.739598\n",
       "gravity**(1/2)_300  1 high        0.023450\n",
       "                    2 medium      0.010536\n",
       "                    3 low         0.028201\n",
       "                    4 no          0.937813\n",
       "gravity**(1/2)_600  1 high        0.026192\n",
       "                    2 medium      0.039883\n",
       "                    3 low         0.075343\n",
       "                    4 no          0.858582\n",
       "gravity**(1/2)_1000 1 high        0.052881\n",
       "                    2 medium      0.077931\n",
       "                    3 low         0.136987\n",
       "                    4 no          0.732201\n",
       "gravity**(1/3)_300  1 high        0.023435\n",
       "                    2 medium      0.009185\n",
       "                    3 low         0.025046\n",
       "                    4 no          0.942334\n",
       "gravity**(1/3)_600  1 high        0.026190\n",
       "                    2 medium      0.035334\n",
       "                    3 low         0.068071\n",
       "                    4 no          0.870405\n",
       "gravity**(1/3)_1000 1 high        0.052249\n",
       "                    2 medium      0.067883\n",
       "                    3 low         0.113313\n",
       "                    4 no          0.766554\n",
       "gravity**(1/5)_300  1 high        0.023475\n",
       "                    2 medium      0.009082\n",
       "                    3 low         0.024925\n",
       "                    4 no          0.942517\n",
       "gravity**(1/5)_600  1 high        0.027660\n",
       "                    2 medium      0.034707\n",
       "                    3 low         0.068490\n",
       "                    4 no          0.869142\n",
       "gravity**(1/5)_1000 1 high        0.055687\n",
       "                    2 medium      0.065740\n",
       "                    3 low         0.108727\n",
       "                    4 no          0.769845"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Johannesburg'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Johannesburg', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ada3d4d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['India']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/1260c38ac9a886bd9baa59c0c638f489-898dbdb56af955c30b12aa5c6e5e9c59:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\IND_Bangalore_2020.tif\n",
      "get road networks from OSM\n",
      "Bangalore done 5.07 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Bangalore done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Bangalore 10.12 mns\n",
      "\n",
      "Bangalore\n",
      "0.0 % fake entry points done 0.02  mns\n",
      "6.4 % fake entry points done 1.89  mns\n",
      "12.7 % fake entry points done 3.79  mns\n",
      "19.1 % fake entry points done 5.69  mns\n",
      "25.4 % fake entry points done 7.56  mns\n",
      "31.8 % fake entry points done 9.53  mns\n",
      "38.2 % fake entry points done 11.51  mns\n",
      "44.5 % fake entry points done 13.3  mns\n",
      "50.9 % fake entry points done 15.2  mns\n",
      "57.3 % fake entry points done 17.11  mns\n",
      "63.6 % fake entry points done 19.12  mns\n",
      "70.0 % fake entry points done 21.32  mns\n",
      "76.3 % fake entry points done 23.27  mns\n",
      "82.7 % fake entry points done 25.2  mns\n",
      "89.1 % fake entry points done 27.02  mns\n",
      "95.4 % fake entry points done 28.79  mns\n",
      "Bangalore 100 % fake entry points done 30.46  mns\n",
      "Bangalore 100 % done 31.36  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "12.37 % 1.79 mns\n",
      "24.74 % 3.01 mns\n",
      "37.11 % 4.23 mns\n",
      "49.48 % 5.42 mns\n",
      "61.85 % 6.62 mns\n",
      "74.22 % 7.84 mns\n",
      "86.59 % 9.07 mns\n",
      "98.96 % 10.3 mns\n",
      "100 % finding combinations done\n",
      "Bangalore 364413 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Bangalore\n",
      "0.0 % done 9.83 mns\n",
      "12.37 % done 13.57 mns\n",
      "24.74 % done 14.43 mns\n",
      "37.11 % done 15.56 mns\n",
      "49.48 % done 16.5 mns\n",
      "61.85 % done 17.2 mns\n",
      "74.22 % done 18.04 mns\n",
      "86.59 % done 18.77 mns\n",
      "98.96 % done 19.59 mns\n",
      "100 % done 20.11 mns\n",
      "\n",
      "Bangalore\n",
      "0.0 % 0.81 mns\n",
      "8.44 % 1.44 mns\n",
      "16.88 % 2.05 mns\n",
      "25.32 % 2.66 mns\n",
      "33.76 % 3.3 mns\n",
      "42.2 % 3.93 mns\n",
      "50.64 % 4.58 mns\n",
      "59.08 % 5.22 mns\n",
      "67.52 % 5.85 mns\n",
      "75.97 % 6.47 mns\n",
      "84.41 % 7.11 mns\n",
      "92.85 % 7.72 mns\n",
      "100 % done 8.23 mns\n",
      "\n",
      "Bangalore\n",
      "entrance 0.14 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 4.88 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 10.35 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 16.12 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\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(['Bangalore'])]\n",
    "cities_adj = cities_adj.reset_index()\n",
    "\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,'D:/Dumps/GEE_city_grids/')\n",
    "\n",
    "# Files are downloaded automatically to the specified path. Files are also stored in Google with a downloadlink:\n",
    "\n",
    "# 2. Information extraction\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_network(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\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)\n",
    "\n",
    "print(' ')\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     road_networks['nodes'],\n",
    "                                     UGS,\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print('')\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           road_networks['graphs'],\n",
    "                           cities_adj['City'],\n",
    "                           population_grids,\n",
    "                           thresholds,\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",
    "suitible_enh = suitible_enhanced(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], \n",
    "                                 cities_adj['City'], \n",
    "                                 thresholds)\n",
    "print('')\n",
    "subgraphs = obtaining_subgraphs(road_networks['graphs'],\n",
    "                                population_grids,\n",
    "                                UGS_entry,\n",
    "                                road_networks['nodes'],\n",
    "                                cities_adj['City'],\n",
    "                                thresholds)\n",
    "print('')\n",
    "Dir_Routes = direct_routing (suitible_enh,\n",
    "                             subgraphs['graphs'],\n",
    "                             road_networks['edges'],\n",
    "                             cities_adj['City'],\n",
    "                            time_sleep = 30)\n",
    "\n",
    "print('')\n",
    "grid_scores = grid_score_summary (Dir_Routes, # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '_Bangalore', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'D:/Dumps/GEE-WP Scores/Gravity_adj/',\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c0bde4ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gee_worldpop_extract (city_file, iso, save_path = None):\n",
    "    \n",
    "    cities = city_file\n",
    "    iso['name'] = np.where(iso['name'] == 'Macedonia','North Macedonia',iso['name'])\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",
    "    print(obj_displ)\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()\n",
    "    \n",
    "    # Block 2 Road networks\n",
    "def road_network (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})\n",
    "# Block 3 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)\n",
    "# Block 4 population grids extraction\n",
    "def city_grids_format(city_grids, cities_area, road_nodes, UGS, 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",
    "        try:\n",
    "            entry_index = [int(road_nodes[i]['geometry'].sindex.nearest(grid)[1])\\\n",
    "                                 for grid in popgrid['centroid']]\n",
    "        except:\n",
    "            entry_index = [int(road_nodes[i]['geometry'].sindex.nearest(grid)[1][0])\\\n",
    "                                 for grid in popgrid['centroid']]\n",
    "            \n",
    "        nearest_index = road_nodes[i].iloc[entry_index]\n",
    "        popgrid['grid_osm'] = nearest_index.reset_index(drop = True)['osmid_var']\n",
    "        popgrid['node_geom'] = nearest_index.reset_index(drop = True)['geometry']\n",
    "        popgrid['node_geom_m'] = nearest_index.reset_index(drop = True)['geometry_m']\n",
    "        popgrid['G-entry cost'] = popgrid['node_geom_m'].distance(popgrid['centroid_m'])\n",
    "        \n",
    "        UGS_all = UGS[i].dissolve().geometry[0]\n",
    "        popgrid['in_out_UGS'] = popgrid.intersection(UGS_all).is_empty == False\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": 4,
   "id": "59bebc3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5 park entry points\n",
    "def UGS_fake_entry(UGS, road_nodes, graphs, cities, pop_grids,\n",
    "                   thresholds, UGS_entry_buf = 25, walk_radius = 500, entry_point_merge = 0):\n",
    "    \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",
    "        print(cities[j].rsplit(',')[0])\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(round(i/len(UGS[j])*100,1),'% fake entry points done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "                \n",
    "        print(cities[j].rsplit(',')[0],'100 % fake entry points done', round((time.time() - start_time) / 60,2),' mns')\n",
    "        \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",
    "        ParkRoad['raw'] = 1\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",
    "        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)\n",
    "# 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": 5,
   "id": "371c0f0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def suitible_enhanced (UGS_entry, pop_grids, road_nodes, cities, thresholds):\n",
    "    start_time = time.time()\n",
    "    suits_all = []\n",
    "    for j in range(len(cities)):\n",
    "        print('get (Euclidean) suitible combinations')\n",
    "        print('0.0 %', round((time.time() - start_time) / 60,2),'mns')\n",
    "        UGSe = UGS_entry[j]\n",
    "        entry_geoms = UGSe.geometry_m\n",
    "        pop = pop_grids[j]\n",
    "        road_node = road_nodes[j]\n",
    "\n",
    "        suits = pd.DataFrame()\n",
    "        cols = ['osmid','Park_No','walk_area','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5']\n",
    "        for i in range(len(entry_geoms)):\n",
    "            max_infl = np.max(UGSe[['raw','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5']], axis = 1)[i]\n",
    "            suit_df = pop[pop.node_geom_m.distance(entry_geoms.iloc[i]) < (max_infl*np.max(thresholds))]\n",
    "        \n",
    "            suit_df['UGSe_osmid_m'] = entry_geoms.iloc[i]\n",
    "            suit_df['Grid_No'] = suit_df.index\n",
    "            suit_df = suit_df[['Grid_No','grid_osm','G-entry cost','in_out_UGS','node_geom_m','UGSe_osmid_m']].reset_index(drop = True)\n",
    "            suit_df['Park_entry_No'] = UGSe.index[i]\n",
    "            #suit_df = pd.merge(suit_df, UGSe[cols], left_on = 'Park_entry_No',right_index = True, how = 'left')\n",
    "            suits = pd.concat([suits,suit_df])\n",
    "            if (i+1) % 500 == 0: print(round((i+1) / len(entry_geoms)*100,2),'%',\n",
    "                                       round((time.time() - start_time) / 60,2),'mns')\n",
    "            \n",
    "        suits = pd.merge(suits, UGSe[cols], left_on = 'Park_entry_No',right_index = True, how = 'left')\n",
    "        suits = suits.reset_index(drop = True)\n",
    "        suits = suits.rename(columns = {'osmid':'Parkroad_osmid','walk_area':'walk_area_m2'})\n",
    "        suits['gridpark_no'] = suits['Grid_No'].astype(str)+'-'+suits['Park_No'].astype(str)\n",
    "        suits['graph_key'] = suits['grid_osm'].astype(str)+'-'+suits['Parkroad_osmid'].astype(str)\n",
    "        suits_all.append(suits)\n",
    "        print('100 % finding combinations done')\n",
    "        print(cities[j],len(suits),'suitible combinations')\n",
    "    return(suits_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6553bf33",
   "metadata": {},
   "outputs": [],
   "source": [
    "def obtaining_subgraphs(graphs, pop_grids, UGS_entry, nodes, cities, thresholds, time_sleep = 30):\n",
    "    print('obtain local graphs')\n",
    "    start_time = time.time()\n",
    "    subgraphs_all = []\n",
    "    suits_all = []\n",
    "    for j in range(len(cities)):\n",
    "        print(cities[j])\n",
    "        Graph = graphs[j]\n",
    "        pop = pop_grids[j]\n",
    "        UGSe = UGS_entry[j].sort_values('osmid')\n",
    "        road_node = nodes[j]\n",
    "        node_geoms = road_node.geometry_m\n",
    "        entry_geoms = UGSe.geometry_m\n",
    "        osmid = UGSe['osmid']\n",
    "        max_infl = np.max(UGSe[['raw','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5']], axis = 1)*(np.max(thresholds))\n",
    "\n",
    "        dist = [node_geoms.distance(Point(i)) for i in entry_geoms]\n",
    "\n",
    "        print('0.0 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        subgraphs = []\n",
    "        UGSe_ids = []\n",
    "        suits = pd.DataFrame()\n",
    "        for i in range(len(entry_geoms)):      \n",
    "            suit = road_node[['geometry_m']]\n",
    "            suit['UGSe_osmid_m'] = entry_geoms.iloc[i]\n",
    "            suit_df = dist[i]\n",
    "            suit_in = suit_df[suit_df <= max_infl.iloc[i]]\n",
    "            UGSe_ids.append(osmid.iloc[i])\n",
    "            suit_in = pd.DataFrame(suit_in).join(node_geoms)\n",
    "            suit_in['Parkroad_osmid'] = osmid.iloc[i]\n",
    "            subgraphs.append(Graph.subgraph(suit_in.index))\n",
    "            suits = pd.concat([suits, suit_in])\n",
    "\n",
    "            if (i+1) % 500 == 0: \n",
    "                print(round((i+1) / len(entry_geoms)*100,2),'% done',\n",
    "                                        round((time.time() - start_time) / 60,2),'mns')\n",
    "                time.sleep(time_sleep)\n",
    "        print('100 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        subgraphs_all.append(pd.Series(subgraphs, index = UGSe_ids))\n",
    "        suits_all.append(suits)\n",
    "    return({'graphs':subgraphs_all,'graph nodes':suits_all})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "aeab83a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_fast (Geo_1, Geo_2):\n",
    "    return((abs(Geo_1.x - Geo_2.x)**2 + abs(Geo_1.y - Geo_2.y)**2).apply(math.sqrt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c670d949",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_routing (suitible_comb, graphs, edges, cities, chunk = 20000, time_sleep = 15):\n",
    "    start_time = time.time()\n",
    "    Routes = []\n",
    "    Lines = []\n",
    "    for j in enumerate(cities):\n",
    "        print(j[1])\n",
    "        \n",
    "        suitible = suitible_comb[j[0]].sort_values('Parkroad_osmid').reset_index()\n",
    "        grouped = suitible[suitible['in_out_UGS'] == False].groupby(['Parkroad_osmid'])['grid_osm'].apply(list)\n",
    "        sets = grouped.apply(np.unique)\n",
    "\n",
    "        Conn = edges[j[0]]\n",
    "        SG = graphs[j[0]]\n",
    "        \n",
    "        SGr = SG.reset_index()\n",
    "        SG = SGr.iloc[pd.Series(SGr['index'].drop_duplicates()).index].set_index('index')[0]\n",
    "\n",
    "        num = int(np.ceil(chunk / sets.apply(len).mean()))\n",
    "        length = int(np.ceil(len(suitible['Parkroad_osmid'].unique())/num))\n",
    "\n",
    "        Routes_df = pd.DataFrame()\n",
    "        Lines_df = pd.DataFrame()\n",
    "        for l in range(length):\n",
    "            comb = suitible[suitible['Parkroad_osmid'].isin(sets.index[l*num:l*num+num])]\n",
    "            sets2 = sets[l*num:l*num+num]\n",
    "\n",
    "            parknode = list(comb['Parkroad_osmid'])\n",
    "            gridnode = list(comb['grid_osm'])\n",
    "            subgraph = SG[sets2.index]\n",
    "\n",
    "            ls = []\n",
    "            ls2 = []\n",
    "            ls3 = []\n",
    "            lod = []\n",
    "            lgk = []\n",
    "            Routes\n",
    "            for i in range(len(sets2)):\n",
    "                path = nx.single_source_dijkstra(subgraph.iloc[i], sets2.index[i], weight = 'length')\n",
    "\n",
    "                incl = np.isin(list(path[0].keys()),sets2.iloc[i])\n",
    "                incl2 = np.isin(list(path[1].keys()),sets2.iloc[i])\n",
    "\n",
    "                # route cost\n",
    "                orig_c = list(np.repeat(sets2.index[i],sum(incl)))\n",
    "                dest_c = list(np.array(list(path[0].keys()))[incl])\n",
    "                cost = list(np.array(list(path[0].values()))[incl])\n",
    "\n",
    "                ls = ls + orig_c\n",
    "                ls2= ls2+ dest_c\n",
    "                ls3= ls3+ cost\n",
    "\n",
    "                # route steps\n",
    "                orig_s = list(np.repeat(sets2.index[i],sum(incl2)))\n",
    "                dest_s = list(np.array(list(path[1].keys()))[incl2])\n",
    "                steps = list(np.array(list(path[1].values()),dtype=object)[incl2])\n",
    "\n",
    "                fr = []\n",
    "                to = []\n",
    "                og = []\n",
    "                de = []\n",
    "                for j in enumerate(steps):\n",
    "                    if len(j[1]) > 1:\n",
    "                        fr.append(j[1][:-1])\n",
    "                        to.append(j[1][1:])\n",
    "                        og.append(list(np.repeat(orig_s[j[0]], len(j[1][:-1]))))\n",
    "                        de.append(list(np.repeat(dest_s[j[0]], len(j[1][:-1]))))\n",
    "                    else:\n",
    "                        pass\n",
    "\n",
    "                fr = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, fr) for i in b]\n",
    "                to = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, to) for i in b]\n",
    "                og = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, og) for i in b]\n",
    "                de = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, de) for i in b]\n",
    "\n",
    "                gk = [str(fr[k])+'-'+str(to[k]) for k in range(len(to))]\n",
    "                gkr = [str(to[k])+'-'+str(fr[k]) for k in range(len(to))]\n",
    "                od = [str(de[k])+'-'+str(og[k]) for k in range(len(og))]\n",
    "\n",
    "                lgk.append(gk)\n",
    "                lod.append(od)\n",
    "\n",
    "            dist_df = pd.DataFrame({'UGSe_id':ls,'GrE_id':ls2,'route cost':ls3})\n",
    "            dist_df['graph_key'] = dist_df['GrE_id'].astype(str)+'-'+dist_df['UGSe_id'].astype(str)\n",
    "\n",
    "            routes = pd.merge(comb, dist_df, on = 'graph_key', how = 'left')\n",
    "            routes['route cost'] = np.where(routes['in_out_UGS'],0,routes['route cost'])\n",
    "            routes = routes[~routes['route cost'].isna()].reset_index(drop = True)\n",
    "\n",
    "            routes['G-entry cost'] = np.where(routes['in_out_UGS'],0,routes['G-entry cost'])\n",
    "\n",
    "            routes['raw_Tcost'] = routes['route cost']+routes['G-entry cost']\n",
    "            routes['grav2_Tcost'] = routes['raw_Tcost'] / routes['size_infl_sqr2']\n",
    "            routes['grav3_Tcost'] = routes['raw_Tcost'] / routes['size_infl_sqr3']\n",
    "            routes['grav5_Tcost'] = routes['raw_Tcost'] / routes['size_infl_sqr5']\n",
    "\n",
    "            lgk = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, lgk) for i in b]\n",
    "            lod = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, lod) for i in b]\n",
    "            \n",
    "            linestr = pd.DataFrame({'route no':lod,'route step':lgk})\n",
    "            \n",
    "            linestr = pd.merge(linestr, Conn.geometry, left_on = 'route step', right_index = True, how = 'left')\n",
    "            linestr = linestr[['route no','geometry']]\n",
    "            linestr = gpd.GeoDataFrame(linestr[['route no','geometry']], crs = 4326)\n",
    "            \n",
    "            linestr = linestr.dissolve('route no')\n",
    "            routes2 = pd.merge(routes, linestr, left_on = 'graph_key', right_index = True, how = 'left')\n",
    "            \n",
    "            Lines_df = pd.concat([Lines_df, linestr])\n",
    "            Routes_df = pd.concat([Routes_df, routes2])\n",
    "            \n",
    "            print(round(l*num / len(sets)*100,2),'%', \n",
    "                  round((time.time() - start_time) / 60,2),'mns')\n",
    "            time.sleep(time_sleep)\n",
    "        Routes_df = Routes_df.sort_values('index')\n",
    "        Routes_df = Routes_df.set_index('index')\n",
    "        Routes_df = Routes_df.reset_index(drop = True)\n",
    "        \n",
    "        Routes_df = Routes_df[Routes_df.columns[~Routes_df.columns.isin(['UGSe_id', 'GrE_id','size_infl_sqr2',\\\n",
    "                                                                  'size_infl_sqr3', 'size_infl_sqr5'])]]\n",
    "        \n",
    "        print('100 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        \n",
    "        Routes.append(Routes_df)\n",
    "        Lines.append(Lines_df)\n",
    "    return(Routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "93760a0f",
   "metadata": {},
   "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']], axis = 1)\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": 10,
   "id": "50bdbdeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_gridpark_comb (routes, var_abbr, pop_grid):\n",
    "    str1 = 'gridpark_' + var_abbr\n",
    "    locals()[str1] = routes.iloc[routes.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": 11,
   "id": "28006c20",
   "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": "64a07b5a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea1e1b7a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2da28c8c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75f67806",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "146ae178",
   "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}