Evaluating Public Transit Effectiveness with GIS

import arcpy
import os

# Function to convert GTFS shapes to features and calculate geodesic length
def process_gtfs_shapes(gtfs_folder, feature_dataset):
	try:
    	# Convert GTFS shapes to features
    	gtfs_shapes_file = os.path.join(gtfs_folder, 'shapes.txt')
    	city_name = os.path.basename(gtfs_folder)
    	out_feature_class_shapes = f'shapes_{city_name}'

    	# Check if the output feature class already exists, delete it if it does
    	if arcpy.Exists(out_feature_class_shapes):
        	arcpy.Delete_management(out_feature_class_shapes)

    	# Convert GTFS shapes to features
    	arcpy.transit.GTFSShapesToFeatures(in_gtfs_shapes_file=gtfs_shapes_file, out_feature_class=out_feature_class_shapes)

    	print(f'Conversion of GTFS shapes to features completed for {city_name}. Output: {out_feature_class_shapes}')

    	# Set the input feature class for geodesic length calculation (for shapes)
    	input_fc_shapes = out_feature_class_shapes

    	# Add a new field called LENGTH in miles (for shapes)
    	field_name_length_shapes = 'LENGTH'
    	arcpy.AddField_management(input_fc_shapes, field_name_length_shapes, 'DOUBLE')

    	# Calculate geodesic length in statute miles (for shapes)
    	expression_length_shapes = '!shape.geodesicLength@MILES!'
    	arcpy.CalculateField_management(input_fc_shapes, field_name_length_shapes, expression_length_shapes, 'PYTHON3')

    	print(f'Geodesic length calculation completed for {city_name}. Field added: {field_name_length_shapes}')

    	# Define projection for the feature class to NAD83 (EPSG:4269)
    	sr = arcpy.SpatialReference(4269)  # NAD83
    	arcpy.management.DefineProjection(out_feature_class_shapes, sr)

	except arcpy.ExecuteError:
    	print(arcpy.GetMessages(2))
	except Exception as e:
    	print(f'Failed to process GTFS shapes for {city_name} with error: {str(e)}')


# Function to convert GTFS stops to features
def process_gtfs_stops(gtfs_folder, feature_dataset):
	try:
    	# Convert GTFS stops to features
    	gtfs_stops_file = os.path.join(gtfs_folder, 'stops.txt')
    	city_name = os.path.basename(gtfs_folder)
    	out_feature_class_stops = f'stops_{city_name}'

    	# Check if the output feature class already exists, delete it if it does
    	if arcpy.Exists(out_feature_class_stops):
        	arcpy.Delete_management(out_feature_class_stops)

    	# Convert GTFS stops to features
    	arcpy.transit.GTFSStopsToFeatures(in_gtfs_stops_file=gtfs_stops_file, out_feature_class=out_feature_class_stops)

    	print(f'Conversion of GTFS stops to features completed for {city_name}. Output: {out_feature_class_stops}')

    	# Define projection for the feature class to NAD83 (EPSG:4269)
    	sr = arcpy.SpatialReference(4269)  # NAD83
    	arcpy.management.DefineProjection(out_feature_class_stops, sr)

	except arcpy.ExecuteError:
    	print(arcpy.GetMessages(2))
	except Exception as e:
    	print(f'Failed to process GTFS stops for {city_name} with error: {str(e)}')


# Function to run GTFSToPublicTransitDataModel
def run_gtfs_to_transit_data_model(gtfs_folder, target_feature_dataset):
	try:
    	arcpy.transit.GTFSToPublicTransitDataModel(
        	in_gtfs_folders=gtfs_folder,
        	target_feature_dataset=target_feature_dataset,
        	interpolate="INTERPOLATE",
        	append="NO_APPEND"
    	)
    	print(f'GTFSToPublicTransitDataModel completed for {os.path.basename(gtfs_folder)}.')
	except arcpy.ExecuteError:
    	# Print detailed error messages
    	print(arcpy.GetMessages(2))
	except Exception as e:
    	print(f'Failed to process {os.path.basename(gtfs_folder)} with error: {str(e)}')


