<jupyter_start><jupyter_text>Advanced Spatial Analysis: Spatial Accessibility 2 OverviewIn this lecture, we will advance the 2SFCA method, which was covered in the previous lecture, by integrating **travel time and distance decay functions**. Compared to the original 2SFCA method, this method is called **Enhanced 2SFCA (E2SFCA)** method. The method is defined as follows: First step:$$\huge R_j = \frac{S_j}{\sum_{k\in {\left\{\color{blue}{t_{kj}} \le \color{blue}{t_0} \right\}}}^{}{P_k}\color{blue}{W_k}}$$where$R_j$: the supply-to-demand ratio of location $j$. $S_j$: the degree of supply (e.g., number of doctors) at location $j$. $P_k$: the degree of demand (e.g., population) at location $k$. $\color{blue}{t_{kj}}$: the travel time between locations $k$ and $j$. $\color{blue}{t_0}$: the threshold travel time of the analysis. $\color{blue}{W_k}$: Weight based on a distance decay function Second step:$$\huge A_i = \sum_{j\in {\left\{\color{blue}{t_{ij}} \le \color{blue}{t_0} \right\}}} R_j\color{blue}{W_j}$$where$A_i$: the accessibility measures at location $i$. $R_j$: the supply-to-demand ratio of location $j$. $\color{blue}{W_j}$: Weight based on a distance decay function<jupyter_code>import geopandas as gpd
import pandas as pd
import osmnx as ox
import time
from tqdm import tqdm, trange
from shapely.geometry import Point, MultiPoint
import networkx as nx
import matplotlib.pyplot as plt
from scipy.stats import pearsonr
from shapely.ops import cascaded_union, unary_union
import utils  # Local file
import warnings
warnings.filterwarnings("ignore")<jupyter_output><empty_output><jupyter_text>Let's see the result firstIn the maps below, E2SFCA method shows the result measured with threshold travel time **15 minutes**. The original 2SFCA method shows the result measured with threshold travel **distance 20Km** (=50 mph (80 km/h) * 0.25 hr (15 minutes)).<jupyter_code>E_step2 = gpd.read_file('./data/result_E2SFCA.shp')
step2 = gpd.read_file('./data/result_2SFCA.shp')
# Plotting accessibility measurement result. 
fig, ax = plt.subplots(1, 2, figsize=(15, 10))

# Enhanced 2SFCA method
E_step2.plot('access', ax=ax[0], figsize=(10,10), legend=True, cmap='Blues', scheme='FisherJenks')
E_step2.loc[E_step2['access'] == 0].plot(ax=ax[0], color='grey', zorder=1)
E_step2.boundary.plot(ax=ax[0], linestyle='dotted', lw=0.5, color='black', zorder=1)

# Original 2SFCA method
step2.plot('access', ax=ax[1], figsize=(10,10), legend=True, cmap='Blues', scheme='FisherJenks')
step2.loc[step2['access'] == 0].plot(ax=ax[1], color='grey', zorder=1)
step2.boundary.plot(ax=ax[1], linestyle='dotted', lw=0.5, color='black', zorder=1)<jupyter_output><empty_output><jupyter_text>Import Data - same data as the previous lecture<jupyter_code># Supply: hospitals in the city of Chicago
hospitals = gpd.read_file('./data/Chicago_Hospital_Info.shp')
hospitals.head(1)
# Demand: population per census tract
tracts = gpd.read_file('./data/Chicago_Tract.shp')
tracts.head(1)
fig, ax = plt.subplots(figsize=(10, 10))

tracts.plot('TotalPop', ax=ax, scheme='FisherJenks', cmap='Blues')
hospitals.plot(markersize='Total_Bed', ax=ax, color='black')
# Mobility: Chicago Road Network
G = ox.io.load_graphml('./data/chicago_road.graphml')
ox.plot_graph(G)
# This function helps you to find the nearest OSM node from a given GeoDataFrame
# If geom type is point, it will take it without modification, but 
# IF geom type is polygon or multipolygon, it will take its centroid to calculate the nearest element. 

