Spatial Epi Public Health GIS

Calculating the Two-Step Floating Catchment Area Index

The Two-Step Floating Catchment Area (2SFCA) index remains the operational standard for quantifying spatial healthcare accessibility in public health surveillance. Unlike static provider-to-population ratios, 2SFCA dynamically weights supply capacity against demand within realistic travel thresholds, applying distance-decay functions to reflect diminishing utilization probability. Deploying this metric at scale requires a deterministic, auditable pipeline that bridges network routing, capacity modeling, and spatial aggregation. Production systems must enforce strict coordinate reference system (CRS) alignment, implement compliant data minimization, and handle routing failures without silent data loss. This guide details a Python-native workflow optimized for agency-scale spatial epidemiology pipelines, aligning with established frameworks for Healthcare Access & Network Analysis Automation.

Phase 1: Compliant Data Ingestion & CRS Enforcement

Healthcare accessibility modeling begins with authoritative supply (clinics, hospitals, FQHCs) and demand (census tracts, block groups, or hexagonal grids) layers. Prior to any spatial operation, strip personally identifiable information (PII), hash facility identifiers, and retain only aggregated capacity metrics (e.g., FTE clinicians, licensed beds, weekly appointment slots). This satisfies HIPAA Safe Harbor and GDPR pseudonymization mandates by design.

CRS misalignment introduces systematic bias in distance calculations and area-weighted aggregations. All inputs must be projected to a locally optimized metric CRS (e.g., UTM zone or state plane) before network routing or Euclidean fallbacks are applied. The following routine enforces projection consistency, repairs invalid geometries, and logs non-compliant records for audit trails:

import geopandas as gpd
import logging
from pyproj import CRS
from shapely.validation import make_valid

logging.basicConfig(level=logging.INFO, format="%(levelname)s: %(message)s")

def standardize_crs_and_validate(
    gdf: gpd.GeoDataFrame, target_epsg: int, layer_name: str
) -> gpd.GeoDataFrame:
    target_crs = CRS.from_epsg(target_epsg)
    if gdf.crs != target_crs:
        logging.info(f"[{layer_name}] Transforming to EPSG:{target_epsg}")
        gdf = gdf.to_crs(target_crs)
    
    # Repair topological errors and drop empty geometries
    gdf["geometry"] = gdf["geometry"].apply(make_valid)
    invalid_mask = gdf.geometry.is_empty
    if invalid_mask.any():
        logging.warning(f"[{layer_name}] Dropped {invalid_mask.sum()} empty/invalid geometries")
    
    return gdf[~invalid_mask].copy()

Reference the official GeoPandas projection documentation for region-specific EPSG selection and transformation accuracy thresholds.

Phase 2: Network-Constrained Travel Matrices & Decay Functions

Euclidean buffers fail to capture real-world mobility constraints, particularly in topographically complex or transit-dependent regions. Production pipelines should generate origin-destination (OD) travel-time matrices using osmnx, pandana, or commercial routing APIs. When network routing fails due to disconnected components or restricted road segments, implement a configurable Euclidean fallback with an explicit warning flag to maintain pipeline continuity.

To prevent O(n²) memory exhaustion in regional deployments, pre-filter the OD matrix using a spatial index (e.g., scipy.spatial.cKDTree or geopandas.sjoin_nearest) before applying the distance-decay function. A Gaussian decay function is standard for healthcare utilization modeling:

import numpy as np
import pandas as pd

def gaussian_decay(dist_meters: np.ndarray, threshold_meters: float, sigma_ratio: float = 0.5) -> np.ndarray:
    """Returns decay weights [0, 1] based on Gaussian kernel."""
    sigma = threshold_meters * sigma_ratio
    weights = np.exp(-0.5 * (dist_meters / sigma) ** 2)
    weights[dist_meters > threshold_meters] = 0.0
    return weights

Phase 3: Deterministic 2SFCA Computation Pipeline

The 2SFCA algorithm executes in two deterministic passes:

  1. Step 1 (Supply Perspective): Calculate the provider-to-population ratio at each facility j, weighted by demand k within the catchment.
  2. Step 2 (Demand Perspective): Aggregate the ratios from all accessible facilities j to each demand location i, again applying distance decay.