# Main function to process GTFS data for each city
def process_gtfs_data(base_workspace, gtfs_directory):
	# Iterate over each city folder in the directory
	for city_folder in os.listdir(gtfs_directory):
    	gtfs_folder = os.path.join(gtfs_directory, city_folder)

    	# Check if it's a directory and contains necessary GTFS files (shapes.txt and stops.txt)
    	if os.path.isdir(gtfs_folder):
        	if 'shapes.txt' in os.listdir(gtfs_folder) and 'stops.txt' in os.listdir(gtfs_folder):
            	# Create a new geodatabase for the city if it doesn't exist
            	gdb_path = os.path.join(base_workspace, f'{city_folder}.gdb')
            	if not arcpy.Exists(gdb_path):
                	arcpy.CreateFileGDB_management(base_workspace, f'{city_folder}.gdb')
                	print(f'Created geodatabase: {gdb_path}')

            	# Set the workspace to the new geodatabase
            	arcpy.env.workspace = gdb_path

            	# Create a feature dataset for GTFS data if it doesn't exist
            	feature_dataset_name = 'GTFS'
            	feature_dataset_path = os.path.join(gdb_path, feature_dataset_name)
            	if not arcpy.Exists(feature_dataset_path):
                	# Create feature dataset with NAD83 (EPSG:4269) spatial reference
                	sr = arcpy.SpatialReference(4269)  # NAD83
                	arcpy.CreateFeatureDataset_management(gdb_path, feature_dataset_name, sr)
                	print(f'Created feature dataset: {feature_dataset_path}')

            	# Process GTFS shapes and stops for the city
            	process_gtfs_shapes(gtfs_folder, feature_dataset_path)
            	process_gtfs_stops(gtfs_folder, feature_dataset_path)

            	# Run GTFSToPublicTransitDataModel for the city
            	run_gtfs_to_transit_data_model(gtfs_folder, feature_dataset_path)

        	else:
            	print(f'Skipping {city_folder} - Missing shapes.txt or stops.txt.')
    	else:
        	print(f'Skipping {city_folder} - Not a valid directory.')

	print('GTFS data processing completed.')

# Execute the main processing function
if __name__ == "__main__":
	# Set environment settings
	base_workspace = r'M:\Documents\ArcGIS\Projects\ALL'  # Base workspace where all geodatabases will reside

	# Path to the directory containing GTFS folders for each city
	gtfs_directory = r'M:\Documents\ArcGIS\Projects\ALL\GTFS'

	# Process GTFS data
	process_gtfs_data(base_workspace, gtfs_directory)

import requests
import zipfile
import io

#Transitland API key
transitland_api_key = ""

# Feed keys for WMATA rail and bus GTFS
feed_keys = {
    "rail": "f-dqc-wmata~rail",
    "bus": "f-dqc-wmata~bus"
}

# Base URL for downloading the latest feed version
base_url = "https://transit.land/api/v2/rest/feeds/{feed_key}/download_latest_feed_version"

# Function to download and extract GTFS feed
def download_and_extract_gtfs(mode, feed_key):
    url = base_url.format(feed_key=feed_key)
    headers = {"apikey": transitland_api_key}

    # Make the request
    response = requests.get(url, headers=headers)

    # Check if the request was successful
    if response.status_code == 200:
        # Extract the contents of the ZIP file
        z = zipfile.ZipFile(io.BytesIO(response.content))
        extract_path = f"path_to_extract_folder/{mode}"  # Replace with your desired extraction path
        z.extractall(extract_path)
        print(f"{mode.capitalize()} GTFS feed downloaded and extracted successfully to {extract_path}.")
    else:
        print(f"Failed to download {mode} GTFS feed. Status code: {response.status_code}")
        try:
            error_message = response.json().get('error')
            print(f"Error message: {error_message}")
        except Exception as e:
            print(f"Failed to parse error response: {str(e)}")

# Download each GTFS feed
for mode, feed_key in feed_keys.items():
    download_and_extract_gtfs(mode, feed_key)