def find_nearest_osm(network, gdf):
    for idx, row in tqdm(gdf.iterrows(), total=gdf.shape[0]):
        if row.geometry.geom_type == 'Point':
            nearest_osm = ox.distance.nearest_nodes(network, 
                                                    X=row.geometry.x, 
                                                    Y=row.geometry.y
                                                   )
        elif row.geometry.geom_type == 'Polygon' or row.geometry.geom_type == 'MultiPolygon':
            nearest_osm = ox.distance.nearest_nodes(network, 
                                        X=row.geometry.centroid.x, 
                                        Y=row.geometry.centroid.y
                                       )
        else:
            print(row.geometry.geom_type)
            continue

        gdf.at[idx, 'nearest_osm'] = nearest_osm

    return gdf

supply = find_nearest_osm(G, hospitals)
demand = find_nearest_osm(G, tracts)<jupyter_output><empty_output><jupyter_text>Advancement 1: Caculate the estimated travel time for each edgeTo calculate the catchment area based on threshold travel time, we need to calculate how long it would take to travel each network edge. OSM network has two attributes that are helpful to calculate the estimated travel time for each edge: `length` and `maxspeed`. If we divide `length` into `maxspeed`, we will get the minimum travel time necessary to travel to the edge.<jupyter_code># Extract the nodes and edges of the network dataset for the future analysis. 
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True, node_geometry=True)
edges.head()<jupyter_output><empty_output><jupyter_text>You can iterate through the edges in the graph (`G`) with `G.edges()` method. This will return the entry of each edge as the form of dictionary.<jupyter_code>for u, v, data in G.edges(data=True):
    print(type(data), data.keys())<jupyter_output><empty_output><jupyter_text>Check the data in the `maxspeed` column. It has the forms of either **string or list**. But, we need to have it as a **numerical form** to do the calculation.<jupyter_code>for u, v, data in G.edges(data=True):
    if 'maxspeed' in data.keys():
        print(data['maxspeed'])
str_test = '55 mph'
str_test.split(' ')
list_test = ['35 mph', '30 mph']
list_test[0].split(' ')<jupyter_output><empty_output><jupyter_text>By splitting either list or string, we can obtain the numerical value of max speed, as shown below.<jupyter_code>for u, v, data in G.edges(data=True):
    if 'maxspeed' in data.keys():
        if type(data['maxspeed']) == list:
            temp_speed = data['maxspeed'][0]  # extract only the first entry if there are many 
        else: 
            temp_speed = data['maxspeed'] 

        temp_speed = temp_speed.split(' ')[0]  # Extract only the number
        
        data['maxspeed'] = temp_speed # Assign back to the original entry<jupyter_output><empty_output><jupyter_text>Examine the replaced values in `maxspeed` column. You will notice that some rows have `NaN` value.<jupyter_code>nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True, node_geometry=True)
edges.head()<jupyter_output><empty_output><jupyter_text>If the `maxspeed` column is empty, we can assign maximum travel speed based on their road type. The type is stored in `highway` column. This  website  shows the kinds of attributes and their meanings.<jupyter_code>def assign_max_speed_with_highway_type(row_):
    """
    Assign the maximum speed of an edge based on its attribute 'highway'
    # https://wiki.openstreetmap.org/wiki/Key:highway
    
    Args:
        row_: (dict) a row of OSMnx network data
    
    Returns:
        temp_speed_: (int) the maximum speed of an edge
    
    """
    
    max_speed_per_type = {'motorway': 50, 
              'motorway_link': 30, 
              'trunk': 50,
              'trunk_link': 30, 
              'primary': 40,
              'primary_link': 30, 
              'secondary': 40,
              'secondary_link': 30, 
              'tertiary': 40, 
              'tertiary_link': 20,
              'residential': 30,
              'living_street': 20,
              'unclassified': 20
             }
    
    # if the variable is a list, grab just the first one. 
    if type(row_['highway']) == list: 
        road_type = row_['highway'][0]
    else:
        road_type = row_['highway']
    
    # If the maximum speed of the road_type is predefined.
    if road_type in max_speed_per_type.keys():
        temp_speed_ = max_speed_per_type[road_type]
    else: # If not defined, just use 20 mph. 
        temp_speed_ = 20
    
    return temp_speed_
