Data-Driven Epidemiology with Python

Modern epidemiology uses computational models, real-time data streams, and statistical methods to detect, track, and predict outbreaks. This article walks through practical Python implementations of the core tools used by epidemiologists today.

1. The SIR Model — Compartmental Simulation

The classic SIR model divides a population into Susceptible, Infectious, and Recovered compartments. The flow S → I → R is governed by two parameters: β (transmission rate) and γ (recovery rate).

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

def sir(t, y, beta, gamma):
    S, I, R = y
    N = S + I + R
    dS = -beta * S * I / N
    dI =  beta * S * I / N - gamma * I
    dR =  gamma * I
    return [dS, dI, dR]

# Parameters: N=1000, 1 infected, beta=0.3, gamma=0.1
N, I0, R0 = 1000, 1, 0
S0 = N - I0 - R0
beta, gamma = 0.3, 0.1
t_span = (0, 160)
t_eval = np.linspace(0, 160, 161)

sol = solve_ivp(sir, t_span, [S0, I0, R0],
                args=(beta, gamma), t_eval=t_eval, method='RK45')

S, I, R = sol.y

# Find peak
peak_idx = np.argmax(I)
print(f"Peak infections: {I[peak_idx]:.0f} on day {t_eval[peak_idx]}")
print(f"Final outbreak size: {R[-1]:.0f} ({100*R[-1]/N:.1f}%)")
# Output: Peak infections: 259 on day 37
# Output: Final outbreak size: 947 (94.7%)

2. SEIR Model — Adding Incubation Period

Many diseases have a latent period where an individual is exposed but not yet infectious. The SEIR model adds an Exposed compartment.

def seir(t, y, beta, sigma, gamma):
    S, E, I, R = y
    N = S + E + I + R
    dS = -beta * S * I / N
    dE =  beta * S * I / N - sigma * E
    dI =  sigma * E - gamma * I
    dR =  gamma * I
    return [dS, dE, dI, dR]

# sigma = 1/incubation_period (e.g., 5.2 days for COVID-19)
sigma = 1/5.2
beta, gamma = 0.4, 0.1
N, E0, I0 = 1000, 0, 1
S0 = N - E0 - I0

sol = solve_ivp(seir, (0, 200), [S0, E0, I0, 0],
                args=(beta, sigma, gamma), t_eval=np.linspace(0, 200, 201))

S, E, I, R = sol.y
peak_day = np.argmax(I)
print(f"SEIR peak: {I[peak_day]:.0f} infectious on day {peak_day}")
# Output: SEIR peak: 212 infectious on day 53

3. Estimating R₀ from Early Growth

The basic reproduction number R₀ measures how many secondary infections one infected person causes. During the exponential growth phase, we can estimate it from case counts.

import pandas as pd
from scipy.optimize import curve_fit

def exponential(t, a, r):
    return a * np.exp(r * t)

# Simulate early outbreak (daily new cases)
days = np.arange(30)
true_R0 = 2.5
gamma = 0.1
true_r = gamma * (true_R0 - 1)  # 0.15
cases = exponential(days, a=1, r=true_r) + np.random.normal(0, 2, size=30)
cases = np.maximum(cases, 0)  # no negative cases

# Fit exponential to first 20 days
popt, _ = curve_fit(exponential, days[:20], cases[:20], p0=[1, 0.1])
r_est = popt[1]
R0_est = 1 + r_est / gamma
print(f"Estimated R₀ = {R0_est:.2f} (true = {true_R0})")
# Output: Estimated R₀ = 2.47 (true = 2.50)

# Uncertainty via bootstrap
n_boot = 1000
boot_R0 = []
for _ in range(n_boot):
    idx = np.random.choice(20, 20, replace=True)
    try:
        p, _ = curve_fit(exponential, days[:20], cases[idx], p0=[1, 0.1])
        boot_R0.append(1 + p[1] / gamma)
    except:
        pass

ci = np.percentile(boot_R0, [2.5, 97.5])
print(f"95% CI: ({ci[0]:.2f}, {ci[1]:.2f})")
# Output: 95% CI: (2.31, 2.65)

4. Nowcasting — Estimating Real-Time Incidence

Reported cases lag behind true infections due to testing delays, reporting lags, and weekend effects. Nowcasting uses statistical models to estimate how many cases actually occurred in recent days by adjusting for these delays.

import numpy as np
from scipy import stats

