Gravity GEE-OSM Model-Func 2.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "675eb1c0",
   "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": 3,
   "id": "d8720c72",
   "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=6Nyq081WQ2fxQk4vf_G50hV-El-gqb9XQ7gCdDos9Zs&tc=zRf1u4bSDil4SW7BAcLzjPIQLU3eRyU-IKx1CYwyy9U&cc=w8N3s9CfLVmyb5vRLFu-42NMmY-mOXNhNyKWWW4c8P4>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=6Nyq081WQ2fxQk4vf_G50hV-El-gqb9XQ7gCdDos9Zs&tc=zRf1u4bSDil4SW7BAcLzjPIQLU3eRyU-IKx1CYwyy9U&cc=w8N3s9CfLVmyb5vRLFu-42NMmY-mOXNhNyKWWW4c8P4</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/1AVHEtk575uHD8-D1tg494NBnIX6eXe9QnW-xLoeWVMhBNniS2D26aLNfXSg\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": 4,
   "id": "674831cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Thresholds and cities\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(['Addis Ababa','Dhaka','Shijiazhuang','Santo Domingo'])]\n",
    "cities_adj = cities_adj.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "7de9fb99",
   "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/f0de31b02bb07953d66a210523e93d01-76a258af2bfa11770a4695c13a1f5b2c:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\ETH_Addis Ababa_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/e374c5e555e31457cffe53390382762f-45c7cc7aa94fed3f50d21c2f19e15a2e:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\BGD_Dhaka_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/aa1d5bb09020fa6804d00e3e57c83081-aba785176e22ca563efe5ea4e8d70269:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\DOM_Santo Domingo_2020.tif\n",
      "Generating URL ...\n",
      "Downloading data from https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/thumbnails/c00d2dec9576c7941a5faf8cf3e6eb04-879651c71c4de2a1d0bce11577d82895:getPixels\n",
      "Please wait ...\n",
      "Data downloaded to C:\\Dumps\\GEE_city_grids\\CHN_Shijiazhuang_2020.tif\n",
      "Wall time: 15.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\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,'C:/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": 35,
   "id": "9e396ef0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Addis Ababa 0.42 mns\n",
      "Dhaka 0.68 mns\n",
      "Santo Domingo 1.78 mns\n",
      "Shijiazhuang 2.14 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Addis Ababa done 1.66 mns\n",
      "Dhaka done 2.68 mns\n",
      "Santo Domingo done 4.66 mns\n",
      "Shijiazhuang done 5.04 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Addis Ababa done\n",
      "Dhaka done\n",
      "Santo Domingo done\n",
      "Shijiazhuang done\n",
      "Wall time: 7min 14s\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": "eb027bb5",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get fake UGS entry points\n",
      "Addis Ababa 0.0 % done 0.01  mns\n",
      "Addis Ababa 73.5 % done 0.42  mns\n",
      "Addis Ababa 100 % done 0.56  mns\n",
      "Dhaka 0.0 % done 0.57  mns\n",
      "Dhaka 19.4 % done 0.81  mns\n",
      "Dhaka 38.8 % done 1.04  mns\n",
      "Dhaka 58.1 % done 1.3  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(' ')\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'], \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(' ')\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": "40c30628",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "# 4. Finding shortest routes.\n",
    "\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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08decce9",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "# 5. summarize scores\n",
    "grid_scores = grid_score_summary (Routes['route summary'], # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n",
    "                                  cities_adj['City'], \n",
    "                                  population_grids, \n",
    "                                  ext = '', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  save_path = 'C:/Dumps/GEE-WP Scores/Gravity/\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "60b25137",
   "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": 11,
   "id": "a44b205a",
   "metadata": {
    "scrolled": false
   },
   "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": 12,
   "id": "5bc1aa68",
   "metadata": {
    "scrolled": false
   },
   "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": 13,
   "id": "3d3ceef5",
   "metadata": {
    "scrolled": false
   },
   "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": 14,
   "id": "cdc51e59",
   "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']], left_on = 'Park_No', right_index = True)\n",
    "\n",
    "        # Get the walkable park size\n",
    "        ParkRoad['park_size_walkable'] = ParkRoad['geometry_m'].buffer(walk_radius).to_crs(4326).intersection(ParkRoad['geometry_y'])\n",
    "        ParkRoad['walk_area'] = ParkRoad['park_size_walkable'].to_crs(3043).area\n",
    "        ParkRoad['park_area'] = ParkRoad['geometry_y'].to_crs(3043).area\n",
    "        ParkRoad['share_walked'] = ParkRoad['walk_area'] / ParkRoad['park_area']\n",
    "        \n",
    "        # Get size inflation factors for the gravity model\n",
    "        ParkRoad['size_infl_factor'] = ParkRoad['walk_area'] / ParkRoad['walk_area'].median()\n",
    "        ParkRoad['size_infl_sqr2'] = ParkRoad['size_infl_factor']**(1/2)\n",
    "        ParkRoad['size_infl_sqr3'] = ParkRoad['size_infl_factor']**(1/3)\n",
    "        ParkRoad['size_infl_sqr5'] = ParkRoad['size_infl_factor']**(1/5)\n",
    "                \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)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "af76feed",
   "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": 16,
   "id": "19711d72",
   "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",
    "        #blockA = block_combinations\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",
    "            #print(blockC, 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",
    "            # 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",
    "                                       'size_infl_sqr2','size_infl_sqr3','size_infl_sqr5','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",
    "            infl2 = list(node_merged['size_infl_sqr2'])\n",
    "            infl3 = list(node_merged['size_infl_sqr3'])\n",
    "            infl5 = list(node_merged['size_infl_sqr5'])\n",
    "\n",
    "            # Get the euclidean distances\n",
    "            mat = list()\n",
    "            mat2 = list()\n",
    "            mat3 = list()\n",
    "            mat4 = 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",
    "                mat2.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl2[j])\n",
    "                mat3.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl3[j])\n",
    "                mat4.append(math.sqrt(abs(plon[j] - glon[j])**2 + abs(plat[j] - glat[j])**2) / infl5[j])\n",
    "\n",
    "            # Check if distances are within 1000m and join remaining info and concat in master df per 1000.\n",
    "            mat_df = pd.DataFrame(mat3)[(np.array(mat) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat2) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat3) <= np.max(thresholds)) | \n",
    "                                        (np.array(mat4) <= np.max(thresholds))]\n",
    "\n",
    "            # join the other gravity euclidean scores and other information\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat), lsuffix='_infl', rsuffix='_entr', how = 'left')\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat2), lsuffix='_entry', rsuffix='_pwr', how = 'left')\n",
    "            mat_df = mat_df.join(pd.DataFrame(mat4), lsuffix='_pwr', rsuffix='_root', how = 'left')\n",
    "            mat_df.columns = ['size_infl_eucl2','raw euclidean','size_infl_eucl3','size_infl_eucl5']    \n",
    "            mat_df = mat_df.join(node_merged)\n",
    "\n",
    "            output = pd.concat([output, mat_df])\n",
    "\n",
    "            print('chunk',(i+1),'/',len4,len(mat_df),'suitible comb.')\n",
    "        # Renaming columns\n",
    "        print('total combinations within distance',len(output))\n",
    "\n",
    "        output.columns = ['size_infl_eucl3','raw euclidean','size_infl_eucl2','size_infl_eucl5',\n",
    "                          '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',\n",
    "                          'park_lon','park_lat','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5',\n",
    "                          'parkshare_walked','park_area','walk_area_m2']\n",
    "\n",
    "        output = output[['raw euclidean','size_infl_eucl2','size_infl_eucl3','size_infl_eucl5',\n",
    "                         'Grid_No','Park_entry_No','Grid_coords_centroid','Grid_m_centroid',\n",
    "                          'Park_No','Parkroad_osmid','Park_geom','Parkroad_coords_centroid','Parkroad_m_centroid',\n",
    "                         'walk_area_m2','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5']]\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": 17,
   "id": "e9d2c955",
   "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": 18,
   "id": "a8f10468",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check grids in or out of UGS\n",
    "def grids_in_UGS (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']\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",
    "\n",
    "        dfGrUGS = pd.DataFrame(pd.DataFrame(np.array(lst)).unstack())\n",
    "        dfGrUGS.columns = ['in_out_UGS']\n",
    "        merged = pd.merge(suitible_comb[i], dfGrUGS, left_on = ['Grid_No','Park_No'], right_index = True, how = 'left')\n",
    "        RoadInOut.append(merged)\n",
    "    return(RoadInOut)    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "07b96b77",
   "metadata": {
    "scrolled": false
   },
   "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",
    "    warnings.filterwarnings(\"ignore\")\n",
    "    start_time = time.time()\n",
    "    \n",
    "    print('comb. by city')\n",
    "    for n in enumerate(cities): # Know how much comb. need to be calculcated.\n",
    "        print(n[1], len(combinations[n[0]]))\n",
    "    print(' ')\n",
    "    \n",
    "    Routes = list()\n",
    "    Routes_detail = list()\n",
    "    for j in range(len(cities)):\n",
    "        suit_raw = combinations[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",
    "        \n",
    "        len2 = int(np.ceil(len(suitible)/block_size)) # get number of blocks (chunks)\n",
    "        Route_parts = pd.DataFrame()\n",
    "        Route_dparts = pd.DataFrame()\n",
    "\n",
    "        # Divide in chunks of block for computational load\n",
    "        for k in range(len2):    \n",
    "            suitible_chunk = suitible.iloc[k*block_size:k*block_size+block_size] # Get block ids\n",
    "\n",
    "            parknode = list(suitible_chunk['Parkroad_osmid']) # UGS road entry ids\n",
    "            gridnode = list(suitible_chunk['grid_osm']) # grid centroid road entry ids\n",
    "\n",
    "            s_mat = list([]) # osmid from\n",
    "            s_mat1 = list([]) # osmid to\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 = [] # sums number of routes containing nearest nodes.\n",
    "            len1 = len(suitible_chunk)\n",
    "\n",
    "            print(cities[j].rsplit(',')[0], k+1,'/',len2, \n",
    "                  'range',k*block_size,'-',k*block_size+np.where(k*block_size+block_size >= len1,len1,block_size))\n",
    "            \n",
    "            for i in range(len(suitible_chunk)):\n",
    "                try:\n",
    "                    shortest = nx.shortest_path(graphs[j], 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_size*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(graphs[j], 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_size*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 the nearest node is taken, which is iterated X times at max, check assumptions, block #0 \n",
    "                        nn_route_finding(graphs[j], suitible_chunk, road_nodes[j],\n",
    "                                   s_mat, s_mat1, s_mat2, s_mat3, s_mat4, s_mat5, mat_nn, # matrice info see above\n",
    "                                   it = i, block = k, block_size = block_size, \n",
    "                                         nn_iter = 10) # max nearest nodes to be found\n",
    "                        \n",
    "                if i % 10000 == 0: print(round((i+block_size*k)/len(suitible)*100,2),'% done',\n",
    "                                         round((time.time() - start_time) / 60,2),'mns')\n",
    "            print(len(mat_nn),'nearest nodes found')\n",
    "\n",
    "            print(round((i+block_size*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])\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)\n",
    "            print('dissolving done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "            \n",
    "            # Concats chunk with others already calculated\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','size_infl_sqr2','size_infl_sqr3','size_infl_sqr5',\n",
    "                                   'raw 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','raw_Tcost','grav2_Tcost','grav3_Tcost','grav5_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({'route summary':Routes,'route detail':Routes_detail})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "460e1c8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_route_finding (Graph, comb, nodes, \n",
    "                mat_from, mat_to, mat_route, mat_step, mat_way, mat_wbin, mat_nn, \n",
    "                it, block, block_size = 250000, nn_iter = 10):\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",
    "    len3 = 0\n",
    "    alt_route = list([])\n",
    "    \n",
    "    gosm = comb['grid_osm'] # grid osmids (origin)\n",
    "    posm = comb['Parkroad_osmid'] # UGS osmids (destination)\n",
    "    node = nodes['geometry'] # road node geoms\n",
    "    node_osm = nodes['osmid_var'] # road node osmids\n",
    "    \n",
    "    while len3 < nn_iter and len(alt_route) < 1: # continue if no more than 10 nearest nodes or if a route is found\n",
    "        \n",
    "        len3 = len3 +1\n",
    "        # Finds nearest node per iteration.\n",
    "        nn = nn_finding(gosm, posm, node, node_osm, it, len3)\n",
    "        \n",
    "         # routing within graph and current and found nearest nodes of grids and UGS\n",
    "        nn_routing(Graph, nn['curr_park'], nn['near_park'], nn['curr_grid'], nn['near_grid'],\n",
    "                        mat_way, mat_wbin, alt_route, len3)\n",
    "        \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 append\n",
    "        mat_nn.append(it+block_size*block)\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(it+block_size*block,len4)))\n",
    "        mat_step.append(list(np.arange(0, len4)))\n",
    "    else: # if no route is found fill values.\n",
    "        mat_from.append(-1)\n",
    "        mat_to.append(-1)\n",
    "        mat_route.append(it+block_size*block)\n",
    "        mat_step.append(-1)\n",
    "        mat_way.append('no way')\n",
    "        mat_wbin.append(2)\n",
    "        print('index',it+block_size*block,'No route')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "b7b34fa8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_finding (grid_osmid, UGS_osmid, node_geom, node_osmid, it, nn_i):\n",
    "    # Grid nearest\n",
    "    g_geom = node_geom[node_osmid == int(grid_osmid[it:it+1])] # Get current grid road entry geometry\n",
    "    g_nearest = pd.DataFrame((abs(float(g_geom.x) - node_geom.x)**2 # Find nearest.\n",
    "                              +abs(float(g_geom.y) - node_geom.y)**2)**(1/2)\n",
    "                            ).join(node_osmid).sort_values(0)\n",
    "\n",
    "    g_grid = g_nearest.iloc[nn_i,1] # Take '1' because 0 will get the current node with distance 0.\n",
    "    g_park = list(UGS_osmid)[it]\n",
    "\n",
    "    p_geom = node_geom[node_osmid == int(UGS_osmid[it:it+1])] # Get current UGS raod entry geometry\n",
    "    p_nearest = pd.DataFrame((abs(float(p_geom.x) - node_geom.x)**2 # Find nearest\n",
    "                              +abs(float(p_geom.y) - node_geom.y)**2)**(1/2)\n",
    "                            ).join(node_osmid).sort_values(0)\n",
    "\n",
    "    p_grid = list(grid_osmid)[it]\n",
    "    p_park = p_nearest.iloc[nn_i,1] # Take '1' because 0 will get the current node with distance 0.\n",
    "    \n",
    "    return({'curr_park':p_grid, 'near_park':p_park, 'curr_grid':g_park, 'near_grid':g_grid}) # return as dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "79bfe490",
   "metadata": {},
   "outputs": [],
   "source": [
    "def nn_routing (Graph, curr_park, near_park, curr_grid, near_grid, mat_way, mat_wbin, found_route, nn_i):\n",
    "    try: # First try from current grid to nearest UGS id.\n",
    "        found_route.append(nx.shortest_path(Graph, curr_park, near_park, \n",
    "                                          'travel_dist', method = 'dijkstra'))\n",
    "        mat_way.append(str(nn_i)+'grid > n-park')\n",
    "        mat_wbin.append(1)\n",
    "    except:\n",
    "        try: # Else try the reverse.\n",
    "            found_route.append(nx.shortest_path(Graph, near_park, curr_park, \n",
    "                                              'travel_dist', method = 'dijkstra'))\n",
    "            mat_way.append(str(nn_i)+'n-park > grid')\n",
    "            mat_wbin.append(0)\n",
    "        except:\n",
    "            try: # If no success try from current UGS id to nearest grid id\n",
    "                found_route.append(nx.shortest_path(Graph, near_grid, curr_grid, \n",
    "                                                  'travel_dist', method = 'dijkstra'))\n",
    "                mat_way.append(str(nn_i)+'n-grid > park')\n",
    "                mat_wbin.append(1)\n",
    "            except:\n",
    "                try: # Else try the reverse\n",
    "                    found_route.append(nx.shortest_path(Graph, curr_grid, near_grid, \n",
    "                                                      'travel_dist', method = 'dijkstra'))\n",
    "                    mat_way.append(str(nn_i)+'park > n-grid')\n",
    "                    mat_wbin.append(0)\n",
    "                except: # if no routes are found pass.\n",
    "                    pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f4984b15",
   "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)):\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')\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": 24,
   "id": "90524c2b",
   "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','size_infl_sqr2','size_infl_sqr3',\n",
    "                                          'size_infl_sqr5','raw 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['raw_Tcost'] = routes2['route_cost'] + routes2['real_G-entry']\n",
    "    routes2['grav2_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr2']\n",
    "    routes2['grav3_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr3']\n",
    "    routes2['grav5_Tcost'] = (routes2['route_cost'] + routes2['real_G-entry']) / routes2['size_infl_sqr5']\n",
    "    return(routes2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "7ffd4567",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 8 determine best parkentry points from each grid, then calculate grid scores\n",
    "# and finally aggregate city access in categories (high, medium, low and no access)\n",
    "def grid_score_summary (routes, cities, pop_grids, ext = '', grid_size = 100, save_path = 'C:/Dumps/GEE-WP Scores/Gravity/'):\n",
    "    start_time = time.time()\n",
    "    popg_acc = pd.DataFrame()\n",
    "    grid_scores = list([])\n",
    "    gridpark = list([])\n",
    "    for n in range(len(cities)):    \n",
    "        print(cities[n])\n",
    "\n",
    "        # For the four distance decay variants regarding park size.\n",
    "        l1 = list(['raw','grav2','grav3','grav5'])\n",
    "        m1 = list(['entrance','gravity**(1/2)','gravity**(1/3)','gravity**(1/5)'])\n",
    "        grid_score = list([])\n",
    "        gridparks = list([])\n",
    "        gridpark.append(gridparks)\n",
    "        popgrid_access = pd.DataFrame()\n",
    "        for i in range(len(l1)):\n",
    "            # Get the lowest indices grouped by a key consisting of grid no and park no (best entry point from a grid to a park)\n",
    "            var_best_routes = best_gridpark_comb (routes[n], l1[i], pop_grids[n])\n",
    "\n",
    "            grdsc = pd.DataFrame()\n",
    "            gridsc = pd.DataFrame()\n",
    "            print(m1[i], round((time.time() - start_time) / 60,2), 'mns')\n",
    "\n",
    "            # For each threshold given, calculate a score\n",
    "            for k in range(len(thresholds)):\n",
    "                \n",
    "                t = thresholds[k]\n",
    "                score = 'tr_'+ str(t)\n",
    "                scores = determine_scores(var_best_routes, pop_grids[n], thresholds[k], l1[i], cities[n], \n",
    "                                          save_path, grid_size = 100)\n",
    "                \n",
    "                grdsc = pd.concat([grdsc, scores['score_w_route']], axis = 1)\n",
    "                gridsc = pd.concat([gridsc, scores['grid_score']])\n",
    "                                \n",
    "                # Group according to the categories just created and sum the populations living in those grids\n",
    "                popgacc = pd.DataFrame()\n",
    "                popgacc[m1[i]+'_'+str(t)] = scores['score_w_route'].groupby(score+'_access')['population'].sum()\n",
    "                popgrid_access = pd.concat([popgrid_access, popgacc],axis=1)   \n",
    "\n",
    "                print('grid ',t)\n",
    "\n",
    "            grid_score.append(grdsc)\n",
    "\n",
    "            gridsc = gridsc.join(pop_grids[n]['geometry'])\n",
    "            gridsc = gpd.GeoDataFrame(gridsc, geometry = 'geometry', crs = 4326)\n",
    "\n",
    "            if not os.path.exists(save_path+str(grid_size)+'m grids/Grid_geoms/'):\n",
    "                os.makedirs(save_path+str(grid_size)+'m grids/Grid_geoms/')\n",
    "\n",
    "            gridsc.to_file(save_path+str(grid_size)+'m grids/Grid_geoms/gridscore_'+ l1[i] + '_' + cities[n] + '.gpkg')\n",
    "\n",
    "            # Detailed scores to files number of cities * ways to measure = number of files.\n",
    "            # Different threshold-scores are in the same dataframe\n",
    "            gridsc = gridsc.loc[:, gridsc.columns!='geometry']\n",
    "\n",
    "            if not os.path.exists(save_path+str(grid_size)+'m grids/Grid_csv/'):\n",
    "                os.makedirs(save_path+str(grid_size)+'m grids/Grid_csv/')\n",
    "\n",
    "            gridsc.to_csv(save_path+str(grid_size)+'m grids/Grid_csv/gridscore_'+ l1[i] + '_' + cities[n] + '.csv')\n",
    "            gridparks.append(var_best_routes)\n",
    "\n",
    "        grid_scores.append(grid_score)\n",
    "\n",
    "        # For each city, divide the population access by group by the total to get its share.\n",
    "        popgrid_access = popgrid_access / popgrid_access.sum()\n",
    "        popgrid_access = pd.DataFrame(popgrid_access.unstack())\n",
    "        popg_acc = pd.concat([popg_acc, popgrid_access], axis = 1)\n",
    "\n",
    "        print(cities[n],'done', round((time.time() - start_time) / 60,2), 'mns')\n",
    "    popg_acc.columns = cities\n",
    "    popg_acc.to_csv(save_path+str(grid_size)+'m grids/popgrid_access.csv')\n",
    "    return(popg_acc)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "8c00432c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def best_gridpark_comb (routes, var_abbr, pop_grid):\n",
    "    Rclean = routes[routes['way_calc'] != 'no way'].reset_index()\n",
    "    str1 = 'gridpark_' + var_abbr\n",
    "    locals()[str1] = Rclean.iloc[Rclean.groupby('gridpark_no')[(str(var_abbr) +'_Tcost')].idxmin()]  \n",
    "\n",
    "    # Get grid information\n",
    "    locals()[str1] = pd.merge(locals()[str1], pop_grid[['population','geometry']],\n",
    "                            left_on = 'Grid_No', right_index = True, how = 'outer')\n",
    "    locals()[str1] = locals()[str1].reset_index()\n",
    "\n",
    "    # formatting\n",
    "    locals()[str1]['Park_No'] = locals()[str1]['Park_No'].fillna(-1)\n",
    "    locals()[str1]['Park_No'] = locals()[str1]['Park_No'].astype(int)\n",
    "    locals()[str1]['Park_entry_No'] = locals()[str1]['Park_entry_No'].fillna(-1)\n",
    "    locals()[str1]['Park_entry_No'] = locals()[str1]['Park_entry_No'].astype(int)\n",
    "    return(locals()[str1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "d3f1ff5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def determine_scores(var_df, pop_grid, thresholds, var_abbr, city, save_path, grid_size = 100):\n",
    "    t = thresholds\n",
    "    str2 = str(t)\n",
    "    score = 'tr_'+ str2\n",
    "\n",
    "    #Only get routes within the threshold given (it loops over every threshold) and calculate the scores\n",
    "    thold = var_df[var_df[var_abbr + '_Tcost'] <= t]\n",
    "    thold[score] = t - thold[var_abbr + '_Tcost']\n",
    "    thold['pop' + score] = thold[score] * thold['population']\n",
    "    thold['walk_area_ha' + str2] = var_df['walk_area_m2'] /10000\n",
    "    thold['walkha_person' + str2] = thold['population'] / thold['walk_area_ha' + str2]\n",
    "\n",
    "    # Join the gridpark information from before.\n",
    "    var_df = var_df.join(thold[[score,'pop' + score,'walk_area_ha' + str2, 'walkha_person' + str2]])\n",
    "    # get the grid_scores\n",
    "    gs = pd.DataFrame()\n",
    "    gs[[score,'pop_' + score,'walkha_' + str2]] = var_df.groupby(\n",
    "            'Grid_No')[score,'pop' + score, 'walk_area_ha' + str2].sum()\n",
    "\n",
    "    gs['walkha_person_' + score] = var_df.groupby('Grid_No')['walkha_person' + str2].mean()\n",
    "\n",
    "    trstr = var_df[var_df[score] > 0]\n",
    "    gs[score + '_parks'] = trstr.groupby('Grid_No')['gridpark_no'].count()\n",
    "\n",
    "    # Add the routes as a dissolved line_geom\n",
    "    gs[score + '_routes'] = gpd.GeoDataFrame(trstr[['Grid_No','geometry_x']],\n",
    "                                                  geometry = 'geometry_x', crs = 4326).dissolve('Grid_No')\n",
    "\n",
    "    # Add parks which grids have access to with its closest access point\n",
    "    gs[score+'Park:entry'] = trstr[trstr['Park_No'] >=0].groupby('Grid_No')['Park_No'].apply(list).astype(str\n",
    "    ) + ':' + trstr[trstr['Park_entry_No'] >=0].groupby('Grid_No')['Park_entry_No'].apply(list).astype(str)\n",
    "                \n",
    "    # determine the thresholds category-score. \n",
    "    # High >= threshold (perfect score to one park), medium is above half perfect, \n",
    "    # low is below this and no is no access to a park for a certain grid within the threshold given\n",
    "    gs[score+'_access'] = np.select([gs[score] >= t, (gs[score] < t) & (\n",
    "    gs[score]>= t/2), (gs[score] < t/2) & (gs[score]> 0), gs[score] <= 0],\n",
    "          ['1 high','2 medium','3 low','4 no'])\n",
    "    gs = gs.join(pop_grid['population'], how = 'outer')\n",
    "            \n",
    "    gs = gpd.GeoDataFrame(gs, geometry = score + '_routes', crs = 4326)\n",
    "            \n",
    "    if not os.path.exists(save_path+str(grid_size)+'m grids/Grid_lines/'):\n",
    "        os.makedirs(save_path+str(grid_size)+'m grids/Grid_lines/')\n",
    "                \n",
    "    gs.to_file(save_path+str(grid_size)+'m grids/Grid_lines/gridscore_'+ var_abbr + '_' + str2 + '_' + city + '.gpkg')\n",
    "            \n",
    "    gsc = gs.loc[:,~gs.columns.isin([score + '_routes'])]\n",
    "\n",
    "    return({'grid_score':gsc,'score_w_route':gs})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a28fb0a1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e0d37c8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2fd114c1",
   "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}