for u, v, data in G.edges(data=True):
    if 'maxspeed' in data.keys():
        if type(data['maxspeed']) == list:
            temp_speed = data['maxspeed'][0]  # extract only numbers         
        else: 
            temp_speed = data['maxspeed']

        temp_speed = temp_speed.split(' ')[0]
        
    else:
        temp_speed = assign_max_speed_with_highway_type(data)
        
    data['maxspeed'] = temp_speed<jupyter_output><empty_output><jupyter_text>Check the `maxspeed` column one more time. You will see all the rows are populated with numerical values.<jupyter_code># Extract the nodes and edges of the network dataset for the future analysis. 
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True, node_geometry=True)
edges.head()<jupyter_output><empty_output><jupyter_text>Since we have `maxspeed` and `length` data ready for every edge, now we can calculate the estimated travel time per edge. You can simply create/add column like the below.<jupyter_code>for u, v, data in G.edges(data=True):
    data['maxspeed_meters'] = int(data['maxspeed']) * 26.8223 # MPH * 1.6 * 1000 / 60; meter per minute
    data['time'] = float(data['length'] / data['maxspeed_meters'])  # Unit: minutes
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True, node_geometry=True)
edges.head()<jupyter_output><empty_output><jupyter_text>In summary, the following codes are what I have explained. ```pythondef network_settings(network):    for u, v, data in network.edges(data=True):        if 'maxspeed' in data.keys():            if type(data['maxspeed']) == list:                temp_speed = data['maxspeed'][0]   extract only numbers            else:                temp_speed = data['maxspeed']            temp_speed = temp_speed.split(' ')[0]        else:            temp_speed = assign_max_speed_with_highway_type(data)        data['maxspeed'] = temp_speed        data['maxspeed_meters'] = int(data['maxspeed']) * 26.8223   MPH * 1.6 * 1000 / 60; meter per minute        data['time'] = float(data['length'] / data['maxspeed_meters'])     create point geometries for the entire graph    for node, data in network.nodes(data=True):        data['geometry'] = Point(data['x'], data['y'])    return networkdef assign_max_speed_with_highway_type(row_):    max_speed_per_type = {'motorway': 50,                          'motorway_link': 30,                          'trunk': 50,                          'trunk_link': 30,                          'primary': 40,                          'primary_link': 30,                          'secondary': 40,                          'secondary_link': 30,                          'tertiary': 40,                          'tertiary_link': 20,                          'residential': 30,                          'living_street': 20,                          'unclassified': 20                          }     if the variable is a list, obtain just the first one.    if type(row_['highway']) == list:        road_type = row_['highway'][0]    else:        road_type = row_['highway']     If the maximum speed of the road_type is predefined.    if road_type in max_speed_per_type.keys():        temp_speed_ = max_speed_per_type[road_type]    else:   If not defined, just use 20 mph.        temp_speed_ = 20    return temp_speed_```<jupyter_code>def network_settings(network):
    for u, v, data in network.edges(data=True):
        if 'maxspeed' in data.keys():
            if type(data['maxspeed']) == list:
                temp_speed = data['maxspeed'][0]  # extract only numbers
            else:
                temp_speed = data['maxspeed']

            temp_speed = temp_speed.split(' ')[0]

        else:
            temp_speed = assign_max_speed_with_highway_type(data)

        data['maxspeed'] = temp_speed
        data['maxspeed_meters'] = int(data['maxspeed']) * 26.8223  # MPH * 1.6 * 1000 / 60; meter per minute
        data['time'] = float(data['length'] / data['maxspeed_meters'])

    # create point geometries for the entire graph
    for node, data in network.nodes(data=True):
        data['geometry'] = Point(data['x'], data['y'])

    return network


def assign_max_speed_with_highway_type(row_):
    max_speed_per_type = {'motorway': 50,
                          'motorway_link': 30,
                          'trunk': 50,
                          'trunk_link': 30,
                          'primary': 40,
                          'primary_link': 30,
                          'secondary': 40,
                          'secondary_link': 30,
                          'tertiary': 40,
                          'tertiary_link': 20,
                          'residential': 30,
                          'living_street': 20,
                          'unclassified': 20
                          }

    # if the variable is a list, obtain just the first one.
    if type(row_['highway']) == list:
        road_type = row_['highway'][0]
    else:
        road_type = row_['highway']

    # If the maximum speed of the road_type is predefined.
    if road_type in max_speed_per_type.keys():
        temp_speed_ = max_speed_per_type[road_type]
    else:  # If not defined, just use 20 mph.
        temp_speed_ = 20

    return temp_speed_
