code_for_a_single_city.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 94,
   "id": "ee0484c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import geopandas as gpd\n",
    "import matplotlib.pyplot as plt\n",
    "import networkx as nx\n",
    "import numpy as np\n",
    "import osmnx as ox\n",
    "import pandas as pd\n",
    "from shapely import geometry\n",
    "import libpysal\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "id": "a1f88eee",
   "metadata": {},
   "outputs": [],
   "source": [
    "distance_threshold = [300,600,1000]\n",
    "distance_threshold_L = max(distance_threshold)\n",
    "crs = 3034\n",
    "\n",
    "# get the network of the city\n",
    "place_name = \"Ghent, Belgium\"\n",
    "G = ox.graph_from_place(place_name, network_type='all', buffer_dist = 2*distance_threshold_L)\n",
    "G_proj = ox.project_graph(G,to_crs='epsg:3043')\n",
    "G_nodes_proj, G_edges_proj = ox.graph_to_gdfs(G_proj, nodes=True, edges=True)\n",
    "edges = G_edges_proj.unary_union\n",
    "\n",
    "# read other data\n",
    "pop = gpd.read_file('file_path') # path for population grid cell\n",
    "green = gpd.read_file('file_path') # path for green space\n",
    "boundry= gpd.read_file('file_path') # path for boundry"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c09a927b",
   "metadata": {},
   "source": [
    "## Calculate the accessibility score using M2sfca"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8740959",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data preprocessing for green space; use the 'green_data_prepocessing' defined below\n",
    "green_access = green_data_prepocessing(park_poly, boundry, distance_threshold_L, destination = 'Entrance')\n",
    "\n",
    "# data preprocessing for population; use the 'population_data_prepocessing' defined below\n",
    "pop_remove_bound = population_data_prepocessing(population, boundry, distance_threshold_L)\n",
    " \n",
    "    \n",
    "# get pairs of pop grid cell and green space that euclidian distance between them is < 1000\n",
    "df_index = get_the_index_within_dist(pop_remove_bound,green_access )\n",
    "    \n",
    "# get the origin and target node ID from OSM using defined function below\n",
    "orig_id_df = origin_node(df_index,pop_remove_bound)\n",
    "target_id_df =  target_node(df_index, green_access)\n",
    "pop_merge = pd.merge(df_index, orig_id_df, how = 'left', on = 'pop_index' )\n",
    "merge = pd.merge(pop_merge, target_id_df, how = 'left', on = 'green_index' )\n",
    "    \n",
    "# get the shortest distance between the origin and target nodes using function definded below\n",
    "merge_dist = OD_Matrix(merge)\n",
    "  \n",
    "# select distance smaller than threshold\n",
    "OD_mile = merge_dist[merge_dist['distance'] < distance_threshold_L]\n",
    "    \n",
    "# if the green space could reach multiple entrance of the same green space\n",
    "# then use the defined fucntion to select the nearest entrance \n",
    "OD_nearest_E = neareat_entrance(OD_mile)\n",
    "   \n",
    "# apply the M2sfca model to get the score\n",
    "pop_score_df = M2sfca(OD_nearest_E)\n",
    "    \n",
    "# merge the result with the whole population grid cell \n",
    "pop_remove_bound['pop_index'] = list(range(len(pop_remove_bound)))\n",
    "final = pd.merge(pop_remove_bound, pop_score_df, how = 'left', on='pop_index')\n",
    "final = final[['PoP2015_Number','geometry','access_score','mean_dist','mean_area','mean_supply']]\n",
    "final['normalize'] = final['access_score']/final['PoP2015_Number']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f32fac65",
   "metadata": {},
   "source": [
    "## The following include the detailed defined function used above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "d36214f1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def green_data_prepocessing(green, boundry, buffer, destination = 'centroids'):\n",
    "    \"\"\"\n",
    "    green space data preprocessing:\n",
    "    green is the green space geodataframe;\n",
    "    boundry is the city's administrative boundry;\n",
    "    destination represents the fake gate of green space, 'centroids' means using centroids of \n",
    "    green space as access point. 'Entrance' means using the intersections points;\n",
    "    buffer is added around the green space to do the intersection with road network, usually 20m\n",
    "\n",
    "    \"\"\"\n",
    "     \n",
    "    # transfer the CRS\n",
    "    boundry = boundry.to_crs(epsg = 3043) \n",
    "    green =  green.to_crs(epsg = 3043)\n",
    "    \n",
    "    # clip green space within the buffered city boundry\n",
    "    # first to add a buffer around the boundry to ensure people \n",
    "    # live around the edge of the boundry could access the green space outside the city but still within distance threshold\n",
    "    boundry_buffer = boundry.buffer(buffer)\n",
    "    green_clip = green.clip(boundry_buffer,keep_geom_type=False)\n",
    "    \n",
    "    # merge the overlap green space into one polygon\n",
    "    s_ = gpd.GeoDataFrame(geometry=[green_clip.unary_union]).explode(\n",
    "        index_parts=False).reset_index(drop=True).set_crs(epsg = 3043)\n",
    "    s_ = gpd.sjoin(s_, green_clip, how='left')\n",
    "    green_dis = s_.dissolve(s_.index, aggfunc='first')\n",
    "    # reset the index \n",
    "    green_dis =green_dis.reset_index(drop=True)\n",
    "    #create the buffer of green space, ususally 20m, prepare for intersetion to get the accesses\n",
    "    green_dis['buffer_20m'] = green_dis['geometry'].buffer(20)\n",
    "    \n",
    "    # define the green areas that overlap with each other within 20m as the same component\n",
    "    W = libpysal.weights.fuzzy_contiguity(green_dis['buffer_20m'])\n",
    "    green_dis['green_components']= W.component_labels\n",
    "    green_space = green_dis.dissolve('green_components')\n",
    "\n",
    "    \n",
    "    #create the centroid of the green space\n",
    "    green_space['centroids'] = green_space['geometry'].centroid\n",
    "    #calculate the area of the green space\n",
    "    green_space['area'] = green_space['geometry'].area\n",
    "    # select area >400m2\n",
    "    green_space = green_space[green_space['area'] > 400]\n",
    "    green_space = green_space.reset_index(drop=True)\n",
    "    \n",
    "    entrance = []  \n",
    "    # choose different fake gate model\n",
    "    # use centroid of green space \n",
    "    if destination == 'centroids':\n",
    "        green_access = green_space[['centroids','area']]\n",
    "        green_access.rename(columns = {'centroids': 'geometry'}, inplace = True)\n",
    "        green_access['park_id'] = list(range(len(green_space)))\n",
    "        \n",
    "    \n",
    "    # use the intersection points: the road network and buffered green space polygon\n",
    "    elif destination == 'Entrance':\n",
    "        # add the buffer of 20m for each components\n",
    "        green_space['buffer'] = green_space['geometry'].buffer(20)\n",
    "        for i in range(len(green_space)):\n",
    "            green_area = green_space.loc[i,'geometry'].area\n",
    "            intersection = green_space['buffer'][i].boundary.intersection(edges)\n",
    "            try:\n",
    "                for point in intersection:\n",
    "                    dic = {'geometry': point, 'area': green_area, 'park_id': i}\n",
    "                    entrance.append(dic)\n",
    "            except: continue\n",
    "    green_access = gpd.GeoDataFrame(entrance)\n",
    "    # if the distance of two fake entrances is less than 50m, then keep the first entrance\n",
    "    green_access['buffer'] = green_access['geometry'].buffer(25)\n",
    "    w_e = libpysal.weights.fuzzy_contiguity(green_access['buffer'])\n",
    "    green_access['entrance_components']= w_e.component_labels\n",
    "    # Delete duplicate rows based on specific columns \n",
    "    green_access =  green_access.drop_duplicates(subset=[\"entrance_components\"], keep='first')\n",
    "    green_access = green_access.reset_index(drop=True)\n",
    "    \n",
    "    return green_access   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "ae9f2434",
   "metadata": {},
   "outputs": [],
   "source": [
    "# preprocessing population data\n",
    "def population_data_prepocessing(pop, boundry, distance_threshold):\n",
    "    \n",
    "    \"\"\"\n",
    "    population data prepocessing:\n",
    "    pop is the population geodataframe;\n",
    "    boundry is the city's administrative boundry; \n",
    "    distance_threshold is the walking or driving distance; \n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    # transfer the CRS\n",
    "    boundry = boundry.to_crs(epsg = 3043) \n",
    "    pop =  pop.to_crs(epsg = 3043)\n",
    "    \n",
    "    # clip population within the city boundry\n",
    "    boundry_2_buffer = boundry.buffer(2*distance_threshold)\n",
    "    pop_clip = pop.clip(boundry_2_buffer,keep_geom_type=False)\n",
    "\n",
    "    # remove grid cell if its population is less than 10 in case of extreme value\n",
    "    pop_remove = pop_clip[(pop_clip['PoP2015_Number'] < 10) == False]\n",
    "    pop_remove['area'] = pop_remove['geometry'].area\n",
    "    \n",
    "    # remove boundry population if the area of grid cell is less than half of the regular area\n",
    "    pop_remove_bound = pop_remove[pop_remove['area'] > pop_remove['area'].max()*0.99]\n",
    "    \n",
    "    #create the centroid of the population grid cell\n",
    "    pop_remove_bound['centroids'] = pop_remove_bound['geometry'].centroid\n",
    "    #create the buffer of each grid cell, which equals the distance threshold\n",
    "    pop_remove_bound['buffer'] = pop_remove_bound['geometry'].buffer(distance_threshold)\n",
    "    # reset the index \n",
    "    pop_remove_bound  = pop_remove_bound.reset_index(drop=True)\n",
    "    \n",
    "    return pop_remove_bound\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "id": "d99a68b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_the_index_within_dist(pop_clip,green_access ):\n",
    "    \n",
    "    \"\"\"\n",
    "    this step is to create a dataframe with pairs of population grid cell and green space\n",
    "    whose euclidian distance is smaller than distance threshold\n",
    "    \n",
    "    \"\"\"\n",
    "  \n",
    "    k= []\n",
    "    for i in range(len(pop_clip)):\n",
    "        for j in range(len(green_access)):\n",
    "            if green_access['geometry'][j].within(pop_clip['buffer'][i]):\n",
    "            #distance. = pop_clip['centroids'][i].distance(access_point['geometry'][j])\n",
    "                index = {'pop_index': i, 'green_index': j,\n",
    "                         'green_area': green_access['area'][j],\n",
    "                         'park_id': green_access['park_id'][j],\n",
    "                         'pop_num': pop_clip.loc[i, 'PoP2015_Number']}\n",
    "                k.append(index)\n",
    "    df_index = pd.DataFrame(k)\n",
    "    return df_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "id": "1b814e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "def origin_node(df_index,pop_clip):\n",
    "    # get the unique value of pop index\n",
    "    pop_unique = df_index['pop_index'].unique()\n",
    "    \n",
    "    orig_id = []\n",
    "    # get the nearest node id of each centroid of population grid cell\n",
    "    for i in range(len(pop_unique)):\n",
    "        pop_index = pop_unique[i]\n",
    "        orig_node = ox.distance.nearest_nodes(G_proj,pop_clip.loc[pop_index,'centroids'].x,\n",
    "                                                pop_clip.loc[pop_index,'centroids'].y)\n",
    "        dic = {'pop_index':pop_unique[i],'orig_node':orig_node}\n",
    "        orig_id.append(dic)\n",
    "    orig_id_df = pd.DataFrame(orig_id)\n",
    "    return orig_id_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "id": "0ecc31d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def target_node(df_index, green_access):\n",
    "    # get the unique value of green index\n",
    "    green_unique = df_index['green_index'].unique()\n",
    "    target_id = []\n",
    "    # get the nearest node id of each green space\n",
    "    for i in range(len(green_unique)):\n",
    "        green_index = green_unique[i]\n",
    "        target_node = ox.distance.nearest_nodes(G_proj,green_access.loc[green_index,'geometry'].x,\n",
    "                                                green_access.loc[green_index,'geometry'].y)\n",
    "        dic = {'green_index':green_unique[i],'target_node': target_node}\n",
    "        target_id.append(dic)\n",
    "    target_id_df = pd.DataFrame(target_id)\n",
    "    return target_id_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "id": "228c1a75",
   "metadata": {},
   "outputs": [],
   "source": [
    "def OD_Matrix(merge):\n",
    "    \"\"\"\n",
    "    get the real distance between each pair of population grid cell and green space \n",
    "    \"\"\"\n",
    "    distance = []\n",
    "    for i in range(len(merge)):\n",
    "        try: \n",
    "            orig_node = merge.loc[i, 'orig_node']\n",
    "            target_node = merge.loc[i, 'target_node']\n",
    "            dist = nx.shortest_path_length(G_proj, source=orig_node, \n",
    "                                           target=target_node, weight='length')\n",
    "            distance.append(dist)\n",
    "        except:\n",
    "            distance.append(np.nan)\n",
    "    merge['distance'] = distance\n",
    "    return merge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "id": "6960af02",
   "metadata": {},
   "outputs": [],
   "source": [
    "def neareat_entrance(OD_mile): # OD_mile is the OD Matrix after selection, e.g real_distance < 1km\n",
    "    unique = OD_mile['pop_index'].unique()\n",
    "    OD_nearest_E = []\n",
    "    for i in unique:\n",
    "        df = OD_mile.loc[OD_mile['pop_index'] == i]\n",
    "        df1 = df.groupby('park_id').min('distance')\n",
    "        for j in df1.index:\n",
    "            dic = {'pop_index': i,'park_id': j,'distance':df1.loc[j,'distance'],\n",
    "              'green_area': df1.loc[j,'green_area'], 'pop_num':df1.loc[j,'pop_num']}\n",
    "            OD_nearest_E.append(dic)\n",
    "    OD_nearest_E = pd.DataFrame(OD_nearest_E)\n",
    "    return OD_nearest_E"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "7fbac8c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "def M2sfca (OD_nearest_E, distance_threshold):\n",
    "    v = -distance_threshold**2/np.log(0.5)\n",
    "    # use gussian distribution: let v= 923325, then the weight for 800m is 0.5\n",
    "    # add a column of weight: apply the decay function on distance\n",
    "    OD_nearest_E['weight'] = np.exp(-(OD_nearest_E['distance'])**2/v).astype(float) \n",
    "    OD_nearest_E['weighted_pop'] = OD_nearest_E['weight'] * OD_nearest_E['pop_num']\n",
    "\n",
    "    # get the sum of weighted population each green space has to serve.\n",
    "    s_w_p = pd.DataFrame(OD_nearest_E.groupby('park_id').sum('weighted_pop')['weighted_pop'])\n",
    "    # delete other columns, because they are useless after groupby\n",
    "    s_w_p = s_w_p.rename({'weighted_pop':'sum_weighted_pop'},axis = 1)\n",
    "    middle = pd.merge(OD_nearest_E,s_w_p, how = 'left', on = 'park_id' )\n",
    "    # calculate the supply-demand ratio for each green space\n",
    "    middle['green_supply'] = middle['green_area']/middle['sum_weighted_pop']\n",
    "    \n",
    "    # caculate the accessbility score for each green space that each population grid cell could reach\n",
    "    middle['access_score'] = middle['weight'] * middle['green_supply']\n",
    "    # add the scores for each population grid cell\n",
    "    pop_score_df = pd.DataFrame(middle.groupby('pop_index').sum('access_score')['access_score'])\n",
    "\n",
    "    # calculate the mean distance of all the green space each population grid cell could reach\n",
    "    mean_dist = middle.groupby('pop_index').mean('distance')['distance']\n",
    "    pop_score_df['mean_dist'] = mean_dist\n",
    "\n",
    "    # calculate the mean area of all the green space each population grid cell could reach\n",
    "    mean_area = middle.groupby('pop_index').mean('green_area')['green_area']\n",
    "    pop_score_df['mean_area'] = mean_area\n",
    "\n",
    "    # calculate the mean supply_demand ratio of all the green space each population grid cell could reach\n",
    "    mean_supply = middle.groupby('pop_index').mean('green_supply')['green_supply']\n",
    "    pop_score_df['mean_supply'] = mean_supply\n",
    "\n",
    "    return pop_score_df\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e030f07",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (ssml)",
   "language": "python",
   "name": "ssml"
  },
  "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
}