Python for Transportation & Urban Planning

import pandas as pd

# Reading data from various sources
df = pd.read_csv('traffic_counts.csv')
df = pd.read_excel('population_data.xlsx', sheet_name='Summary')
df = pd.read_sql('SELECT * FROM parcels', connection)

# Basic data exploration
df.head()          # First 5 rows
df.info()          # Data types and null values
df.describe()      # Statistical summary
df.shape           # Dimensions (rows, columns)

# Handle missing values (common in planning datasets)
df.dropna(subset=['population'])           # Remove rows with missing population
df.fillna({'income': df['income'].median()})  # Fill with median income

# Data type conversions
df['census_tract'] = df['census_tract'].astype(str)
df['survey_date'] = pd.to_datetime(df['survey_date'])

# Remove duplicates
df.drop_duplicates(subset=['parcel_id'])

# Group by transportation analysis zones
zone_summary = df.groupby('taz_id').agg({
    'population': 'sum',
    'employment': 'sum',
    'trip_generation': 'mean',
    'income': 'median'
}).round(2)

# Multiple grouping levels
mode_zone = df.groupby(['taz_id', 'mode_choice']).size().unstack(fill_value=0)

# Time-based aggregation
daily_counts = df.groupby(df['timestamp'].dt.date)['vehicle_count'].sum()

import geopandas as gpd
from shapely.geometry import Point, Polygon

# Read various spatial formats
parcels = gpd.read_file('parcels.shp')
transit_stops = gpd.read_file('transit_stops.geojson')
census_tracts = gpd.read_file('census_tracts.gpkg')

# Create GeoDataFrame from coordinates
transit_df = pd.DataFrame({
    'stop_id': [1, 2, 3],
    'ridership': [150, 200, 300],
    'longitude': [-122.4, -122.3, -122.2],
    'latitude': [37.8, 37.7, 37.6]
})

geometry = [Point(xy) for xy in zip(transit_df.longitude, transit_df.latitude)]
transit_gdf = gpd.GeoDataFrame(transit_df, geometry=geometry, crs='EPSG:4326')

# Buffer analysis - critical for accessibility studies
transit_buffer = transit_stops.buffer(400)  # 400-meter walking distance
parcels_near_transit = parcels[parcels.within(transit_buffer.unary_union)]

# Spatial joins
parcels_with_zones = gpd.sjoin(
    parcels, 
    zoning_districts, 
    how='left', 
    predicate='within'
)

# Overlay operations
intersection = gpd.overlay(
    flood_zones, 
    development_parcels, 
    how='intersection'
)

# Calculate areas and distances
parcels['area_sqft'] = parcels.geometry.area
parcels['perimeter'] = parcels.geometry.length

# Centroid calculations
parcels['centroid'] = parcels.geometry.centroid
parcels['centroid_x'] = parcels.centroid.x
parcels['centroid_y'] = parcels.centroid.y

# Distance calculations
parcels['dist_to_transit'] = parcels.geometry.distance(
    transit_stops.geometry.unary_union
)

import networkx as nx

# Create network from edge list
road_network = nx.from_pandas_edgelist(
    edges_df, 
    source='from_node', 
    target='to_node', 
    edge_attr=['length', 'speed_limit', 'capacity']
)

# Add node attributes
node_attrs = dict(zip(nodes_df.node_id, nodes_df.signal_type))
nx.set_node_attributes(road_network, node_attrs, 'signal_type')

# Single shortest path
shortest_path = nx.shortest_path(
    road_network, 
    source=origin_node, 
    target=destination_node, 
    weight='travel_time'
)

# All shortest paths from one origin
distances = nx.single_source_shortest_path_length(
    road_network, 
    origin_node, 
    weight='travel_time'
)

# Travel time matrix
travel_times = dict(nx.all_pairs_shortest_path_length(
    road_network, 
    weight='travel_time'
))

# Betweenness centrality - identifies critical corridors
betweenness = nx.betweenness_centrality(road_network, weight='length')

# Closeness centrality - measures accessibility
closeness = nx.closeness_centrality(road_network, distance='travel_time')

# Degree centrality - connectivity measure
degree = nx.degree_centrality(road_network)

# PageRank - importance in network
pagerank = nx.pagerank(road_network, weight='traffic_volume')

import matplotlib.pyplot as plt
import seaborn as sns

# Set planning-friendly style
plt.style.use('seaborn-v0_8')
sns.set_palette("husl")

# Mode share visualization
mode_counts = df['mode_choice'].value_counts()
plt.figure(figsize=(10, 6))
plt.pie(mode_counts.values, labels=mode_counts.index, autopct='%1.1f%%')
plt.title('Transportation Mode Share')
plt.show()

# Trip length distribution
plt.figure(figsize=(12, 6))
plt.hist(df['trip_distance'], bins=50, alpha=0.7, edgecolor='black')
plt.xlabel('Trip Distance (miles)')
plt.ylabel('Frequency')
plt.title('Distribution of Trip Distances')
plt.axvline(df['trip_distance'].mean(), color='red', linestyle='--', 
           label=f'Mean: {df["trip_distance"].mean():.1f} miles')
plt.legend()
plt.show()

import folium
from folium.plugins import HeatMap, MarkerCluster

# Create base map
m = folium.Map(
    location=[city_lat, city_lon], 
    zoom_start=12, 
    tiles='OpenStreetMap'
)