# Mobility: Chicago Road Network
G = ox.io.load_graphml('./data/chicago_road.graphml')

G = network_settings(G)

# Extract the nodes and edges of the network dataset for the future analysis. 
nodes, edges = ox.graph_to_gdfs(G, nodes=True, edges=True, node_geometry=True)
edges.head()<jupyter_output><empty_output><jupyter_text>--- *Exercise*Now we will investigate how the catchment area differs if we utilize threshold travel distance or travel time.  Assuming the overall travel speed in the study area is 50 mph, we will compare the catchment areas drawn by 20 km (15 minutes driving distance with 50 mph) and 15 minutes. Change the value of `supply_idx` from 0 to 33, and investigate how the catchment looks different at different supply locations.  ```python In summary, the following is the necessary code to create a catchment area from a given location. threshold_dist = 20000   15 minute driving distance with 50mph. threshold_time = 15supply_idx = 0  Range can be 0 - 33 1. Calculate accessible nodes in the network dataset from a given location temp_nodes_time = nx.single_source_dijkstra_path_length(G, supply.loc[supply_idx, 'nearest_osm'], threshold_time, weight='time')temp_nodes_dist = nx.single_source_dijkstra_path_length(G, supply.loc[supply_idx, 'nearest_osm'], threshold_dist, weight='length') 2. Extract the locations (or coordinates) of accessible nodes based on the OSMID.access_nodes_time = nodes.loc[nodes.index.isin(temp_nodes_time.keys()), 'geometry']access_nodes_dist = nodes.loc[nodes.index.isin(temp_nodes_dist.keys()), 'geometry'] 3. Create a convex hull with the locations of the nodes.access_nodes_time = gpd.GeoSeries(access_nodes_time.unary_union.convex_hull, crs="EPSG:4326")access_nodes_dist = gpd.GeoSeries(access_nodes_dist.unary_union.convex_hull, crs="EPSG:4326") Result. demand_time = demand.loc[demand['geometry'].centroid.within(access_nodes_time[0])]demand_dist = demand.loc[demand['geometry'].centroid.within(access_nodes_dist[0])]print(f"threshold by time: {demand_time.shape[0]}")print(f"threshold by distance: {demand_dist.shape[0]}") Plot graphsfig, ax = plt.subplots(figsize=(10, 10))tracts.plot('TotalPop', ax=ax, scheme='FisherJenks', cmap='Blues')access_nodes_time.boundary.plot(ax=ax, color='red', linewidth=4)access_nodes_dist.boundary.plot(ax=ax, color='blue', linewidth=4)edges.plot(ax=ax, color='black', lw=0.5)supply_loc = supply.loc[supply.index==supply_idx]supply_loc.plot(markersize='Total_Bed', ax=ax, color='black')```---<jupyter_code># In summary, the following is the necessary code to create a catchment area from a given location. 
threshold_dist = 20000  # 15 minute driving distance with 50mph. 
threshold_time = 15

supply_idx = 0 # Range can be 0 - 33

# 1. Calculate accessible nodes in the network dataset from a given location 
temp_nodes_time = nx.single_source_dijkstra_path_length(G, supply.loc[supply_idx, 'nearest_osm'], threshold_time, weight='time')
temp_nodes_dist = nx.single_source_dijkstra_path_length(G, supply.loc[supply_idx, 'nearest_osm'], threshold_dist, weight='length')

# 2. Extract the locations (or coordinates) of accessible nodes based on the OSMID.
access_nodes_time = nodes.loc[nodes.index.isin(temp_nodes_time.keys()), 'geometry']
access_nodes_dist = nodes.loc[nodes.index.isin(temp_nodes_dist.keys()), 'geometry']

# 3. Create a convex hull with the locations of the nodes.
access_nodes_time = gpd.GeoSeries(access_nodes_time.unary_union.convex_hull, crs="EPSG:4326")
access_nodes_dist = gpd.GeoSeries(access_nodes_dist.unary_union.convex_hull, crs="EPSG:4326")