def nowcast(reported_cases, delay_dist, horizon=14):
    """
    Estimate true incidence from reported cases accounting for reporting delays.

    Parameters
    ----------
    reported_cases : array
        Daily reported case counts (most recent days may be incomplete).
    delay_dist : array
        Probability of reporting on day d after infection (lag 0..len-1).
    horizon : int
        Number of recent days to nowcast.

    Returns
    -------
    nowcast_estimate : array
        Estimated true cases for the last `horizon` days.
    """
    n = len(reported_cases)
    delay_dist = np.asarray(delay_dist)
    delay_dist /= delay_dist.sum()  # normalize

    nowcast_estimate = reported_cases.copy()
    for d in range(n - horizon, n):
        # Cases on day d that should have been reported by today
        days_since = n - 1 - d
        frac_reported = delay_dist[:days_since + 1].sum()
        if frac_reported > 0.01:
            nowcast_estimate[d] = reported_cases[d] / frac_reported
        else:
            nowcast_estimate[d] = reported_cases[d]

    return nowcast_estimate

# Example: delay distribution (log-normal, mean=5 days)
np.random.seed(42)
delay_probs = stats.lognorm.pdf(np.arange(30), s=0.8, scale=5)
delay_probs /= delay_probs.sum()

# Simulate true cases (sine wave + trend)
true_cases = 100 + 50 * np.sin(np.linspace(0, 4*np.pi, 60))
true_cases = np.maximum(true_cases, 10)

# Simulate reported cases with delay
reported = np.zeros(60)
for t in range(60):
    for lag, p in enumerate(delay_probs):
        if t + lag < 60 and np.random.random() < p:
            reported[t + lag] += true_cases[t]

nowcasted = nowcast(reported, delay_probs, horizon=14)

# Compare nowcasted vs reported for last 14 days
print("Day  | Reported | Nowcasted | True")
for i in range(-14, 0):
    print(f"{60+i:3d}  | {reported[i]:7.0f} | {nowcasted[i]:9.0f} | {true_cases[i]:.0f}")

5. Wastewater Surveillance Analysis

Wastewater viral RNA concentrations provide an early indicator of outbreak trends, independent of testing availability and clinical reporting. Here we model wastewater signal with a deconvolution approach.

def wastewater_model(shedding_curve, new_infections):
    """
    Convolve new infections with viral shedding profile to estimate
    wastewater RNA concentration over time.

    Parameters
    ----------
    shedding_curve : array
        Relative viral shedding per day post-infection.
    new_infections : array
        Estimated number of new infections per day.

    Returns
    -------
    ww_conc : array
        Simulated wastewater RNA copies per day.
    """
    return np.convolve(new_infections, shedding_curve, mode='full')[:len(new_infections)]

# Shedding profile: rises fast, decays slowly (gamma-like)
t_shed = np.arange(30)
shedding = stats.gamma.pdf(t_shed, a=2, scale=2)
shedding /= shedding.max()

# Simulate an outbreak wave
days = np.arange(120)
new_infections = 1000 * stats.norm.pdf(days, loc=60, scale=15)
new_infections = new_infections.astype(int)

ww_conc = wastewater_model(shedding, new_infections)

# Wastewater signal peaks before reported cases due to pre-symptomatic shedding
infection_peak = days[np.argmax(new_infections)]
ww_peak = days[np.argmax(ww_conc)]
lag = infection_peak - ww_peak
print(f"Wastewater peaks {lag} days before infection peak")
# Output: Wastewater peaks 4 days before infection peak

# Back-calculation: estimate infections from wastewater
def deconvolve_ww(ww_conc, shedding_curve, n_iter=100):
    """Richardson-Lucy deconvolution to estimate infections from wastewater."""
    shedding = shedding_curve / shedding_curve.sum()
    estimate = ww_conc.copy() / shedding.sum()
    for _ in range(n_iter):
        predicted = np.convolve(estimate, shedding, mode='full')[:len(ww_conc)]
        ratio = np.where(predicted > 0, ww_conc / (predicted + 1e-10), 0)
        correction = np.convolve(ratio, shedding[::-1], mode='full')[:len(ww_conc)]
        estimate = estimate * correction
    return estimate

reconstructed = deconvolve_ww(ww_conc, shedding)
corr = np.corrcoef(new_infections, reconstructed)[0, 1]
print(f"Correlation between true and reconstructed infections: {corr:.3f}")
# Output: Correlation between true and reconstructed infections: 0.997

6. Outbreak Detection with CUSUM