# Add transit stops with clustering
marker_cluster = MarkerCluster().add_to(m)
for idx, stop in transit_stops.iterrows():
    folium.Marker(
        location=[stop.geometry.y, stop.geometry.x],
        popup=f"Stop: {stop.stop_name}
Ridership: {stop.daily_boardings}",
        icon=folium.Icon(color='blue', icon='bus')
    ).add_to(marker_cluster)

# Add choropleth map for demographics
folium.Choropleth(
    geo_data=census_tracts,
    data=demographic_data,
    columns=['census_tract', 'median_income'],
    key_on='feature.properties.tract_id',
    fill_color='YlOrRd',
    fill_opacity=0.7,
    line_opacity=0.2,
    legend_name='Median Household Income'
).add_to(m)

m.save('planning_map.html')

from scipy import stats
import scipy.spatial.distance as distance

# Hypothesis testing for planning studies
# Compare before/after traffic volumes
before_traffic = df[df['period'] == 'before']['volume']
after_traffic = df[df['period'] == 'after']['volume']

# Paired t-test for before/after comparison
t_stat, p_value = stats.ttest_rel(before_traffic, after_traffic)
print(f"T-statistic: {t_stat:.3f}, p-value: {p_value:.3f}")

# Chi-square test for mode choice independence
contingency_table = pd.crosstab(df['income_group'], df['mode_choice'])
chi2, p_val, dof, expected = stats.chi2_contingency(contingency_table)

# Correlation analysis
correlation, p_value = stats.pearsonr(df['density'], df['transit_ridership'])
spearman_corr, p_val_spearman = stats.spearmanr(df['density'], df['transit_ridership'])

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import RandomForestRegressor
from sklearn.cluster import KMeans
from sklearn.metrics import mean_absolute_error, r2_score

# Predict transit ridership based on built environment
features = ['population_density', 'job_density', 'walk_score', 
           'parking_spaces', 'distance_to_cbd']
X = df[features]
y = df['transit_ridership']

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Linear regression model
lr_model = LinearRegression()
lr_model.fit(X_train, y_train)
y_pred = lr_model.predict(X_test)

print(f"R² Score: {r2_score(y_test, y_pred):.3f}")
print(f"Mean Absolute Error: {mean_absolute_error(y_test, y_pred):.0f}")

import requests
import json

# Census API for demographic data
census_key = 'your_census_api_key'
census_url = f"https://api.census.gov/data/2021/acs/acs5"
params = {
    'get': 'B01003_001E,B19013_001E,B25001_001E',  # Population, income, housing
    'for': 'tract:*',
    'in': 'state:06 county:075',  # San Francisco County
    'key': census_key
}

response = requests.get(census_url, params=params)
census_data = response.json()
census_df = pd.DataFrame(census_data[1:], columns=census_data[0])

# GTFS Real-time transit data
gtfs_url = "https://api.511.org/transit/vehiclepositions"
headers = {'api_key': 'your_511_api_key'}
transit_response = requests.get(gtfs_url, headers=headers)
transit_data = transit_response.json()

# Process JSON responses into DataFrames
def process_census_data(json_data):
    df = pd.DataFrame(json_data[1:], columns=json_data[0])
    df['population'] = pd.to_numeric(df['B01003_001E'], errors='coerce')
    df['median_income'] = pd.to_numeric(df['B19013_001E'], errors='coerce')
    df['housing_units'] = pd.to_numeric(df['B25001_001E'], errors='coerce')
    df['fips'] = df['state'] + df['county'] + df['tract']
    return df[['fips', 'population', 'median_income', 'housing_units']]

# Batch API requests with error handling
def batch_api_requests(urls, max_retries=3):
    results = []
    for url in urls:
        for attempt in range(max_retries):
            try:
                response = requests.get(url, timeout=10)
                response.raise_for_status()
                results.append(response.json())
                break
            except requests.exceptions.RequestException as e:
                if attempt == max_retries - 1:
                    print(f"Failed to fetch {url}: {e}")
                    results.append(None)
    return results

# Filter and transform data efficiently
high_density_zones = [zone for zone in zones if zone['density'] > 10000]
trip_distances = [calculate_distance(trip['origin'], trip['destination']) 
                 for trip in trip_data]

# Conditional operations
mode_categories = ['sustainable' if mode in ['walk', 'bike', 'transit'] 
                  else 'auto' for mode in df['mode_choice']]

# Nested comprehensions for complex operations
zone_accessibility = {
    zone: [stop for stop in transit_stops 
           if distance(zone['centroid'], stop['location']) < 400]
    for zone in analysis_zones
}

def calculate_gravity_model(origins, destinations, impedance_function):
    """
    Calculate gravity model for trip distribution
    """
    interactions = np.zeros((len(origins), len(destinations)))
    
    for i, origin in enumerate(origins):
        for j, destination in enumerate(destinations):
            distance = haversine_distance(origin['coords'], destination['coords'])
            impedance = impedance_function(distance)
            interactions[i, j] = (origin['production'] * 
                                destination['attraction'] * impedance)
    
    return interactions

def accessibility_index(zone_centroid, opportunities, decay_function):
    """
    Calculate accessibility index for a zone
    """
    accessibility = 0
    for opportunity in opportunities:
        distance = haversine_distance(zone_centroid, opportunity['location'])
        accessibility += opportunity['size'] * decay_function(distance)
    return accessibility

def mode_split_logit(utilities):
    """
    Calculate mode split using logit model
    """
    exp_utilities = np.exp(utilities)
    probabilities = exp_utilities / np.sum(exp_utilities)
    return probabilities