# Installing ALL the relevant libraries for this unified notebook

import os
import pandas as pd
import geopandas as gpd
import contextily as ctx
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib import lines as mlines
from matplotlib.colors import TwoSlopeNorm, LinearSegmentedColormap
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import numpy as np
import osmnx
import networkx as nx
from shapely.geometry import LineString, MultiLineString, Point, Polygon
from shapely.ops import nearest_points
import r5py
from datetime import datetime, timedelta
import zipfile
import shutil
import gtfs_kit as gk

zip_file_path = "OpenData_TTC_Schedules.zip"

# Directory where you want to extract all files
output_directory = "new-gtfs"

# Create the output directory if it doesn't exist
os.makedirs(output_directory, exist_ok=True)

# Open the zip file and extract all files
with zipfile.ZipFile(zip_file_path, 'r') as zip_ref:
    # Extract all files to the output directory
    zip_ref.extractall(output_directory)

# Pulling the Metrolinx rail shapefile, now under a different variable
mr = gpd.read_file(os.path.join(sq2, 'RTP_TRANSIT_NETWORK.shp'))

# Converting it to the UTM 17N NAD 1983 projection to facilitate accurate distance calculation                   
mr = mr.to_crs('EPSG:2958')

# Singling out the Ontario Line in the Metrolinx shapefile
ontario = mr[mr['NAME'] == 'Ontario Line']

# Getting the geometry of the Ontario Line
line = ontario.geometry.iloc[0]

# Returning Ontario Line back to WGS84 for visualisation purposes later
ontario = ontario.to_crs(4326)

# Since line is a MultiLineString, need to convert it to a LineString instead
all_coords = [coord for geom in line.geoms for coord in geom.coords]
ontario_linestring = LineString(all_coords)

ontario_linestring

# Error in the above LineString could be due to inclusion of erroneous coordinates at the start of the LineString
all_coords = all_coords[2:]  # Remove the first two coordinates
ontario_linestring = LineString(all_coords)

ontario_linestring

# Making a function to convert the LineString into a set of points with cumulative distances calculated

def line_to_shapes_with_distance(line, shape_id='a'):
    """
    Convert a LineString geometry into a DataFrame compatible with GTFS shapes.txt format,
    including cumulative distance.
    """
    
    data = []  # This is the empty table data
    sequence = 1  # shapes.txt in GTFS assigns '1' as the starting point of any new route
    cumulative_distance = 0

    # Process points in the LineString
    points = list(line.coords)
    for i in range(len(points)):
        x, y = points[i]
        coord = Point(x, y)

        if i > 0:  # Calculate Euclidean distance from the previous point
            prev_x, prev_y = points[i - 1]
            segment_distance = ((x - prev_x)**2 + (y - prev_y)**2)**0.5 / 1000  # Convert to km
            cumulative_distance += segment_distance

        data.append({
            "shape_id": shape_id,
            "shape_pt_sequence": sequence,
            "shape_dist_traveled": cumulative_distance,  # In kilometers
            'geometry': coord
        })
        sequence += 1

    return gpd.GeoDataFrame(data, crs='EPSG:2958')

# Converting the Ontario Line shape into a suitable format for shapes.txt file
ol_shapes = line_to_shapes_with_distance(ontario_linestring, shape_id = '10000')

# Re-converting it to WGS84 because the file needs to provide point coordinates in lat and lon
ol_shapes = ol_shapes.to_crs(4326)
ol_shapes['shape_pt_lat'] = ol_shapes.geometry.y
ol_shapes['shape_pt_lon'] = ol_shapes.geometry.x

# Limiting the columns present in the cleaned table
ol_shapes = ol_shapes[['shape_id', 'shape_pt_lat', 'shape_pt_lon', 'shape_pt_sequence', 'shape_dist_traveled']]
ol_shapes

# Locating the shapes.txt in the unzipped folder
shape_path = os.path.join(output_directory, 'shapes.txt')

# Reading the shapes.txt file
ttcshape = pd.read_csv(shape_path)

# Concatenating the Ontario Line shape data to the existing shapes.txt
totalshape = pd.concat([ttcshape, ol_shapes])

# Overwriting the shapes.txt file
totalshape.to_csv(shape_path, index = False)
totalshape

# Pulling the Metrolinx stations shapefile, now under a different variable
ms = gpd.read_file(os.path.join(sq2, 'RTP_POINTS.shp'))

# Setting the CRS
ms = ms.to_crs(4326)

# Singling out stations that are part of Ontario Line
ol_stops = ms[ms['NAME'] == 'Ontario Line']

# Creating columns to make it compatible with stops.txt file
ol_stops['stop_name'] = ol_stops['LOCATION_N']
ol_stops['stop_id'] = 'TBA'
ol_stops['stop_code'] = 'N.A.'
ol_stops['wheelchair_boarding'] = 1
ol_stops['stop_lat'] = ol_stops.geometry.y
ol_stops['stop_lon'] = ol_stops.geometry.x

C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)
C:\Users\Khalis\anaconda3\Lib\site-packages\geopandas\geodataframe.py:1819: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame.
Try using .loc[row_indexer,col_indexer] = value instead

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  super().__setitem__(key, value)

# Stations are ordered by manually adding in stop_ids
ol_stops.loc[ol_stops['stop_name'] == 'Exhibition', 'stop_id'] = 101
ol_stops.loc[ol_stops['stop_name'] == 'King / Bathurst', 'stop_id'] = 102
ol_stops.loc[ol_stops['stop_name'] == 'Queen / Spadina', 'stop_id'] = 103
ol_stops.loc[ol_stops['stop_name'] == 'Osgoode', 'stop_id'] = 104
ol_stops.loc[ol_stops['stop_name'] == 'Queen', 'stop_id'] = 105
ol_stops.loc[ol_stops['stop_name'] == 'Moss Park', 'stop_id'] = 106
ol_stops.loc[ol_stops['stop_name'] == 'Corktown', 'stop_id'] = 107
ol_stops.loc[ol_stops['stop_name'] == 'East Harbour', 'stop_id'] = 108
ol_stops.loc[ol_stops['stop_name'] == 'Riverside/Leslieville', 'stop_id'] = 109
ol_stops.loc[ol_stops['stop_name'] == 'Gerrard', 'stop_id'] = 110
ol_stops.loc[ol_stops['stop_name'] == 'Pape', 'stop_id'] = 111
ol_stops.loc[ol_stops['stop_name'] == 'Cosburn', 'stop_id'] = 112
ol_stops.loc[ol_stops['stop_name'] == 'Thorncliffe Park', 'stop_id'] = 113
ol_stops.loc[ol_stops['stop_name'] == 'Flemingdon Park', 'stop_id'] = 114
ol_stops.loc[ol_stops['stop_name'] == 'Science Centre', 'stop_id'] = 115