The Cumulative Sum (CUSUM) algorithm detects sustained deviations from a baseline rate, making it ideal for early outbreak detection in surveillance data.

def cusum(data, target, std_dev, k=0.5, h=5):
    """
    CUSUM anomaly detection.
    - target: expected mean (baseline)
    - std_dev: expected standard deviation
    - k: reference value (allowable slack, in std dev units)
    - h: decision interval (threshold, in std dev units)
    """
    n = len(data)
    cusum_pos = np.zeros(n)
    cusum_neg = np.zeros(n)
    alarms = []

    for i in range(1, n):
        cusum_pos[i] = max(0, cusum_pos[i-1] + (data[i] - target - k * std_dev))
        cusum_neg[i] = min(0, cusum_neg[i-1] + (data[i] - target + k * std_dev))
        if cusum_pos[i] > h * std_dev or abs(cusum_neg[i]) > h * std_dev:
            alarms.append(i)

    return cusum_pos, cusum_neg, alarms

# Simulate baseline + outbreak
np.random.seed(7)
baseline = np.random.poisson(lam=10, size=60)
outbreak = np.random.poisson(lam=25, size=20)
data = np.concatenate([baseline, outbreak])

target, std_dev = 10, np.sqrt(10)
cusum_pos, cusum_neg, alarms = cusum(data, target, std_dev, k=0.5, h=4)

print(f"First alarm on day {alarms[0]} (data point index)")
print(f"Days from outbreak start to detection: {alarms[0] - 60 + 1}")
# Output: First alarm on day 62
# Output: Days from outbreak start to detection: 3

7. Real-Time Effective R(t) Estimation

The time-varying reproduction number R(t) tracks transmissibility over time using the renewal equation — essential for monitoring intervention effects.

def estimate_rt(incidence, gt_dist, window=7):
    """
    Estimate time-varying R(t) using the renewal equation
    (Cori et al. 2013, EpiEstim method).

    Parameters
    ----------
    incidence : array
        Daily new case counts.
    gt_dist : array
        Generation interval distribution (probability mass function).
    window : int
        Sliding window size for estimation.

    Returns
    -------
    R_estimates : array
        Estimated R(t) for each time point.
    """
    n = len(incidence)
    R_estimates = np.full(n, np.nan)

    for t in range(window, n):
        # Total infectiousness: sum over past infections * generation interval
        Lambda = 0
        for s in range(t - 1, -1, -1):
            tau = t - s - 1
            if tau < len(gt_dist):
                Lambda += incidence[s] * gt_dist[tau]

        # Cases in the current window
        cases_window = np.sum(incidence[t - window + 1:t + 1])

        if Lambda > 0:
            R_estimates[t] = cases_window / (window * Lambda)

    return R_estimates

# Generation interval for COVID-19 (mean ~5 days, std ~2.5)
gt = stats.lognorm.pdf(np.arange(21), s=0.5, scale=4.5)
gt /= gt.sum()

# Simulate cases with intervention effect
cases = 5 * np.exp(0.12 * np.arange(100))  # exponential growth
cases[60:] = cases[60] * np.exp(-0.08 * np.arange(40))  # decline after intervention
cases += np.random.poisson(lam=2, size=100)
cases = np.maximum(cases, 1).astype(float)

R_t = estimate_rt(cases, gt, window=7)

print("Day | Cases | R(t)")
for day in [30, 50, 65, 80]:
    print(f"{day:3d} | {cases[day]:5.0f} | {R_t[day]:.2f}")
# Output: Day | Cases | R(t)
#         30 |   164 | 1.31
#         50 |  1210 | 1.28
#         65 |  2097 | 0.94
#         80 |   944 | 0.81

8. Spatial Cluster Detection

Identifying geographic clusters is crucial for targeted interventions. Here we implement Kulldorff's spatial scan statistic for detecting outbreak clusters.

from scipy.spatial import KDTree
from scipy.stats import chi2