# Result. 
demand_time = demand.loc[demand['geometry'].centroid.within(access_nodes_time[0])]
demand_dist = demand.loc[demand['geometry'].centroid.within(access_nodes_dist[0])]

print(f"threshold by time: {demand_time.shape[0]}")
print(f"threshold by distance: {demand_dist.shape[0]}")

# Plot graphs
fig, ax = plt.subplots(figsize=(10, 10))

tracts.plot('TotalPop', ax=ax, scheme='FisherJenks', cmap='Blues')

access_nodes_time.boundary.plot(ax=ax, color='red', linewidth=4)
access_nodes_dist.boundary.plot(ax=ax, color='blue', linewidth=4)

edges.plot(ax=ax, color='black', lw=0.5)

supply_loc = supply.loc[supply.index==supply_idx]
supply_loc.plot(markersize='Total_Bed', ax=ax, color='yellow', zorder=2)<jupyter_output><empty_output><jupyter_text>Here, the main difference is created from `nx.single_source_dijkstra_path_length()`, which is the function calculating the nodes that are accessible within a certain threshold.<jupyter_code>nx.single_source_dijkstra_path_length(G=G,
                                      source=supply.loc[0, 'nearest_osm'], 
                                      cutoff=20000,
                                      weight='length'
                                     )
nx.single_source_dijkstra_path_length(G=G,
                                      source=supply.loc[0, 'nearest_osm'], 
                                      cutoff=15,
                                      weight='time'
                                     )
nx.single_source_dijkstra_path_length(G=G,
                                      source=supply.loc[0, 'nearest_osm'], 
                                      cutoff=15,
#                                       weight='time'
                                     )<jupyter_output><empty_output><jupyter_text>Advancement 2: Apply distance decay functions for catchment areasHere, we will also start with the result first.<jupyter_code>def calculate_catchment_area(network, nearest_osm, minutes, distance_unit='time'):
    polygons = gpd.GeoDataFrame(crs="EPSG:4326")

    # Create convex hull for each travel time (minutes), respectively.
    for minute in minutes:
        access_nodes = nx.single_source_dijkstra_path_length(network, nearest_osm, minute, weight=distance_unit)
        convex_hull = gpd.GeoSeries(nx.get_node_attributes(network.subgraph(access_nodes), 'geometry')).unary_union.convex_hull
        polygon = gpd.GeoDataFrame({'minutes': [minute], 'geometry': [convex_hull]}, crs="EPSG:4326")
        polygon = polygon.set_index('minutes')
        polygons = polygons.append(polygon)

    # Calculate the differences between convex hulls which created in the previous section.
    polygons_ = polygons.copy(deep=True)
    for idx, minute in enumerate(minutes):
        if idx != 0:
            current_polygon = polygons.loc[[minute]]
            previous_polygons = polygons.loc[[minutes[idx-1]]]
            diff_polygon = gpd.overlay(current_polygon, previous_polygons, how="difference")
            if diff_polygon.shape[0] != 0:
                polygons_.at[minute, 'geometry'] = diff_polygon['geometry'].values[0]

    return polygons_.copy(deep=True)<jupyter_output><empty_output><jupyter_text>Demonstration of a catchment area drawn from the supply location 0 for 5, 10 15 minutes threshold travel time.<jupyter_code>supply_idx = 0 # Range can be 0 - 33

# Calculate catchment areas
areas = calculate_catchment_area(G, supply.loc[supply_idx, 'nearest_osm'], [5, 10,15], distance_unit='time')
areas['val'] = areas.index.astype(str)

# Plot graphs
fig, ax = plt.subplots(figsize=(10, 10))

tracts.plot('TotalPop', ax=ax, scheme='FisherJenks', cmap='Blues')
areas.plot('val', categorical=True, alpha=0.7, ax=ax)
areas.boundary.plot(ax=ax, color='black')
edges.plot(ax=ax, color='black', lw=0.5)

