E2SFCA GEE-OSM Model-Func 2.1.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "442b3921",
   "metadata": {},
   "outputs": [],
   "source": [
    "# system packages\n",
    "import sys\n",
    "import time\n",
    "import warnings\n",
    "import os\n",
    "\n",
    "# non-geo numeric packages\n",
    "import numpy as np\n",
    "import math\n",
    "from itertools import product, combinations\n",
    "import pandas as pd\n",
    "\n",
    "# network and OSM packages\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "city_geo = ox.geocoder.geocode_to_gdf\n",
    "\n",
    "# Earth engine packages\n",
    "import ee\n",
    "import geemap\n",
    "\n",
    "# General geo-packages\n",
    "import libpysal\n",
    "import rasterio\n",
    "import geopandas as gpd\n",
    "import shapely\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "07277cdb",
   "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=-yXf2x-oYzjj-AVDnnwX8tUZGY-TxBhGXlbUP3HH5V8&tc=ybVo6Qlr0ZS5a2PN30VWi96zi1ZwLlx20Zy8f3qR-Ww&cc=ZepHQLeOdjdy0U8EfSNkK4BDd6wqpWAsMSxnn3fdq2w>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=-yXf2x-oYzjj-AVDnnwX8tUZGY-TxBhGXlbUP3HH5V8&tc=ybVo6Qlr0ZS5a2PN30VWi96zi1ZwLlx20Zy8f3qR-Ww&cc=ZepHQLeOdjdy0U8EfSNkK4BDd6wqpWAsMSxnn3fdq2w</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/1AVHEtk54-OVQYSEUi5UwQ43Qi96mWaIh3o6c8JPh0abGhJKU391VvjVfm1w\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": 21,
   "id": "47f4b61b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 0 cities and thresholds\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m\n",
    "\n",
    "# Extract iso-3166 country codes\n",
    "iso = pd.read_excel('iso_countries.xlsx')\n",
    "\n",
    "# Extract cities list\n",
    "cities = pd.read_excel('cities.xlsx') # all cities\n",
    "\n",
    "# 'cities_adj' serves by default as city-input for functions\n",
    "# cities_adj = cities\n",
    "# cities_adj = cities[cities['Included (Y/N)'] == 'Y']\n",
    "cities_adj = cities[cities['City'].isin(['Indore','Toronto','Lagos','Belo Horizonte','Kuala Lumpur',\\\n",
    "                                         'Philadelphia','Casablanca','Bologna','Agra','Hanoi'])]\n",
    "cities_adj = cities_adj.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "fa5e18f7",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/9639064b3ba36889c3d41f642545c48b-465ea5e5d59bc2878890c813e42b91ad:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\IND_Agra_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/df7170e30e9c2a8c44d3a0aa48d87e91-d2408cdb856a6fb1d4305f58242638fa:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\BRA_Belo Horizonte_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/5d0f17f70b75ee205c22dac890f7dd5f-06d47fe65fdcc50c4bcfcd7a84043c2e:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\ITA_Bologna_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/78a904632c412a448ef0e05e43131273-32caedddf26eaf2a5629534ace89cac1:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\MAR_Casablanca_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/e2dcaa7fc8ef07556ab7052aba2ca3a2-c973fc79ec1bf83ea24bb19b34c0b1e9:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\VNM_Hanoi_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/c3f5f28a194d2dd3f6d9013a230bfb38-ba249ae1b5789ecfd518ceae0b5fef53:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\IND_Indore_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/a6ebc706049fdb30294d61c63e9050de-7ef3df24d1a0f8354c33fc48a513386b:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\PAK_Islamabad_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/7e0cbe6efebd974c1e109a103daac86d-1c935046545d1521ffaee0f4549a033d:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\NGA_Lagos_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/daee5df0df5eba2b0539dcc485cbfbf5-6b3f08d855cabc7345ea5d9336f59e02:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\USA_Philadelphia_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/aec17c02897877f7969726425d40f20b-536698a5078c484588e9076767981e83:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to D:\\Dumps\\GEE_city_grids\\CAN_Toronto_2020.tif\n",
      "CPU times: total: 1.66 s\n",
      "Wall time: 28.4 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 1. Required preprocess for information extraction\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# In essence, we use Google Earth Engine to extract a country's grid raster and clip it with the city's preferred OSM area\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:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "df72d9e2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Agra 0.47 mns\n",
      "Belo Horizonte 1.12 mns\n",
      "Bologna 1.32 mns\n",
      "Casablanca 1.6 mns\n",
      "Hanoi 1.88 mns\n",
      "Indore 2.82 mns\n",
      "Islamabad 3.79 mns\n",
      "Lagos 6.8 mns\n",
      "Philadelphia 7.52 mns\n",
      "Toronto 8.3 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Agra done 0.73 mns\n",
      "Belo Horizonte done 1.63 mns\n",
      "Bologna done 2.15 mns\n",
      "Casablanca done 2.97 mns\n",
      "Hanoi done 6.86 mns\n",
      "Indore done 8.39 mns\n",
      "Islamabad done 9.33 mns\n",
      "Lagos done 11.51 mns\n",
      "Philadelphia done 13.43 mns\n",
      "Toronto done 15.9 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Agra done\n",
      "Belo Horizonte done\n",
      "Bologna done\n",
      "Casablanca done\n",
      "Hanoi done\n",
      "Indore done\n",
      "Islamabad done\n",
      "Lagos done\n",
      "Philadelphia done\n",
      "Toronto done\n",
      "CPU times: total: 24min 10s\n",
      "Wall time: 24min 42s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 2. Information extraction\n",
    "\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(files,\n",
    "                                     cities_adj['OSM_area'],\n",
    "                                     grid_size = 100) # aggregating upwards to i.e. 200m, 300m etc. is possible\n",
    "print(' ')\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_networks(cities_adj, # Get 'all' (drive,walk,bike) network\n",
    "                              thresholds,\n",
    "                              undirected = True)\n",
    "print(' ')\n",
    "\n",
    "# Extract urban greenspace (UGS)\n",
    "UGS = urban_greenspace(cities_adj, \n",
    "                       thresholds,\n",
    "                       one_UGS_buf = 25, # buffer at which UGS is seen as one\n",
    "                       min_UGS_size = 400) # WHO sees this as minimum UGS size (400m2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a673017",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get fake UGS entry points\n",
      "Agra 0.0 % done 0.02  mns\n",
      "Agra 100 % done 0.88  mns\n",
      "Belo Horizonte 0.0 % done 1.0  mns\n",
      "Belo Horizonte 16.6 % done 1.63  mns\n",
      "Belo Horizonte 33.2 % done 2.25  mns\n",
      "Belo Horizonte 49.8 % done 2.71  mns\n",
      "Belo Horizonte 66.3 % done 3.07  mns\n",
      "Belo Horizonte 82.9 % done 3.53  mns\n",
      "Belo Horizonte 99.5 % done 3.93  mns\n",
      "Belo Horizonte 100 % done 3.97  mns\n",
      "Bologna 0.0 % done 3.98  mns\n",
      "Bologna 29.3 % done 4.3  mns\n",
      "Bologna 58.7 % done 4.56  mns\n",
      "Bologna 88.0 % done 4.82  mns\n",
      "Bologna 100 % done 4.92  mns\n",
      "Casablanca 0.0 % done 4.94  mns\n",
      "Casablanca 52.9 % done 6.09  mns\n",
      "Casablanca 100 % done 6.41  mns\n",
      "Hanoi 0.0 % done 6.47  mns\n",
      "Hanoi 29.4 % done 7.16  mns\n",
      "Hanoi 58.8 % done 7.88  mns\n",
      "Hanoi 88.2 % done 8.45  mns\n",
      "Hanoi 100 % done 8.72  mns\n",
      "Indore 0.0 % done 8.72  mns\n",
      "Indore 19.3 % done 9.15  mns\n",
      "Indore 38.5 % done 9.51  mns\n",
      "Indore 57.8 % done 9.93  mns\n",
      "Indore 77.1 % done 10.35  mns\n",
      "Indore 96.3 % done 10.7  mns\n",
      "Indore 100 % done 10.78  mns\n",
      "Islamabad 0.0 % done 11.06  mns\n",
      "Islamabad 30.9 % done 11.75  mns\n",
      "Islamabad 61.7 % done 12.14  mns\n",
      "Islamabad 92.6 % done 12.57  mns\n",
      "Islamabad 100 % done 12.68  mns\n",
      "Lagos 0.0 % done 12.78  mns\n",
      "Lagos 85.5 % done 13.73  mns\n",
      "Lagos 100 % done 13.94  mns\n",
      "Philadelphia 0.0 % done 13.97  mns\n",
      "Philadelphia 12.3 % done 14.9  mns\n",
      "Philadelphia 24.6 % done 15.8  mns\n",
      "Philadelphia 36.9 % done 16.64  mns\n",
      "Philadelphia 49.3 % done 17.53  mns\n",
      "Philadelphia 61.6 % done 18.3  mns\n",
      "Philadelphia 73.9 % done 19.12  mns\n",
      "Philadelphia 86.2 % done 19.96  mns\n",
      "Philadelphia 98.5 % done 20.74  mns\n",
      "Philadelphia 100 % done 21.02  mns\n",
      "Toronto 0.0 % done 21.04  mns\n",
      "Toronto 8.9 % done 22.2  mns\n",
      "Toronto 17.9 % done 23.39  mns\n",
      "Toronto 26.8 % done 24.49  mns\n",
      "Toronto 35.7 % done 25.54  mns\n",
      "Toronto 44.7 % done 26.65  mns\n",
      "Toronto 53.6 % done 27.63  mns\n",
      "Toronto 62.6 % done 28.71  mns\n",
      "Toronto 71.5 % done 29.75  mns\n",
      "Toronto 80.4 % done 30.81  mns\n",
      "Toronto 89.4 % done 31.76  mns\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# 3. Preprocess information for route finding\n",
    "\n",
    "# Get fake entry points (between UGS and buffer limits)\n",
    "UGS_entry = UGS_fake_entry(UGS, \n",
    "                           road_networks['nodes'], \n",
    "                           cities_adj['City'],\n",
    "                           UGS_entry_buf = 25, # road nodes within 25 meters are seen as fake entry points\n",
    "                           walk_radius = 500, # assume that the average person only views a UGS up to 500m in radius\n",
    "                                                # more attractive\n",
    "                           entry_point_merge = 0) # merges closeby fake UGS entry points within X meters \n",
    "                                                    # what may be done for performance\n",
    "print(' ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5a4a064",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Checks all potential suitible combinations (points that fall within max threshold Euclidean distance from the ego)\n",
    "suitible = suitible_combinations(UGS_entry, \n",
    "                                 population_grids, \n",
    "                                 road_networks['nodes'], # For finding nearest grid entry points\n",
    "                                 thresholds,\n",
    "                                 cities_adj['City'],\n",
    "                                 chunk_size = 10000000) # calculating per chunk of num UGS entry points * num pop_grids\n",
    "                                                        # Preventing normal PC meltdown, set lower if PC gets stuck\n",
    "print(' ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c16479f2",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# Checks if grids are already in a UGS\n",
    "suitible_InOut_UGS = grids_in_UGS (suitible, UGS, population_grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b98ff692",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# 4. Finding shortest routes.\n",
    "Routes = route_finding (road_networks['graphs'], # graphs of the road networks\n",
    "               suitible_InOut_UGS, # potential suitible routes with grid-UGS comb. separated in or out UGS.\n",
    "               road_networks['nodes'], \n",
    "               road_networks['edges'], \n",
    "               cities_adj['City'], \n",
    "               block_size = 250000, # Chunk to spread dataload.\n",
    "               nn_iter = 10) # max amount of nearest nodes to be found (both for UGS entry and grid-centroid road entries)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a49890ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# 5. summarize scores\n",
    "min_gridUGS = min_gridUGS_comb (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/', \n",
    "                             grid_size = 100,\n",
    "                             ext = '_4_5')\n",
    "\n",
    "E2SCFA_score['score summary']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "313b0bd2",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gee_worldpop_extract (city_file, iso, save_path = None):\n",
    "    \n",
    "    cities = city_file\n",
    "    \n",
    "    # Get included city areas\n",
    "    OSM_incl = [cities[cities['City'] == city]['OSM_area'].tolist()[0].rsplit(', ') for city in cities['City'].tolist()]\n",
    "\n",
    "    # Get the city geoms\n",
    "    obj = [city_geo(city).dissolve()['geometry'].tolist()[0] for city in OSM_incl]\n",
    "\n",
    "    # Get the city countries\n",
    "    obj_displ = [city_geo(city).dissolve()['display_name'].tolist()[0].rsplit(', ')[-1]for city in OSM_incl]\n",
    "    obj_displ = np.where(pd.Series(obj_displ).str.contains(\"Ivoire\"),\"CIte dIvoire\",obj_displ)\n",
    "\n",
    "    # Get the country's iso-code\n",
    "    iso_list = [iso[iso['name'] == ob]['alpha3'].tolist()[0] for ob in obj_displ]\n",
    "\n",
    "    # Based on the iso-code return the worldpop 2020\n",
    "    ee_worldpop = [ee.ImageCollection(\"WorldPop/GP/100m/pop\")\\\n",
    "        .filter(ee.Filter.date('2020'))\\\n",
    "        .filter(ee.Filter.inList('country', [io])).first() for io in iso_list]\n",
    "\n",
    "    # Clip the countries with the city geoms.\n",
    "    clipped = [ee_worldpop[i].clip(shapely.geometry.mapping(obj[i])) for i in range(0,len(obj))]\n",
    "\n",
    "    # Create path if non-existent\n",
    "    if save_path == None:\n",
    "        path = ''\n",
    "    else:\n",
    "        path = save_path\n",
    "        if not os.path.exists(path):\n",
    "                    os.makedirs(path)\n",
    "\n",
    "    # Export as TIFF file.\n",
    "    # Stored in form path + USA_Los Angeles_2020.tif\n",
    "    filenames = [path+iso_list[i]+'_'+cities['City'][i]+'_2020.tif' for i in range(len(obj))]\n",
    "    [geemap.ee_export_image(clipped[i], filename = filenames[i]) for i in range(0,len(obj))]\n",
    "    return(filenames)\n",
    "    sys.stdout.flush()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e5a4c548",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 2 population grids extraction\n",
    "def city_grids_format(city_grids, cities_area, grid_size = 100):\n",
    "    start_time = time.time()\n",
    "    grids = []\n",
    "    print(str(grid_size) + 'm resolution grids extraction')\n",
    "    for i in range(len(city_grids)):\n",
    "        \n",
    "        # Open the raster file\n",
    "        with rasterio.open(city_grids[i]) as src:\n",
    "            band= src.read() # the population values\n",
    "            aff = src.transform # the raster bounds and size (affine)\n",
    "        \n",
    "        # Get the rowwise arrays, get a 2D dataframe\n",
    "        grid = pd.DataFrame()\n",
    "        for b in enumerate(band[0]):\n",
    "            grid = pd.concat([grid, pd.Series(b[1],name=b[0])],axis=1)\n",
    "        grid= grid.unstack().reset_index()\n",
    "        \n",
    "        # Unstack df to columns\n",
    "        grid.columns = ['row','col','value']\n",
    "        grid['minx'] = aff[2]+aff[0]*grid['col']\n",
    "        grid['miny'] = aff[5]+aff[4]*grid['row']\n",
    "        grid['maxx'] = aff[2]+aff[0]*grid['col']+aff[0]\n",
    "        grid['maxy'] = aff[5]+aff[4]*grid['row']+aff[4]\n",
    "        \n",
    "        # Create polygon from affine bounds and row/col indices\n",
    "        grid['geometry'] = [Polygon([(grid.minx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.miny[i]),\n",
    "                                   (grid.maxx[i],grid.maxy[i]),\n",
    "                                   (grid.minx[i],grid.maxy[i])])\\\n",
    "                          for i in range(len(grid))]\n",
    "        \n",
    "        # Set the df as geo-df\n",
    "        grid = gpd.GeoDataFrame(grid, crs = 4326) \n",
    "\n",
    "        # Get dissolvement_key for dissolvement. \n",
    "        grid['row3'] = np.floor(grid['row']/(grid_size/100)).astype(int)\n",
    "        grid['col3'] = np.floor(grid['col']/(grid_size/100)).astype(int)\n",
    "        grid['dissolve_key'] = grid['row3'].astype(str) +'-'+ grid['col3'].astype(str)\n",
    "        \n",
    "        # Define a city's OSM area as Polygon.\n",
    "        geo_ls = gpd.GeoSeries(city_geo(cities_area[i].split(', ')).dissolve().geometry)\n",
    "        \n",
    "        # Intersect grids with the city boundary Polygon.\n",
    "        insec = grid.intersection(geo_ls.tolist()[0])\n",
    "        \n",
    "        # Exclude grids outside the specified city boundaries\n",
    "        insec = insec[insec.area > 0]\n",
    "        \n",
    "        # Join in other information.\n",
    "        insec = gpd.GeoDataFrame(geometry = insec, crs = 4326).join(grid.loc[:, grid.columns != 'geometry'])\n",
    "        \n",
    "        # Dissolve into block by block grids\n",
    "        popgrid = insec[['dissolve_key','geometry','row3','col3']].dissolve('dissolve_key')\n",
    "        \n",
    "        # Get those grids populations and area. Only blocks with population and full blocks\n",
    "        popgrid['population'] = round(insec.groupby('dissolve_key')['value'].sum()).astype(int)\n",
    "        popgrid['area_m'] = round(gpd.GeoSeries(popgrid['geometry'], crs = 4326).to_crs(3043).area).astype(int)\n",
    "        popgrid = popgrid[popgrid['population'] > 0]\n",
    "        popgrid = popgrid[popgrid['area_m'] / popgrid['area_m'].max() > 0.95]\n",
    "\n",
    "        # Get centroids and coords\n",
    "        popgrid['centroid'] = popgrid['geometry'].centroid\n",
    "        popgrid['centroid_m'] = gpd.GeoSeries(popgrid['centroid'], crs = 4326).to_crs(3043)\n",
    "        popgrid['grid_lon'] = popgrid['centroid_m'].x\n",
    "        popgrid['grid_lat'] = popgrid['centroid_m'].y\n",
    "        popgrid = popgrid.reset_index()\n",
    "\n",
    "        minx = popgrid.bounds['minx']\n",
    "        maxx = popgrid.bounds['maxx']\n",
    "        miny = popgrid.bounds['miny']\n",
    "        maxy = popgrid.bounds['maxy']\n",
    "\n",
    "        # Some geometries result in a multipolygon when dissolving (like i.e. 0.05 meters), coords error.\n",
    "        # Therefore recreate the polygon.\n",
    "        Poly = []\n",
    "        for k in range(len(popgrid)):\n",
    "            Poly.append(Polygon([(minx[k],maxy[k]),(maxx[k],maxy[k]),(maxx[k],miny[k]),(minx[k],miny[k])]))\n",
    "        popgrid['geometry'] = Poly\n",
    "\n",
    "        grids.append(popgrid)\n",
    "\n",
    "        print(city_grids[i].rsplit('_')[3], round((time.time() - start_time)/60,2),'mns')\n",
    "    return(grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "5a8007ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 3 Road networks\n",
    "def road_networks (cities, thresholds, undirected = False):\n",
    "    print('get road networks from OSM')\n",
    "    start_time = time.time()\n",
    "    graphs = list()\n",
    "    road_nodes = list()\n",
    "    road_edges = list()\n",
    "    road_conn = list()\n",
    "\n",
    "    for i in enumerate(cities['OSM_area']):\n",
    "        # Get graph, road nodes and edges\n",
    "        road_node = pd.DataFrame()\n",
    "        roads = pd.DataFrame()\n",
    "        \n",
    "        # For each included OSM_area get the roads\n",
    "        for district in i[1].rsplit(', '):\n",
    "            graph = ox.graph_from_place(district, network_type = \"all\", buffer_dist = (np.max(thresholds)+1000))\n",
    "            node, edge = ox.graph_to_gdfs(graph)\n",
    "            road_node = pd.concat([road_node, node], axis = 0)\n",
    "            roads = pd.concat([roads, edge], axis = 0)\n",
    "        \n",
    "        # Eliminate lists in the df which prevents drop of duplicate columns\n",
    "        road_edge = pd.DataFrame([[c[0] if isinstance(c,list) else c for c in roads[col]]\\\n",
    "                              for col in roads]).transpose()\n",
    "        road_edge.columns = roads.columns\n",
    "        road_edge.index = roads.index\n",
    "        road_edge = gpd.GeoDataFrame(road_edge, crs = 4326)\n",
    "        \n",
    "        # Return the unique nodes and edges of the (often) adjacent OSM_areas.\n",
    "        road_node = road_node.drop_duplicates()\n",
    "        road_edge = road_edge.drop_duplicates()\n",
    "        \n",
    "        # Road nodes format\n",
    "        road_node = road_node.to_crs(4326)\n",
    "        road_node['geometry_m'] = gpd.GeoSeries(road_node['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_node['osmid_var'] = road_node.index\n",
    "        road_node = gpd.GeoDataFrame(road_node, geometry = 'geometry', crs = 4326)\n",
    "\n",
    "        # format road edges\n",
    "        road_edge['geometry_m'] = gpd.GeoSeries(road_edge['geometry'], crs = 4326).to_crs(3043)\n",
    "        road_edge = road_edge.reset_index()\n",
    "        road_edge.rename(columns={'u':'from', 'v':'to', 'key':'keys'}, inplace=True)\n",
    "        road_edge['key'] = road_edge['from'].astype(str) + '-' + road_edge['to'].astype(str)\n",
    "        \n",
    "        if undirected == True:\n",
    "            # Apply one-directional to both for walking\n",
    "            both = road_edge[road_edge['oneway'] == False]\n",
    "            one = road_edge[road_edge['oneway'] == True]\n",
    "            rev = pd.DataFrame()\n",
    "            rev[['from','to']] = one[['to','from']]\n",
    "            rev = pd.concat([rev,one.iloc[:,2:]],axis = 1)\n",
    "            edge_bidir = pd.concat([both, one, rev])\n",
    "            edge_bidir = edge_bidir.reset_index()\n",
    "            edge_bidir['oneway'] = False\n",
    "        else:\n",
    "            edge_bidir = road_edge\n",
    "\n",
    "        # Exclude highways and ramps on edges    \n",
    "        edge_filter = edge_bidir[(edge_bidir['highway'].str.contains('motorway') | \n",
    "              (edge_bidir['highway'].str.contains('trunk') & \n",
    "               edge_bidir['maxspeed'].astype(str).str.contains(\n",
    "                   '40 mph|45 mph|50 mph|55 mph|60 mph|65|70|75|80|85|90|95|100|110|120|130|140'))) == False]\n",
    "        road_edges.append(edge_filter)\n",
    "\n",
    "        # Exclude isolated nodes\n",
    "        fltrnodes = pd.Series(list(edge_filter['from']) + list(edge_filter['to'])).unique()\n",
    "        newnodes = road_node[road_node['osmid_var'].isin(fltrnodes)]\n",
    "        road_nodes.append(newnodes)\n",
    "\n",
    "        # Get only necessary road connections columns for network performance\n",
    "        road_con = edge_filter[['osmid','key','length','geometry']]\n",
    "        road_con = road_con.set_index('key')\n",
    "\n",
    "        road_conn.append(road_con)\n",
    "\n",
    "        # formatting to graph again.\n",
    "        newnodes = newnodes.loc[:, ~newnodes.columns.isin(['geometry_m', 'osmid_var'])]\n",
    "        edge_filter = edge_filter.set_index(['from','to','keys'])\n",
    "        edge_filter = edge_filter.loc[:, ~edge_filter.columns.isin(['geometry_m', 'key'])]\n",
    "\n",
    "        graph2 = ox.graph_from_gdfs(newnodes, edge_filter)\n",
    "\n",
    "        graphs.append(graph2)\n",
    "        print(cities['City'][i[0]].rsplit(',')[0], 'done', round((time.time() - start_time) / 60,2),'mns')\n",
    "    return({'graphs':graphs,'nodes':road_nodes,'edges':road_conn,'edges long':road_edges})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "de012d0d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 4 city greenspace\n",
    "def urban_greenspace (cities, thresholds, one_UGS_buf = 25, min_UGS_size = 400):\n",
    "    print('get urban greenspaces from OSM')\n",
    "    parks_in_range = list()\n",
    "    for i in enumerate(cities['OSM_area']):\n",
    "        # Tags seen as Urban Greenspace (UGS) require the following:\n",
    "        # 1. Tag represent an area\n",
    "        # 2. The area is outdoor\n",
    "        # 3. The area is (semi-)publically available\n",
    "        # 4. The area is likely to contain trees, grass and/or greenery\n",
    "        # 5. The area can reasonable be used for walking or recreational activities\n",
    "        tags = {'landuse':['allotments','forest','greenfield','village_green'],\\\n",
    "                'leisure':['garden','fitness_station','nature_reserve','park','playground'],\\\n",
    "                'natural':'grassland'}\n",
    "        gdf = ox.geometries_from_place(i[1].rsplit(', '),tags = tags,buffer_dist = np.max(thresholds))\n",
    "        gdf = gdf[(gdf.geom_type == 'Polygon') | (gdf.geom_type == 'MultiPolygon')]\n",
    "        greenspace = gdf.reset_index()    \n",
    "        warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "        green_buffer = gpd.GeoDataFrame(geometry = greenspace.to_crs(3043).buffer(one_UGS_buf).to_crs(4326))\n",
    "        greenspace['geometry_w_buffer'] = green_buffer\n",
    "        greenspace['geometry_w_buffer'] = gpd.GeoSeries(greenspace['geometry_w_buffer'], crs = 4326)\n",
    "        greenspace['geom buffer diff'] = greenspace['geometry_w_buffer'].difference(greenspace['geometry'])\n",
    "\n",
    "        # This function group components in itself that overlap (with the buffer set of 25 metres)\n",
    "        # https://stackoverflow.com/questions/68036051/geopandas-self-intersection-grouping\n",
    "        W = libpysal.weights.fuzzy_contiguity(greenspace['geometry_w_buffer'])\n",
    "        greenspace['components'] = W.component_labels\n",
    "        parks = greenspace.dissolve('components')\n",
    "\n",
    "        # Exclude parks below 0.04 ha.\n",
    "        parks = parks[parks.to_crs(3043).area > min_UGS_size]\n",
    "        print(cities['City'][i[0]], 'done')\n",
    "        parks = parks.reset_index()\n",
    "        parks['geometry_m'] = parks['geometry'].to_crs(3043)\n",
    "        parks['park_area'] = parks['geometry_m'].area\n",
    "        parks_in_range.append(parks)\n",
    "    return(parks_in_range)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "06ac19c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5 park entry points\n",
    "def UGS_fake_entry(UGS, road_nodes, cities, UGS_entry_buf = 25, walk_radius = 500, entry_point_merge = 0):\n",
    "    print('get fake UGS entry points')\n",
    "    start_time = time.time()\n",
    "    ParkRoads = list()\n",
    "    for j in range(len(cities)):\n",
    "        ParkRoad = pd.DataFrame()\n",
    "        mat = list()\n",
    "        # For all\n",
    "        for i in range(len(UGS[j])):\n",
    "            dist = road_nodes[j]['geometry'].to_crs(3043).distance(UGS[j]['geometry'].to_crs(\n",
    "                3043)[i])\n",
    "            buf_nodes = road_nodes[j][(dist < UGS_entry_buf) & (dist > 0)]\n",
    "            mat.append(list(np.repeat(i, len(buf_nodes))))\n",
    "            ParkRoad = pd.concat([ParkRoad, buf_nodes])\n",
    "            if i % 100 == 0: print(cities[j].rsplit(',')[0], round(i/len(UGS[j])*100,1),'% done', \n",
    "                                  round((time.time() - start_time) / 60,2),' mns')\n",
    "        # Park no list conversion\n",
    "        mat_u = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat) for i in b]\n",
    "\n",
    "        # Format\n",
    "        ParkRoad['Park_No'] = mat_u\n",
    "        ParkRoad = ParkRoad.reset_index()\n",
    "        ParkRoad['park_lon'] = ParkRoad['geometry_m'].x\n",
    "        ParkRoad['park_lat'] = ParkRoad['geometry_m'].y\n",
    "        \n",
    "        # Get the road nodes intersecting with the parks' buffer\n",
    "        ParkRoad = pd.merge(ParkRoad, UGS[j][['geometry','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",
    "    return(ParkRoads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0c537f03",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 5.5 (not in use, buffer is 0, thus retains all the park entry points as is)\n",
    "def simplify_UGS_entry(fake_UGS_entry, entry_point_merge = 0):\n",
    "    # Get buffer of nodes close to each other.\n",
    "    # Get the buffer\n",
    "    ParkComb = fake_UGS_entry\n",
    "    ParkComb['geometry_m_buffer'] = ParkComb['geometry_m'].buffer(entry_point_merge)\n",
    "\n",
    "    # Get and merge components\n",
    "    M = libpysal.weights.fuzzy_contiguity(ParkComb['geometry_m_buffer'])\n",
    "    ParkComb['components'] = M.component_labels\n",
    "\n",
    "    # Take centroid of merged components\n",
    "    centr = gpd.GeoDataFrame(ParkComb, geometry = 'geometry_x', crs = 4326).dissolve('components')['geometry_x'].centroid\n",
    "    centr = gpd.GeoDataFrame(centr)\n",
    "    centr.columns = ['comp_centroid']\n",
    "\n",
    "    # Get node closest to the centroid of all merged nodes, which accesses the road network.\n",
    "    ParkComb = pd.merge(ParkComb, centr, left_on = 'components', right_index = True)\n",
    "    ParkComb['centr_dist'] = ParkComb['geometry_x'].distance(ParkComb['comp_centroid'])\n",
    "    ParkComb = ParkComb.iloc[ParkComb.groupby('components')['centr_dist'].idxmin()]\n",
    "    return(ParkComb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e354607f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 6 grid-parkentry combinations within euclidean threshold distance\n",
    "def suitible_combinations(UGS_entry, pop_grids, road_nodes, thresholds, cities, chunk_size = 10000000):\n",
    "    print('get potential (Euclidean) suitible combinations')\n",
    "    start_time = time.time()\n",
    "    RoadComb = list()\n",
    "    for l in range(len(cities)):\n",
    "        print(cities[l])\n",
    "        len1 = len(pop_grids[l])\n",
    "        len2 = len(UGS_entry[l])\n",
    "\n",
    "        # Reduce the size of combinations per iteration\n",
    "        len4 = 1\n",
    "        len5 = len1 * len2\n",
    "        blockC = len5\n",
    "        while blockC > chunk_size:\n",
    "            blockC = len5 / len4\n",
    "            len4 = len4+1\n",
    "\n",
    "        # Amount of grids taken per iteration block\n",
    "        block = round(len1 / len4)\n",
    "\n",
    "        output = pd.DataFrame()\n",
    "        # Checking all the combinations at once is too performance intensive, it is broken down per 1000 (or what you want)\n",
    "        for i in range(len4):\n",
    "            try:\n",
    "                # Check all grid-park combinations per block\n",
    "                l1, l2 = range(i*block,(i+1)*block), range(0,len2)\n",
    "                listed = pd.DataFrame(list(product(l1, l2)))\n",
    "\n",
    "                # Merge grid and park information\n",
    "                grid_merged = pd.merge(listed, \n",
    "                                       pop_grids[l][['grid_lon','grid_lat','centroid','centroid_m']],\n",
    "                                       left_on = 0, right_index = True)\n",
    "                node_merged = pd.merge(grid_merged, \n",
    "                                       UGS_entry[l][['Park_No','osmid','geometry_x','geometry_y','geometry_m','park_lon','park_lat',\n",
    "                                           'share_walked','park_area','walk_area']], \n",
    "                                       left_on = 1, right_index = True)\n",
    "\n",
    "                # Preset index for merging\n",
    "                node_merged['key'] = range(0,len(node_merged))\n",
    "                node_merged = node_merged.set_index('key')\n",
    "                node_merged = node_merged.loc[:, ~node_merged.columns.isin(['index'])]\n",
    "\n",
    "                # Create lists for better computational performance\n",
    "                glon = list(node_merged['grid_lon'])\n",
    "                glat = list(node_merged['grid_lat'])\n",
    "                plon = list(node_merged['park_lon'])\n",
    "                plat = list(node_merged['park_lat'])\n",
    "\n",
    "                # Get the euclidean distances\n",
    "                mat = list()\n",
    "                for j in range(len(node_merged)):\n",
    "                    mat.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2))\n",
    "\n",
    "                # Check if distances are within 1000m and join remaining info and concat in master df per 1000.\n",
    "                mat_df = pd.DataFrame(mat)[(np.array(mat) <= np.max(thresholds))]\n",
    "\n",
    "                # join the other gravity euclidean scores and other information\n",
    "                mat_df.columns = ['Euclidean']    \n",
    "                mat_df = mat_df.join(node_merged)\n",
    "\n",
    "                output = pd.concat([output, mat_df])\n",
    "                print('in chunk',(i+1),'/',len4,len(mat_df),'suitible comb.')\n",
    "            except:\n",
    "                print('chunk',(i+1),'/',len4,'out of range')\n",
    "                pass\n",
    "            \n",
    "            \n",
    "        # Renaming columns\n",
    "        print('total combinations within distance',len(output))\n",
    "\n",
    "        output.columns = ['Euclidean','Grid_No','Park_entry_No','grid_lon','grid_lat','Grid_coords_centroid','Grid_m_centroid',\n",
    "                      'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid','park_lon',\n",
    "                      'park_lat','parkshare_walked','park_area','walk_area_m2']\n",
    "\n",
    "        output = output[['Euclidean','Grid_No','Park_entry_No','Grid_coords_centroid','Grid_m_centroid','walk_area_m2',\n",
    "                     'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid','park_area']]\n",
    "\n",
    "        # Reinstate geographic elements\n",
    "        output = gpd.GeoDataFrame(output, geometry = 'Grid_coords_centroid', crs = 4326)\n",
    "        output['Grid_m_centroid'] = gpd.GeoSeries(output['Grid_m_centroid'], crs = 3043)\n",
    "        output['Parkroad_coords_centroid'] = gpd.GeoSeries(output['Parkroad_coords_centroid'], crs = 4326)\n",
    "        output['Parkroad_m_centroid'] = gpd.GeoSeries(output['Parkroad_m_centroid'], crs = 3043)\n",
    "\n",
    "        # Get the nearest entrance point for the grid centroids\n",
    "        output = gridroad_entry(output, road_nodes[l])\n",
    "\n",
    "        print('100 % gridentry done', round((time.time() - start_time) / 60,2),' mns')\n",
    "        RoadComb.append(output)\n",
    "    return (RoadComb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6058fb80",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gridroad_entry (suitible_comb, road_nodes):    \n",
    "    start_time = time.time()\n",
    "    mat5 = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        try:\n",
    "            nearest = int(road_nodes['geometry'].sindex.nearest(suitible_comb['Grid_coords_centroid'].iloc[i])[1])\n",
    "            mat5.append(road_nodes['osmid_var'].iloc[nearest])\n",
    "        except: \n",
    "            # sometimes two nodes are the exact same distance, then the first in the list is taken.\n",
    "            nearest = int(road_nodes['geometry'].sindex.nearest(suitible_comb['Grid_coords_centroid'].iloc[i])[1][0])\n",
    "            mat5.append(road_nodes['osmid_var'].iloc[nearest])\n",
    "        if i % 250000 == 0: print(round(i/len(suitible_comb)*100,1),'% gridentry done', round((time.time() - start_time) / 60,2),' mns')\n",
    "    # format resulting dataframe\n",
    "    suitible_comb['grid_osm'] = mat5\n",
    "    suitible_comb = pd.merge(suitible_comb, road_nodes['geometry'], left_on = 'grid_osm', right_index = True)\n",
    "    suitible_comb['geometry_m'] = gpd.GeoSeries(suitible_comb['geometry'], crs = 4326).to_crs(3043)\n",
    "    suitible_comb = suitible_comb.reset_index()\n",
    "    return(suitible_comb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "6af9d8e8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def grids_in_UGS (suitible_comb, UGS, pop_grid): \n",
    "    print('grids in UGS')\n",
    "    start_time = time.time()\n",
    "    RoadInOut = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        gridUGS = pop_grid[i]['centroid'].intersection(UGS[i].dissolve().geometry[0]).is_empty == False\n",
    "        gridUGS.name = 'in_out_UGS'\n",
    "        merged = pd.merge(suitible_comb[i], gridUGS, left_on = 'Grid_No', right_index = True)\n",
    "        RoadInOut.append(merged)\n",
    "        print(i)\n",
    "    return(RoadInOut) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "cb7b21c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check grids in or out of UGS\n",
    "def GIU (suitible_comb, UGS, pop_grid): \n",
    "    start_time = time.time()\n",
    "    RoadInOut = list()\n",
    "    for i in range(len(suitible_comb)):\n",
    "        UGS_geoms = UGS[i]['geometry'].to_crs(4326)\n",
    "        grid = pop_grid[i]['centroid']\n",
    "        lst = list()\n",
    "        print('Check grids within UGS')\n",
    "        for l in enumerate(UGS_geoms):\n",
    "            lst.append(grid.intersection(l[1]).is_empty == False)\n",
    "            if l[0] % 100 == 0: print(l[0], round((time.time() - start_time) / 60,2),' mns')\n",
    "        print(np.array(lst).shape)\n",
    "        dfGrUGS = pd.DataFrame(pd.DataFrame(np.array(lst)).unstack())\n",
    "        dfGrUGS = dfGrUGS.reset_index()\n",
    "        dfGrUGS.columns = ['Grid_No','Park_No', 'in_out_UGS']\n",
    "        dfGrUGS['Gridpark_No'] = dfGrUGS['Grid_No'].astype(str)+\"-\"+dfGrUGS['Park_No'].astype(str)\n",
    "        dfGrUGS = dfGrUGS[['Gridpark_No','in_out_UGS']]\n",
    "        dfGrUGS = dfGrUGS.set_index('Gridpark_No')\n",
    "        suitible_comb[i]['Gridpark_No'] = suitible_comb[i]['Grid_No'].astype(str)+\"-\"+suitible_comb[i]['Park_No'].astype(str)\n",
    "        suitible_comb[i] = suitible_comb[i].set_index('Gridpark_No')\n",
    "        print(dfGrUGS.columns)\n",
    "        suitible_comb[i].join(dfGrUGS, how = 'left')\n",
    "        merged = pd.merge(suitible_comb[i], dfGrUGS, on = ['Gridpark_No'], how = 'left')\n",
    "        RoadInOut.append(merged)\n",
    "    return(RoadInOut)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "6d5b61e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 7 calculate route networks of all grid-parkentry combinations within euclidean threshold distance\n",
    "def route_finding (graphs, combinations, road_nodes, road_edges, cities, block_size = 250000, nn_iter = 10):\n",
    "\n",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    start_time = time.time()\n",
    "\n",
    "    Routes = list()\n",
    "    Routes_detail = list()\n",
    "    for j in range(len(cities)):\n",
    "        Graph = graphs[j]\n",
    "        suit_raw = combinations[j] # iloc to test the iteration speed.\n",
    "        nodes = road_nodes[j]\n",
    "\n",
    "        In_UGS = suit_raw[suit_raw['in_out_UGS'] == True] # Check if a grid centroid is in an UGS\n",
    "        suitible = suit_raw[suit_raw['in_out_UGS'] == False].reset_index(drop = True) # recreate a subsequential index\n",
    "                                                                                      # for the other grids outside UGS\n",
    "        block = block_size # Execute with chunks for performance improvement.\n",
    "\n",
    "        Route_parts = pd.DataFrame()\n",
    "        Route_dparts = pd.DataFrame()\n",
    "        len2 = int(np.ceil(len(suitible)/block))\n",
    "        # Divide in chunks of block for computational load\n",
    "        for k in range(len2):    \n",
    "            suitible_chunk = suitible.iloc[k*block:k*block+block] # Select chunk\n",
    "\n",
    "            parknode = list(suitible_chunk['Parkroad_osmid'])\n",
    "            gridnode = list(suitible_chunk['grid_osm'])\n",
    "\n",
    "            s_mat = list([]) # origin (normally grid) osmid\n",
    "            s_mat1 = list([]) # destination (normally UGS) osmid\n",
    "            s_mat2 = list([]) # route id\n",
    "            s_mat3 = list([]) # step id\n",
    "            s_mat4 = list([]) # way calculated\n",
    "            s_mat5 = list([]) # way calculated id\n",
    "            mat_nn = [] # found nearest nodes by block\n",
    "            len1 = len(suitible_chunk)\n",
    "\n",
    "            print(cities[j].rsplit(',')[0], k+1,'/',len2,'range',k*block,'-',k*block+np.where(k*block+block >= len1,len1,block))\n",
    "            for i in range(len(suitible_chunk)):\n",
    "                try: \n",
    "                    # from grid to UGS.\n",
    "                    shortest = nx.shortest_path(Graph, gridnode[i], parknode[i], 'travel_dist', method = 'dijkstra')\n",
    "                    s_mat.append(shortest)\n",
    "                    shortest_to = list(shortest[1:len(shortest)])\n",
    "                    shortest_to.append(-1)\n",
    "                    s_mat1.append(shortest_to)\n",
    "                    s_mat2.append(list(np.repeat(i+block*k, len(shortest))))\n",
    "                    s_mat3.append(list(np.arange(0, len(shortest))))\n",
    "                    s_mat4.append('normal way')\n",
    "                    s_mat5.append(1)\n",
    "                except:\n",
    "                    try:\n",
    "                        # Check the reverse\n",
    "                        shortest = nx.shortest_path(Graph, parknode[i], gridnode[i], 'travel_dist', method = 'dijkstra')\n",
    "                        s_mat.append(shortest)\n",
    "                        shortest_to = list(shortest[1:len(shortest)])\n",
    "                        shortest_to.append(-1)\n",
    "                        s_mat1.append(shortest_to)\n",
    "                        s_mat2.append(list(np.repeat(i+block*k, len(shortest))))\n",
    "                        s_mat3.append(list(np.arange(0, len(shortest))))\n",
    "                        s_mat4.append('reverse way')\n",
    "                        s_mat5.append(0)\n",
    "                    except:\n",
    "                        # Otherwise find nearest nodes (grid and UGS) and try to find routes between them\n",
    "                        nn_route_finding(Graph, suitible_chunk, nodes, s_mat, s_mat1, s_mat2, s_mat3,\n",
    "                                             s_mat4, s_mat5, mat_nn, i, block, k, nn_iter)\n",
    "                        \n",
    "                if i % 10000 == 0: print(round((i+block*k)/len(suitible)*100,2),'% done',\n",
    "                                         round((time.time() - start_time) / 60,2),'mns')\n",
    "            print('for', len(mat_nn),'routes nearest nodes found')\n",
    "\n",
    "            print(round((i+block*k)/len(suitible)*100,2),'% pathfinding done', round((time.time() - start_time) / 60,2),'mns')\n",
    "\n",
    "            # Formats route information by route and step (detailed)\n",
    "            routes = route_formatting(s_mat, s_mat1, s_mat2, s_mat3, road_edges[j]) # Formats lists to routes detail.\n",
    "            print('formatting done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            # Summarizes information by route\n",
    "            routes2 = route_summarization(routes, suitible_chunk, road_nodes[j], s_mat4, s_mat5) # formats routes to summary\n",
    "            print('dissolving done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            Route_parts = pd.concat([Route_parts, routes2])\n",
    "            Route_dparts = pd.concat([Route_dparts, routes])\n",
    "\n",
    "        # Format grids in UGS to enable smooth df concat\n",
    "        In_UGS = In_UGS.set_geometry(In_UGS['Grid_coords_centroid'])\n",
    "        In_UGS = In_UGS[['geometry','Grid_No','grid_osm','Park_No','Park_entry_No','Parkroad_osmid',\n",
    "                                   'Grid_m_centroid','walk_area_m2',\n",
    "                                   'Euclidean','geometry_m']]\n",
    "\n",
    "        In_UGS['realG_osmid'] = suit_raw['Parkroad_osmid']\n",
    "        In_UGS['realP_osmid'] = suit_raw['grid_osm']\n",
    "        In_UGS['way_calc'] = 'grid in UGS'\n",
    "\n",
    "        Route_parts = pd.concat([Route_parts,In_UGS])\n",
    "        Route_parts = Route_parts.reset_index(drop = True)\n",
    "\n",
    "        Route_parts['gridpark_no'] = Route_parts['Grid_No'].astype(str) +'-'+ Route_parts['Park_No'].astype(str)\n",
    "\n",
    "        # All fill value 0 because no routes are calculated for grid centroids in UGSs\n",
    "        to_fill = ['way-id','route_cost','steps','real_G-entry','Tcost']                                   \n",
    "        Route_parts[to_fill] = Route_parts[to_fill].fillna(0)  \n",
    "            \n",
    "        Routes.append(Route_parts)\n",
    "        Routes_detail.append(Route_dparts)\n",
    "    return(Routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "854c5949",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_route_finding(graph, suitible_chunk, nodes, mat_from, mat_to, mat_route, mat_step,\n",
    "                                             mat_way, mat_wbin, mat_nn, i, block, k, nn_iter):\n",
    "                        \n",
    "    # Order in route for nearest node:\n",
    "    # 1. gridnode to nearest to the original failed parknode\n",
    "    # 2. The reverse of 1.\n",
    "    # 3. nearest gridnode to the failed one and route to park\n",
    "    # 4. The reverse of 3.\n",
    "                        \n",
    "    gridosm = suitible_chunk['grid_osm'] # grid osmid\n",
    "    UGSosm = suitible_chunk['Parkroad_osmid'] # UGS osmid\n",
    "    nodeosm = nodes['osmid_var'] # road node osmid\n",
    "    nodegeom = nodes['geometry'] # road node geometry\n",
    "                        \n",
    "    len3 = 0\n",
    "    alt_route = list([])\n",
    "    while len3 < nn_iter and len(alt_route) < 1: # If a route is found (alt_route == 1) or until max iterations\n",
    "\n",
    "        len3 = len3 +1\n",
    "                            \n",
    "        nn = nn_finding(gridosm, UGSosm, nodeosm, nodegeom, nodes, i, len3) # finds nearest node.\n",
    "\n",
    "        nn_routing (graph, nn['currUGS'], nn['nearUGS'], nn['currgrid'], nn['neargrid'], \n",
    "                                        mat_way, mat_wbin, len3, alt_route) # executes route finding in try order.\n",
    "    if len(alt_route) == 0: \n",
    "        alt = alt_route \n",
    "    else: \n",
    "        alt = alt_route[0]\n",
    "    len4 = len(alt)\n",
    "    if len4 > 0: # If a route is found\n",
    "        mat_nn.append(i+block*k)\n",
    "        mat_from.append(alt)\n",
    "        shortest_to = list(alt[1:len(alt)])\n",
    "        shortest_to.append(-1)\n",
    "        mat_to.append(shortest_to)\n",
    "        mat_route.append(list(np.repeat(i+block*k,len4)))\n",
    "        mat_step.append(list(np.arange(0, len4)))\n",
    "    else: # If a route is not found\n",
    "        mat_from.append(-1)\n",
    "        mat_to.append(-1)\n",
    "        mat_route.append(i+block*k)\n",
    "        mat_step.append(-1)\n",
    "        mat_way.append('no way')\n",
    "        mat_wbin.append(2)\n",
    "        print('index',i+block*k,'No route')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "3ecd5703",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_finding (gridosm, UGSosm, nodeosm, nodegeom, nodes, i, nn_i): \n",
    "    # Grid nearest\n",
    "    g_geom = nodegeom[nodeosm == int(gridosm[i:i+1])] # Get geom of current node UGS\n",
    "    g_nearest = pd.DataFrame((abs(float(g_geom.x) - nodegeom.x)**2 # Check distance UGS\n",
    "    +abs(float(g_geom.y) - nodegeom.y)**2)**(1/2)\n",
    "                            ).join(nodeosm).sort_values(0) # sort by distance ascending UGS\n",
    "\n",
    "    g_grid = g_nearest.iloc[nn_i,1] # get the nearest node according to the nn_iter UGS entry\n",
    "    g_park = list(UGSosm)[i] # current node\n",
    "        \n",
    "    p_geom = nodegeom[nodeosm == int(UGSosm[i:i+1])] # get the geom of the current node grid\n",
    "    p_nearest = pd.DataFrame((abs(float(p_geom.x) - nodegeom.x)**2 # Check distance grid\n",
    "    +abs(float(p_geom.y) - nodegeom.y)**2)**(1/2)\n",
    "                            ).join(nodeosm).sort_values(0) # sort by distance ascending grid\n",
    "\n",
    "    p_grid = list(gridosm)[i] # current node\n",
    "    p_park = p_nearest.iloc[nn_i,1] # get the nearest node to the nn_iter grid\n",
    "    return({'currUGS':p_grid, 'nearUGS':p_park,'currgrid':g_park, 'neargrid':g_grid})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "50b89e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Improve: 2-to-2 instead of 1-to-all.\n",
    "\n",
    "def nn_routing (graph, curr_UGS, near_UGS, curr_grid, near_grid, mat_way, mat_wbin, nn_i, found_route):\n",
    "    try:\n",
    "        found_route.append(nx.shortest_path(graph, curr_UGS, near_UGS, \n",
    "                                          'travel_dist', method = 'dijkstra'))\n",
    "        mat_way.append(str(nn_i)+'grid > n-park') # grid to nearest unseen UGS node\n",
    "        mat_wbin.append(1)\n",
    "    except:\n",
    "        try:\n",
    "            found_route.append(nx.shortest_path(graph, near_UGS, curr_UGS, \n",
    "                                              'travel_dist', method = 'dijkstra'))\n",
    "            mat_way.append(str(nn_i)+'n-park > grid') # nearest unseen UGS node to grid\n",
    "            mat_wbin.append(0)\n",
    "        except:\n",
    "            try:\n",
    "                found_route.append(nx.shortest_path(graph, curr_grid, near_grid, \n",
    "                                                  'travel_dist', method = 'dijkstra'))\n",
    "                mat_way.append(str(nn_i)+'n-grid > park') # nearest grid node to UGS\n",
    "                mat_wbin.append(1)\n",
    "            except:\n",
    "                try:\n",
    "                    found_route.append(nx.shortest_path(graph, near_grid, curr_grid, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                    mat_way.append(str(nn_i)+'park > n-grid') # UGS to nearest grid node\n",
    "                    mat_wbin.append(0)\n",
    "                except:\n",
    "                    try:\n",
    "                        found_route.append(nx.shortest_path(graph, near_grid, near_UGS, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                        mat_way.append(str(nn_i)+'park > n-grid') # UGS to nearest grid node\n",
    "                        mat_wbin.append(0)\n",
    "                    except:\n",
    "                        try:\n",
    "                            found_route.append(nx.shortest_path(graph, near_UGS, near_grid, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                            mat_way.append(str(nn_i)+'park > n-grid') # UGS to nearest grid node\n",
    "                            mat_wbin.append(1)\n",
    "                        except:\n",
    "                            pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1cf68b13",
   "metadata": {},
   "outputs": [],
   "source": [
    "def route_formatting(mat_from, mat_to, mat_route, mat_step, road_edges):\n",
    "    # Unpack lists\n",
    "    s_mat_u = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_from) for i in b]\n",
    "    s_mat_u1 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_to) for i in b]\n",
    "    s_mat_u2 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_route) for i in b]\n",
    "    s_mat_u3 = [i for b in map(lambda x:[x] if not isinstance(x, list) else x, mat_step) for i in b]\n",
    "\n",
    "    # Format df\n",
    "    routes = pd.DataFrame([s_mat_u,s_mat_u1,s_mat_u2,s_mat_u3]).transpose()\n",
    "    routes.columns = ['from','to','route','step']\n",
    "    mat_key = list([])\n",
    "    for n in range(len(routes)): # get key of origin and destination\n",
    "        mat_key.append(str(int(s_mat_u[n])) + '-' + str(int(s_mat_u1[n])))\n",
    "    routes['key'] = mat_key\n",
    "    routes = routes.set_index('key')\n",
    "\n",
    "    # Add route information\n",
    "    routes = routes.join(road_edges, how = 'left') # to add road node information\n",
    "    routes = gpd.GeoDataFrame(routes, geometry = 'geometry', crs = 4326)\n",
    "    routes = routes.sort_values(by = ['route','step'])\n",
    "    return(routes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "4e347555",
   "metadata": {},
   "outputs": [],
   "source": [
    "def route_summarization(routes, suitible_comb, road_nodes, mat_way, mat_wbin):\n",
    "    # dissolve route\n",
    "    routes2 = routes[['route','geometry']].dissolve('route')\n",
    "\n",
    "    # get used grid- and parkosm. Differs at NN-route.\n",
    "    route_reset = routes.reset_index()\n",
    "    origin = route_reset['from'].iloc[list(route_reset.groupby('route')['step'].idxmin()),]\n",
    "    origin = origin.reset_index().iloc[:,-1]\n",
    "    dest = route_reset['from'].iloc[list(route_reset.groupby('route')['step'].idxmax()),]\n",
    "    dest = dest.reset_index().iloc[:,-1]\n",
    "\n",
    "    # grid > park = 1, park > grid = 0, no way = 2, detailed way in way_calc.\n",
    "    routes2['way-id'] = mat_wbin\n",
    "    routes2['realG_osmid'] = np.where(routes2['way-id'] == 1, origin, dest)\n",
    "    routes2['realP_osmid'] = np.where(routes2['way-id'] == 1, dest, origin)\n",
    "    routes2['way_calc'] = mat_way\n",
    "\n",
    "    # get route cost, steps, additional information.\n",
    "    routes2['route_cost'] = routes.groupby('route')['length'].sum()\n",
    "    routes2['steps'] = routes.groupby('route')['step'].max()\n",
    "    routes2['index'] = suitible_comb.index\n",
    "    routes2 = routes2.set_index(['index'])\n",
    "    routes2.index = routes2.index.astype(int)\n",
    "    routes2 = pd.merge(routes2, suitible_comb[['Grid_No','grid_osm','Park_No','Park_entry_No','Parkroad_osmid',\n",
    "                                          'Grid_m_centroid','walk_area_m2','Euclidean']],\n",
    "                                            left_index = True, right_index = True)\n",
    "    routes2 = pd.merge(routes2, road_nodes['geometry_m'], how = 'left', left_on = 'realG_osmid', right_index = True)\n",
    "    # calculate distance of used road-entry for grid-centroid.\n",
    "    routes2['real_G-entry'] = round(gpd.GeoSeries(routes2['Grid_m_centroid'], crs = 3043).distance(routes2['geometry_m']),3)\n",
    "                                    \n",
    "    # Calculcate total route cost for the four gravity variants\n",
    "    routes2['Tcost'] = routes2['route_cost'] + routes2['real_G-entry']\n",
    "    return(routes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "1ce9ee7b",
   "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": 20,
   "id": "3467e686",
   "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": "fc2146e7",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0710956a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c4704676",
   "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
}