# Some interchanges are renamed to prevent possible confusion with the existing stops in stops.txt
ol_stops.loc[ol_stops['stop_id'] == 101, 'stop_name'] = 'Exhibition - Ontario Line'
ol_stops.loc[ol_stops['stop_id'] == 104, 'stop_name'] = 'Osgoode - Ontario Line'
ol_stops.loc[ol_stops['stop_id'] == 105, 'stop_name'] = 'Queen - Ontario Line'
ol_stops.loc[ol_stops['stop_id'] == 111, 'stop_name'] = 'Pape - Ontario Line'
ol_stops.loc[ol_stops['stop_id'] == 115, 'stop_name'] = 'Science Centre - Ontario Line'

# Ordering based on the stop_ids manually inserted earlier
stops_df = ol_stops.sort_values('stop_id', ascending = True)

# Limiting the columns present in the cleaned table
stops_df = stops_df[['stop_id', 'stop_code', 'stop_name', 'stop_lat', 'stop_lon', 'wheelchair_boarding']]
stops_df

# Locating the stops.txt in the unzipped folder
stop_path = os.path.join(output_directory, 'stops.txt')

# Reading the stops.txt file
ttcstops = pd.read_csv(stop_path)

# Concatenating the Ontario Line stop data to the existing stops.txt
totalstops = pd.concat([ttcstops, stops_df])

# Overwriting the stops.txt file
totalstops.to_csv(stop_path, index = False)
totalstops

# Creating a function to format stop times data, because TTC operating hours straddle between two calendar days