supply_loc = supply.loc[supply.index==supply_idx]
supply_loc.plot(markersize='Total_Bed', ax=ax, color='black')
areas = calculate_catchment_area(G, supply.loc[0, 'nearest_osm'], [5, 10,15], distance_unit='time')
areas
areas.loc[5, 'geometry']
areas.loc[10, 'geometry']
areas.loc[15, 'geometry']<jupyter_output><empty_output><jupyter_text>The steps are as follows to create multiple polygons with a hole inside. 1. Create polygons based on each step of threshold travel time (e.g., 5, 10, 15 minutes)2. Calculate the difference between a polygon with a bigger threshold travel time and the one with a smaller threshold travel time (e.g., 15 minute polygon - 10 minute polygon).<jupyter_code>minutes = [5, 10, 15]
polygons = gpd.GeoDataFrame(crs="EPSG:4326")

# Create convex hull for each travel time (minutes), respectively.
for minute in minutes:
    # Get the accessible nodes within a certain threshold travel time from the network
    access_nodes = nx.single_source_dijkstra_path_length(G, 
                                                         supply.loc[supply_idx, 'nearest_osm'],
                                                         minute, 
                                                         weight='time'
                                                        )
    # Create the convex hull of accessible nodes
    convex_hull = gpd.GeoSeries(nx.get_node_attributes(G.subgraph(access_nodes), 'geometry')).unary_union.convex_hull
    
    # `convex_hull` is a Shapely Polygon, so need to convert it to GeoDataFrame `polygon`
    polygon = gpd.GeoDataFrame({'minutes': [minute], 'geometry': [convex_hull]}, crs="EPSG:4326")
    
    # Append a GeoDataFrame to another GeoDataFrame
    polygons = polygons.append(polygon)
    
polygons = polygons.set_index('minutes')
polygons<jupyter_output><empty_output><jupyter_text>The result will have multiple polygons but they don't have a hole in it.<jupyter_code>polygons.loc[5, 'geometry']
polygons.loc[10, 'geometry']
polygons.loc[15, 'geometry']<jupyter_output><empty_output><jupyter_text>We can take advantage of `gpd.overlay()` to calculate the difference between two polygons.<jupyter_code>gpd.overlay(polygons.loc[[15]], polygons.loc[[10]], how="difference")
# gpd.overlay(polygons.loc[[15]], polygons.loc[[10]], how="difference").plot()
# Key here is that you are subtracting a polygon at the second location from the one at the first location
gpd.overlay(polygons.loc[[10]], polygons.loc[[15]], how="difference")<jupyter_output><empty_output><jupyter_text>If you make a loop to automate this process, it will be as below.<jupyter_code># Calculate the differences between convex hulls which created in the previous section.

minutes = [5, 10, 15]

polygons_ = polygons.copy(deep=True)
for idx, minute in enumerate(minutes):
    print(f'The index of {minute} is {idx}')
    
    current_idx = idx
    previous_idx = idx-1
    print(f'In the loop, the current index is {current_idx}, and previous index is {previous_idx}')
    
    current_threshold = minutes[current_idx]
    previous_threshold = minutes[previous_idx]
    print(f'In the loop, the current threshold time is {current_threshold}, and previous threshold time is {previous_threshold}')
    
    print('#-----------#')
#     if idx != 0:
#         current_polygon = polygons.loc[[minute]]
#         previous_polygons = polygons.loc[[minutes[idx-1]]]
#         diff_polygon = gpd.overlay(current_polygon, previous_polygons, how="difference")
#         if diff_polygon.shape[0] != 0:
#             polygons_.at[minute, 'geometry'] = diff_polygon['geometry'].values[0]<jupyter_output><empty_output><jupyter_text>Again, the following summarizes steps that create multiple catchment areas from a single origin.<jupyter_code>def calculate_catchment_area(network, nearest_osm, minutes, distance_unit='time'):
    polygons = gpd.GeoDataFrame(crs="EPSG:4326")

    # Create convex hull for each travel time (minutes), respectively.
    for minute in minutes:
        access_nodes = nx.single_source_dijkstra_path_length(network, nearest_osm, minute, weight=distance_unit)
        convex_hull = gpd.GeoSeries(nx.get_node_attributes(network.subgraph(access_nodes), 'geometry')).unary_union.convex_hull
        polygon = gpd.GeoDataFrame({'minutes': [minute], 'geometry': [convex_hull]}, crs="EPSG:4326")
        polygon = polygon.set_index('minutes')
        polygons = polygons.append(polygon)

    # Calculate the differences between convex hulls which created in the previous section.
    polygons_ = polygons.copy(deep=True)
    for idx, minute in enumerate(minutes):
        if idx != 0:
            current_polygon = polygons.loc[[minute]]
            previous_polygons = polygons.loc[[minutes[idx-1]]]
            diff_polygon = gpd.overlay(current_polygon, previous_polygons, how="difference")
            if diff_polygon.shape[0] != 0:
                polygons_.at[minute, 'geometry'] = diff_polygon['geometry'].values[0]

    return polygons_.copy(deep=True)<jupyter_output><empty_output><jupyter_text>Implementation of the advancements to the accessibility measurements The original 2SFCA method calculates supply-to-demand ratio (Step 1) as shown below.<jupyter_code># Calculate supply-to-demand ratio of supply location 0
