Gravity GEE-OSM Model 3.0.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "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": 7,
   "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=m64lvihVXahIHqiE0UaCd4yLymz4mo8NLSsiFZkPsaw&tc=q_IuREh-VxUCjwkbwE_90m4hwESgHRtgjkzNn3lWXQ8&cc=ePYRkaTg-dec9-nlnWGAcbpFkGouLzFuZYXuKHGJUY4>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=m64lvihVXahIHqiE0UaCd4yLymz4mo8NLSsiFZkPsaw&tc=q_IuREh-VxUCjwkbwE_90m4hwESgHRtgjkzNn3lWXQ8&cc=ePYRkaTg-dec9-nlnWGAcbpFkGouLzFuZYXuKHGJUY4</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/1AbUR2VP3gvYyzJv5UPWJqog4nBJBs8UESTB5KOjVPME0JltoNAH2NXkjCUU\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": 42,
   "id": "8049e46b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Vietnam']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/e2dcaa7fc8ef07556ab7052aba2ca3a2-2501388843005501fd75142997a66968:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\VNM_Hanoi_2020.tif\n",
      "get road networks from OSM\n",
      "Hanoi done 3.6 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Hanoi done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Hanoi 0.47 mns\n",
      "\n",
      "Hanoi\n",
      "0.0 % fake entry points done 0.0  mns\n",
      "29.4 % fake entry points done 0.64  mns\n",
      "58.8 % fake entry points done 1.3  mns\n",
      "88.2 % fake entry points done 1.88  mns\n",
      "Hanoi 100 % fake entry points done 2.12  mns\n",
      "Hanoi 100 % done 2.13  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "33.9 % 0.18 mns\n",
      "67.8 % 0.37 mns\n",
      "100 % finding combinations done\n",
      "\n",
      "obtain local graphs\n",
      "Hanoi\n",
      "0.0 % done 1.22 mns\n",
      "33.9 % done 1.52 mns\n",
      "67.8 % done 1.84 mns\n",
      "100 % done 2.17 mns\n",
      "\n",
      "Hanoi\n",
      "1.93 % done 0.26 mns\n",
      "3.86 % done 0.53 mns\n",
      "5.8 % done 0.78 mns\n",
      "7.73 % done 1.04 mns\n",
      "9.66 % done 1.29 mns\n",
      "11.59 % done 1.54 mns\n",
      "13.52 % done 1.8 mns\n",
      "15.46 % done 2.08 mns\n",
      "17.39 % done 2.31 mns\n",
      "19.32 % done 2.58 mns\n",
      "21.25 % done 2.84 mns\n",
      "23.18 % done 3.09 mns\n",
      "25.12 % done 3.34 mns\n",
      "27.05 % done 3.6 mns\n",
      "28.98 % done 3.87 mns\n",
      "30.91 % done 4.12 mns\n",
      "32.84 % done 4.38 mns\n",
      "34.78 % done 4.64 mns\n",
      "36.71 % done 4.91 mns\n",
      "38.64 % done 5.17 mns\n",
      "40.57 % done 5.44 mns\n",
      "42.5 % done 5.69 mns\n",
      "44.44 % done 5.96 mns\n",
      "46.37 % done 6.21 mns\n",
      "48.3 % done 6.47 mns\n",
      "50.23 % done 6.74 mns\n",
      "52.16 % done 7.01 mns\n",
      "54.1 % done 7.29 mns\n",
      "56.03 % done 7.56 mns\n",
      "57.96 % done 7.81 mns\n",
      "59.89 % done 8.09 mns\n",
      "61.82 % done 8.37 mns\n",
      "63.76 % done 8.63 mns\n",
      "65.69 % done 8.94 mns\n",
      "67.62 % done 9.22 mns\n",
      "69.55 % done 9.48 mns\n",
      "71.48 % done 9.73 mns\n",
      "73.42 % done 10.0 mns\n",
      "75.35 % done 10.26 mns\n",
      "77.28 % done 10.53 mns\n",
      "79.21 % done 10.79 mns\n",
      "81.14 % done 11.06 mns\n",
      "83.08 % done 11.32 mns\n",
      "85.01 % done 11.59 mns\n",
      "86.94 % done 11.86 mns\n",
      "88.87 % done 12.13 mns\n",
      "90.8 % done 12.42 mns\n",
      "92.74 % done 12.68 mns\n",
      "94.67 % done 12.95 mns\n",
      "96.6 % done 13.22 mns\n",
      "98.53 % done 13.5 mns\n",
      "missing suitible comb in graph 4.2186 %\n",
      "100 % done 13.77 mns\n",
      "\n",
      "Hanoi\n",
      "entrance 0.04 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/2) 0.85 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/3) 1.8 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "gravity**(1/5) 2.65 mns\n",
      "grid  300\n",
      "grid  600\n",
      "grid  1000\n",
      "Hanoi done 3.38 mns\n",
      "CPU times: total: 26min\n",
      "Wall time: 26min 10s\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>Hanoi</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.077236</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.016875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.083527</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.822362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.120740</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.099747</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.203736</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.575777</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">entrance_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.237582</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.162705</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.198612</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.401101</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.078192</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.042870</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.158614</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.720323</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.157872</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.201049</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.284090</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.356989</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/2)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.323774</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.268241</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.108636</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.299349</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.077138</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.025224</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.101957</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.795681</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.111969</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.141721</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.271441</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.474869</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/3)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.278883</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.223186</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.197362</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.300569</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_300</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.076773</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.017639</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.084419</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.821169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_600</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.107272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.112444</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.233493</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.546792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">gravity**(1/5)_1000</th>\n",
       "      <th>1 high</th>\n",
       "      <td>0.252231</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2 medium</th>\n",
       "      <td>0.185326</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3 low</th>\n",
       "      <td>0.240855</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4 no</th>\n",
       "      <td>0.321588</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                             Hanoi\n",
       "entrance_300        1 high    0.077236\n",
       "                    2 medium  0.016875\n",
       "                    3 low     0.083527\n",
       "                    4 no      0.822362\n",
       "entrance_600        1 high    0.120740\n",
       "                    2 medium  0.099747\n",
       "                    3 low     0.203736\n",
       "                    4 no      0.575777\n",
       "entrance_1000       1 high    0.237582\n",
       "                    2 medium  0.162705\n",
       "                    3 low     0.198612\n",
       "                    4 no      0.401101\n",
       "gravity**(1/2)_300  1 high    0.078192\n",
       "                    2 medium  0.042870\n",
       "                    3 low     0.158614\n",
       "                    4 no      0.720323\n",
       "gravity**(1/2)_600  1 high    0.157872\n",
       "                    2 medium  0.201049\n",
       "                    3 low     0.284090\n",
       "                    4 no      0.356989\n",
       "gravity**(1/2)_1000 1 high    0.323774\n",
       "                    2 medium  0.268241\n",
       "                    3 low     0.108636\n",
       "                    4 no      0.299349\n",
       "gravity**(1/3)_300  1 high    0.077138\n",
       "                    2 medium  0.025224\n",
       "                    3 low     0.101957\n",
       "                    4 no      0.795681\n",
       "gravity**(1/3)_600  1 high    0.111969\n",
       "                    2 medium  0.141721\n",
       "                    3 low     0.271441\n",
       "                    4 no      0.474869\n",
       "gravity**(1/3)_1000 1 high    0.278883\n",
       "                    2 medium  0.223186\n",
       "                    3 low     0.197362\n",
       "                    4 no      0.300569\n",
       "gravity**(1/5)_300  1 high    0.076773\n",
       "                    2 medium  0.017639\n",
       "                    3 low     0.084419\n",
       "                    4 no      0.821169\n",
       "gravity**(1/5)_600  1 high    0.107272\n",
       "                    2 medium  0.112444\n",
       "                    3 low     0.233493\n",
       "                    4 no      0.546792\n",
       "gravity**(1/5)_1000 1 high    0.252231\n",
       "                    2 medium  0.185326\n",
       "                    3 low     0.240855\n",
       "                    4 no      0.321588"
      ]
     },
     "execution_count": 42,
     "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(['Hanoi'])]\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",
    "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 = '_Hanoi', # 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": 9,
   "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",
    "        entry_index = [int(road_nodes[i]['geometry'].sindex.nearest(grid)[1])\\\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": 10,
   "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": 11,
   "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",
    "    return(suits_all)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "bd28337b",
   "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",
    "        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: 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": 13,
   "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": null,
   "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['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",
    "        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": 18,
   "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']])\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": 19,
   "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": 20,
   "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": "1bbb2699",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9439f540",
   "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
}