E2SFCA GEE-OSM Model 3.0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 4,
   "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": 5,
   "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=AdYpXcSmNdWtjhoLdw98ofc71Z1rJXefqUyIX8oap6U&tc=YIXSOxvVWZdLwx5ya7RsnJ5zlZrhslgCQKKfsn9Bg4g&cc=yAqKv2tyjQ-49oPcFuyb9ZC2jReeSDal7XrK6zns3z4>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=AdYpXcSmNdWtjhoLdw98ofc71Z1rJXefqUyIX8oap6U&tc=YIXSOxvVWZdLwx5ya7RsnJ5zlZrhslgCQKKfsn9Bg4g&cc=yAqKv2tyjQ-49oPcFuyb9ZC2jReeSDal7XrK6zns3z4</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/1AbUR2VNhidX5cZ-qMKUGvd6fFOgvwZRdvOFZuHqAwbrW4xNj4F2mTqfWZFk\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": 18,
   "id": "2bfd480a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['United Kingdom']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/5a7391877736faaea03fabc69abfd803-8e887c9784636a4c8031381791c5adb7:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\GBR_Cambridge_2020.tif\n",
      "get road networks from OSM\n",
      "Cambridge done 0.45 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Cambridge done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Cambridge 0.18 mns\n",
      "\n",
      "get fake UGS entry points\n",
      "Cambridge 0.0 % done 0.0  mns\n",
      "Cambridge 48.3 % done 0.26  mns\n",
      "Cambridge 96.6 % done 0.49  mns\n",
      "Cambridge 100 % done 0.55  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "18.28 % 0.11 mns\n",
      "36.56 % 0.27 mns\n",
      "54.84 % 0.47 mns\n",
      "73.13 % 0.71 mns\n",
      "91.41 % 0.99 mns\n",
      "100 % finding combinations done\n",
      "Cambridge 1156686 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Cambridge\n",
      "0.0 % done 0.63 mns\n",
      "18.28 % done 0.71 mns\n",
      "36.56 % done 0.89 mns\n",
      "54.84 % done 1.18 mns\n",
      "73.13 % done 1.54 mns\n",
      "91.41 % done 1.99 mns\n",
      "100 % done 2.25 mns\n",
      "\n",
      "Cambridge\n",
      "1.88 % done 0.3 mns\n",
      "3.76 % done 0.55 mns\n",
      "5.63 % done 0.83 mns\n",
      "7.51 % done 1.15 mns\n",
      "9.39 % done 1.43 mns\n",
      "11.27 % done 1.71 mns\n",
      "13.14 % done 2.07 mns\n",
      "15.02 % done 2.43 mns\n",
      "16.9 % done 2.8 mns\n",
      "18.78 % done 3.11 mns\n",
      "20.65 % done 3.42 mns\n",
      "22.53 % done 3.7 mns\n",
      "24.41 % done 4.03 mns\n",
      "26.29 % done 4.34 mns\n",
      "28.16 % done 4.69 mns\n",
      "30.04 % done 4.99 mns\n",
      "31.92 % done 5.31 mns\n",
      "33.8 % done 5.69 mns\n",
      "35.67 % done 6.06 mns\n",
      "37.55 % done 6.46 mns\n",
      "39.43 % done 6.84 mns\n",
      "41.31 % done 7.18 mns\n",
      "43.18 % done 7.59 mns\n",
      "45.06 % done 7.91 mns\n",
      "46.94 % done 8.23 mns\n",
      "48.82 % done 8.58 mns\n",
      "50.69 % done 8.93 mns\n",
      "52.57 % done 9.28 mns\n",
      "54.45 % done 9.63 mns\n",
      "56.33 % done 9.94 mns\n",
      "58.21 % done 10.25 mns\n",
      "60.08 % done 10.61 mns\n",
      "61.96 % done 10.98 mns\n",
      "63.84 % done 11.34 mns\n",
      "65.72 % done 11.69 mns\n",
      "67.59 % done 11.98 mns\n",
      "69.47 % done 12.3 mns\n",
      "71.35 % done 12.64 mns\n",
      "73.23 % done 13.0 mns\n",
      "75.1 % done 13.38 mns\n",
      "76.98 % done 13.67 mns\n",
      "78.86 % done 14.03 mns\n",
      "80.74 % done 14.45 mns\n",
      "82.61 % done 14.89 mns\n",
      "84.49 % done 15.38 mns\n",
      "86.37 % done 15.9 mns\n",
      "88.25 % done 16.36 mns\n",
      "90.12 % done 16.74 mns\n",
      "92.0 % done 17.17 mns\n",
      "93.88 % done 17.52 mns\n",
      "95.76 % done 17.95 mns\n",
      "97.63 % done 18.4 mns\n",
      "99.51 % done 18.88 mns\n",
      "100 % done 19.41 mns\n",
      "\n",
      "300 Cambridge\n",
      "600 Cambridge\n",
      "1000 Cambridge\n",
      "CPU times: total: 24min 13s\n",
      "Wall time: 24min 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>City</th>\n",
       "      <th>Cambridge</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>137,073.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 300</th>\n",
       "      <td>231.98</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 300</th>\n",
       "      <td>65.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 300</th>\n",
       "      <td>53,626.32</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 300</th>\n",
       "      <td>85.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 300</th>\n",
       "      <td>147.14</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 600</th>\n",
       "      <td>227.26</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 600</th>\n",
       "      <td>240.29</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 600</th>\n",
       "      <td>86,428.86</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 600</th>\n",
       "      <td>99.74</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 600</th>\n",
       "      <td>141.28</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 1000</th>\n",
       "      <td>186.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 1000</th>\n",
       "      <td>481.96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 1000</th>\n",
       "      <td>77,711.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 1000</th>\n",
       "      <td>66.52</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 1000</th>\n",
       "      <td>99.83</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                      Cambridge\n",
       "0                        137,073.00\n",
       "Sc-access 300                231.98\n",
       "M-dist 300                    65.14\n",
       "M-area 300                53,626.32\n",
       "M-supply 300                  85.50\n",
       "Sc-norm 300                  147.14\n",
       "Sc-access 600                227.26\n",
       "M-dist 600                   240.29\n",
       "M-area 600                86,428.86\n",
       "M-supply 600                  99.74\n",
       "Sc-norm 600                  141.28\n",
       "Sc-access 1000               186.42\n",
       "M-dist 1000                  481.96\n",
       "M-area 1000               77,711.03\n",
       "M-supply 1000                 66.52\n",
       "Sc-norm 1000                  99.83"
      ]
     },
     "execution_count": 18,
     "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(['Cambridge'])]\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",
    "print('')\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Dir_Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, \n",
    "                             population_grids, \n",
    "                             thresholds, \n",
    "                             cities_adj['City'], \n",
    "                             save_path = 'D:/Dumps/GEE-WP Scores/E2SFCA_adj/', \n",
    "                             grid_size = 100,\n",
    "                             ext = '_Cambridge')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b4ccf2f",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['United Kingdom']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/5273964e6b7f04c3f26a4c716e76327d-f204b51c07ad3cae06c91097c22d98b7:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\GBR_Liverpool_2020.tif\n",
      "get road networks from OSM\n",
      "Liverpool done 0.94 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Liverpool done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Liverpool 0.58 mns\n",
      "\n",
      "get fake UGS entry points\n",
      "Liverpool 0.0 % done 0.0  mns\n",
      "Liverpool 38.5 % done 0.53  mns\n",
      "Liverpool 76.9 % done 1.03  mns\n",
      "Liverpool 100 % done 1.32  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "16.74 % 0.19 mns\n",
      "33.48 % 0.41 mns\n",
      "50.22 % 0.69 mns\n",
      "66.96 % 0.99 mns\n",
      "83.7 % 1.33 mns\n",
      "100 % finding combinations done\n",
      "Liverpool 1401443 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Liverpool\n",
      "0.0 % done 1.38 mns\n",
      "16.74 % done 1.46 mns\n",
      "33.48 % done 1.6 mns\n",
      "50.22 % done 1.8 mns\n",
      "66.96 % done 2.07 mns\n",
      "83.7 % done 2.39 mns\n",
      "100 % done 2.75 mns\n",
      "\n",
      "Liverpool\n",
      "1.85 % done 0.42 mns\n",
      "3.69 % done 0.83 mns\n",
      "5.54 % done 1.2 mns\n",
      "7.39 % done 1.57 mns\n",
      "9.24 % done 1.99 mns\n",
      "11.08 % done 2.39 mns\n",
      "12.93 % done 2.88 mns\n",
      "14.78 % done 3.4 mns\n",
      "16.62 % done 3.9 mns\n",
      "18.47 % done 4.36 mns\n",
      "20.32 % done 4.88 mns\n",
      "22.16 % done 5.3 mns\n",
      "24.01 % done 5.81 mns\n",
      "25.86 % done 6.32 mns\n",
      "27.71 % done 6.72 mns\n",
      "29.55 % done 7.15 mns\n",
      "31.4 % done 7.64 mns\n",
      "33.25 % done 8.12 mns\n",
      "35.09 % done 8.58 mns\n",
      "36.94 % done 9.1 mns\n",
      "38.79 % done 9.68 mns\n",
      "40.64 % done 10.19 mns\n",
      "42.48 % done 10.83 mns\n",
      "44.33 % done 11.49 mns\n",
      "46.18 % done 12.14 mns\n",
      "48.02 % done 12.73 mns\n",
      "49.87 % done 13.44 mns\n",
      "51.72 % done 13.93 mns\n",
      "53.56 % done 14.5 mns\n",
      "55.41 % done 15.05 mns\n",
      "57.26 % done 15.58 mns\n",
      "59.11 % done 16.16 mns\n",
      "60.95 % done 16.72 mns\n",
      "62.8 % done 17.29 mns\n",
      "64.65 % done 17.86 mns\n",
      "66.49 % done 18.48 mns\n",
      "68.34 % done 19.04 mns\n",
      "70.19 % done 19.63 mns\n",
      "72.04 % done 20.33 mns\n",
      "73.88 % done 20.94 mns\n",
      "75.73 % done 21.52 mns\n",
      "77.58 % done 22.11 mns\n",
      "79.42 % done 22.62 mns\n",
      "81.27 % done 23.18 mns\n",
      "83.12 % done 23.78 mns\n",
      "84.96 % done 24.47 mns\n",
      "86.81 % done 25.13 mns\n",
      "88.66 % done 25.77 mns\n",
      "90.51 % done 26.41 mns\n",
      "92.35 % done 27.01 mns\n",
      "94.2 % done 27.61 mns\n",
      "96.05 % done 28.22 mns\n",
      "97.89 % done 28.79 mns\n",
      "99.74 % done 29.37 mns\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(['Liverpool'])]\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",
    "print('')\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Dir_Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, \n",
    "                             population_grids, \n",
    "                             thresholds, \n",
    "                             cities_adj['City'], \n",
    "                             save_path = 'D:/Dumps/GEE-WP Scores/E2SFCA_adj/', \n",
    "                             grid_size = 100,\n",
    "                             ext = '_Liverpool')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e8582b76",
   "metadata": {},
   "outputs": [],
   "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(['Oslo'])]\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",
    "print('')\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Dir_Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, \n",
    "                             population_grids, \n",
    "                             thresholds, \n",
    "                             cities_adj['City'], \n",
    "                             save_path = 'D:/Dumps/GEE-WP Scores/E2SFCA_adj/', \n",
    "                             grid_size = 100,\n",
    "                             ext = '_Oslo')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "feb07c0a",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['India']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/dbd35b028bdbef09ee8ed57193db8381-ae6730b87eb597e718cd229c3458afed:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\IND_Chandigarh_2020.tif\n",
      "get road networks from OSM\n",
      "Chandigarh done 0.64 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Chandigarh done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Chandigarh 0.59 mns\n",
      "\n",
      "get fake UGS entry points\n",
      "Chandigarh 0.0 % done 0.0  mns\n",
      "Chandigarh 14.1 % done 0.36  mns\n",
      "Chandigarh 28.2 % done 0.73  mns\n",
      "Chandigarh 42.3 % done 1.09  mns\n",
      "Chandigarh 56.4 % done 1.45  mns\n",
      "Chandigarh 70.5 % done 1.81  mns\n",
      "Chandigarh 84.6 % done 2.19  mns\n",
      "Chandigarh 98.7 % done 2.55  mns\n",
      "Chandigarh 100 % done 2.61  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "16.09 % 0.17 mns\n",
      "32.18 % 0.37 mns\n",
      "48.26 % 0.56 mns\n",
      "64.35 % 0.76 mns\n",
      "80.44 % 0.97 mns\n",
      "96.53 % 1.19 mns\n",
      "100 % finding combinations done\n",
      "Chandigarh 345485 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Chandigarh\n",
      "0.0 % done 0.92 mns\n",
      "16.09 % done 1.02 mns\n",
      "32.18 % done 1.1 mns\n",
      "48.26 % done 1.18 mns\n",
      "64.35 % done 1.27 mns\n",
      "80.44 % done 1.37 mns\n",
      "96.53 % done 1.48 mns\n",
      "100 % done 1.5 mns\n",
      "\n",
      "Chandigarh\n",
      "1.78 % done 0.24 mns\n",
      "3.57 % done 0.44 mns\n",
      "5.35 % done 0.67 mns\n",
      "7.13 % done 0.92 mns\n",
      "8.92 % done 1.16 mns\n",
      "10.7 % done 1.41 mns\n",
      "12.48 % done 1.66 mns\n",
      "14.27 % done 1.9 mns\n",
      "16.05 % done 2.15 mns\n",
      "17.83 % done 2.39 mns\n",
      "19.61 % done 2.64 mns\n",
      "21.4 % done 2.9 mns\n",
      "23.18 % done 3.15 mns\n",
      "24.96 % done 3.41 mns\n",
      "26.75 % done 3.64 mns\n",
      "28.53 % done 3.9 mns\n",
      "30.31 % done 4.23 mns\n",
      "32.1 % done 4.51 mns\n",
      "33.88 % done 4.79 mns\n",
      "35.66 % done 5.81 mns\n",
      "37.45 % done 6.16 mns\n",
      "39.23 % done 6.52 mns\n",
      "41.01 % done 6.84 mns\n",
      "42.8 % done 7.12 mns\n",
      "44.58 % done 7.37 mns\n",
      "46.36 % done 7.63 mns\n",
      "48.15 % done 7.86 mns\n",
      "49.93 % done 8.08 mns\n",
      "51.71 % done 8.31 mns\n",
      "53.5 % done 8.53 mns\n",
      "55.28 % done 8.75 mns\n",
      "57.06 % done 8.97 mns\n",
      "58.84 % done 9.19 mns\n",
      "60.63 % done 9.42 mns\n",
      "62.41 % done 9.66 mns\n",
      "64.19 % done 9.89 mns\n",
      "65.98 % done 10.13 mns\n",
      "67.76 % done 10.36 mns\n",
      "69.54 % done 10.6 mns\n",
      "71.33 % done 10.83 mns\n",
      "73.11 % done 11.05 mns\n",
      "74.89 % done 11.28 mns\n",
      "76.68 % done 11.53 mns\n",
      "78.46 % done 11.79 mns\n",
      "80.24 % done 12.03 mns\n",
      "82.03 % done 12.26 mns\n",
      "83.81 % done 12.51 mns\n",
      "85.59 % done 12.77 mns\n",
      "87.38 % done 13.03 mns\n",
      "89.16 % done 13.27 mns\n",
      "90.94 % done 13.51 mns\n",
      "92.72 % done 13.76 mns\n",
      "94.51 % done 13.99 mns\n",
      "96.29 % done 14.22 mns\n",
      "98.07 % done 14.45 mns\n",
      "99.86 % done 14.68 mns\n",
      "100 % done 14.81 mns\n",
      "\n",
      "300 Chandigarh\n",
      "600 Chandigarh\n",
      "1000 Chandigarh\n",
      "CPU times: total: 20min 40s\n",
      "Wall time: 23min 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>City</th>\n",
       "      <th>Chandigarh</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>893,833.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 300</th>\n",
       "      <td>54.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 300</th>\n",
       "      <td>36.33</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 300</th>\n",
       "      <td>497,217.88</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 300</th>\n",
       "      <td>36.11</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 300</th>\n",
       "      <td>1.90</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 600</th>\n",
       "      <td>52.95</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 600</th>\n",
       "      <td>142.07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 600</th>\n",
       "      <td>651,607.44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 600</th>\n",
       "      <td>37.46</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 600</th>\n",
       "      <td>1.77</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 1000</th>\n",
       "      <td>51.55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 1000</th>\n",
       "      <td>282.71</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 1000</th>\n",
       "      <td>753,911.21</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 1000</th>\n",
       "      <td>34.55</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 1000</th>\n",
       "      <td>1.68</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                     Chandigarh\n",
       "0                        893,833.00\n",
       "Sc-access 300                 54.62\n",
       "M-dist 300                    36.33\n",
       "M-area 300               497,217.88\n",
       "M-supply 300                  36.11\n",
       "Sc-norm 300                    1.90\n",
       "Sc-access 600                 52.95\n",
       "M-dist 600                   142.07\n",
       "M-area 600               651,607.44\n",
       "M-supply 600                  37.46\n",
       "Sc-norm 600                    1.77\n",
       "Sc-access 1000                51.55\n",
       "M-dist 1000                  282.71\n",
       "M-area 1000              753,911.21\n",
       "M-supply 1000                 34.55\n",
       "Sc-norm 1000                   1.68"
      ]
     },
     "execution_count": 21,
     "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(['Chandigarh'])]\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",
    "print('')\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (Dir_Routes, population_grids, UGS)\n",
    "\n",
    "E2SCFA_score = E2SCFA_scores(min_gridUGS, \n",
    "                             population_grids, \n",
    "                             thresholds, \n",
    "                             cities_adj['City'], \n",
    "                             save_path = 'D:/Dumps/GEE-WP Scores/E2SFCA_adj/', \n",
    "                             grid_size = 100,\n",
    "                             ext = '_Chandigarh')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c0bde4ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gee_worldpop_extract (city_file, iso, save_path = None):\n",
    "    \n",
    "    cities = city_file\n",
    "    \n",
    "    # Get included city areas\n",
    "    OSM_incl = [cities[cities['City'] == city]['OSM_area'].tolist()[0].rsplit(', ') for city in cities['City'].tolist()]\n",
    "\n",
    "    # Get the city geoms\n",
    "    obj = [city_geo(city).dissolve()['geometry'].tolist()[0] for city in OSM_incl]\n",
    "\n",
    "    # Get the city countries\n",
    "    obj_displ = [city_geo(city).dissolve()['display_name'].tolist()[0].rsplit(', ')[-1]for city in OSM_incl]\n",
    "    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": 7,
   "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",
    "    print('get fake UGS entry points')\n",
    "    start_time = time.time()\n",
    "    ParkRoads = list()\n",
    "    for j in range(len(cities)):\n",
    "        ParkRoad = pd.DataFrame()\n",
    "        mat = list()\n",
    "        # For all\n",
    "        for i in range(len(UGS[j])):\n",
    "            dist = road_nodes[j]['geometry'].to_crs(3043).distance(UGS[j]['geometry'].to_crs(\n",
    "                3043)[i])\n",
    "            buf_nodes = road_nodes[j][(dist < UGS_entry_buf) & (dist > 0)]\n",
    "            mat.append(list(np.repeat(i, len(buf_nodes))))\n",
    "            ParkRoad = pd.concat([ParkRoad, buf_nodes])\n",
    "            if i % 100 == 0: print(cities[j].rsplit(',')[0], round(i/len(UGS[j])*100,1),'% done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "        # Park no list conversion\n",
    "        mat_u = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat) for i in b]\n",
    "\n",
    "        # Format\n",
    "        ParkRoad['Park_No'] = mat_u\n",
    "        ParkRoad = ParkRoad.reset_index()\n",
    "        ParkRoad['park_lon'] = ParkRoad['geometry_m'].x\n",
    "        ParkRoad['park_lat'] = ParkRoad['geometry_m'].y\n",
    "        \n",
    "        # Get the road nodes intersecting with the parks' buffer\n",
    "        ParkRoad = pd.merge(ParkRoad, UGS[j][['geometry','park_area']], left_on = 'Park_No', right_index = True)\n",
    "\n",
    "        # Get the walkable park size\n",
    "        ParkRoad['park_size_walkable'] = ParkRoad['geometry_m'].buffer(walk_radius).to_crs(4326).intersection(ParkRoad['geometry_y'].to_crs(4326))\n",
    "        ParkRoad['walk_area'] = ParkRoad['park_size_walkable'].to_crs(3043).area\n",
    "        ParkRoad['park_area'] = ParkRoad['geometry_y'].to_crs(3043).area\n",
    "        ParkRoad['share_walked'] = ParkRoad['walk_area'] / ParkRoad['park_area']\n",
    "                \n",
    "        # Merge fake UGS entry points if within X meters of each other for better system performance\n",
    "        # Standard no merging\n",
    "        ParkRoad = simplify_UGS_entry(ParkRoad, entry_point_merge = 0)\n",
    "                \n",
    "        ParkRoads.append(ParkRoad)\n",
    "\n",
    "        print(cities[j].rsplit(',')[0],'100 % done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "        \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": 8,
   "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','park_area']\n",
    "        for i in range(len(entry_geoms)):\n",
    "            suit_df = pop[pop.node_geom_m.distance(entry_geoms.iloc[i]) < 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",
    "            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','park_area':'park_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": 9,
   "id": "6553bf33",
   "metadata": {},
   "outputs": [],
   "source": [
    "def obtaining_subgraphs(graphs, pop_grids, UGS_entry, nodes, cities, thresholds):\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",
    "\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 <= np.max(thresholds)]\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: print(round((i+1) / len(entry_geoms)*100,2),'% done',\n",
    "                                        round((time.time() - start_time) / 60,2),'mns')\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": 10,
   "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": 11,
   "id": "d879c1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_routing (suitible_comb, graphs, edges, cities):\n",
    "    start_time = time.time()\n",
    "    Routes = []\n",
    "    for j in enumerate(cities):\n",
    "        print(j[1])\n",
    "        comb = suitible_comb[j[0]]\n",
    "        grouped = comb[comb['in_out_UGS'] == False].groupby(['Parkroad_osmid'])['grid_osm'].apply(list)\n",
    "        sets = grouped.apply(np.unique)\n",
    "        \n",
    "        parknode = list(comb['Parkroad_osmid'])\n",
    "        gridnode = list(comb['grid_osm'])\n",
    "        Conn = edges[j[0]]\n",
    "        subgraph = graphs[j[0]]\n",
    "        subgraph = subgraph[sets.index]\n",
    "\n",
    "        ls = []\n",
    "        ls2 = []\n",
    "        ls3 = []\n",
    "        linestr = pd.Series()\n",
    "        for i in range(len(sets)):\n",
    "            path = nx.single_source_dijkstra(subgraph.iloc[i], sets.index[i], weight = 'length')\n",
    "\n",
    "            # Only include routes that were prespecified, order depends on route cost, low to high, and steps low to high.\n",
    "            # Grid destinations from UGS entry points therefore may be ranked differently, and subsetted separately.\n",
    "            incl = np.isin(list(path[0].keys()),sets.iloc[i])\n",
    "            incl2 = np.isin(list(path[1].keys()),sets.iloc[i])\n",
    "\n",
    "            # route cost\n",
    "            orig_c = list(np.repeat(sets.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(sets.index[i],sum(incl2)))\n",
    "            dest_s = list(np.array(list(path[1].keys()))[incl2])\n",
    "            steps = list(np.array(list(path[1].values()))[incl2])\n",
    "\n",
    "            fr = []\n",
    "            to = []\n",
    "            og = []\n",
    "            de = []\n",
    "            for j in enumerate(steps):\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",
    "                \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",
    "            if len(od) > 0:\n",
    "                ser_gk = pd.DataFrame(gkr, index = gk)\n",
    "                ser_gk['od'] = od\n",
    "                ser_gk = ser_gk.join(Conn.geometry, how = 'left')\n",
    "                ser_gk = ser_gk.set_index(ser_gk.columns[0], drop = True)\n",
    "                ser_gk = ser_gk.join(Conn.geometry, how = 'left', rsuffix = '_r')\n",
    "                ser_gk['geom'] = np.where(ser_gk.geometry.isna(),ser_gk.geometry_r,ser_gk.geometry)\n",
    "                diss = gpd.GeoDataFrame(ser_gk[['od','geom']], geometry = 'geom', crs = 4326).dissolve(by = 'od')\n",
    "                diss = diss.geom\n",
    "                linestr = pd.concat([linestr, diss], axis = 0)\n",
    "\n",
    "            if (i+1) % 50 == 0: print(round((i+1) / len(sets)*100,2),'% done',round((time.time() - start_time) / 60,2),'mns')\n",
    "\n",
    "        linestr = linestr.rename('geometry')\n",
    "\n",
    "        dist_df = pd.DataFrame([ls, ls2, ls3]).transpose()\n",
    "        dist_df.columns = ['UGSe_id','GrE_id','route cost']\n",
    "        dist_df['UGSe_id'] = [int(i) for i in dist_df['UGSe_id']]\n",
    "        dist_df['GrE_id'] = [int(i) for i in dist_df['GrE_id']]\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['G-entry cost'] = np.where(routes['in_out_UGS'],0,routes['G-entry cost'])\n",
    "        \n",
    "        routes['Tcost'] = routes['route cost']+routes['G-entry cost']\n",
    "        \n",
    "        routes = pd.merge(routes, linestr, left_on = 'graph_key', right_index = True, how = 'left')\n",
    "        \n",
    "        routes2 = routes.iloc[routes['route cost'].dropna().index].reset_index(drop = True)\n",
    "        \n",
    "        print('100 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        \n",
    "        Routes.append(routes2)\n",
    "    return(Routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "50bdbdeb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_gridUGS_comb (routes, grids, UGS):\n",
    "    gp_nearest = []\n",
    "    for i in range(len(routes)):\n",
    "        gp_nn = routes[i][routes[i]['Tcost'] <= max(thresholds)]\n",
    "        gp_nn = pd.merge(gp_nn, grids[i]['population'], left_on='Grid_No', right_index = True)\n",
    "        gp_nn = pd.merge(gp_nn, UGS[i]['park_area'], left_on = 'Park_No', right_index = True)\n",
    "        gp_nn = gp_nn.reset_index()\n",
    "\n",
    "        gp_nn = gp_nn.iloc[gp_nn.groupby('gridpark_no')['Tcost'].idxmin()]\n",
    "        gp_nn.index.name = 'idx'\n",
    "        gp_nn = gp_nn.sort_values('idx')\n",
    "        gp_nn = gp_nn.reset_index()\n",
    "        gp_nearest.append(gp_nn)\n",
    "    gp_nearest[0].sort_values('Grid_No')\n",
    "    return(gp_nearest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "93760a0f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def E2SCFA_scores(min_gridUGS_comb, grids, thresholds, cities, \n",
    "                  save_path = 'D:/Dumps/GEE-WP Scores/E2SFCA/', grid_size = 100, ext = ''):\n",
    "    pd.options.display.float_format = '{:20,.2f}'.format\n",
    "    E2SFCA_cities = []\n",
    "    E2SFCA_summary = pd.DataFrame()\n",
    "    for i in range(len(cities)):\n",
    "        E2SFCA_score = grids[i][['population','geometry']]\n",
    "        for j in range(len(thresholds)):\n",
    "            subset = min_gridUGS_comb[i][min_gridUGS_comb[i]['Tcost'] <= thresholds[j]]\n",
    "\n",
    "            # use gussian distribution: let v= 923325, then the weight for 800m is 0.5\n",
    "            v = -thresholds[j]**2/np.log(0.5)\n",
    "\n",
    "            # add a column of weight: apply the decay function on distance\n",
    "            subset['weight'] = np.exp(-(subset['Tcost']**2/v)).astype(float)\n",
    "            subset['pop_weight'] = subset['weight'] * subset['population']\n",
    "\n",
    "            # get the sum of weighted population each green space has to serve.\n",
    "            s_w_p = pd.DataFrame(subset.groupby('Park_No').sum('pop_weight')['pop_weight'])\n",
    "\n",
    "            # delete other columns, because they are useless after groupby\n",
    "            s_w_p = s_w_p.rename({'pop_weight':'pop_weight_sum'},axis = 1)\n",
    "            middle = pd.merge(subset,s_w_p, how = 'left', on = 'Park_No' )\n",
    "\n",
    "            # calculate the supply-demand ratio for each green space\n",
    "            middle['green_supply'] = middle['park_area']/middle['pop_weight_sum']\n",
    "\n",
    "            # caculate the accessbility score for each green space that each population grid cell could reach\n",
    "            middle['Sc-access'] = middle['weight'] * middle['green_supply']\n",
    "            # add the scores for each population grid cell\n",
    "            pop_score_df = pd.DataFrame(middle.groupby('Grid_No').sum('Sc-access')['Sc-access'])\n",
    "\n",
    "            # calculate the mean distance of all the green space each population grid cell could reach\n",
    "            mean_dist = middle.groupby('Grid_No').mean('Tcost')['Tcost']\n",
    "            pop_score_df['M-dist'] = mean_dist\n",
    "\n",
    "            # calculate the mean area of all the green space each population grid cell could reach\n",
    "            mean_area = middle.groupby('Grid_No').mean('park_area')['park_area']\n",
    "            pop_score_df['M-area'] = mean_area\n",
    "\n",
    "            # calculate the mean supply_demand ratio of all the green space each population grid cell could reach\n",
    "            mean_supply = middle.groupby('Grid_No').mean('green_supply')['green_supply']\n",
    "            pop_score_df['M-supply'] = mean_supply\n",
    "\n",
    "            pop_score = pop_score_df\n",
    "\n",
    "            pop_score_df = pop_score_df.join(grids[i]['population'], how = 'right')\n",
    "            pop_score_df['Sc-norm'] = pop_score_df['Sc-access'] / pop_score_df['population']\n",
    "\n",
    "            pop_score_df = pop_score_df.loc[:, pop_score_df.columns != 'population']\n",
    "            pop_score_df = pop_score_df.add_suffix(' '+str(thresholds[j]))\n",
    "            E2SFCA_score = E2SFCA_score.join(pop_score_df, how = 'left')\n",
    "\n",
    "            print(thresholds[j], cities[i])\n",
    "\n",
    "        E2SFCA_score = E2SFCA_score.fillna(0)\n",
    "        \n",
    "        if not os.path.exists(save_path+str(grid_size)+'m grids'+'/grid_geoms/'):\n",
    "            os.makedirs(save_path+str(grid_size)+'m grids'+'/grid_geoms/')\n",
    "        \n",
    "        E2SFCA_score.to_file(save_path+str(grid_size)+'m grids'+'/grid_geoms/'+cities[i]+'.gpkg') # Detailed scores\n",
    "        pop_sum = pd.Series(E2SFCA_score['population'].sum()).astype(int)\n",
    "        mean_metrics = E2SFCA_score.loc[:, ~E2SFCA_score.columns.isin(['population','geometry'])].mean()\n",
    "        E2SFCA_sum = pd.concat([pop_sum, mean_metrics])\n",
    "        E2SFCA_summary = pd.concat([E2SFCA_summary, E2SFCA_sum], axis = 1) # summarized results\n",
    "        E2SFCA_cities.append(E2SFCA_score)\n",
    "        \n",
    "        if not os.path.exists(save_path):\n",
    "            os.makedirs(save_path)\n",
    "        \n",
    "        E2SFCA_score.loc[:, E2SFCA_score.columns != 'geometry'].to_csv(save_path+cities[i]+'.csv')\n",
    "    E2SFCA_summary.columns = cities\n",
    "    \n",
    "    if not os.path.exists(save_path):\n",
    "        os.makedirs(save_path)\n",
    "    \n",
    "    E2SFCA_summary.to_csv(save_path+str(grid_size)+'m grids'+'all_cities'+ext+'.csv')\n",
    "    E2SFCA_summary\n",
    "    return({'score summary':E2SFCA_summary,'score detail':E2SFCA_cities})"
   ]
  },
  {
   "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": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}