Gravity WorldPoP-OSM Model 2.1-Func.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "675eb1c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import packages\n",
    "import geopandas as gpd\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import libpysal\n",
    "import networkx as nx\n",
    "import osmnx as ox\n",
    "import time\n",
    "import os\n",
    "from shapely import geometry\n",
    "from shapely.geometry import Point, MultiLineString, LineString, Polygon, MultiPolygon\n",
    "from shapely.ops import nearest_points, polygonize\n",
    "import shapely\n",
    "from itertools import product, combinations\n",
    "import math\n",
    "import warnings\n",
    "import socket\n",
    "from wpgpDownload.utils.dl import wpFtp\n",
    "from wpgpDownload.utils.isos import Countries\n",
    "from wpgpDownload.utils.convenience_functions import download_country_covariates as dl\n",
    "from wpgpDownload.utils.wpcsv import Product\n",
    "import georasters as gr\n",
    "from wpgpDownload.utils.convenience_functions import refresh_csv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "50e92f97",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Block 0 cities and assumptions\n",
    "\n",
    "start = time.time()\n",
    "\n",
    "cities = ['Dhaka Metropolitan']\n",
    "\n",
    "# idea to convert to dask-pandas and dask-geopandas\n",
    "# https://towardsdatascience.com/pandas-with-dask-for-an-ultra-fast-notebook-e2621c3769f\n",
    "# Or with Koalas (Spark-like pandas)\n",
    "\n",
    "# Assumptions\n",
    "thresholds = [300, 600, 1000] # route threshold in metres. WHO guideline speaks of access within 300m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d841a73e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "if prefer dwnl from terminal: \n",
      "wpgpDownload download -i BGD --id 5004\n",
      " \n",
      "downloaded:\n",
      "BGD downloaded 0.01 mns\n"
     ]
    }
   ],
   "source": [
    "# 1. Required preprocess for information extraction\n",
    "\n",
    "# Let's ignore depreciation warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "# Get the city boundaries\n",
    "bound_df = ox.geocoder.geocode_to_gdf(cities) # gets city boundaries from OSM\n",
    "\n",
    "# Get unique iso-codes of selected cities (only load country raster once)\n",
    "unique_iso = iso_countries(bound_df, # Finding the country of the bounded city\n",
    "                           cities)\n",
    "print(' ')\n",
    "\n",
    "print('downloaded:')\n",
    "# Get raster of countries (if automatic download is preferred (standard))\n",
    "raster = countries_grids(unique_iso,\n",
    "                         r'D:\\Dumps\\WorldPoP_Grids') # custom path, where grid files can be stored without downloading"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "9e396ef0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100m resolution grids extraction\n",
      "Dhaka Metropolitan 0.13 mns\n",
      " \n",
      "get road networks from OSM\n",
      "Dhaka Metropolitan done 0.86 mns\n",
      " \n",
      "get urban greenspaces from OSM\n",
      "Dhaka Metropolitan done\n"
     ]
    }
   ],
   "source": [
    "# 2. Information extraction\n",
    "\n",
    "# Clip cities from countries, format population grids\n",
    "population_grids = city_grids_format(bound_df, # city boundaries\n",
    "                                     unique_iso,\n",
    "                                     raster, # country raster\n",
    "                                     cities, \n",
    "                                     grid_size = 100)\n",
    "print(' ')\n",
    "\n",
    "# Get road networks\n",
    "road_networks = road_networks(cities, # Get 'all' (drive,walk,bike) network\n",
    "                                 thresholds,\n",
    "                                 undirected = True)\n",
    "\n",
    "# Road network returns a dict of keys consisting of graphs, nodes, edges and edges_full\n",
    "road_networks.keys()\n",
    "\n",
    "print(' ')\n",
    "# Extracting UGS\n",
    "UGS = urban_greenspace(cities, \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": 32,
   "id": "eb027bb5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "get fake UGS entry points\n",
      "Dhaka Metropolitan 0.0 % done 0.0  mns\n",
      "Dhaka Metropolitan 59.5 % done 0.21  mns\n",
      "Dhaka Metropolitan 100 % done 0.36  mns\n",
      " \n",
      "get potential (Euclidean) suitible combinations\n",
      "Dhaka Metropolitan\n",
      "chunk 1 / 3 21890 suitible comb.\n",
      "chunk 2 / 3 25070 suitible comb.\n",
      "chunk 3 / 3 5533 suitible comb.\n",
      "total combinations within distance 52493\n",
      "0.0 % gridentry done 0.0  mns\n",
      "100 % gridentry done 1.1  mns\n",
      " \n",
      "Check grids within UGS\n",
      "0 0.0  mns\n",
      "100 0.25  mns\n"
     ]
    }
   ],
   "source": [
    "# 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, \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,\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": 36,
   "id": "971ac326",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "77999.2338892479"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "UGS_entry[0]['park_size_walkable'].to_crs(3043).area.mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "40c30628",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "comb. by city\n",
      "Dhaka Metropolitan 52493\n",
      " \n",
      "Dhaka Metropolitan 1 / 1 range 0 - 50265\n",
      "0.0 % done 0.0 mns\n",
      "19.89 % done 0.06 mns\n",
      "39.79 % done 0.13 mns\n",
      "59.68 % done 0.21 mns\n",
      "79.58 % done 0.4 mns\n",
      "99.47 % done 0.47 mns\n",
      "0 nearest nodes found\n",
      "100.0 % pathfinding done 0.48 mns\n",
      "formatting done 0.71 mns\n",
      "dissolving done 0.98 mns\n"
     ]
    }
   ],
   "source": [
    "# 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, \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": 37,
   "id": "08decce9",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dhaka Metropolitan\n",
      "entrance 0.13 mns\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_20096\\1291535314.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m# 5. summarize scores\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m grid_scores = grid_score_summary (Routes['route summary'], # Shortest routes by the Dijkstra algorithm, with gravity variant distance adj.\n\u001b[0m\u001b[0;32m      3\u001b[0m                                   \u001b[0mcities\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m                                   \u001b[0mpopulation_grids\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m                                   \u001b[0mext\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m''\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;31m# At multiple runs, the extention prevents the summarized file to be overwritten.\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_20096\\1112879558.py\u001b[0m in \u001b[0;36mgrid_score_summary\u001b[1;34m(routes, cities, pop_grids, ext, grid_size)\u001b[0m\n\u001b[0;32m     29\u001b[0m                 \u001b[0mt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mthresholds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     30\u001b[0m                 \u001b[0mscore\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'tr_'\u001b[0m\u001b[1;33m+\u001b[0m \u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 31\u001b[1;33m                 \u001b[0mscores\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdetermine_scores\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvar_best_routes\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpop_grids\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mthresholds\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0ml1\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcities\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mn\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgrid_size\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m100\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     32\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     33\u001b[0m                 \u001b[0mgrdsc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpd\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconcat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mgrdsc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mscores\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'score_w_route'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_20096\\1729266355.py\u001b[0m in \u001b[0;36mdetermine_scores\u001b[1;34m(var_df, pop_grid, thresholds, var_abbr, city, grid_size)\u001b[0m\n\u001b[0;32m     44\u001b[0m         \u001b[0mos\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmakedirs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'D:Dumps/Scores output WP2-OSM/'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrid_size\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'm grids/Grid_lines/'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     45\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 46\u001b[1;33m     \u001b[0mgs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mto_file\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'D:Dumps/Scores output WP2-OSM/'\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mgrid_size\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m+\u001b[0m\u001b[1;34m'm grids/Grid_lines/gridscore_'\u001b[0m\u001b[1;33m+\u001b[0m \u001b[0mvar_abbr\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'_'\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mstr2\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'_'\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0mcity\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'.gpkg'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     47\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     48\u001b[0m     \u001b[0mgsc\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mloc\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m~\u001b[0m\u001b[0mgs\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcolumns\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misin\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mscore\u001b[0m \u001b[1;33m+\u001b[0m \u001b[1;34m'_routes'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\ssml\\lib\\site-packages\\geopandas\\geodataframe.py\u001b[0m in \u001b[0;36mto_file\u001b[1;34m(self, filename, driver, schema, index, **kwargs)\u001b[0m\n\u001b[0;32m   1112\u001b[0m         \u001b[1;32mfrom\u001b[0m \u001b[0mgeopandas\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mio\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mfile\u001b[0m \u001b[1;32mimport\u001b[0m \u001b[0m_to_file\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1113\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1114\u001b[1;33m         \u001b[0m_to_file\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfilename\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdriver\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mschema\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mindex\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m   1115\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m   1116\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mset_crs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcrs\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mepsg\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0minplace\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mallow_override\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\ssml\\lib\\site-packages\\geopandas\\io\\file.py\u001b[0m in \u001b[0;36m_to_file\u001b[1;34m(df, filename, driver, schema, index, mode, crs, **kwargs)\u001b[0m\n\u001b[0;32m    394\u001b[0m             \u001b[0mfilename\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdriver\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdriver\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcrs_wkt\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcrs_wkt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mschema\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mschema\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    395\u001b[0m         ) as colxn:\n\u001b[1;32m--> 396\u001b[1;33m             \u001b[0mcolxn\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwriterecords\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0miterfeatures\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    397\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    398\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m~\\miniconda3\\envs\\ssml\\lib\\site-packages\\fiona\\collection.py\u001b[0m in \u001b[0;36mwriterecords\u001b[1;34m(self, records)\u001b[0m\n\u001b[0;32m    359\u001b[0m         \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmode\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32min\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;34m'a'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m'w'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    360\u001b[0m             \u001b[1;32mraise\u001b[0m \u001b[0mIOError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"collection not open for writing\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 361\u001b[1;33m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msession\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwriterecs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mrecords\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m    362\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_len\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msession\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mget_length\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m    363\u001b[0m         \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_bounds\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32mfiona\\ogrext.pyx\u001b[0m in \u001b[0;36mfiona.ogrext.WritingSession.writerecs\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;32mfiona\\ogrext.pyx\u001b[0m in \u001b[0;36mfiona.ogrext.OGRFeatureBuilder.build\u001b[1;34m()\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "# 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, \n",
    "                                  population_grids, \n",
    "                                  ext = '', # At multiple runs, the extention prevents the summarized file to be overwritten.\n",
    "                                  grid_size = 100) # Size of the grid in meters\n",
    "grid_scores"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3e7c8c5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def iso_countries(bounds, cities):\n",
    "    # bound_df = ox.geocoder.geocode_to_gdf(cities)\n",
    "    # The 'Countries' is a list of iso-countries and descriptions from the package wpgpDownload.utils.isos\n",
    "    C = pd.DataFrame(Countries)\n",
    "    start_time = time.time()\n",
    "    iso_countries = []\n",
    "    print('if prefer dwnl from terminal: ')\n",
    "    \n",
    "    # Check the display name in the city boundaries to get the country name (enabling only specifying city in front)\n",
    "    for i in bounds['display_name']:\n",
    "        country = i.rsplit(',')[-1][1:]\n",
    "        iso = C[C['name'] == country].iloc[0,1]\n",
    "        # Get unique ISO countries, so all country-grids are only loaded once\n",
    "        if iso not in iso_countries:\n",
    "            iso_countries.append(iso)\n",
    "            \n",
    "            # List data and extract raster file download string with 2020 population (if download manually is preferred)\n",
    "            products = Product(iso)\n",
    "            Results = products.description_contains('people per grid-cell 2020')\n",
    "            list1 = []\n",
    "            for p in Results:\n",
    "                prints = '%s/%s\\t%s\\t%s' % (p.idx, p.country_name,p.dataset_name,p.path)\n",
    "                list1.append(prints)\n",
    "            print('wpgpDownload download -i',iso,'--id',list1[0].split(\"\\t\")[0].split('/')[0])\n",
    "    \n",
    "    return(iso_countries)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c68fbf7c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def countries_grids(iso_countries, download_dir = ' '):\n",
    "    start_time = time.time()\n",
    "    blocks = []\n",
    "    for iso in iso_countries:\n",
    "        # Check if raster files already exist on the system path or a manually specified path\n",
    "        path1 = os.getcwd() +'\\\\'+ iso.lower() + '_ppp_2020.tif'\n",
    "        path2 = download_dir +'\\\\'+ iso.lower() + '_ppp_2020.tif'\n",
    "        # First check the manual path\n",
    "        if os.path.exists(path2): \n",
    "            block = gr.from_file(path2)\n",
    "            blocks.append(block)\n",
    "        else:\n",
    "            # Then the system path\n",
    "            if os.path.exists(path1): \n",
    "                block = gr.from_file(path1)\n",
    "                blocks.append(block)\n",
    "            else:\n",
    "                # Otherwise run a suprocess (spr.run) command to download via the terminal in notebook.\n",
    "                runstr = 'wpgpDownload download -i '+ iso+ ' -f people --datasets'\n",
    "                p1 = spr.run('wpgpDownload download -i '+ iso+ ' -f people --datasets', \n",
    "                                    shell = True, \n",
    "                                    capture_output = True)\n",
    "                # decode the output to a list of available datasets from WorldPoP\n",
    "                datasets = p1.stdout.decode().rsplit('\\n')\n",
    "\n",
    "                # The first population raster grid (id-sorted) is the general one, without specifying to demographic groups\n",
    "                for i in enumerate(datasets):\n",
    "                    if '2020' in i[1]:\n",
    "                        ds = datasets[i[0]].rsplit('\\t')[0]\n",
    "                        print(ds)\n",
    "                        # if we found the file, we can stop the loop (we don't need the demograhically specified files)\n",
    "                        break\n",
    "                # Construct the download string\n",
    "                dwnl = 'wpgpDownload download -i '+iso+' --id '+str(ds)\n",
    "                # Get the specified file (terminal)\n",
    "                spr.run(dwnl, shell = True)\n",
    "                # Extract the file\n",
    "                block = gr.from_file(path1)\n",
    "                blocks.append(block)\n",
    "        print(iso,'downloaded', round((time.time() - start_time)/60,2),'mns')\n",
    "    return(blocks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "a44b205a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Block 2 population grids extraction\n",
    "def city_grids_format(bounds, iso_countries, country_grids, cities, 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(cities)):\n",
    "        C = pd.DataFrame(Countries)\n",
    "        iso = C[bounds['display_name'][i].rsplit(',')[-1][1:] == C['name']].iloc[0,1]\n",
    "        contains = [j for j, x in enumerate(iso_countries) if x == iso][0]\n",
    "\n",
    "        # Clip the city from the country\n",
    "        clipped = country_grids[contains].clip(bounds['geometry'][i])\n",
    "        clipped = clipped[0].to_geopandas()\n",
    "\n",
    "        # Get dissolvement_key for dissolvement. \n",
    "        clipped['row3'] = np.floor(clipped['row']/(grid_size/100)).astype(int)\n",
    "        clipped['col3'] = np.floor(clipped['col']/(grid_size/100)).astype(int)\n",
    "        clipped['dissolve_key'] = clipped['row3'].astype(str) +'-'+ clipped['col3'].astype(str)\n",
    "\n",
    "        # Dissolve into block by block grids\n",
    "        popgrid = clipped[['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(clipped.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) which is in my mind an coords error\n",
    "        # I therefore create one 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(cities[i].rsplit(',')[0], round((time.time() - start_time)/60,2),'mns')\n",
    "    return(grids)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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 cities:\n",
    "        # Get graph, road nodes and edges\n",
    "        graph = ox.graph_from_place(i, network_type = \"all\", buffer_dist = (np.max(thresholds)+1000))\n",
    "        #graphs.append(graph)\n",
    "\n",
    "        road_node, road_edge = ox.graph_to_gdfs(graph)\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 = road_edge.to_crs(4326)\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(i.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 cities:\n",
    "        gdf = ox.geometries_from_place(i, tags={'leisure':'park'}, 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(i, 'done')\n",
    "        parks = parks.reset_index()\n",
    "        parks['geometry_m'] = parks['geometry'].to_crs(3043)\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(it+block_size*block,'No route between grid and park-entry and their both',nn_iter,'alternatives')"
   ]
  },
  {
   "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):\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], 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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_geoms/'):\n",
    "                os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_geoms/')\n",
    "\n",
    "            gridsc.to_file('D:Dumps/Scores output WP2-OSM/'+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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_csv/'):\n",
    "                os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_csv/')\n",
    "\n",
    "            gridsc.to_csv('D:/Dumps/Scores output WP2-OSM/'+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('D:/Dumps/Scores output WP2-OSM/'+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, 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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_lines/'):\n",
    "        os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_lines/')\n",
    "                \n",
    "    gs.to_file('D:Dumps/Scores output WP2-OSM/'+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": 28,
   "id": "7d91c262",
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPTIONAL\n",
    "# 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_ranks (routes, UGS_entry, cities, pop_grids, ext = '', grid_size = 100, ranks = 1):\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",
    "            best_gridpark_comb_ranks = best_gridpark_comb (routes[n], UGS_entry[n], l1[i], pop_grids[n], ranks = 1)\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], 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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_geoms/'):\n",
    "                os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_geoms/')\n",
    "\n",
    "            gridsc.to_file('D:Dumps/Scores output WP2-OSM/'+str(grid_block_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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_csv/'):\n",
    "                os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids/Grid_csv/')\n",
    "\n",
    "            gridsc.to_csv('D:/Dumps/Scores output WP2-OSM/'+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('D:/Dumps/Scores output WP2-OSM/'+str(grid_block_size)+'m grids/popgrid_access.csv')\n",
    "    return(popg_acc)    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "8c5881ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "# OPTIONAL\n",
    "def best_gridpark_comb_ranks (routes, UGS_entry, var_abbr, pop_grid, ranks = 1):\n",
    "    str1 = 'gridpark_' + var_abbr\n",
    "    locals()[str1] = ranked_route_scoring (routes, UGS_entry, cities, var_abbr, ranks) \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": null,
   "id": "a758cd87",
   "metadata": {},
   "outputs": [],
   "source": [
    "def determine_scores(var_df, pop_grid, thresholds, var_abbr, city, 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('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids local/Grid_lines/'):\n",
    "        os.makedirs('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids local/Grid_lines/')\n",
    "                \n",
    "    gs.to_file('D:Dumps/Scores output WP2-OSM/'+str(grid_size)+'m grids local/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": 32,
   "id": "d90c64ad",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "99.65 mns\n"
     ]
    }
   ],
   "source": [
    "print(round((time.time() - start) / 60,2),'mns')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "2ea78431",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "665a9824",
   "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
}