def format_gtfs_time(dt, start_of_service):
    """
    Convert a datetime object into GTFS-compliant time (HH:MM:SS).
    """
    elapsed = dt - start_of_service
    total_seconds = int(elapsed.total_seconds())
    hours = (total_seconds // 3600) + 6  # This code only records stop times as x hours into service, 
                                         # thus adding 6 to reflect the actual time since TTC starts at 0600hrs 
    minutes = (total_seconds % 3600) // 60
    seconds = total_seconds % 60
    
    return f"{hours:02}:{minutes:02}:{seconds:02}"

# Creating the stop times for all eastbound journeys from Exhibition to Science Centre

# Pulling the formatted stops data from earlier
stops = ol_stops.sort_values('stop_id', ascending = True)

# Converting it to UTM 17N NAD83 to calculate cumulative distance travelled at each stop
stops = stops.to_crs(2958)

# Snapping the station coordinates to the LineString in Step 2
stops['nearest_point'] = stops['geometry'].apply(
    lambda station: nearest_points(station, ontario_linestring)[1]
)

# Calculating cumulative distances travelled
stops['shape_dist_traveled'] = stops['nearest_point'].apply(
    lambda pt: ontario_linestring.project(pt) / 1000  
)

# Normalize shape_dist_traveled
total_distance = stops['shape_dist_traveled'].max()  # Total length of the route
stops['distance_ratio'] = stops['shape_dist_traveled'] / total_distance

# Parameters
journey_time_minutes = 45  # End-to-end journey time
boarding_time_per_stop = timedelta(seconds=30)  # 30 seconds boarding time
headway_minutes = 2  # Frequency of trains every 2 minutes

# Operating hours
start_time = datetime.strptime("06:00:00", "%H:%M:%S")  # Start of service
end_time = start_time + timedelta(hours=20)  # End of service (2 AM next day)

# Generate trip start times
trip_start_times = []
current_time = start_time
while current_time < end_time:
    trip_start_times.append(current_time)
    current_time += timedelta(minutes=headway_minutes)

# Adjust travel times based on distance ratios
stops['adjusted_travel_time'] = stops['distance_ratio'] * journey_time_minutes

# Generate stop_times.txt
stop_times = []
starting_trip_id = 50000  # Starting trip_id

for trip_index, trip_start_time in enumerate(trip_start_times):
    cumulative_time = trip_start_time
    trip_id = starting_trip_id + trip_index  # Increment trip ID for each trip
    
    for i, stop in enumerate(stops.itertuples()):
        if i == 0:
            # First stop: arrival time is the trip start time
            arrival_time = cumulative_time
        else:
            # Subsequent stops: calculate arrival time using adjusted travel time
            arrival_time = trip_start_time + timedelta(minutes=stop.adjusted_travel_time)
        
        departure_time = arrival_time + boarding_time_per_stop  # Add boarding time
        
        stop_times.append({
            "trip_id": trip_id,
            "arrival_time": format_gtfs_time(arrival_time, start_time),
            "departure_time": format_gtfs_time(departure_time, start_time),
            "stop_id": stop.stop_id,
            "stop_sequence": i + 1,
            "pickup_type": 0,
            "drop_off_type": 0,
            "shape_dist_traveled": stop.shape_dist_traveled
        })

# Saving the stop times as a DataFrame
stop_times_df = pd.DataFrame(stop_times)

# Creating the stop times for all westbound journeys from Science Centre to Exhibition

# Pulling the formatted stops data from earlier
rev_stops = ol_stops.sort_values('stop_id', ascending = False)

# Converting it to UTM 17N NAD83 to calculate cumulative distance travelled at each stop
rev_stops = rev_stops.to_crs(2958)

# Snapping the station coordinates to the LineString in Step 2
rev_stops['nearest_point'] = rev_stops['geometry'].apply(
    lambda station: nearest_points(station, ontario_linestring)[1]
)

# Calculating cumulative distances travelled
rev_stops['shape_dist_traveled'] = rev_stops['nearest_point'].apply(
    lambda pt: (ontario_linestring.length / 1000) - (ontario_linestring.project(pt) / 1000)  
)

# Normalize shape_dist_traveled
total_distance = rev_stops['shape_dist_traveled'].max()  # Total length of the route
rev_stops['distance_ratio'] = rev_stops['shape_dist_traveled'] / total_distance

# Parameters
journey_time_minutes = 45  # End-to-end journey time
boarding_time_per_stop = timedelta(seconds=30)  # 30 seconds boarding time
headway_minutes = 2  # Frequency of trains every 2 minutes

# Operating hours
start_time = datetime.strptime("06:00:00", "%H:%M:%S")  # Start of service
end_time = start_time + timedelta(hours=20)  # End of service (2 AM next day)

# Generate trip start times
trip_start_times = []
current_time = start_time
while current_time < end_time:
    trip_start_times.append(current_time)
    current_time += timedelta(minutes=headway_minutes)

# Adjust travel times based on distance ratios
rev_stops['adjusted_travel_time'] = rev_stops['distance_ratio'] * journey_time_minutes

# Generate stop_times.txt
rev_stop_times = []
starting_trip_id = 50600  # Starting trip_id

for trip_index, trip_start_time in enumerate(trip_start_times):
    cumulative_time = trip_start_time
    trip_id = starting_trip_id + trip_index  # Increment trip ID for each trip
    
    for i, stop in enumerate(rev_stops.itertuples()):
        if i == 0:
            # First stop: arrival time is the trip start time
            arrival_time = cumulative_time
        else:
            # Subsequent stops: calculate arrival time using adjusted travel time
            arrival_time = trip_start_time + timedelta(minutes=stop.adjusted_travel_time)
        
        departure_time = arrival_time + boarding_time_per_stop  # Add boarding time
        
        rev_stop_times.append({
            "trip_id": trip_id,
            "arrival_time": format_gtfs_time(arrival_time, start_time),
            "departure_time": format_gtfs_time(departure_time, start_time),
            "stop_id": stop.stop_id,
            "stop_sequence": i + 1,
            "pickup_type": 0,
            "drop_off_type": 0,
            "shape_dist_traveled": stop.shape_dist_traveled
        })

# Saving the stop times as a DataFrame
rev_stop_times_df = pd.DataFrame(rev_stop_times)

# Putting the eastward and westward stop times into a combined stop time table
ol_stop_times = pd.concat([stop_times_df, rev_stop_times_df])
ol_stop_times

# Locating the stop_times.txt in the unzipped folder
stoptime_path = os.path.join(output_directory, 'stop_times.txt')

# Reading the stop_times.txt file
ttcstoptimes = pd.read_csv(stoptime_path)

# Concatenating the Ontario Line stop times data to the existing stop_times.txt
totalstoptimes = pd.concat([ttcstoptimes, ol_stop_times])

# Overwriting the stop_times.txt file
totalstoptimes.to_csv(stoptime_path, index = False)
totalstoptimes

# Extract unique trip_ids from ol_stop_times
unique_ids = ol_stop_times['trip_id'].unique()

# Define additional fields for trips.txt
service_id = "1"  # Assuming daily service
route_id = "3000"  # Matches the route_id from your route DataFrame

# Create trips data based on unique trip_ids
trips = []
for trip_id in unique_ids:
    if trip_id <= 50599:
        trips.append({
            "route_id": route_id,
            "service_id": service_id,
            "trip_id": trip_id,
            "trip_headsign": "ONTARIO LINE towards SCIENCE CENTRE",  # Customize as needed
            "direction_id": 0,
            "shape_id": 10000,
            "wheelchair_accessible": 1,
            "bikes_allowed": 1
        })
    else:
        trips.append({
            "route_id": route_id,
            "service_id": service_id,
            "trip_id": trip_id,
            "trip_headsign": "ONTARIO LINE towards EXHIBITION",  # Customize as needed
            "direction_id": 1,
            "shape_id": 10000,
            "wheelchair_accessible": 1,
            "bikes_allowed": 1
        })

# Convert to a DataFrame
trips_df = pd.DataFrame(trips)
trips_df

# Locating the trips.txt in the unzipped folder
trip_path = os.path.join(output_directory, 'trips.txt')

# Reading the trips.txt file
ttctrips = pd.read_csv(trip_path)

# Concatenating the Ontario Line trips data to the existing trips.txt
totaltrips = pd.concat([ttctrips, trips_df])

# Overwriting the trips.txt file
totaltrips.to_csv(trip_path, index = False)
totaltrips

# Making the Ontario Line Route details
ol_route = pd.DataFrame({'route_id':['3000'],
                         'agency_id': ['1'],
                         'route_short_name': ['3'],
                         'route_long_name': ['ONTARIO LINE'],
                         'route_type': ['1'],
                         'route_color': ['008000'],
                         'route_text_color': ['FFFFFF']
                        }
                       )

# Locating the routes.txt in the unzipped folder
route_path = os.path.join(output_directory, 'routes.txt')

# Reading the routes.txt file
ttcroutes = pd.read_csv(route_path)

# Concatenating the Ontario Line route details to the existing routes.txt
totalroutes = pd.concat([ttcroutes, ol_route])

# Overwriting the routes.txt file
totalroutes.to_csv(route_path, index = False)
totalroutes

# Path for the output zip file
output_directory_zipped = "new-gtfs.zip"

# Create a zip file of the folder
shutil.make_archive(output_directory_zipped.replace('.zip', ''), 'zip', output_directory)

'C:\\Users\\Khalis\\Documents\\UofT Matters\\2024-2025\\GIS and Programming\\FINAL PROJECT\\GGR375-Kasandra-Muhammad-William-CodeEnvironment\\new-gtfs.zip'

old_feed = gk.read_feed('OpenData_TTC_Schedules.zip', dist_units = 'km')
old_feed.describe()

new_feed = gk.read_feed('new-gtfs.zip', dist_units = 'km')
new_feed.describe()

# Loading the Metrolinx stations shapefile data under a different variable
new_stn = gpd.read_file(os.path.join(sq2, 'RTP_POINTS.shp'))

# Filtering for Line 1 stations along Yonge and Line 2 Stations along Bloor, intersecting at Bloor-Yonge
new_stn = new_stn[
    new_stn['LOCATION_N'].isin([
        'King', 'Queen', 'Dundas', 'College', 'Wellesley', 
        'Bloor-Yonge', 'Sherbourne', 'Castle Frank',
        'Moss Park', 'Corktown', 'East Harbour'
    ]) & 
    new_stn['NAME'].isin(['Line 1: Yonge-University Subway', 
                          'Line 2: Bloor-Danforth Subway',
                          'Ontario Line'
                         ]
                        )
]

new_stn = new_stn.to_crs("EPSG:4326") # Changing the crs

# Getting the nearest centre nodes to subway stops in walking network

nearest_nodesTTC = []
for geom in new_stn.geometry :  # Loop through the subway stops
    nearest_node = osmnx.distance.nearest_nodes(G, X=geom.x, Y=geom.y) # Extracting the xy coords to get the nearest nodes on the OSM network
    nearest_nodesTTC.append(nearest_node)

center_nodes = nearest_nodesTTC

# Making the 5min-isochrones from TTC subway stops

TTCstop_5min_iso = {}  

for center_node in center_nodes:
    subgraph = nx.ego_graph(G, 
                            center_node, 
                            radius=fivemin*60, 
                            distance='walk_time') 
    
    node_points = [Point((data['x'], data['y'])) \
                   for node, data in subgraph.nodes(data=True)]
    
    polygon = Polygon(gpd.GeoSeries(node_points).union_all().convex_hull)
    
    TTCstop_5min_iso[center_node] = polygon  

TTCstop_5min_union = gpd.GeoSeries(TTCstop_5min_iso.values()).union_all()
TTCstop_5min_union = gpd.GeoDataFrame(geometry=[TTCstop_5min_union], crs='EPSG:4326')

# Making the 10min-isochrones from TTC subway stops

TTCstop_10min_iso = {}

for center_node in center_nodes:
    subgraph = nx.ego_graph(G, 
                            center_node, 
                            radius=tenmin*60, 
                            distance='walk_time')
    
    node_points = [Point((data['x'], data['y'])) \
                   for node, data in subgraph.nodes(data=True)]
    
    polygon = Polygon(gpd.GeoSeries(node_points).union_all().convex_hull)
    
    TTCstop_10min_iso[center_node] = polygon

TTCstop_10min_union = gpd.GeoSeries(TTCstop_10min_iso.values()).union_all()
TTCstop_10min_union = gpd.GeoDataFrame(geometry=[TTCstop_10min_union], crs='EPSG:4326')

# Making the 20min-isochrones from TTC subway stops

TTCstop_20min_iso = {}

for center_node in center_nodes:
    subgraph = nx.ego_graph(G, 
                            center_node, 
                            radius=twentymin*60, 
                            distance='walk_time')
    
    node_points = [Point((data['x'], data['y'])) \
                   for node, data in subgraph.nodes(data=True)]
    
    polygon = Polygon(gpd.GeoSeries(node_points).union_all().convex_hull)
    
    TTCstop_20min_iso[center_node] = polygon

TTCstop_20min_union = gpd.GeoSeries(TTCstop_20min_iso.values()).union_all()
TTCstop_20min_union = gpd.GeoDataFrame(geometry=[TTCstop_20min_union], crs='EPSG:4326')

# Visualising all walking isochrones

NEW_TTC_walk_unified_fig, ax = plt.subplots(figsize = (15,15))

TTCstop_20min_union.plot(ax=ax, color= 'coral', alpha=0.7, edgecolor="white")
TTCstop_10min_union.plot(ax=ax, color= 'yellow', alpha=0.7, edgecolor="white")
TTCstop_5min_union.plot(ax=ax, color= 'green', alpha=0.7, edgecolor="white") 
osmnx.plot_graph(G, ax=ax, node_size=0, edge_linewidth=0.5, show=False)
moss_park.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=4)
ttc.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax, edgecolor='blue', facecolor='none', linewidth=2)
new_stn.plot(ax=ax, color="red", edgecolor="black", markersize=200)