i= 0
dist = 20000
supply_ = supply.copy(deep=True)
supply_['ratio'] = 0

# Create a catchment area from a given location
temp_nodes = nx.single_source_dijkstra_path_length(G, supply.loc[i, 'nearest_osm'], dist, weight='length')
access_nodes = nodes.loc[nodes.index.isin(temp_nodes.keys()), 'geometry']
access_nodes = gpd.GeoSeries(access_nodes.unary_union.convex_hull, crs="EPSG:4326")

# Calculate the population within the catchment area
temp_demand = demand.loc[demand['geometry'].centroid.within(access_nodes[0]), 'TotalPop'].sum()

# Calculate the number of hospital beds in each hospital
temp_supply = supply.loc[i, 'Total_Bed']

# Calculate the number of hospital beds available for 100,000 people 
supply_.at[i, 'ratio'] = temp_supply / temp_demand * 100000
supply_.at[i, 'ratio']<jupyter_output><empty_output><jupyter_text>The Enhanced 2SFCA method calculates supply-to-demand ratio (Step 1) as shown below.<jupyter_code>minutes = [5, 10, 15]
weights = {5: 1, 10: 0.68, 15: 0.22}

i= 0
supply_ = supply.copy(deep=True)
supply_['ratio'] = 0

# Create multiple catchment areas from a given location
ctmt_area = calculate_catchment_area(G, supply.loc[i, 'nearest_osm'], minutes)

# Calculate the population within each catchment area
ctmt_area_pops = 0
for c_idx, c_row in ctmt_area.iterrows():
    temp_pop = demand.loc[demand['geometry'].centroid.within(c_row['geometry']), 'TotalPop'].sum()
    print(f'Catchment area within {c_idx} minutes has {temp_pop} people and its weight is {weights[c_idx]}')
    ctmt_area_pops += temp_pop * weights[c_idx]
    print(f'Accumulated pop is {ctmt_area_pops}')

# Calculate the number of hospital beds in each hospital
temp_supply = supply.loc[i, 'Total_Bed']