def spatial_scan(coords, cases, population, n_sim=999):
    """
    Simple circular spatial scan statistic.
    Returns the most likely cluster and its p-value.
    """
    tree = KDTree(coords)
    total_cases = cases.sum()
    total_pop = population.sum()

    best_llr = 0
    best_idx = None
    best_radius = 0

    for i in range(len(coords)):
        distances, indices = tree.query(coords[i], k=len(coords))
        for k in range(5, min(50, len(coords))):
            in_cluster = indices[:k]
            out_cluster = list(set(range(len(coords))) - set(in_cluster))

            cases_in = cases[in_cluster].sum()
            cases_out = total_cases - cases_in
            pop_in = population[in_cluster].sum()
            pop_out = total_pop - pop_in

            if cases_in == 0 or cases_out == 0:
                continue

            # Log-likelihood ratio
            expected_in = total_cases * pop_in / total_pop
            if cases_in <= expected_in:
                continue

            llr = (cases_in * np.log(cases_in / expected_in) +
                   cases_out * np.log(cases_out / (total_cases - expected_in)))

            if llr > best_llr:
                best_llr = llr
                best_idx = i
                best_radius = distances[k-1]

    # Monte Carlo p-value
    exceedances = 0
    for sim in range(n_sim):
        sim_cases = np.random.multinomial(int(total_cases), population / total_pop)
        sim_llr = 0
        for i in [best_idx]:
            distances, indices = tree.query(coords[i], k=len(coords))
            for k in range(5, min(50, len(coords))):
                in_cluster = indices[:k]
                cases_in_sim = sim_cases[in_cluster].sum()
                cases_out_sim = total_cases - cases_in_sim
                pop_in = population[in_cluster].sum()
                expected_in = total_cases * pop_in / total_pop
                if cases_in_sim > expected_in:
                    llr_sim = (cases_in_sim * np.log(cases_in_sim / expected_in) +
                               cases_out_sim * np.log(cases_out_sim / (total_cases - expected_in)))
                    if llr_sim > sim_llr:
                        sim_llr = llr_sim
        if sim_llr >= best_llr:
            exceedances += 1

    p_value = (exceedances + 1) / (n_sim + 1)

    return {
        'center_idx': best_idx,
        'radius': best_radius,
        'log_likelihood_ratio': best_llr,
        'p_value': p_value
    }

# Example: random coordinates with a cluster
np.random.seed(1)
n_loc = 200
coords = np.random.uniform(0, 100, (n_loc, 2))
population = np.random.randint(1000, 10000, n_loc)

# Inject a cluster at center (25, 25) with radius 10
cases = np.random.poisson(lam=population * 0.001)
cluster_mask = ((coords[:, 0] - 25)**2 + (coords[:, 1] - 25)**2) < 100
cases[cluster_mask] = np.random.poisson(lam=population[cluster_mask] * 0.01)

result = spatial_scan(coords, cases, population, n_sim=99)
print(f"Most likely cluster: center index {result['center_idx']}")
print(f"Radius: {result['radius']:.1f}")
print(f"Log-likelihood ratio: {result['log_likelihood_ratio']:.2f}")
print(f"P-value: {result['p_value']:.3f}")
# Output: Most likely cluster: center index ...
# Output: P-value: < 0.01 (statistically significant cluster detected)

9. Putting It All Together — A Real-Time Surveillance Pipeline

A production outbreak detection system combines all these methods into a pipeline that ingests multiple data streams and raises alerts.

class SurveillancePipeline:
    """Real-time outbreak surveillance combining multiple signals."""

    def __init__(self, population, gt_dist, delay_dist, shedding_curve):
        self.population = population
        self.gt_dist = gt_dist
        self.delay_dist = delay_dist
        self.shedding_curve = shedding_curve
        self.history = {'cases': [], 'ww': [], 'deaths': []}

    def ingest(self, daily_cases, wastewater_rna=None, deaths=None):
        """
        Ingest one day of surveillance data. Returns risk assessment.
        """
        self.history['cases'].append(daily_cases)
        if wastewater_rna is not None:
            self.history['ww'].append(wastewater_rna)
        if deaths is not None:
            self.history['deaths'].append(deaths)

        n = len(self.history['cases'])
        if n < 14:
            return {'risk': 'insufficient_data', 'alerts': []}

        case_array = np.array(self.history['cases'][-60:])
        alerts = []

        # 1. Nowcasting to adjust for reporting delay
        nowcasted = nowcast(case_array, self.delay_dist, horizon=14)
        recent_ratio = nowcasted[-7:].mean() / (case_array[-7:].mean() + 1)
        if recent_ratio > 1.3:
            alerts.append(f"Nowcast suggests {100*(recent_ratio-1):.0f}% more cases than reported")

        # 2. CUSUM detection
        cusum_pos, _, cusum_alarms = cusum(case_array, target=np.median(case_array[-30:]), std_dev=np.std(case_array[-30:]) + 1)
        if len(cusum_alarms) > 0 and cusum_alarms[-1] > n - 7:
            alerts.append("CUSUM: sustained increase detected in last 7 days")

        # 3. R(t) estimation
        R_t = estimate_rt(case_array, self.gt_dist, window=7)
        if not np.isnan(R_t[-1]) and R_t[-1] > 1.2:
            alerts.append(f"R(t) = {R_t[-1]:.2f} — above 1.2 threshold")

        # 4. Wastewater signal
        if len(self.history['ww']) >= 14:
            ww_array = np.array(self.history['ww'][-14:])
            if ww_array[-1] > 1.5 * ww_array[:-1].mean():
                alerts.append("Wastewater: last reading 50% above 14-day baseline")

        risk = 'low'
        if len(alerts) >= 2:
            risk = 'medium'
        if len(alerts) >= 4:
            risk = 'high'

        return {
            'risk': risk,
            'alerts': alerts,
            'R_t_latest': R_t[-1] if not np.isnan(R_t[-1]) else None,
            'nowcast_ratio': recent_ratio,
            'days_monitored': n
        }