Rail GTFS feed downloaded and extracted successfully to path_to_extract_folder/rail.
Bus GTFS feed downloaded and extracted successfully to path_to_extract_folder/bus.

# Station data with lines serviced
station_data = {
	'Addison Road-Seat Pleasant': ['Silver Line', 'Blue Line'],
	'Anacostia': ['Green Line'],
	'Archives-Navy Memorial-Penn Quarter': ['Green Line', 'Yellow Line'],
	'Arlington Cemetery': ['Blue Line'],
	'Ashburn': ['Silver Line'],
	'Ballston-MU': ['Orange Line', 'Silver Line'],
	'Benning Road': ['Silver Line', 'Blue Line'],
	'Bethesda': ['Red Line'],
	'Braddock Road': ['Blue Line', 'Yellow Line'],
	'Branch Ave': ['Green Line'],
	'Brookland-CUA': ['Red Line'],
	'Capitol Heights': ['Silver Line', 'Blue Line'],
	'Capitol South': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Cheverly': ['Orange Line'],
	'Clarendon': ['Orange Line', 'Silver Line'],
	'Cleveland Park': ['Red Line'],
	'College Park-U of Md': ['Green Line'],
	'Columbia Heights': ['Green Line'],
	'Congress Heights': ['Green Line'],
	'Court House': ['Orange Line', 'Silver Line'],
	'Crystal City': ['Blue Line', 'Yellow Line'],
	'Deanwood': ['Orange Line'],
	'Downtown Largo': ['Silver Line', 'Blue Line'],
	'Dunn Loring-Merrifield': ['Orange Line'],
	'Dupont Circle': ['Red Line'],
	'East Falls Church': ['Orange Line', 'Silver Line'],
	'Eastern Market': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Eisenhower Avenue': ['Yellow Line'],
	'Farragut North': ['Red Line'],
	'Farragut West': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Federal Center SW': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Federal Triangle': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Foggy Bottom-GWU': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Forest Glen': ['Red Line'],
	'Fort Totten': ['Green Line', 'Red Line'],
	'Franconia-Springfield': ['Blue Line'],
	'Friendship Heights': ['Red Line'],
	'Gallery Pl-Chinatown': ['Green Line', 'Yellow Line', 'Red Line'],
	'Georgia Ave-Petworth': ['Green Line'],
	'Glenmont': ['Red Line'],
	'Greenbelt': ['Green Line'],
	'Greensboro': ['Silver Line'],
	'Grosvenor-Strathmore': ['Red Line'],
	'Herndon': ['Silver Line'],
	'Huntington': ['Yellow Line'],
	'Hyattsville Crossing': ['Green Line'],
	'Innovation Center': ['Silver Line'],
	'Judiciary Square': ['Red Line'],
	'King St-Old Town': ['Blue Line', 'Yellow Line'],
	'L\'Enfant Plaza': ['Orange Line', 'Silver Line', 'Blue Line', 'Yellow Line', 'Green Line'],
	'Landover': ['Orange Line'],
	'Loudoun Gateway': ['Silver Line'],
	'McLean': ['Silver Line'],
	'McPherson Square': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Medical Center': ['Red Line'],
	'Metro Center': ['Red Line', 'Orange Line', 'Silver Line', 'Blue Line'],
	'Minnesota Ave': ['Orange Line'],
	'Morgan Boulevard': ['Silver Line', 'Blue Line'],
	'Mt Vernon Sq 7th St-Convention Center': ['Green Line', 'Yellow Line'],
	'Navy Yard-Ballpark': ['Green Line'],
	'Naylor Road': ['Green Line'],
	'New Carrollton': ['Orange Line'],
	'NoMa-Gallaudet U': ['Red Line'],
	'North Bethesda': ['Red Line'],
	'Pentagon': ['Blue Line', 'Yellow Line'],
	'Pentagon City': ['Blue Line', 'Yellow Line'],
	'Potomac Ave': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Potomac Yard': ['Blue Line', 'Yellow Line'],
	'Reston Town Center': ['Silver Line'],
	'Rhode Island Ave-Brentwood': ['Red Line'],
	'Rockville': ['Red Line'],
	'Ronald Reagan Washington National Airport': ['Blue Line', 'Yellow Line'],
	'Rosslyn': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Shady Grove': ['Red Line'],
	'Shaw-Howard U': ['Green Line'],
	'Silver Spring': ['Red Line'],
	'Smithsonian': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Southern Avenue': ['Green Line'],
	'Spring Hill': ['Silver Line'],
	'Stadium-Armory': ['Orange Line', 'Silver Line', 'Blue Line'],
	'Suitland': ['Green Line'],
	'Takoma': ['Red Line'],
	'Tenleytown-AU': ['Red Line'],
	'Twinbrook': ['Red Line'],
	'Tysons': ['Silver Line'],
	'U Street/African-Amer Civil War Memorial/Cardozo': ['Green Line'],
	'Union Station': ['Red Line'],
	'Van Dorn Street': ['Blue Line'],
	'Van Ness-UDC': ['Red Line'],
	'Vienna/Fairfax-GMU': ['Orange Line'],
	'Virginia Square-GMU': ['Orange Line', 'Silver Line'],
	'Washington Dulles International Airport': ['Silver Line'],
	'Waterfront': ['Green Line'],
	'West Falls Church': ['Orange Line'],
	'West Hyattsville': ['Green Line'],
	'Wheaton': ['Red Line'],
	'Wiehle-Reston East': ['Silver Line'],
	'Woodley Park-Zoo/Adams Morgan': ['Red Line']
}

