M.Sc. Statistics · Vadodara Municipal Corporation · 2024–2026

Urban
Inequality in Vadodara

A thesis-grade geo-spatial analysis combining accessibility modelling, Bayesian spatial statistics, hydrodynamic flood simulation, and graph neural networks to measure and map opportunity across 19 urban wards.

19Urban Wards

0.73Moran's I

4hSWMM Simulation

BYM2Spatial Model

GCNGraph Neural Net

Explore the Network → View Repository

Scroll to explore

Context

What is this project?

Vadodara (Baroda), Gujarat's third-largest city, exhibits sharp spatial disparities in access to healthcare, education, transit, and flood safety. This project builds a rigorous, reproducible measurement framework.

Urban inequality is rarely a single-dimension problem. A ward may have excellent schools but be chronically flood-prone. Another may have strong transit coverage but poor healthcare access. The Urban Opportunity Index (UOI) synthesises these four dimensions into a single, interpretable ward-level score using normalised, weighted components.

The spatial distribution of the UOI is not random — it exhibits strong spatial autocorrelation (Global Moran's I = 0.73), meaning high-opportunity and low-opportunity wards cluster together. This clustering motivates the use of Bayesian spatial models (BYM2) rather than standard regression.

The flood risk component is computed dynamically using SWMM DYNWAVE hydrodynamic simulation — not a static raster overlay — making it physically grounded and sensitive to road-network topology.

🗺️

Open Data

OSMnx road networks, Open Vadodara Smart City Bus GTFS, SRTM 30m DEM, VMC ward boundaries.

📐

Rigorous Statistics

Bayesian BYM2 with PC priors, Moran's I, LISA clusters, Kolm–Pollak EDE inequality adjustment.

🌊

Physics-Based Flood

SWMM St Venant DYNWAVE engine simulates 4-hour cloudburst. Road failures computed per ward.

🧠

Graph Neural Networks

3-layer GCN extracts spatially-contextualised ward embeddings, fed into BYM2 as covariates.

Urban Opportunity Index

The UOI Formula

Four normalised components, weighted by policy relevance, produce a 0–100 score where higher means better access and lower risk.

// formula_v2.4 · balanced_weights

UOI_i = 0.30 · Ĥᵢ + 0.30 · Êᵢ + 0.20 · T̂ᵢ + 0.20 · (1 − F̂ᵢ)

All components min-max normalised to [0,1] across all 19 wards. F̂ᵢ is inverted so high flood exposure reduces the score. Final scale is 0–100. Kolm–Pollak EDE adjustment with κ ∈ [0.1, 2.0] validates inequality-adjusted rankings.

30%

Healthcare (Ĥ)

Drive-time isochrone to nearest hospital, clinic, or health post via OSMnx road network.

OSMnx · ISOCHRONE

30%

Education (Ê)

Walk-time to nearest school, college, or university. Captures local educational access.

OSMnx · WALK

20%

Transit (T̂)

Bus stop density + route coverage from Open Vadodara Smart City Bus GTFS data.

GTFS · city2graph

20%

Flood Risk (1−F̂)

Inverted SWMM hydraulic flood road failure percentage. Dynamic, not a static raster.

SWMM · DYNWAVE

Ward UOI Scores — All 19 Wards

Methodology

Five analytical stages

Each stage is reproducible and runs in order. The pipeline processes open data into the final Bayesian posterior.

01 // DATA COLLECTION

Spatial data ingestion

Road networks downloaded via OSMnx. Bus stops from Open Vadodara GTFS. SRTM 30m DEM for elevations. VMC ward boundaries as the analytical unit.

OSMnxGTFSSRTMGeoPandas

02 // ACCESSIBILITY

Travel-time isochrones

Network-based drive and walk times from each ward centroid to nearest facilities. Transit stop density computed from GTFS via city2graph.travel_summary_graph().

city2graphisochroneGTFS

03 // COMPOSITE INDEX

UOI construction

Min-max normalise all four components. Apply weights. Compute Gini coefficient and Kolm–Pollak EDE with κ-sweep from 0.1 to 2.0 for inequality sensitivity.

NormalisationKolm-PollakGini

04 // FLOOD SIMULATION

SWMM DYNWAVE

Major road network (~700 nodes) filtered from OSMnx. 4-hour triangular storm hyetograph. St Venant DYNWAVE engine. Ward-level flood road percentage aggregated.

SWMM 5.2DYNWAVECFL dt=5s

05 // BAYESIAN INFERENCE

BYM2 + GCN

BYM2 with ICAR prior via PyMC 5 NUTS. city2graph Queen contiguity adjacency. Optional GCN penultimate-layer embeddings as covariates. 4 chains × 2000 draws.

PyMC 5BYM2ICARGCN

06 // SENSITIVITY

Robustness checks

OAT (±5pp), Monte Carlo Dirichlet N=10,000, Sobol' first-order indices, Kolm–Pollak κ-sweep, spatial CI mapping. 9/19 wards achieve full classification robustness.

OATMonte CarloSobol'

Spatial Statistics

Measuring inequality patterns

Global and local spatial autocorrelation metrics reveal how opportunity clusters — and where inequality concentrates.

Global Moran's I

Measures spatial autocorrelation of UOI across all ward pairs under the Queen contiguity weight matrix (city2graph).

I = 0.73

Global Moran's I · Queen contiguity · p < 0.001

Strong positive spatial autocorrelation. High-opportunity wards cluster with high-opportunity neighbours. Justifies BYM2 spatial random effect over standard OLS.

HH: Opportunity clusters LL: Deprivation clusters HL: High outlier LH: Low outlier

Inequality Metrics

Gini coefficient and Kolm–Pollak EDE measure the distribution of opportunity across wards.

Gini Coefficient0.218

0 (perfect equality)1 (maximum inequality)

Kolm–Pollak EDE (κ-sweep)

κ = 0.152.8

κ = 0.551.4

κ = 1.050.6

κ = 2.049.8

Spearman ρ ≥ 0.99 across full κ range → rankings are robust to inequality aversion parameter choice.

Bayesian Spatial Model

BYM2 — Besag York Mollié 2

The primary inferential model. Decomposes spatial variation into structured clustering (ICAR) and unstructured noise, with interpretable Penalised Complexity hyperpriors.

// likelihood

yᵢ | μᵢ, σε ~ Normal(μᵢ, σε²)

// linear predictor

μᵢ = α + Σⱼ βⱼ xᵢⱼ + bᵢ

// BYM2 spatial effect

bᵢ = σb · (√φ · u*ᵢ + √(1−φ) · v*ᵢ)

// ICAR structured effect

u* ~ ICAR(W) [scaled, sum-to-zero]

// IID unstructured noise

v*ᵢ ~ Normal(0, 1) [IID]

// PC priors (Riebler et al. 2016)

σb ~ HalfNormal(0.5) P(σb > 1.0) = 0.01
φ ~ Beta(0.5, 0.5) [neutral]

⚠ Critical: do NOT use pm.CAR

pm.CAR(alpha=0.99, tau=1.0) approximates ICAR but hardcodes precision and produces biased σb and φ estimates. Use the explicit ICAR log-potential: pm.Potential("icar", icar_logp(u_raw, node1, node2)) with node arrays from city2graph Queen contiguity.

Parameter	Prior	Interpretation
α	Normal(0,1)	Baseline UOI (weakly informative)
βⱼ	Normal(0,1)	Component j coefficient
σb	HalfNormal(0.5) [PC]	Total spatial SD
φ	Beta(0.5,0.5) [PC]	Spatial fraction of σb
σε	HalfNormal(1.0)	Residual noise SD

MCMC Configuration

4 chains × 2,000 draws
tune = 1,000 steps
target_accept = 0.90
adjacency: city2graph Queen

Interpretation of φ

φ → 1: spatial variation is entirely structured (neighbouring wards are similar — ICAR dominates).

φ → 0: variation is ward-specific noise with no spatial pattern.

Expected range for UOI: φ ∈ [0.55, 0.85], confirming strong neighbourhood effects.

Hydrodynamic Simulation

SWMM — Why it hung for 17 minutes

The original build_swm.py caused a 17-minute execution freeze due to four compounding failures — all diagnosed empirically and fixed.

Network scale

✕ network_type='drive' → ~18,000 nodes

✓ custom_filter='primary|secondary|tertiary' → ~700 nodes (20× reduction)

Catastrophic constant inflow

✕ baseline=0.05 CMS × 18k nodes = 900 CMS

✓ TIMESERIES storm hyetograph, peak 0.005 CMS/node (100× less)

Random elevations

✕ np.random.uniform(30,40) → 50% uphill conduits, closed loops

✓ Synthetic DEM slope toward Vishwamitri river (real topography)

Zero junction storage

✕ MaxDepth=0 → instant surcharge at every node

✓ MaxDepth=2.0m, Aponded=200m², 5 boundary outfalls

Result: 17+ minutes → <30 seconds

Storm Timeline

0h30m1h1h30 ⚡2h2h303h4h

CFL Condition

Δt ≤ Δx / c
c = √(gA/T) ≈ 3.1 m/s
Δx_min ≈ 50m
→ max stable dt = 16s
→ routing_step = 5s ✓

Storm Design

Vadodara 10-yr return: ~50mm/hr peak
Catchment per node: ~500m²
Runoff coefficient C = 0.70
Q_peak = 50/3600 × 500 × 0.70 = 0.005 CMS

Reproduce this

Clean pipeline — 11 steps

Run these scripts in order. Each step feeds the next.

03_accessibility_engine.py

Accessibility modelling

OSMnx network download → drive/walk isochrones to hospitals, schools.

03b_transit_accessibility.py

Transit graph

city2graph.load_gtfs() + travel_summary_graph(). Requires vadodara_gtfs.zip.

04_risk_assessment.py

Flood risk raster

SRTM 30m DEM → flood exposure % per ward. Place DEM at data/raw/.

05_calculate_uoi.py

UOI composite

Min-max normalise, apply weights, compute Gini + EDE. Output: vadodara_final_uoi_balanced.gpkg

05b_prepare_bym2_inputs.py

BYM2 data prep

Standardise features, build adjacency arrays. Output: bym2_inputs.npz

10_city2graph_integration.py

city2graph graphs

Queen/Rook contiguity, KNN, Delaunay, POI group_nodes. Exports bym2_node1/2.npy.

06_bym2_pymc_model.py

BYM2 model

ICAR Potential, PC priors, NUTS 4 chains × 2000. Output: wards_bym2_results.gpkg + bym2_idata.nc

build_swm_fixed.py

SWMM network builder

Major roads only, storm timeseries, DEM slope, 5 outfalls. Output: vadodara_overland.inp

run_swmm_fixed.py

SWMM simulation

DYNWAVE 4-hour run, RPT parser, ward aggregation. Adds swmm_flood_road_pct to .gpkg

GNN.py (fix 5 bugs first)

GCN + BYM2 hybrid

3-layer GCN embeddings → BYM2 covariate. Fix: replace pm.CAR with ICAR Potential; use LOO-CV predictions.

streamlit run dashboard.py

Streamlit dashboard

5-KPI dark-theme dashboard at localhost:8501. Reads gnn_bym2_results layer from .gpkg.

Theory

Theoretical foundations

Key citations and theoretical extracts supporting each analytical choice.

ANSELIN, L. (1988)

Spatial Econometrics

OLS fails when observations are spatially autocorrelated — breaking the independence assumption inflates Type I error rates and produces biased coefficients. The LM tests (LM-Lag, LM-Error, RLM versions) provide a data-driven procedure to select between SLM and SEM, which is the gateway to BYM2.

RIEBLER ET AL. (2016)

BYM2 Penalised Complexity

The BYM2 reparameterisation makes σb and φ directly interpretable. PC priors penalise deviation from the base model (constant risk), placing probability mass on simpler models unless the data demands complexity. The mixing parameter φ ∈ [0,1] quantifies the proportion of spatial variation attributable to neighbourhood structure.

P(σb > U) = a where U=0.5/0.31, a=0.01

ANSELIN (1995)

LISA Cluster Maps

Local Indicators of Spatial Association decompose global Moran's I into four typologies: High-High (opportunity clusters — policy reinforcement), Low-Low (deprivation clusters — priority intervention), High-Low (positive outlier), Low-High (isolated deprived ward). The LISA map is the primary policy communication output.

KOLM (1976), POLLAK (1971)

Kolm–Pollak EDE

The Equally Distributed Equivalent measures what constant UOI would deliver the same social welfare as the observed distribution. Unlike Atkinson's measure, it uses a negative exponential social welfare function making it unit-specific. The κ-sweep from 0.1 to 2.0 tests sensitivity to the degree of inequality aversion.

EDE_κ = −(1/κ)·log[(1/n)Σᵢ exp(−κ·xᵢ)]

EPA / SWMM 5.2

DYNWAVE St Venant

DYNWAVE solves the full 1D Saint-Venant equations: continuity ∂A/∂t + ∂Q/∂x = q and momentum ∂Q/∂t + ∂(Q²/A)/∂x + gA(∂H/∂x + Sf) = 0. The CFL stability condition Δt ≤ Δx/c must be satisfied at all timesteps. Violation causes DYNWAVE to internally substep to milliseconds — the root cause of the 17-minute hang.

SATO (2025) / BOEING (2025)

city2graph + OSMnx

city2graph converts geospatial datasets to graph representations integrated with GeoPandas, NetworkX and PyTorch Geometric. The contiguity_graph() function returns a MultiIndex edge GeoDataFrame where level_0/level_1 are 0-indexed ward positions — directly usable as BYM2 node1/node2 adjacency arrays without manual mapping.

Repository

All deliverables

Python

06_bym2_pymc_model.py

BYM2 ICAR model in PyMC 5. Correct ICAR via pm.Potential, PC priors, city2graph adjacency.

Complete

Python

10_city2graph_integration.py

Contiguity, proximity (KNN/Delaunay/Gabriel), POI group_nodes, transit graph, heterogeneous graph.

Complete

Python

uoi_sensitivity_analysis.py

OAT, Monte Carlo Dirichlet N=10k, Sobol' indices, Kolm–Pollak κ-sweep, spatial CI mapping.

Complete

Python

build_swm_fixed.py

Major roads only, storm TIMESERIES, DEM slope, 5 outfalls. Fixed all 4 root causes of 17-min hang.

Complete

Python

run_swmm_fixed.py

Regex edge ID parser, column-discovery RPT parser, road-type thresholds, ward aggregation.

Complete

Python

GNN.py

GCN + BYM2 hybrid. 5 bugs documented: replace pm.CAR, use LOO-CV, fix priors, target_accept, spatial CV.

Needs fixes

Python

dashboard.py

Streamlit dark-theme dashboard. 5 KPI cards, pydeck 3D map, layer toggles, correlation matrix.

Complete

HTML

vadodara_network.html

Self-contained 142KB interactive flood network. Pre/Post/Diff modes, storm timeline, pan+zoom.

Complete

SWMM .inp

vadodara_overland_fixed.inp

Patched .inp file. 5 missing junction nodes (OSMnx outfall candidates) added back. Passes swmm_open().

Complete

UrbanInequality in Vadodara

What is this project?

The UOI Formula

Five analytical stages

Road network & flood simulation

Measuring inequality patterns

BYM2 — Besag York Mollié 2

SWMM — Why it hung for 17 minutes

Clean pipeline — 11 steps

Theoretical foundations

All deliverables

Urban
Inequality in Vadodara