# Manually create legend handles
ttc_handle = mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line')
ontario_handle = mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line')
stops_handle = mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=10, label='Subway Stop', linewidth = 0)
fivemin_handle = mpatches.Patch(color='green', label='5 min walking isochrone')
tenmin_handle = mpatches.Patch(color='yellow', label='10 min walking isochrone')
twentymin_handle = mpatches.Patch(color='coral', label='20 min walking isochrone')
handles = [fivemin_handle, tenmin_handle, twentymin_handle, ttc_handle, ontario_handle, stops_handle]

NEW_TTC_walk_unified_fig.suptitle('Walking Isochrones from TTC Subway Stops with Ontario Line', fontsize=20, fontweight="bold")
ax.legend(handles=handles, loc='lower right', fontsize=16)
ax.axis('off')

plt.tight_layout(rect=[0, 0, 1, 0.95])  # Adjust layout to fit suptitle and legend
plt.savefig("FUTURE_TTC_walking_maps_unified.png", dpi=300, bbox_inches="tight")
plt.show()

# First, get the nearest centre nodes to subway stops in cycling network

nearest_nodesTTCbike = []
for geom in new_stn.geometry :  # Loop through the subway stops
    nearest_node = osmnx.distance.nearest_nodes(D, X=geom.x, Y=geom.y) # Extracting the xy coords to get the nearest nodes on the OSM network
    nearest_nodesTTCbike.append(nearest_node)

center_nodes = nearest_nodesTTCbike

# Making the 5min-isochrones from TTC subway stops

TTCstopbike_5min_iso = {}

for center_node in center_nodes:
    subgraph = nx.ego_graph(D, 
                            center_node, 
                            radius=fivemin*60, 
                            distance='cycle_time')
    
    node_points = [Point((data['x'], data['y'])) \
                   for node, data in subgraph.nodes(data=True)]
    
    polygon = Polygon(gpd.GeoSeries(node_points).union_all().convex_hull)
    
    TTCstopbike_5min_iso[center_node] = polygon

TTCstopbike_5min_union = gpd.GeoSeries(TTCstopbike_5min_iso.values()).union_all()
TTCstopbike_5min_union = gpd.GeoDataFrame(geometry=[TTCstopbike_5min_union], crs='EPSG:4326')

# Visualising all cycling isochrones

NEW_TTC_cycle_unified_fig, ax = plt.subplots(figsize = (15,15))

TTCstopbike_5min_union.plot(ax=ax, color= 'green', alpha=0.7, edgecolor="white") 
osmnx.plot_graph(D, ax=ax, node_size=0, edge_linewidth=0.5, show=False)
moss_park.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=4)
ttc.plot(ax=ax, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax, edgecolor='blue', facecolor='none', linewidth=2)
new_stn.plot(ax=ax, color="red", edgecolor="black", markersize=200)

# Manually create legend handles
ttc_handle = mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line')
ontario_handle = mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line')
stops_handle = mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=10, label='Subway Stop', linewidth = 0)
fivemin_handle = mpatches.Patch(color='green', label='5 min cycling isochrone')
handles = [fivemin_handle, ttc_handle, ontario_handle, stops_handle]

NEW_TTC_cycle_unified_fig.suptitle('Cycling Isochrones from TTC Subway Stops with Ontario Line', fontsize=20, fontweight="bold")
ax.legend(handles=handles, loc='lower right', fontsize=16)
ax.axis('off')

plt.tight_layout(rect=[0, 0, 1, 0.95])  # Adjust layout to fit suptitle and legend
plt.savefig("FUTURE_TTC_cycling_maps_unified.png", dpi=300, bbox_inches="tight")
plt.show()

# This part will again rely on r5py library to do network analysis. Documentation can be accessed via this link: https://r5py.readthedocs.io/en/stable/

# First, set the network data that will be used by the library to conduct travel analysis
new_network = r5py.TransportNetwork('Toronto.osm.pbf', # OSM Walking network in a pbf, obtained from https://download.bbbike.org/osm/bbbike/.
                                    'new-gtfs.zip') # The new GTFS created

# Setting up the environment to measure travel times from Moss Park DA centroids to Toronto Neighbourhood centroids based on new transport network

matrix3 = r5py.TravelTimeMatrixComputer(new_network,
                                        origins = origins,  # Refer to Sub-Question 2, Part I, Step 5
                                        destinations = destinations,  # Refer to Sub-Question 2, Part I, Step 5
                                        departure = datetime(2024, 11, 18, 8, 00), # Assuming that one starts travel on 18 November 2024, at 8 am
                                        transport_modes = [r5py.TransportMode.WALK,
                                                           r5py.TransportMode.TRANSIT] # Assuming that one EITHER walks OR walk and take the TTC
                                       )

# The command to compute travel times based on the above matrix
tt_3 = matrix3.compute_travel_times()

# The result is a table of one origin to one destination with its travel time. 
# Pivoting allows one to compare travel times based on a common origin or destination
tt_3_pivot = tt_3.pivot(index="to_id", columns="from_id", values="travel_time")

# Pivoted table is merged with the Toronto neighbourhood GeoDataFrame to facilitate mapping
new_scenario = tor_gdf.merge(tt_3_pivot,
                             left_on = 'AREA_ID',
                             right_on = 'to_id')

# Using the columns to identify average travel time between Toronto Neighbourhood centroids and all the Moss Park DA centroids
# Refer to Sub-Question 2, Part I, Step 5 for identification of Moss Park DAs
new_scenario['avg_travel_time'] = new_scenario[mosspark_DAs].mean(axis = 1)

# Categorising them into intervals of 15 mins
bins = [0, 15, 30, 45, 60, float('inf')]
labels = ['<15 min', '15-30 min', '30-45 min', '45-60 min', '>60 min']
new_scenario ['classify'] = pd.cut(new_scenario['avg_travel_time'],
                                   bins = bins,
                                   labels = labels,
                                   right = False)