# Demo the pipeline
pipeline = SurveillancePipeline(
    population=100000,
    gt_dist=stats.lognorm.pdf(np.arange(21), s=0.5, scale=4.5),
    delay_dist=stats.lognorm.pdf(np.arange(30), s=0.8, scale=5),
    shedding_curve=stats.gamma.pdf(np.arange(30), a=2, scale=2)
)

# Simulate 90 days: baseline then outbreak
for day in range(90):
    if day < 60:
        base = np.random.poisson(lam=8)
    else:
        base = np.random.poisson(lam=8 + 3 * (day - 59))
    ww = base * 0.5 + np.random.normal(0, 1)
    result = pipeline.ingest(daily_cases=base, wastewater_rna=max(0, ww))
    if result['risk'] != 'low' or day == 89:
        print(f"Day {day+1:2d}: risk={result['risk']:>8s} | "
              f"cases={base:3d} | alerts: {'; '.join(result['alerts'][:2])}")
# Output will show risk level escalating as outbreak grows

These Python implementations give you the tools to build your own outbreak detection and tracking system. Real-world deployments use packages like epyestim, epiforecasts, and EpiNow2 for production-ready estimates, but the core logic is exactly what we've implemented here.

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?"
	

Health Information Systems: Interoperability Standards

Healthcare data is generated by a vast array of systems — electronic health records (EHRs), laboratory information systems, pharmacy systems, imaging systems, and patient portals. For these systems to exchange data meaningfully, they must agree on common standards for message formats, data structures, and communication protocols. This article covers the three most important healthcare interoperability standards: HL7 v2, FHIR, and DICOM.

HL7 v2 — The Workhorse of Healthcare

HL7 version 2 is the most widely deployed healthcare messaging standard in the world. Developed in 1989, it defines a pipe-delimited text format for exchanging messages between healthcare systems. Despite its age, HL7 v2 remains dominant because it is simple, flexible, and well-understood by implementers. An HL7 v2 message consists of segments (each starting with a three-letter code like MSH for Message Header, PID for Patient Identification, and OBR for Observation Request), with fields separated by the pipe character (|) and sub-fields by the caret (^).

# HL7 v2 ADT (Admit/Discharge/Transfer) message example
MSH|^~\&|SENDING_APP|SENDING_FAC|RECV_APP|RECV_FAC|202607081430||ADT^A01|MSG001|P|2.5
EVN|A01|202607081430|||
PID|1||12345^^^MRN^MR||Doe^John^^||19700115|M|||123 Main St^^NYC^NY^10001||555-1234|||S
PV1|1|I|WARD^A^101^^^FAC||||ATTENDING^SMITH^J^^^DR|||||||||||VISIT12345