# Calculate the number of hospital beds available for 100,000 people 
supply_.at[i, 'ratio'] = temp_supply / ctmt_area_pops * 100000
supply_.at[i, 'ratio']
ctmt_area.loc[5, 'geometry']
demand.loc[demand['geometry'].centroid.within(ctmt_area.loc[5, 'geometry'])].plot()
ctmt_area.loc[10, 'geometry']
demand.loc[demand['geometry'].centroid.within(ctmt_area.loc[10, 'geometry'])].plot()
ctmt_area.loc[15, 'geometry']
demand.loc[demand['geometry'].centroid.within(ctmt_area.loc[15, 'geometry'])].plot()<jupyter_output><empty_output><jupyter_text>In summary, we can define the functions for Enhanced 2SFCA method as shown below.<jupyter_code>def step1_E2SFCA(supply, supply_attr, demand, demand_attr, mobility, thresholds, weights):
    """
    Input: 
    - supply (GeoDataFrame): stores locations and attributes of supply 
    - supply_attr (str): the column of `supply` to be used for the analysis
    - demand (GeoDataFrame): stores locations and attributes of demand 
    - demand_attr (str): the column of `demand` to be used for the analysis
    - mobility (NetworkX MultiDiGraph): Network Dataset obtained from OSMnx
    - thresholds (list): the list of threshold travel times e.g., [5, 10, 15]
    - weights (dict): keys: threshold travel time, values: weigths according to the threshold travel times 
                      e.g., [5: 1, 10: 0.68, 15: 0.22]
    
    Output:
    - supply_ (GeoDataFrame): 
      a copy of supply and it stores supply-to-demand ratio of each supply at `ratio` column
    """
    
    # Your code here (Change the name of the variable according to the inputs)
    supply_ = supply.copy(deep=True)
    supply_['ratio'] = 0

    for i in trange(supply.shape[0]):

        # Create multiple catchment areas from a given location
        ctmt_area = calculate_catchment_area(mobility, supply.loc[i, 'nearest_osm'], thresholds)
    
        # Calculate the population within each catchment area
        ctmt_area_pops = 0
        for c_idx, c_row in ctmt_area.iterrows():
            temp_pop = demand.loc[demand['geometry'].centroid.within(c_row['geometry']), demand_attr].sum()
            ctmt_area_pops += temp_pop*weights[c_idx]

        # Calculate the number of hospital beds in each hospital
        temp_supply = supply.loc[i, supply_attr]

        # Calculate the number of hospital beds available for 100,000 people 
        supply_.at[i, 'ratio'] = temp_supply / ctmt_area_pops * 100000

    
    return supply_


def step2_E2SFCA(result_step1, demand, mobility, thresholds, weights):
    """
    Input: 
    - result_step1 (GeoDataFrame): stores locations and 'ratio' attribute that resulted in step1 
    - demand (GeoDataFrame): stores locations and attributes of demand 
    - mobility (NetworkX MultiDiGraph): Network Dataset obtained from OSMnx
    - thresholds (list): the list of threshold travel times e.g., [5, 10, 15]
    - weights (dict): keys: threshold travel time, values: weigths according to the threshold travel times 
                      e.g., [5: 1, 10: 0.68, 15: 0.22]
    
    Output:
    - demand_ (GeoDataFrame): 
      a copy of demand and it stores the final accessibility measures of each demand location at `ratio` column
    """
    
    # Your code here (Change the name of the variable according to the inputs)
    demand_ = demand.copy(deep=True)
    demand_['access'] = 0

    for j in trange(demand.shape[0]):
        ctmt_area = calculate_catchment_area(mobility, demand.loc[j, 'nearest_osm'], thresholds)

        ctmt_area_ratio = 0
        for c_idx, c_row in ctmt_area.iterrows():
            temp_ratio = result_step1.loc[result_step1['geometry'].centroid.within(c_row['geometry']), 'ratio'].sum()
            ctmt_area_ratio += temp_ratio * weights[c_idx]
            
        demand_.at[j, 'access'] = ctmt_area_ratio

    return demand_
minutes = [5, 10, 15]
weights = {5: 1, 10: 0.68, 15: 0.22}

E_step1 = step1_E2SFCA(supply, 'Total_Bed', demand, 'TotalPop', G, minutes, weights)
E_step2 = step2_E2SFCA(E_step1, demand, G, minutes, weights)
# Plotting accessibility measurement result. 
fig, ax = plt.subplots(figsize=(10,10))

hospitals.plot(markersize='Total_Bed', ax=ax, color='black', zorder=2)
E_step2.plot('access', ax=ax, figsize=(10,10), legend=True, cmap='Blues', scheme='FisherJenks')
E_step2.loc[step2['access'] == 0].plot(ax=ax, color='grey', zorder=1)
E_step2.boundary.plot(ax=ax, linestyle='dotted', lw=0.5, color='black', zorder=1)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11436758339302956, 'fraction_numerical': 0.01915470146696982, 'mean_word_length': 3.488118535085267, 'pattern_counts': {'":': 0, '<': 97, '<?xml version=': 0, '>': 97, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 30, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '501605', 'license': 'MIT', 'n_tokens_mistral': 10661, 'n_tokens_neox': 9868, 'n_words': 2711, 'path': 'Week11/Advanced_Spatial_Analysis_Spatial_Accessibility_2.ipynb', 'repo_name': 'jparkgeo/GEOG489', 'size': '32107'}