# Essentially repeating the four cells above, with the difference being the inclusion that one may ALSO cycle along the whole journey, or parts of it
matrix4 = r5py.TravelTimeMatrixComputer(new_network,
                                        origins = origins,
                                        destinations = destinations,
                                        departure = datetime(2024, 11, 18, 8, 00),
                                        transport_modes = [r5py.TransportMode.WALK,
                                                           r5py.TransportMode.BICYCLE,
                                                           r5py.TransportMode.TRANSIT]
                                       )

tt_4 = matrix4.compute_travel_times()

tt_4_pivot = tt_4.pivot(index="to_id", columns="from_id", values="travel_time")

new_scenario_bikes = tor_gdf.merge(tt_4_pivot,
                                   left_on = 'AREA_ID',
                                   right_on = 'to_id')

new_scenario_bikes['avg_travel_time'] = new_scenario_bikes[mosspark_DAs].mean(axis = 1)

bins = [0, 15, 30, 45, 60, float('inf')]
labels = ['<15 min', '15-30 min', '30-45 min', '45-60 min', '>60 min']
new_scenario_bikes['classify'] = pd.cut(new_scenario_bikes['avg_travel_time'],
                                        bins = bins,
                                        labels = labels,
                                        right = False)

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_TTC_acc, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, color='black', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
ax1.set_title('Walking and TTC only')
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, color='black', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
ax2.set_title('Walking, Cycling and TTC')
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min'),
    mpatches.Patch(color='#78c679', label='15-30 min'),
    mpatches.Patch(color='#ffffcc', label='30-45 min'),
    mpatches.Patch(color='#fdae61', label='45-60 min'),
    mpatches.Patch(color='#a50026', label='>60 min'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line')
]

# Placing the custom legend in the middle
new_TTC_acc.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.5),
    fontsize=8,
    title="Average Travel Time",
    title_fontsize=10
)


new_TTC_acc.suptitle('Possible Multi-Modal City Accessibility from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_city_access.png", dpi=300, bbox_inches="tight")
plt.show()

# From Sub-Question 3, Part II, Step 1, filter the walk-TTC network results to show journeys of 30 min and less
new_scenario_30min = new_scenario[new_scenario['classify'].isin(['<15 min', '15-30 min'])]

# Getting the dissolved geometry to facilitate OSM retrieval
new_scen_30min_poly = new_scenario_30min.union_all()

# Conducting the OSM retrieval
new_scen_30min_clinic = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_clinic)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_clinic = gpd.GeoDataFrame({"id": range(len(new_scen_30min_clinic)), 
                                          "geometry": new_scen_30min_clinic.geometry.centroid}, 
                                         crs='EPSG:4326'
                                        )

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\597819928.py:12: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_clinic.geometry.centroid},

# From Sub-Question 3, Part II, Step 1, filter the walk-TTC network results to show journeys of 30 min and less
new_scenario_bikes_30min = new_scenario_bikes[new_scenario_bikes['classify'].isin(['<15 min', '15-30 min'])]

# Getting the dissolved geometry to facilitate OSM retrieval
new_scen_bikes_30min_poly = new_scenario_bikes_30min.union_all()

# Conducting the OSM retrieval
new_scen_bikes_30min_clinic = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_clinic)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_clinic = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_clinic)), 
                                                "geometry": new_scen_bikes_30min_clinic.geometry.centroid}, 
                                               crs='EPSG:4326'
                                              )

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\2298355300.py:12: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_clinic.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_clinic_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_clinic.plot(ax=ax1, color="yellow", edgecolor="black", markersize=100)
old_scen_30min_clinic.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_clinic) - len(old_scen_30min_clinic)} more clinics within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_clinic.plot(ax=ax2, color="yellow", edgecolor="black", markersize=100)
old_scen_bikes_30min_clinic.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_clinic) - len(old_scen_bikes_30min_clinic)} more clinics within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Clinic', linewidth = 0),
    mlines.Line2D([], [], marker='o', color='yellow',markeredgecolor="black", markersize=8, label='New Clinic', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_clinic_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_clinic_30min.suptitle('Difference in Multi-Modal Accessibility to Clinics from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_clinic_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Conducting the OSM retrieval, using new_scen_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_30min_hospital = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_hospital)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_hospital = gpd.GeoDataFrame({"id": range(len(new_scen_30min_hospital)), 
                                          "geometry": new_scen_30min_hospital.geometry.centroid}, 
                                         crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\1160977940.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_hospital.geometry.centroid},

# Conducting the OSM retrieval, using new_scen_bikes_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_hospital = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_hospital)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_hospital = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_hospital)), 
                                                "geometry": new_scen_bikes_30min_hospital.geometry.centroid}, 
                                               crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\3910666539.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_hospital.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_hospital_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_hospital.plot(ax=ax1, color="yellow", edgecolor="black", markersize=100)
old_scen_30min_hospital.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_hospital) - len(old_scen_30min_hospital)} more hospitals within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_hospital.plot(ax=ax2, color="yellow", edgecolor="black", markersize=100)
old_scen_bikes_30min_hospital.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_hospital) - len(old_scen_bikes_30min_hospital)} more hospitals within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Hospital', linewidth = 0),
    mlines.Line2D([], [], marker='o', color='yellow',markeredgecolor="black", markersize=8, label='New Hospital', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_hospital_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_hospital_30min.suptitle('Difference in Multi-Modal Accessibility to Hospitals from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_hospital_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Conducting the OSM retrieval, using new_scen_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_30min_grocery = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_grocery)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_grocery = gpd.GeoDataFrame({"id": range(len(new_scen_30min_grocery)), 
                                          "geometry": new_scen_30min_grocery.geometry.centroid}, 
                                         crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\1977468625.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_grocery.geometry.centroid},

# Conducting the OSM retrieval, using new_scen_bikes_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_grocery = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_grocery)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_grocery = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_grocery)), 
                                                "geometry": new_scen_bikes_30min_grocery.geometry.centroid}, 
                                               crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\3842587446.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_grocery.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_grocery_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_grocery.plot(ax=ax1, color="yellow", edgecolor="black", markersize=100)
old_scen_30min_grocery.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_grocery) - len(old_scen_30min_grocery)} more grocery stores within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_grocery.plot(ax=ax2, color="yellow", edgecolor="black", markersize=100)
old_scen_bikes_30min_grocery.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_grocery) - len(old_scen_bikes_30min_grocery)} more grocery stores within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Grocery Store', linewidth = 0),
    mlines.Line2D([], [], marker='o', color='yellow',markeredgecolor="black", markersize=8, label='New Grocery Store', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_grocery_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_grocery_30min.suptitle('Difference in Multi-Modal Accessibility to Grocery Stores from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_grocery_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Conducting the OSM retrieval, using new_scen_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_30min_pow = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_pow)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_pow = gpd.GeoDataFrame({"id": range(len(new_scen_30min_pow)), 
                                          "geometry": new_scen_30min_pow.geometry.centroid}, 
                                         crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\4001689017.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_pow.geometry.centroid},