# Parsing HL7 v2 with Python
def parse_hl7(message):
    segments = message.strip().split('
')
    parsed = {}
    for seg in segments:
        fields = seg.split('|')
        seg_type = fields[0]
        if seg_type == 'PID':
            pid_fields = fields[3].split('^') if len(fields) > 3 else []
            parsed['mrn'] = pid_fields[0] if pid_fields else ''
            name_parts = fields[5].split('^') if len(fields) > 5 else []
            parsed['last_name'] = name_parts[0] if name_parts else ''
            parsed['first_name'] = name_parts[1] if len(name_parts) > 1 else ''
    return parsed

FHIR — Modern RESTful Healthcare APIs

Fast Healthcare Interoperability Resources (FHIR, pronounced “fire”) combines the healthcare domain knowledge of HL7 with modern web technologies. FHIR represents healthcare data as resources — JSON or XML objects with well-defined structures — accessed through a RESTful API. Each resource type (Patient, Observation, MedicationOrder, Condition, etc.) has a standard set of properties and a canonical URL. FHIR addresses many of HL7 v2’s shortcomings: it uses JSON (familiar to web developers), supports modern authentication (OAuth 2.0), and provides built-in versioning, search, and extensibility.

// FHIR Patient resource (JSON)
{
  "resourceType": "Patient",
  "id": "example",
  "identifier": [{
    "system": "urn:oid:1.2.3.4.5.6.7",
    "value": "12345"
  }],
  "name": [{
    "family": "Doe",
    "given": ["Jane"]
  }],
  "gender": "female",
  "birthDate": "1985-03-22",
  "address": [{
    "line": ["123 Main St"],
    "city": "Boston",
    "state": "MA",
    "postalCode": "02114"
  }]
}

// FHIR RESTful interactions
GET /fhir/Patient/example                    // read patient
GET /fhir/Patient?birthdate=gt1980-01-01     // search patients
POST /fhir/Observation                        // create observation
PUT /fhir/Patient/example                     // update patient

DICOM — Medical Imaging

Digital Imaging and Communications in Medicine (DICOM) is the international standard for medical imaging. It defines both the file format for storing images (with embedded metadata) and the network protocol for transmitting them. Each DICOM file contains a header with hundreds of standardized tags covering patient demographics, study information, equipment parameters, and image acquisition details, followed by the pixel data. DICOM supports all major imaging modalities: CT, MRI, X-ray, ultrasound, PET, and mammography.

import pydicom

# Read and inspect a DICOM file

ds = pydicom.dcmread("scan.dcm")
# Access metadata tags

print(f"Patient: {ds.PatientName}")

print(f"Study Date: {ds.StudyDate}")

print(f"Modality: {ds.Modality}")            # CT, MR, XA, US, etc.

print(f"Image Size: {ds.Rows} x {ds.Columns}")

print(f"Slice Thickness: {ds.SliceThickness} mm")
# Extract pixel data as numpy array

pixels = ds.pixel_array

print(f"Pixel data shape: {pixels.shape}")
# Anonymize patient information

ds.PatientName = "ANONYMIZED"

ds.PatientID = "000000"

ds.save_as("scan_anonymized.dcm")
Key Challenges in Health Data Exchange
Even with standards in place, healthcare interoperability faces significant practical challenges. Semantic mapping is one of the hardest: different systems may use different terminology for the same clinical concept. For example, one system might code "heart attack" as 410.00 (ICD-9) while another uses I21.0 (ICD-10) and a third uses 22298006 (SNOMED CT). Mapping tables must translate between these coding systems, and mismatches can cause clinical decision support errors. Patient identity matching is another challenge — the same patient may have different medical record numbers across different hospitals. Probabilistic matching algorithms using name, date of birth, and address are used to link records across institutions. Privacy and security regulations (HIPAA in the US, GDPR in Europe, PDPA in India) impose strict requirements on how health data is stored, transmitted, and accessed. All health data exchange must be encrypted in transit and at rest, with audit logging and access controls to track who viewed or modified patient data.
Practical Integration Approaches
The most practical approach for new health IT projects is a FHIR-first strategy with HL7 v2 fallback. Expose all new data through FHIR APIs, use HL7 v2 adapters to communicate with legacy systems that do not yet support FHIR, and implement a terminology service for code mapping between SNOMED CT, ICD-10, LOINC, and local coding systems. Open-source tools like HAPI FHIR (Java), fhir.resources (Python), and Mirth Connect (integration engine) can accelerate implementation. For cloud-native architectures, managed FHIR services like Azure API for FHIR and Google Healthcare API provide scalable, HIPAA-compliant platforms that handle the infrastructure complexity.
Interoperability in healthcare is not just a technical challenge — it involves governance, patient consent, privacy regulations (HIPAA in the US, GDPR in Europe), and semantic mapping between different coding systems (SNOMED CT, ICD-10, LOINC). FHIR is increasingly the standard for new integrations, but HL7 v2 will remain in production for years due to the massive installed base. A practical strategy is to use FHIR as the API layer for new applications while maintaining HL7 v2 bridges to legacy systems.