# Dictionary to count lines per station
lines_per_station = {}

# Count lines per station
for station, lines in station_data.items():
	lines_per_station[station] = len(lines)

# Calculate average number of lines per station
total_stations = len(station_data)
total_lines = sum(lines_per_station.values())
average_lines_per_station = total_lines / total_stations

# Print results
print(f"Average number of lines per station: {average_lines_per_station:.2f}")


—-------------

# Average lines per station and total stations
average_lines_per_station = 1.52
total_stations = 91

# Define operational hours for Tuesday in minutes
operating_hours_tuesday = {
	'Tuesday': 19 * 60  # 5:00 AM to 12:00 AM
}

# Define peak service frequencies in minutes
peak_frequencies = {
	'Red Line': 8,
	'Orange Line': 12,
	'Silver Line': 12,
	'Blue Line': 12,
	'Yellow Line': 12,
	'Green Line': 8
}

peak_frequency_minutes = sum(peak_frequencies.values()) / len(peak_frequencies)

# Print the average peak frequency
print(f"Average peak frequency: {peak_frequency_minutes:.2f} minutes")


# Calculate total trains per day for all stations combined
trains_per_hour = 60 // peak_frequency_minutes
total_trains_per_day = trains_per_hour * 19  # 19 hours operating on tuesday

# Calculate total lines
total_lines = average_lines_per_station * total_stations

# Calculate total trains based on lines and frequency
total_trains_combined = total_lines * total_trains_per_day

# Print the result
print(f"Total trains coming through all DC stations combined per day: {total_trains_combined:.2f}")

import jenkspy

# Given LITA values
lita_values = [
	8.5, 7.5, 11.17, 14.17, 15.17,
	10.17, 9.83, 12.17, 11.17, 10.17
]

# Step 1: Sort data (already sorted for demonstration)
sorted_lita = sorted(lita_values)

# Step 2: Determine number of classes
n_classes = min(len(np.unique(sorted_lita)), 10)  # Limit to 10 classes maximum

# Step 3: Calculate Jenks natural breaks
breaks = jenkspy.jenks_breaks(sorted_lita, n_classes)

# Step 4: Assign classes based on breaks
def assign_class(value, breaks):
	for i in range(1, len(breaks)):
    	if value <= breaks[i]:
        	return i
	return len(breaks) - 1

# Assign classes to each LITA value
classes = [assign_class(value, breaks) for value in lita_values]

# Print results
print("LITA Values\tClass")
for value, cls in zip(lita_values, classes):
	print(f"{value}\t\t{cls}")