# Conducting the OSM retrieval, using new_scen_bikes_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_pow = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_pow)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_pow = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_pow)), 
                                                "geometry": new_scen_bikes_30min_pow.geometry.centroid}, 
                                               crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\2215535961.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_pow.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_pow_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_pow.plot(ax=ax1, color="yellow", edgecolor="black", markersize=100)
old_scen_30min_pow.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_pow) - len(old_scen_30min_pow)} more places of worship within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_pow.plot(ax=ax2, color="yellow", edgecolor="black", markersize=100)
old_scen_bikes_30min_pow.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_pow) - len(old_scen_bikes_30min_pow)} more places of worship within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Place of Worship', linewidth = 0),
    mlines.Line2D([], [], marker='o', color='yellow',markeredgecolor="black", markersize=8, label='New Place of Worship', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_pow_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_pow_30min.suptitle('Difference in Multi-Modal Accessibility to Places of Worship from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_pow_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Conducting the OSM retrieval, using new_scen_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_30min_lib = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_lib)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_lib = gpd.GeoDataFrame({"id": range(len(new_scen_30min_lib)), 
                                          "geometry": new_scen_30min_lib.geometry.centroid}, 
                                         crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\3133954879.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_lib.geometry.centroid},

# Conducting the OSM retrieval, using new_scen_bikes_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_lib = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_lib)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_lib = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_lib)), 
                                                "geometry": new_scen_bikes_30min_lib.geometry.centroid}, 
                                               crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\884825382.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_lib.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_lib_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_lib.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_lib) - len(old_scen_30min_lib)} more libraries within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_lib.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_lib) - len(old_scen_bikes_30min_lib)} more libraries within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Library', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_lib_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_lib_30min.suptitle('Difference in Multi-Modal Accessibility to Libraries from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_library_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Conducting the OSM retrieval, using new_scen_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_30min_sf = osmnx.features.features_from_polygon(new_scen_30min_poly, tag_sf)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_30min_sf = gpd.GeoDataFrame({"id": range(len(new_scen_30min_sf)), 
                                          "geometry": new_scen_30min_sf.geometry.centroid}, 
                                         crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\2716691422.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_30min_sf.geometry.centroid},

# Conducting the OSM retrieval, using new_scen_bikes_30min_poly created in Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_sf = osmnx.features.features_from_polygon(new_scen_bikes_30min_poly, tag_sf)

# Cleaning into a GeoDataFrame of points rather than a mix, for visualisation consistency
new_scen_bikes_30min_sf = gpd.GeoDataFrame({"id": range(len(new_scen_bikes_30min_sf)), 
                                                "geometry": new_scen_bikes_30min_sf.geometry.centroid}, 
                                               crs='EPSG:4326')

