E2SFCA GEE-OSM Model 3.1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "71639573",
   "metadata": {},
   "outputs": [],
   "source": [
    "# system packages\n",
    "import sys\n",
    "import time\n",
    "import warnings\n",
    "import os\n",
    "\n",
    "# non-geo numeric packages\n",
    "import numpy as np\n",
    "import math\n",
    "from itertools import product, combinations\n",
    "import pandas as pd\n",
    "\n",
    "# network and OSM packages\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "city_geo = ox.geocoder.geocode_to_gdf\n",
    "\n",
    "# Earth engine packages\n",
    "import ee\n",
    "import geemap\n",
    "\n",
    "# General geo-packages\n",
    "import libpysal\n",
    "import rasterio\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "4f87c388",
   "metadata": {},
   "outputs": [
    {
     "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>geometry</th>\n",
       "      <th>bbox_north</th>\n",
       "      <th>bbox_south</th>\n",
       "      <th>bbox_east</th>\n",
       "      <th>bbox_west</th>\n",
       "      <th>place_id</th>\n",
       "      <th>osm_type</th>\n",
       "      <th>osm_id</th>\n",
       "      <th>lat</th>\n",
       "      <th>lon</th>\n",
       "      <th>display_name</th>\n",
       "      <th>class</th>\n",
       "      <th>type</th>\n",
       "      <th>importance</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>POLYGON ((130.02576 32.31688, 130.02594 32.316...</td>\n",
       "      <td>32.32</td>\n",
       "      <td>32.32</td>\n",
       "      <td>130.03</td>\n",
       "      <td>130.03</td>\n",
       "      <td>278956318</td>\n",
       "      <td>way</td>\n",
       "      <td>924413802</td>\n",
       "      <td>32.32</td>\n",
       "      <td>130.03</td>\n",
       "      <td>Amakusa City Hall, Kawaura Tomitsu Branch, Sak...</td>\n",
       "      <td>amenity</td>\n",
       "      <td>townhall</td>\n",
       "      <td>0.20</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                            geometry           bbox_north  \\\n",
       "0  POLYGON ((130.02576 32.31688, 130.02594 32.316...                32.32   \n",
       "\n",
       "            bbox_south            bbox_east            bbox_west   place_id  \\\n",
       "0                32.32               130.03               130.03  278956318   \n",
       "\n",
       "  osm_type     osm_id                  lat                  lon  \\\n",
       "0      way  924413802                32.32               130.03   \n",
       "\n",
       "                                        display_name    class      type  \\\n",
       "0  Amakusa City Hall, Kawaura Tomitsu Branch, Sak...  amenity  townhall   \n",
       "\n",
       "            importance  \n",
       "0                 0.20  "
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "city_geo('Kumamoto City')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c23309aa",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<p>To authorize access needed by Earth Engine, open the following\n",
       "        URL in a web browser and follow the instructions:</p>\n",
       "        <p><a href=https://code.earthengine.google.com/client-auth?scopes=https%3A//www.googleapis.com/auth/earthengine%20https%3A//www.googleapis.com/auth/devstorage.full_control&request_id=3PbLqUmBvwdITQxYFqBvhrkv5YLuiyqmTaDLXRIGguE&tc=Np4pwiTJhyqxL6vbtamERM7YoP4FgBKVWmNu_8FnZAQ&cc=VMZq-fMXm4PUb6gYZTrmmWMdbToCuuoXaf7oPfQhOuo>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=3PbLqUmBvwdITQxYFqBvhrkv5YLuiyqmTaDLXRIGguE&tc=Np4pwiTJhyqxL6vbtamERM7YoP4FgBKVWmNu_8FnZAQ&cc=VMZq-fMXm4PUb6gYZTrmmWMdbToCuuoXaf7oPfQhOuo</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/1AbUR2VPNtosClEEvAYVUM6aovF9ZgWy9FW-Y61YmsBL7yno94EEnbtHAYP4\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": 13,
   "id": "885b9c5b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['Japan']\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/6f1a7b1cf67cfcaf3e199e84fc208293-a6b3054675f927bde2e13995aecda08d:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\JPN_Kumamoto_2020.tif\n",
      "get road networks from OSM\n",
      "Kumamoto done 2.94 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Kumamoto done\n",
      " \n",
      "100m resolution grids extraction\n",
      "Kumamoto 2.21 mns\n",
      "\n",
      "get fake UGS entry points\n",
      "Kumamoto 0.0 % done 0.0  mns\n",
      "Kumamoto 26.8 % done 0.48  mns\n",
      "Kumamoto 53.6 % done 0.88  mns\n",
      "Kumamoto 80.4 % done 1.28  mns\n",
      "Kumamoto 100 % done 1.57  mns\n",
      "\n",
      "get (Euclidean) suitible combinations\n",
      "0.0 % 0.0 mns\n",
      "34.18 % 0.38 mns\n",
      "68.35 % 0.81 mns\n",
      "100 % finding combinations done\n",
      "Kumamoto 341618 suitible combinations\n",
      "\n",
      "obtain local graphs\n",
      "Kumamoto\n",
      "0.0 % done 0.78 mns\n",
      "34.18 % done 0.82 mns\n",
      "68.35 % done 1.39 mns\n",
      "100 % done 45.61 mns\n",
      "\n",
      "Kumamoto\n",
      "0.0 % 0.14 mns\n",
      "9.68 % 0.52 mns\n",
      "19.36 % 0.89 mns\n",
      "29.04 % 1.27 mns\n",
      "38.72 % 1.65 mns\n",
      "48.4 % 2.03 mns\n",
      "58.08 % 2.39 mns\n",
      "67.76 % 2.77 mns\n",
      "77.44 % 3.15 mns\n",
      "87.12 % 3.5 mns\n",
      "96.8 % 3.77 mns\n",
      "100 % done 4.03 mns\n",
      "\n",
      "300 Kumamoto\n",
      "600 Kumamoto\n",
      "1000 Kumamoto\n",
      "CPU times: total: 10min 33s\n",
      "Wall time: 58min 23s\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>Kumamoto</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>843,427.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 300</th>\n",
       "      <td>61.42</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 300</th>\n",
       "      <td>9.44</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 300</th>\n",
       "      <td>29,409.06</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 300</th>\n",
       "      <td>45.37</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 300</th>\n",
       "      <td>29.63</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 600</th>\n",
       "      <td>60.36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 600</th>\n",
       "      <td>59.47</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 600</th>\n",
       "      <td>47,604.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 600</th>\n",
       "      <td>53.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 600</th>\n",
       "      <td>28.58</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-access 1000</th>\n",
       "      <td>58.69</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-dist 1000</th>\n",
       "      <td>186.17</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-area 1000</th>\n",
       "      <td>66,218.64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>M-supply 1000</th>\n",
       "      <td>50.99</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>Sc-norm 1000</th>\n",
       "      <td>27.19</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "City                       Kumamoto\n",
       "0                        843,427.00\n",
       "Sc-access 300                 61.42\n",
       "M-dist 300                     9.44\n",
       "M-area 300                29,409.06\n",
       "M-supply 300                  45.37\n",
       "Sc-norm 300                   29.63\n",
       "Sc-access 600                 60.36\n",
       "M-dist 600                    59.47\n",
       "M-area 600                47,604.54\n",
       "M-supply 600                  53.54\n",
       "Sc-norm 600                   28.58\n",
       "Sc-access 1000                58.69\n",
       "M-dist 1000                  186.17\n",
       "M-area 1000               66,218.64\n",
       "M-supply 1000                 50.99\n",
       "Sc-norm 1000                  27.19"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "# Thresholds and cities\n",
    "os.environ['KMP_DUPLICATE_LIB_OK'] = 'True'\n",
    "\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m\n",
    "\n",
    "# Extract cities list\n",
    "iso = pd.read_excel('iso_countries.xlsx')\n",
    "cities = pd.read_excel('cities.xlsx')\n",
    "cities_adj = cities[cities['City'].isin(['Kumamoto'])]\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 = '_Kumamoto')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c0bde4ec",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gee_worldpop_extract (city_file, iso, save_path = None):\n",
    "    \n",
    "    cities = city_file\n",
    "    iso['name'] = np.where(iso['name'] == 'Macedonia','North Macedonia',iso['name'])\n",
    "    \n",
    "    # Get included city areas\n",
    "    OSM_incl = [cities[cities['City'] == city]['OSM_area'].tolist()[0].rsplit(', ') for city in cities['City'].tolist()]\n",
    "\n",
    "    # Get the city geoms\n",
    "    obj = [city_geo(city).dissolve()['geometry'].tolist()[0] for city in OSM_incl]\n",
    "    \n",
    "    # Get the city countries\n",
    "    obj_displ = [city_geo(city).dissolve()['display_name'].tolist()[0].rsplit(', ')[-1]for city in OSM_incl]\n",
    "    print(obj_displ)\n",
    "    obj_displ = np.where(pd.Series(obj_displ).str.contains(\"Ivoire\"),\"CIte dIvoire\",obj_displ)\n",
    "\n",
    "    # Get the country's iso-code\n",
    "    iso_list = [iso[iso['name'] == ob]['alpha3'].tolist()[0] for ob in obj_displ]\n",
    "\n",
    "    # Based on the iso-code return the worldpop 2020\n",
    "    ee_worldpop = [ee.ImageCollection(\"WorldPop/GP/100m/pop\")\\\n",
    "        .filter(ee.Filter.date('2020'))\\\n",
    "        .filter(ee.Filter.inList('country', [io])).first() for io in iso_list]\n",
    "\n",
    "    # Clip the countries with the city geoms.\n",
    "    clipped = [ee_worldpop[i].clip(shapely.geometry.mapping(obj[i])) for i in range(0,len(obj))]\n",
    "\n",
    "    # Create path if non-existent\n",
    "    if save_path == None:\n",
    "        path = ''\n",
    "    else:\n",
    "        path = save_path\n",
    "        if not os.path.exists(path):\n",
    "                    os.makedirs(path)\n",
    "\n",
    "    # Export as TIFF file.\n",
    "    # Stored in form path + USA_Los Angeles_2020.tif\n",
    "    filenames = [path+iso_list[i]+'_'+cities['City'][i]+'_2020.tif' for i in range(len(obj))]\n",
    "    [geemap.ee_export_image(clipped[i], filename = filenames[i]) for i in range(0,len(obj))]\n",
    "    return(filenames)\n",
    "    sys.stdout.flush()\n",
    "    \n",
    "    # Block 2 Road networks\n",
    "def road_network (cities, thresholds, undirected = False):\n",
    "    print('get road networks from OSM')\n",
    "    start_time = time.time()\n",
    "    graphs = list()\n",
    "    road_nodes = list()\n",
    "    road_edges = list()\n",
    "    road_conn = list()\n",
    "\n",
    "    for i in enumerate(cities['OSM_area']):\n",
    "        # Get graph, road nodes and edges\n",
    "        road_node = pd.DataFrame()\n",
    "        roads = pd.DataFrame()\n",
    "        \n",
    "        # For each included OSM_area get the roads\n",
    "        for district in i[1].rsplit(', '):\n",
    "            graph = ox.graph_from_place(district, network_type = \"all\", buffer_dist = (np.max(thresholds)+1000))\n",
    "            node, edge = ox.graph_to_gdfs(graph)\n",
    "            road_node = pd.concat([road_node, node], axis = 0)\n",
    "            roads = pd.concat([roads, edge], axis = 0)\n",
    "        \n",
    "        # Eliminate lists in the df which prevents drop of duplicate columns\n",
    "        road_edge = pd.DataFrame([[c[0] if isinstance(c,list) else c for c in roads[col]]\\\n",
    "                              for col in roads]).transpose()\n",
    "        road_edge.columns = roads.columns\n",
    "        road_edge.index = roads.index\n",
    "        road_edge = gpd.GeoDataFrame(road_edge, crs = 4326)\n",
    "        \n",
    "        # Return the unique nodes and edges of the (often) adjacent OSM_areas.\n",
    "        road_node = road_node.drop_duplicates()\n",
    "        road_edge = road_edge.drop_duplicates()\n",
    "        \n",
    "        # Road nodes format\n",
    "        road_node = road_node.to_crs(4326)\n",
    "        road_node['geometry_m'] = gpd.GeoSeries(road_node['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_node['osmid_var'] = road_node.index\n",
    "        road_node = gpd.GeoDataFrame(road_node, geometry = 'geometry', crs = 4326)\n",
    "\n",
    "        # format road edges\n",
    "        road_edge['geometry_m'] = gpd.GeoSeries(road_edge['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_edge = road_edge.reset_index()\n",
    "        road_edge.rename(columns={'u':'from', 'v':'to', 'key':'keys'}, inplace=True)\n",
    "        road_edge['key'] = road_edge['from'].astype(str) + '-' + road_edge['to'].astype(str)\n",
    "        \n",
    "        if undirected == True:\n",
    "            # Apply one-directional to both for walking\n",
    "            both = road_edge[road_edge['oneway'] == False]\n",
    "            one = road_edge[road_edge['oneway'] == True]\n",
    "            rev = pd.DataFrame()\n",
    "            rev[['from','to']] = one[['to','from']]\n",
    "            rev = pd.concat([rev,one.iloc[:,2:]],axis = 1)\n",
    "            edge_bidir = pd.concat([both, one, rev])\n",
    "            edge_bidir = edge_bidir.reset_index()\n",
    "            edge_bidir['oneway'] = False\n",
    "        else:\n",
    "            edge_bidir = road_edge\n",
    "\n",
    "        # Exclude highways and ramps on edges    \n",
    "        edge_filter = edge_bidir[(edge_bidir['highway'].str.contains('motorway') | \n",
    "              (edge_bidir['highway'].str.contains('trunk') & \n",
    "               edge_bidir['maxspeed'].astype(str).str.contains(\n",
    "                   '40 mph|45 mph|50 mph|55 mph|60 mph|65|70|75|80|85|90|95|100|110|120|130|140'))) == False]\n",
    "        road_edges.append(edge_filter)\n",
    "\n",
    "        # Exclude isolated nodes\n",
    "        fltrnodes = pd.Series(list(edge_filter['from']) + list(edge_filter['to'])).unique()\n",
    "        newnodes = road_node[road_node['osmid_var'].isin(fltrnodes)]\n",
    "        road_nodes.append(newnodes)\n",
    "\n",
    "        # Get only necessary road connections columns for network performance\n",
    "        road_con = edge_filter[['osmid','key','length','geometry']]\n",
    "        road_con = road_con.set_index('key')\n",
    "\n",
    "        road_conn.append(road_con)\n",
    "\n",
    "        # formatting to graph again.\n",
    "        newnodes = newnodes.loc[:, ~newnodes.columns.isin(['geometry_m', 'osmid_var'])]\n",
    "        edge_filter = edge_filter.set_index(['from','to','keys'])\n",
    "        edge_filter = edge_filter.loc[:, ~edge_filter.columns.isin(['geometry_m', 'key'])]\n",
    "\n",
    "        graph2 = ox.graph_from_gdfs(newnodes, edge_filter)\n",
    "\n",
    "        graphs.append(graph2)\n",
    "        print(cities['City'][i[0]].rsplit(',')[0], 'done', round((time.time() - start_time) / 60,2),'mns')\n",
    "    return({'graphs':graphs,'nodes':road_nodes,'edges':road_conn,'edges long':road_edges})\n",
    "# Block 3 city greenspace\n",
    "def urban_greenspace (cities, thresholds, one_UGS_buf = 25, min_UGS_size = 400):\n",
    "    print('get urban greenspaces from OSM')\n",
    "    parks_in_range = list()\n",
    "    for i in enumerate(cities['OSM_area']):\n",
    "        # Tags seen as Urban Greenspace (UGS) require the following:\n",
    "        # 1. Tag represent an area\n",
    "        # 2. The area is outdoor\n",
    "        # 3. The area is (semi-)publically available\n",
    "        # 4. The area is likely to contain trees, grass and/or greenery\n",
    "        # 5. The area can reasonable be used for walking or recreational activities\n",
    "        tags = {'landuse':['allotments','forest','greenfield','village_green'],\\\n",
    "                'leisure':['garden','fitness_station','nature_reserve','park','playground'],\\\n",
    "                'natural':'grassland'}\n",
    "        gdf = ox.geometries_from_place(i[1].rsplit(', '),tags = tags,buffer_dist = np.max(thresholds))\n",
    "        gdf = gdf[(gdf.geom_type == 'Polygon') | (gdf.geom_type == 'MultiPolygon')]\n",
    "        greenspace = gdf.reset_index()    \n",
    "        warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "        green_buffer = gpd.GeoDataFrame(geometry = greenspace.to_crs(3043).buffer(one_UGS_buf).to_crs(4326))\n",
    "        greenspace['geometry_w_buffer'] = green_buffer\n",
    "        greenspace['geometry_w_buffer'] = gpd.GeoSeries(greenspace['geometry_w_buffer'], crs = 4326)\n",
    "        greenspace['geom buffer diff'] = greenspace['geometry_w_buffer'].difference(greenspace['geometry'])\n",
    "\n",
    "        # This function group components in itself that overlap (with the buffer set of 25 metres)\n",
    "        # https://stackoverflow.com/questions/68036051/geopandas-self-intersection-grouping\n",
    "        W = libpysal.weights.fuzzy_contiguity(greenspace['geometry_w_buffer'])\n",
    "        greenspace['components'] = W.component_labels\n",
    "        parks = greenspace.dissolve('components')\n",
    "\n",
    "        # Exclude parks below 0.04 ha.\n",
    "        parks = parks[parks.to_crs(3043).area > min_UGS_size]\n",
    "        print(cities['City'][i[0]], 'done')\n",
    "        parks = parks.reset_index()\n",
    "        parks['geometry_m'] = parks['geometry'].to_crs(3043)\n",
    "        parks['park_area'] = parks['geometry_m'].area\n",
    "        parks_in_range.append(parks)\n",
    "    return(parks_in_range)\n",
    "# Block 4 population grids extraction\n",
    "def city_grids_format(city_grids, cities_area, road_nodes, UGS, grid_size = 100):\n",
    "    start_time = time.time()\n",
    "    grids = []\n",
    "    print(str(grid_size) + 'm resolution grids extraction')\n",
    "    for i in range(len(city_grids)):\n",
    "        \n",
    "        # Open the raster file\n",
    "        with rasterio.open(city_grids[i]) as src:\n",
    "            band= src.read() # the population values\n",
    "            aff = src.transform # the raster bounds and size (affine)\n",
    "        \n",
    "        # Get the rowwise arrays, get a 2D dataframe\n",
    "        grid = pd.DataFrame()\n",
    "        for b in enumerate(band[0]):\n",
    "            grid = pd.concat([grid, pd.Series(b[1],name=b[0])],axis=1)\n",
    "        grid= grid.unstack().reset_index()\n",
    "        \n",
    "        # Unstack df to columns\n",
    "        grid.columns = ['row','col','value']\n",
    "        grid['minx'] = aff[2]+aff[0]*grid['col']\n",
    "        grid['miny'] = aff[5]+aff[4]*grid['row']\n",
    "        grid['maxx'] = aff[2]+aff[0]*grid['col']+aff[0]\n",
    "        grid['maxy'] = aff[5]+aff[4]*grid['row']+aff[4]\n",
    "        \n",
    "        # Create polygon from affine bounds and row/col indices\n",
    "        grid['geometry'] = [Polygon([(grid.minx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.maxy[i]),\n",
    "                                   (grid.minx[i],grid.maxy[i])])\\\n",
    "                          for i in range(len(grid))]\n",
    "        \n",
    "        # Set the df as geo-df\n",
    "        grid = gpd.GeoDataFrame(grid, crs = 4326) \n",
    "\n",
    "        # Get dissolvement_key for dissolvement. \n",
    "        grid['row3'] = np.floor(grid['row']/(grid_size/100)).astype(int)\n",
    "        grid['col3'] = np.floor(grid['col']/(grid_size/100)).astype(int)\n",
    "        grid['dissolve_key'] = grid['row3'].astype(str) +'-'+ grid['col3'].astype(str)\n",
    "        \n",
    "        # Define a city's OSM area as Polygon.\n",
    "        geo_ls = gpd.GeoSeries(city_geo(cities_area[i].split(', ')).dissolve().geometry)\n",
    "        \n",
    "        # Intersect grids with the city boundary Polygon.\n",
    "        insec = grid.intersection(geo_ls.tolist()[0])\n",
    "        \n",
    "        # Exclude grids outside the specified city boundaries\n",
    "        insec = insec[insec.area > 0]\n",
    "        \n",
    "        # Join in other information.\n",
    "        insec = gpd.GeoDataFrame(geometry = insec, crs = 4326).join(grid.loc[:, grid.columns != 'geometry'])\n",
    "        \n",
    "        # Dissolve into block by block grids\n",
    "        popgrid = insec[['dissolve_key','geometry','row3','col3']].dissolve('dissolve_key')\n",
    "        \n",
    "        # Get those grids populations and area. Only blocks with population and full blocks\n",
    "        popgrid['population'] = round(insec.groupby('dissolve_key')['value'].sum()).astype(int)\n",
    "        popgrid['area_m'] = round(gpd.GeoSeries(popgrid['geometry'], crs = 4326).to_crs(3043).area).astype(int)\n",
    "        popgrid = popgrid[popgrid['population'] > 0]\n",
    "        popgrid = popgrid[popgrid['area_m'] / popgrid['area_m'].max() > 0.95]\n",
    "\n",
    "        # Get centroids and coords\n",
    "        popgrid['centroid'] = popgrid['geometry'].centroid\n",
    "        popgrid['centroid_m'] = gpd.GeoSeries(popgrid['centroid'], crs = 4326).to_crs(3043)\n",
    "        popgrid['grid_lon'] = popgrid['centroid_m'].x\n",
    "        popgrid['grid_lat'] = popgrid['centroid_m'].y\n",
    "        popgrid = popgrid.reset_index()\n",
    "\n",
    "        minx = popgrid.bounds['minx']\n",
    "        maxx = popgrid.bounds['maxx']\n",
    "        miny = popgrid.bounds['miny']\n",
    "        maxy = popgrid.bounds['maxy']\n",
    "\n",
    "        # Some geometries result in a multipolygon when dissolving (like i.e. 0.05 meters), coords error.\n",
    "        # Therefore recreate the polygon.\n",
    "        Poly = []\n",
    "        for k in range(len(popgrid)):\n",
    "            Poly.append(Polygon([(minx[k],maxy[k]),(maxx[k],maxy[k]),(maxx[k],miny[k]),(minx[k],miny[k])]))\n",
    "        popgrid['geometry'] = Poly\n",
    "        \n",
    "        try:\n",
    "            entry_index = [int(road_nodes[i]['geometry'].sindex.nearest(grid)[1])\\\n",
    "                                 for grid in popgrid['centroid']]\n",
    "        except:\n",
    "            entry_index = [int(road_nodes[i]['geometry'].sindex.nearest(grid)[1][0])\\\n",
    "                                 for grid in popgrid['centroid']]\n",
    "            \n",
    "        nearest_index = road_nodes[i].iloc[entry_index]\n",
    "        popgrid['grid_osm'] = nearest_index.reset_index(drop = True)['osmid_var']\n",
    "        popgrid['node_geom'] = nearest_index.reset_index(drop = True)['geometry']\n",
    "        popgrid['node_geom_m'] = nearest_index.reset_index(drop = True)['geometry_m']\n",
    "        popgrid['G-entry cost'] = popgrid['node_geom_m'].distance(popgrid['centroid_m'])\n",
    "        \n",
    "        UGS_all = UGS[i].dissolve().geometry[0]\n",
    "        popgrid['in_out_UGS'] = popgrid.intersection(UGS_all).is_empty == False\n",
    "        \n",
    "        grids.append(popgrid)\n",
    "\n",
    "        print(city_grids[i].rsplit('_')[3], round((time.time() - start_time)/60,2),'mns')\n",
    "    return(grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "59bebc3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5 park entry points\n",
    "def UGS_fake_entry(UGS, road_nodes, graphs, cities, pop_grids,\n",
    "                   thresholds, UGS_entry_buf = 25, walk_radius = 500, entry_point_merge = 0):\n",
    "    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": 5,
   "id": "371c0f0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def suitible_enhanced (UGS_entry, pop_grids, road_nodes, cities, thresholds):\n",
    "    start_time = time.time()\n",
    "    suits_all = []\n",
    "    for j in range(len(cities)):\n",
    "        print('get (Euclidean) suitible combinations')\n",
    "        print('0.0 %', round((time.time() - start_time) / 60,2),'mns')\n",
    "        UGSe = UGS_entry[j]\n",
    "        entry_geoms = UGSe.geometry_m\n",
    "        pop = pop_grids[j]\n",
    "        road_node = road_nodes[j]\n",
    "\n",
    "        suits = pd.DataFrame()\n",
    "        cols = ['osmid','Park_No','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": 6,
   "id": "6553bf33",
   "metadata": {},
   "outputs": [],
   "source": [
    "def obtaining_subgraphs(graphs, pop_grids, UGS_entry, nodes, cities, thresholds, time_sleep = 30):\n",
    "    print('obtain local graphs')\n",
    "    start_time = time.time()\n",
    "    subgraphs_all = []\n",
    "    suits_all = []\n",
    "    for j in range(len(cities)):\n",
    "        print(cities[j])\n",
    "        Graph = graphs[j]\n",
    "        pop = pop_grids[j]\n",
    "        UGSe = UGS_entry[j].sort_values('osmid')\n",
    "        road_node = nodes[j]\n",
    "        node_geoms = road_node.geometry_m\n",
    "        entry_geoms = UGSe.geometry_m\n",
    "        osmid = UGSe['osmid']\n",
    "\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(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: \n",
    "                print(round((i+1) / len(entry_geoms)*100,2),'% done',\n",
    "                                        round((time.time() - start_time) / 60,2),'mns')\n",
    "                time.sleep(time_sleep)\n",
    "        print('100 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        subgraphs_all.append(pd.Series(subgraphs, index = UGSe_ids))\n",
    "        suits_all.append(suits)\n",
    "    return({'graphs':subgraphs_all,'graph nodes':suits_all})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "aeab83a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def distance_fast (Geo_1, Geo_2):\n",
    "    return((abs(Geo_1.x - Geo_2.x)**2 + abs(Geo_1.y - Geo_2.y)**2).apply(math.sqrt))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "d879c1ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "def direct_routing (suitible_comb, graphs, edges, cities, chunk = 20000, time_sleep = 15):\n",
    "    start_time = time.time()\n",
    "    Routes = []\n",
    "    Lines = []\n",
    "    for j in enumerate(cities):\n",
    "        print(j[1])\n",
    "        \n",
    "        suitible = suitible_comb[j[0]].sort_values('Parkroad_osmid').reset_index()\n",
    "        grouped = suitible[suitible['in_out_UGS'] == False].groupby(['Parkroad_osmid'])['grid_osm'].apply(list)\n",
    "        sets = grouped.apply(np.unique)\n",
    "\n",
    "        Conn = edges[j[0]]\n",
    "        SG = graphs[j[0]]\n",
    "        SG = SG[sets.index]\n",
    "        \n",
    "        SGr = SG.reset_index()\n",
    "        SG = SGr.iloc[pd.Series(SGr['Parkroad_osmid'].drop_duplicates()).index].set_index('Parkroad_osmid')[0]\n",
    "\n",
    "        num = int(np.ceil(chunk / sets.apply(len).mean()))\n",
    "        length = int(np.ceil(len(suitible['Parkroad_osmid'].unique())/num))\n",
    "\n",
    "        Routes_df = pd.DataFrame()\n",
    "        Lines_df = pd.DataFrame()\n",
    "        for l in range(length):\n",
    "            comb = suitible[suitible['Parkroad_osmid'].isin(sets.index[l*num:l*num+num])]\n",
    "            sets2 = sets[l*num:l*num+num]\n",
    "\n",
    "            parknode = list(comb['Parkroad_osmid'])\n",
    "            gridnode = list(comb['grid_osm'])\n",
    "            subgraph = SG[sets2.index]\n",
    "            #sets2 = sets2[subgraph.index]\n",
    "\n",
    "            ls = []\n",
    "            ls2 = []\n",
    "            ls3 = []\n",
    "            lod = []\n",
    "            lgk = []\n",
    "            Routes\n",
    "            for i in range(len(sets2)):\n",
    "                path = nx.single_source_dijkstra(subgraph.iloc[i], sets2.index[i], weight = 'length')\n",
    "\n",
    "                incl = np.isin(list(path[0].keys()),sets2.iloc[i])\n",
    "                incl2 = np.isin(list(path[1].keys()),sets2.iloc[i])\n",
    "\n",
    "                # route cost\n",
    "                orig_c = list(np.repeat(sets2.index[i],sum(incl)))\n",
    "                dest_c = list(np.array(list(path[0].keys()))[incl])\n",
    "                cost = list(np.array(list(path[0].values()))[incl])\n",
    "\n",
    "                ls = ls + orig_c\n",
    "                ls2= ls2+ dest_c\n",
    "                ls3= ls3+ cost\n",
    "\n",
    "                # route steps\n",
    "                orig_s = list(np.repeat(sets2.index[i],sum(incl2)))\n",
    "                dest_s = list(np.array(list(path[1].keys()))[incl2])\n",
    "                steps = list(np.array(list(path[1].values()),dtype=object)[incl2])\n",
    "\n",
    "                fr = []\n",
    "                to = []\n",
    "                og = []\n",
    "                de = []\n",
    "                for j in enumerate(steps):\n",
    "                    if len(j[1]) > 1:\n",
    "                        fr.append(j[1][:-1])\n",
    "                        to.append(j[1][1:])\n",
    "                        og.append(list(np.repeat(orig_s[j[0]], len(j[1][:-1]))))\n",
    "                        de.append(list(np.repeat(dest_s[j[0]], len(j[1][:-1]))))\n",
    "                    else:\n",
    "                        pass\n",
    "\n",
    "                fr = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, fr) for i in b]\n",
    "                to = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, to) for i in b]\n",
    "                og = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, og) for i in b]\n",
    "                de = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, de) for i in b]\n",
    "\n",
    "                gk = [str(fr[k])+'-'+str(to[k]) for k in range(len(to))]\n",
    "                gkr = [str(to[k])+'-'+str(fr[k]) for k in range(len(to))]\n",
    "                od = [str(de[k])+'-'+str(og[k]) for k in range(len(og))]\n",
    "\n",
    "                lgk.append(gk)\n",
    "                lod.append(od)\n",
    "\n",
    "            dist_df = pd.DataFrame({'UGSe_id':ls,'GrE_id':ls2,'route cost':ls3})\n",
    "            dist_df['graph_key'] = dist_df['GrE_id'].astype(str)+'-'+dist_df['UGSe_id'].astype(str)\n",
    "\n",
    "            routes = pd.merge(comb, dist_df, on = 'graph_key', how = 'left')\n",
    "            routes['route cost'] = np.where(routes['in_out_UGS'],0,routes['route cost'])\n",
    "            routes = routes[~routes['route cost'].isna()].reset_index(drop = True)\n",
    "\n",
    "            routes['G-entry cost'] = np.where(routes['in_out_UGS'],0,routes['G-entry cost'])\n",
    "\n",
    "            routes['Tcost'] = routes['route cost']+routes['G-entry cost']\n",
    "\n",
    "            lgk = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, lgk) for i in b]\n",
    "            lod = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, lod) for i in b]\n",
    "            \n",
    "            linestr = pd.DataFrame({'route no':lod,'route step':lgk})\n",
    "            \n",
    "            linestr = pd.merge(linestr, Conn.geometry, left_on = 'route step', right_index = True, how = 'left')\n",
    "            linestr = linestr[['route no','geometry']]\n",
    "            linestr = gpd.GeoDataFrame(linestr[['route no','geometry']], crs = 4326)\n",
    "            \n",
    "            linestr = linestr.dissolve('route no')\n",
    "            routes2 = pd.merge(routes, linestr, left_on = 'graph_key', right_index = True, how = 'left')\n",
    "            \n",
    "            Lines_df = pd.concat([Lines_df, linestr])\n",
    "            Routes_df = pd.concat([Routes_df, routes2])\n",
    "            \n",
    "            print(round(l*num / len(sets)*100,2),'%', \n",
    "                  round((time.time() - start_time) / 60,2),'mns')\n",
    "            time.sleep(time_sleep)\n",
    "        Routes_df = Routes_df.sort_values('index')\n",
    "        Routes_df = Routes_df.set_index('index')\n",
    "        Routes_df = Routes_df.reset_index(drop = True)\n",
    "        \n",
    "        Routes_df = Routes_df[Routes_df.columns[~Routes_df.columns.isin(['UGSe_id', 'GrE_id'])]]\n",
    "        \n",
    "        print('100 % done',round((time.time() - start_time) / 60,2),'mns')\n",
    "        \n",
    "        Routes.append(Routes_df)\n",
    "        Lines.append(Lines_df)\n",
    "    return(Routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "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": 10,
   "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": "75f67806",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c28d55c7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b05a4e52",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "add37f56",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}