The two floating-catchment passes chain together as follows:

flowchart TD
  subgraph S1["Step 1 — Supply perspective"]
    A["For each facility j"] --> B["Sum decay-weighted demand in catchment"]
    B --> C["Ratio R_j = capacity / weighted demand"]
  end
  subgraph S2["Step 2 — Demand perspective"]
    D["For each demand site i"] --> E["Sum decay-weighted R_j of reachable facilities"]
    E --> F["Accessibility score A_i"]
  end
  C --> D

Vectorized execution using pandas merge operations and sparse filtering ensures reproducibility and sub-linear scaling for county-level deployments:

def calculate_2sfca(
    od_matrix: pd.DataFrame,
    supply_gdf: gpd.GeoDataFrame,
    demand_gdf: gpd.GeoDataFrame,
    capacity_col: str,
    population_col: str,
    threshold_meters: float,
    sigma_ratio: float = 0.5
) -> gpd.GeoDataFrame:
    # od_matrix must contain: 'from_id', 'to_id', 'distance_m'
    # Step 1: Facility-to-Population Ratios
    decay_weights = gaussian_decay(od_matrix["distance_m"].values, threshold_meters, sigma_ratio)
    od_matrix["decay"] = decay_weights
    
    # Weighted demand at each facility
    od_matrix["weighted_pop"] = od_matrix["decay"] * od_matrix[population_col]
    facility_ratios = od_matrix.groupby("to_id")["weighted_pop"].sum()
    
    # Avoid division by zero for facilities with zero accessible population
    facility_ratios = facility_ratios.replace(0, np.nan)
    supply_gdf = supply_gdf.merge(facility_ratios, left_on="facility_id", right_index=True, how="left")
    supply_gdf["r_j"] = supply_gdf[capacity_col] / supply_gdf["weighted_pop"]
    supply_gdf["r_j"] = supply_gdf["r_j"].fillna(0.0)
    
    # Step 2: Demand-to-Facility Aggregation
    od_matrix = od_matrix.merge(supply_gdf[["facility_id", "r_j"]], left_on="to_id", right_on="facility_id", how="left")
    od_matrix["access_score"] = od_matrix["decay"] * od_matrix["r_j"]
    demand_access = od_matrix.groupby("from_id")["access_score"].sum().reset_index()
    
    demand_gdf = demand_gdf.merge(demand_access, left_on="demand_id", right_on="from_id", how="left")
    demand_gdf["access_score"] = demand_gdf["access_score"].fillna(0.0)
    return demand_gdf

Phase 4: Edge Cases, Validation & Audit-Ready Output

Production deployments must explicitly handle zero-capacity facilities, isolated populations, and routing timeouts. Implement fallback logic that flags records where network routing returned NaN distances, and document the proportion of Euclidean fallbacks in the pipeline metadata. Use Python’s native logging module to record execution parameters, CRS transformations, and data lineage, ensuring compliance with ISO 19115 geographic metadata standards.

Validate outputs against known benchmarks:

  • Verify that access_score distributions align with historical utilization patterns.
  • Cross-check catchment boundaries against transit deserts or federally designated Health Professional Shortage Areas (HPSAs).
  • Ensure all spatial joins preserve topology and do not artificially inflate accessibility scores through overlapping service areas.

The resulting accessibility surface feeds directly into downstream equity scoring frameworks. When integrated with demographic stratification layers, the 2SFCA index becomes a foundational input for Spatial Equity Index Calculation, enabling targeted resource allocation and policy intervention.

Production Deployment Notes

  • Parallelization: Partition OD matrices by spatial tiles (e.g., 10km × 10km grids) and process asynchronously using concurrent.futures or Dask.
  • Storage: Export intermediate matrices and final outputs to Apache Parquet or GeoPackage to preserve schema integrity and enable columnar querying.
  • Audit Trail: Embed pipeline version, CRS EPSG codes, decay parameters, and execution timestamps in a companion JSON manifest for regulatory review.
  • Monitoring: Implement automated checks for sudden accessibility score deviations (>15% variance) between quarterly runs, triggering manual routing validation before public release.