C:\Users\Khalis\AppData\Local\Temp\ipykernel_34552\3060287667.py:6: UserWarning: Geometry is in a geographic CRS. Results from 'centroid' are likely incorrect. Use 'GeoSeries.to_crs()' to re-project geometries to a projected CRS before this operation.

  "geometry": new_scen_bikes_30min_sf.geometry.centroid},

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_sf_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_sf.plot(ax=ax1, color="yellow", edgecolor="black", markersize=100)
old_scen_30min_sf.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_sf) - len(old_scen_30min_sf)} more social facilities within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_sf.plot(ax=ax2, color="yellow", edgecolor="black", markersize=100)
old_scen_bikes_30min_sf.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_sf) - len(old_scen_bikes_30min_sf)} more social facilities within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=8, label='Social Facility', linewidth = 0),
    mlines.Line2D([], [], marker='o', color='yellow',markeredgecolor="black", markersize=8, label='New Social Facility', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_sf_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_sf_30min.suptitle('Difference in Multi-Modal Accessibility to Social Facilities from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_sf_access.png", dpi=300, bbox_inches="tight")
plt.show()

# Data is clipped based on journeys of 30 min and less on the walk-TTC network, refer to Sub-Question 3, Part II, Step 4
new_scen_30min_wifi = wifi.clip(new_scenario_30min)

# Data is clipped based on journeys of 30 min and less on the walk-cycle-TTC network, refer to Sub-Question 3, Part II, Step 4
new_scen_bikes_30min_wifi = wifi.clip(new_scenario_bikes_30min)

# Putting the walk-TTC only network analysis and walk-cycle-TTC network analysis side by side

new_multimodal_wifi_30min, (ax1, ax2) = plt.subplots(1, 2, figsize = (15,8))


# The left side will be the walk-TTC network analysis
new_scenario_30min.plot(column='classify', ax=ax1, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax1, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax1, color='blue', linewidth=2)
moss_park.plot(ax=ax1, facecolor='grey')
new_scen_30min_wifi.plot(ax=ax1, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_30min. Needed to prevent the map from displaying the full TTC network map
ax1_bounds = new_scenario_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax1.set_xlim(ax1_bounds[0], ax1_bounds[2])  # min_x to max_x
ax1.set_ylim(ax1_bounds[1], ax1_bounds[3])  # min_y to max_y

ax1.set_title(f"{len(new_scen_30min_wifi) - len(old_scen_30min_wifi)} more public Wifi spots within 30 min of walk-TTC journeys")
ax1.axis('off')


# The right side will be walk-cycle-TTC network analysis
new_scenario_bikes_30min.plot(column='classify', ax=ax2, cmap = 'RdYlGn_r', legend=False)
ttc.plot(ax=ax2, edgecolor='black', facecolor='none', linewidth=2)
ontario.plot(ax=ax2, color='blue', linewidth=2)
moss_park.plot(ax=ax2, facecolor='grey')
new_scen_bikes_30min_wifi.plot(ax=ax2, color="red", edgecolor="black", markersize=100)

# Set axis limits based on the bounds of old_scenario_bikes_30min. Needed to prevent the map from displaying the full TTC network map.
ax2_bounds = new_scenario_bikes_30min.total_bounds  # (min_x, min_y, max_x, max_y)

# Apply bounds to the axis
ax2.set_xlim(ax2_bounds[0], ax2_bounds[2])  # min_x to max_x
ax2.set_ylim(ax2_bounds[1], ax2_bounds[3])  # min_y to max_y

ax2.set_title(f"{len(new_scen_bikes_30min_wifi) - len(old_scen_bikes_30min_wifi)} more public Wifi spots within 30 min of walk-cycle-TTC journeys")
ax2.axis('off')


# Creating custom legend elements because both maps have the same colour scheme and classification
legend_handles = [
    mpatches.Patch(color='#006837', label='<15 min travel'),
    mpatches.Patch(color='#78c679', label='15-30 min travel'),
    mlines.Line2D([], [], color='black', linewidth=3, label='Existing Subway Line'),
    mlines.Line2D([], [], color='blue', linewidth=3, label='Ontario Line'),
    mlines.Line2D([], [], marker='o', color='red',markeredgecolor="black", markersize=10, label='Public Wi-Fi Spot', linewidth = 0)
]

# Placing the custom legend in the middle
new_multimodal_wifi_30min.legend(
    handles=legend_handles,
    loc='center',
    bbox_to_anchor=(0.5, 0.25),
    fontsize=8,
    title="Legend",
    title_fontsize=10
)


new_multimodal_wifi_30min.suptitle('Difference in Multi-Modal Accessibility to Public Wi-Fi Spots from Moss Park (in grey)', fontsize=20, fontweight="bold")
plt.savefig("future_wifi_access.png", dpi=300, bbox_inches="tight")
plt.show()

	shape_id	shape_pt_lat	shape_pt_lon	shape_pt_sequence	shape_dist_traveled
0	10000	43.636378	-79.418197	1	0.000000
1	10000	43.636947	-79.416117	2	0.179373
2	10000	43.636950	-79.416106	3	0.180341
3	10000	43.636953	-79.416094	4	0.181308
4	10000	43.636956	-79.416083	5	0.182276
...	...	...	...	...	...
550	10000	43.720298	-79.338589	551	15.263804
551	10000	43.720384	-79.338611	552	15.273481
552	10000	43.720470	-79.338633	553	15.283157
553	10000	43.720541	-79.338654	554	15.291221
554	10000	43.720569	-79.338662	555	15.294458

	shape_id	shape_pt_lat	shape_pt_lon	shape_pt_sequence	shape_dist_traveled
0	1035570	43.775633	-79.346844	1	0.000000
1	1035570	43.775924	-79.346974	2	0.033500
2	1035570	43.776032	-79.346990	3	0.045700
3	1035570	43.776118	-79.346959	4	0.056100
4	1035570	43.776160	-79.346917	5	0.061100
...	...	...	...	...	...
550	10000	43.720298	-79.338589	551	15.263804
551	10000	43.720384	-79.338611	552	15.273481
552	10000	43.720470	-79.338633	553	15.283157
553	10000	43.720541	-79.338654	554	15.291221
554	10000	43.720569	-79.338662	555	15.294458

	stop_id	stop_code	stop_name	stop_desc	stop_lat	stop_lon	zone_id	stop_url	location_type	parent_station	stop_timezone	wheelchair_boarding
0	262	662	Danforth Rd at Kennedy Rd	NaN	43.714379	-79.260939	NaN	NaN	NaN	NaN	NaN	1
1	263	929	Davenport Rd at Bedford Rd	NaN	43.674448	-79.399659	NaN	NaN	NaN	NaN	NaN	1
2	264	940	Davenport Rd at Dupont St	NaN	43.675511	-79.401938	NaN	NaN	NaN	NaN	NaN	2
3	265	1871	Davisville Ave at Cleveland St	NaN	43.702088	-79.378112	NaN	NaN	NaN	NaN	NaN	1
4	266	11700	Disco Rd at Attwell Dr	NaN	43.701362	-79.594843	NaN	NaN	NaN	NaN	NaN	1
...	...	...	...	...	...	...	...	...	...	...	...	...
340	111	N.A.	Pape - Ontario Line	NaN	43.679727	-79.345215	NaN	NaN	NaN	NaN	NaN	1
350	112	N.A.	Cosburn	NaN	43.689400	-79.348952	NaN	NaN	NaN	NaN	NaN	1
351	113	N.A.	Thorncliffe Park	NaN	43.705226	-79.350054	NaN	NaN	NaN	NaN	NaN	1
352	114	N.A.	Flemingdon Park	NaN	43.715024	-79.337171	NaN	NaN	NaN	NaN	NaN	1
341	115	N.A.	Science Centre - Ontario Line	NaN	43.720569	-79.338662	NaN	NaN	NaN	NaN	NaN	1

	trip_id	arrival_time	departure_time	stop_id	stop_sequence	pickup_type	drop_off_type	shape_dist_traveled
0	50000	06:00:00	06:00:30	101	1	0	0	0.000000
1	50000	06:04:38	06:05:08	102	2	0	0	1.579833
2	50000	06:06:57	06:07:27	103	3	0	0	2.363818
3	50000	06:09:22	06:09:52	104	4	0	0	3.184820
4	50000	06:11:11	06:11:41	105	5	0	0	3.803558
...	...	...	...	...	...	...	...	...
8995	51199	26:31:48	26:32:18	105	11	0	0	11.490900
8996	51199	26:33:37	26:34:07	104	12	0	0	12.109638
8997	51199	26:36:02	26:36:32	103	13	0	0	12.930640
8998	51199	26:38:21	26:38:51	102	14	0	0	13.714625
8999	51199	26:43:00	26:43:30	101	15	0	0	15.294458

	trip_id	arrival_time	departure_time	stop_id	stop_sequence	stop_headsign	pickup_type	drop_off_type	shape_dist_traveled
0	47907132	18:45:00	18:45:00	15182	1	NaN	0	0	NaN
1	47907132	18:45:56	18:45:56	3807	2	NaN	0	0	0.284000
2	47907132	18:46:42	18:46:42	6904	3	NaN	0	0	0.519700
3	47907132	18:47:55	18:47:55	1163	4	NaN	0	0	0.890700
4	47907132	18:48:58	18:48:58	7723	5	NaN	0	0	1.214300
...	...	...	...	...	...	...	...	...	...
8995	51199	26:31:48	26:32:18	105	11	NaN	0	0	11.490900
8996	51199	26:33:37	26:34:07	104	12	NaN	0	0	12.109638
8997	51199	26:36:02	26:36:32	103	13	NaN	0	0	12.930640
8998	51199	26:38:21	26:38:51	102	14	NaN	0	0	13.714625
8999	51199	26:43:00	26:43:30	101	15	NaN	0	0	15.294458

Kasandra Tworzyanski, Muhammad Khalis, William Chan

GGR375 Introduction to Programming in GIS - Final Project

Due: December 11, 2024

Research Question: How would the connectivity for residents of Moss Park, Toronto, change as a result of the Ontario Line?

Sub-Questions¶

Sub-Question 3: How will this change if Ontario Line is open for service today?¶

Part I: Cleaning Ontario Line Data to Make New GTFS¶

1. Extracting Files from Existing TTC GTFS Zipped Folder¶

2. Cleaning Ontario Line Shape Data and Adding it to Relevant .txt File¶

3. Cleaning Ontario Line Station Data and Adding it to Relevant .txt File¶

4. Creating Stop Time Data and Adding it to Relevant .txt File¶

5. Creating Ontario Line Trip Data and Adding it to Relevant .txt File¶

6. Creating Ontario Line Route Details and Adding it to Relevant .txt File¶

7. Zipping and Validating the Updated GTFS File¶

Part II: Measuring Changes in Access¶

1. Changes in Transportation Access (Walking to TTC Subway Stops)¶

2. Changes in Transportation Access (Cycling to TTC Subway Stops)¶

4. Changes to Clinic Access¶

5. Changes to Hospital Access¶

6. Changes to Grocery Store Access¶

7. Changes to Place of Worship Access¶

8. Changes to Library Access¶

10. Changes to Public Wifi Access¶

	stop_id	stop_code	stop_name	stop_lat	stop_lon	wheelchair_boarding
342	101	N.A.	Exhibition - Ontario Line	43.636386	-79.418215	1
344	102	N.A.	King / Bathurst	43.643816	-79.402499	1
345	103	N.A.	Queen / Spadina	43.648679	-79.396341	1
338	104	N.A.	Osgoode - Ontario Line	43.650764	-79.386576	1
339	105	N.A.	Queen - Ontario Line	43.652343	-79.379220	1
346	106	N.A.	Moss Park	43.654276	-79.370209	1
347	107	N.A.	Corktown	43.652149	-79.364385	1
343	108	N.A.	East Harbour	43.655192	-79.347359	1
348	109	N.A.	Riverside/Leslieville	43.659850	-79.345699	1
349	110	N.A.	Gerrard	43.667428	-79.342489	1
340	111	N.A.	Pape - Ontario Line	43.679727	-79.345215	1
350	112	N.A.	Cosburn	43.689400	-79.348952	1
351	113	N.A.	Thorncliffe Park	43.705226	-79.350054	1
352	114	N.A.	Flemingdon Park	43.715024	-79.337171	1
341	115	N.A.	Science Centre - Ontario Line	43.720569	-79.338662	1

	route_id	service_id	trip_id	trip_headsign	direction_id	shape_id	wheelchair_accessible	bikes_allowed
0	3000	1	50000	ONTARIO LINE towards SCIENCE CENTRE	0	10000	1	1
1	3000	1	50001	ONTARIO LINE towards SCIENCE CENTRE	0	10000	1	1
2	3000	1	50002	ONTARIO LINE towards SCIENCE CENTRE	0	10000	1	1
3	3000	1	50003	ONTARIO LINE towards SCIENCE CENTRE	0	10000	1	1
4	3000	1	50004	ONTARIO LINE towards SCIENCE CENTRE	0	10000	1	1
...	...	...	...	...	...	...	...	...
1195	3000	1	51195	ONTARIO LINE towards EXHIBITION	1	10000	1	1
1196	3000	1	51196	ONTARIO LINE towards EXHIBITION	1	10000	1	1
1197	3000	1	51197	ONTARIO LINE towards EXHIBITION	1	10000	1	1
1198	3000	1	51198	ONTARIO LINE towards EXHIBITION	1	10000	1	1
1199	3000	1	51199	ONTARIO LINE towards EXHIBITION	1	10000	1	1

	route_id	service_id	trip_id	trip_headsign	trip_short_name	direction_id	block_id	shape_id	wheelchair_accessible	bikes_allowed
0	72932	1	47907138	EAST - 10 VAN HORNE towards VICTORIA PARK	NaN	0	2085726.0	1035571	1	1
1	72932	1	47907139	EAST - 10 VAN HORNE towards VICTORIA PARK	NaN	0	2085726.0	1035570	1	1
2	72932	1	47907145	EAST - 10 VAN HORNE towards VICTORIA PARK	NaN	0	2085726.0	1035570	1	1
3	72932	1	47907144	EAST - 10 VAN HORNE towards VICTORIA PARK	NaN	0	2085726.0	1035570	1	1
4	72932	1	47907143	EAST - 10 VAN HORNE towards VICTORIA PARK	NaN	0	2085726.0	1035570	1	1
...	...	...	...	...	...	...	...	...	...	...
1195	3000	1	51195	ONTARIO LINE towards EXHIBITION	NaN	1	NaN	10000	1	1
1196	3000	1	51196	ONTARIO LINE towards EXHIBITION	NaN	1	NaN	10000	1	1
1197	3000	1	51197	ONTARIO LINE towards EXHIBITION	NaN	1	NaN	10000	1	1
1198	3000	1	51198	ONTARIO LINE towards EXHIBITION	NaN	1	NaN	10000	1	1
1199	3000	1	51199	ONTARIO LINE towards EXHIBITION	NaN	1	NaN	10000	1	1

	route_id	agency_id	route_short_name	route_long_name	route_desc	route_type	route_url	route_color	route_text_color
0	73152	1	1	LINE 1 (YONGE-UNIVERSITY)	NaN	1	NaN	D5C82B	000000
1	72932	1	10	VAN HORNE	NaN	3	NaN	FF0000	FFFFFF
2	72933	1	100	FLEMINGDON PARK	NaN	3	NaN	FF0000	FFFFFF
3	72934	1	101	DOWNSVIEW PARK	NaN	3	NaN	FF0000	FFFFFF
4	72935	1	102	MARKHAM RD.	NaN	3	NaN	FF0000	FFFFFF
...	...	...	...	...	...	...	...	...	...
212	73143	1	989	WESTON EXPRESS	NaN	3	NaN	008000	FFFFFF
213	73144	1	99	ARROW ROAD	NaN	3	NaN	FF0000	FFFFFF
214	73145	1	995	YORK MILLS EXPRESS	NaN	3	NaN	008000	FFFFFF
215	73146	1	996	WILSON EXPRESS	NaN	3	NaN	008000	FFFFFF
0	3000	1	3	ONTARIO LINE	NaN	1	NaN	008000	FFFFFF

	indicator	value
0	agencies	[TTC]
1	timezone	America/Toronto
2	start_date	20241117
3	end_date	20241221
4	num_routes	216
5	num_trips	99178
6	num_stops	9373
7	num_shapes	1469
8	sample_date	20241121
9	num_routes_active_on_sample_date	213
10	num_trips_active_on_sample_date	37763
11	num_stops_active_on_sample_date	9325

Kasandra Tworzyanski, Muhammad Khalis, William Chan

GGR375 Introduction to Programming in GIS - Final Project

Due: December 11, 2024

Research Question: How would the connectivity for residents of Moss Park, Toronto, change as a result of the Ontario Line?

Sub-Questions¶

Sub-Question 3: How will this change if Ontario Line is open for service today?¶

Part I: Cleaning Ontario Line Data to Make New GTFS¶

1. Extracting Files from Existing TTC GTFS Zipped Folder¶

2. Cleaning Ontario Line Shape Data and Adding it to Relevant .txt File¶

3. Cleaning Ontario Line Station Data and Adding it to Relevant .txt File¶

4. Creating Stop Time Data and Adding it to Relevant .txt File¶

5. Creating Ontario Line Trip Data and Adding it to Relevant .txt File¶

6. Creating Ontario Line Route Details and Adding it to Relevant .txt File¶

7. Zipping and Validating the Updated GTFS File¶

Part II: Measuring Changes in Access¶

1. Changes in Transportation Access (Walking to TTC Subway Stops)¶

2. Changes in Transportation Access (Cycling to TTC Subway Stops)¶

3. Changes in Multi-Modal Accessibility to the City¶

4. Changes to Clinic Access¶

5. Changes to Hospital Access¶

6. Changes to Grocery Store Access¶

7. Changes to Place of Worship Access¶

8. Changes to Library Access¶

9. Changes to Social Facility Access¶

10. Changes to Public Wifi Access¶