Using GIS for Public Health Surveillance

Geographic Information Systems (GIS) are essential tools in public health. By overlaying health data with spatial layers — population density, environmental hazards, healthcare facility locations, and transportation networks — GIS reveals patterns that are invisible in tabular data. This article covers the key spatial analysis techniques used in public health surveillance with practical Python examples.

Core GIS Concepts for Health

Spatial data comes in two main formats: vector data (points for events like disease cases or hospital locations, lines for roads and rivers, polygons for administrative boundaries like districts or census tracts) and raster data (continuous surfaces like temperature, elevation, or population density). In public health, point data representing individual cases is often aggregated to polygon boundaries (counties, states) for analysis and visualization to protect patient privacy. The choice of aggregation level matters — the Modifiable Areal Unit Problem (MAUP) means that different boundary definitions can produce different analysis results from the same underlying data.

import geopandas as gpd
import matplotlib.pyplot as plt

# Load health facility locations
facilities = gpd.read_file("facilities.geojson")
print(facilities.head())
print(facilities.crs)  # coordinate reference system

# Load district boundaries
districts = gpd.read_file("districts.geojson")

# Spatial join: count facilities per district
facility_counts = gpd.sjoin(facilities, districts, how="left", predicate="within")
counts = facility_counts.groupby("district_name").size().reset_index(name="facility_count")

# Merge counts with district geometry
districts = districts.merge(counts, on="district_name", how="left")
districts["facility_count"] = districts["facility_count"].fillna(0)

# Plot
districts.plot(column="facility_count", legend=True,
               legend_kwds={"label": "Healthcare Facilities per District"})
plt.title("Healthcare Facility Distribution")
plt.savefig("facility_map.png")

Spatial Clustering and Hotspot Detection

Identifying disease clusters is a core public health surveillance activity. Two common approaches are Kernel Density Estimation (KDE), which creates a smooth surface of case density, and spatial scan statistics (Kulldorff’s method), which identifies circular or elliptical regions with statistically elevated case counts. Moran’s I measures global spatial autocorrelation — whether cases cluster more than expected by chance across the entire study area — while Getis-Ord Gi* identifies local hotspots where high values cluster together. These methods help epidemiologists detect outbreaks early, target interventions, and allocate resources efficiently.

from sklearn.neighbors import KernelDensity
import numpy as np

# Case coordinates (latitude, longitude)
cases = gpd.read_file("disease_cases.geojson")
coords = np.array([(p.x, p.y) for p in cases.geometry])

# Kernel density estimation
kde = KernelDensity(bandwidth=0.05, metric="haversine")
kde.fit(np.radians(coords))

# Evaluate density on a grid
grid_x, grid_y = np.meshgrid(np.linspace(72, 78, 200), np.linspace(18, 22, 200))
grid_coords = np.radians(np.column_stack([grid_x.ravel(), grid_y.ravel()]))
density = np.exp(kde.score_samples(grid_coords)).reshape(grid_x.shape)

# Visualize hotspot
plt.figure(figsize=(10, 8))
plt.contourf(grid_x, grid_y, density, levels=20, cmap="Reds")
plt.scatter(coords[:, 0], coords[:, 1], alpha=0.3, s=10, c="black")
plt.colorbar(label="Case Density")
plt.title("Disease Case Density — Kernel Density Estimate")
plt.savefig("hotspot_map.png")

Spatial Accessibility Analysis

Access to healthcare is not just about distance — road networks, transportation options, and travel times all matter. A simple approach is buffer analysis (show areas within a certain distance of a facility), but more realistic models use network analysis along road networks. The Enhanced Two-Step Floating Catchment Area (E2SFCA) method accounts for both supply (facility capacity) and demand (population) to measure accessibility. In Python, the osmnx library can fetch road networks from OpenStreetMap and compute travel times along actual roads rather than straight-line distances.

# Simple distance-based accessibility
import geopy.distance

def nearest_facility_distance(case_point, facility_points):
    distances = [geopy.distance.distance(
        (case_point.y, case_point.x),
        (facility.y, facility.x)
    ).km for facility in facility_points.geometry]
    return min(distances)

cases["nearest_km"] = cases.geometry.apply(
    lambda p: nearest_facility_distance(p, facilities)
)

# Percentage of population within 5 km of a facility
within_5km = cases[cases["nearest_km"] <= 5]
print(f"Population within 5 km of nearest facility: "
      f"{len(within_5km)} / {len(cases)} ({100*len(within_5km)/len(cases):.1f}%)")
Creating Interactive Maps with Folium
Interactive web maps are powerful tools for communicating spatial health data to stakeholders. Folium wraps Leaflet.js (a leading open-source mapping library) with a Pythonic API, letting you create zoomable, clickable maps with markers, popups, and choropleth layers. You can overlay disease case locations, health facility catchment areas, and district-level statistics on a single interactive map that can be embedded in dashboards or shared as standalone HTML files.
import folium

# Base map centered on the study area
m = folium.Map(location=[20.0, 75.0], zoom_start=5, tiles="OpenStreetMap")

# Add case points with popups
for _, case in cases.iterrows():
    folium.CircleMarker(
        location=[case.geometry.y, case.geometry.x],
        radius=5,
        color="red",
        fill=True,
        popup=f"Date: {case['date']}, Diagnosis: {case['diagnosis']}"
    ).add_to(m)

# Add health facilities
for _, facility in facilities.iterrows():
    folium.Marker(
        location=[facility.geometry.y, facility.geometry.x],
        icon=folium.Icon(color="green", icon="plus", prefix="fa"),
        popup=f"{facility['name']} - {facility['type']}"
    ).add_to(m)

# Add choropleth layer for district-level rates
folium.Choropleth(
    geo_data=districts.to_json(),
    data=counts,
    columns=["district_name", "rate_per_100k"],
    key_on="feature.properties.district_name",
    fill_color="YlOrRd",
    legend_name="Incidence Rate (per 100,000)"
).add_to(m)

m.save("public_health_dashboard.html")
Spatial analysis in public health is most valuable when it leads to action. A well-designed map that shows a cluster of tuberculosis cases near a specific water source, or a gap in immunization coverage in a particular district, provides evidence that can drive resource allocation and policy decisions. The combination of GeoPandas for analysis and Folium for visualization makes Python a complete platform for public health GIS work.
Tools like QGIS (open-source), GeoPandas (Python), and Folium (interactive web maps) make spatial analysis accessible to public health practitioners. For production surveillance systems, platforms like DHIS2 include built-in GIS modules, and custom solutions can be built with PostGIS for spatial databases and GeoServer for map serving. The key is to combine epidemiological domain expertise with spatial thinking — the question is not just "how many cases?" but "where are the cases, and what spatial factors might explain the pattern?"