Spatial Epi Public Health GIS

Computing Greenness Exposure for Birth Cohorts

This guide, part of Remote Sensing Indices for Health Studies, works through a specific and unusually error-prone assignment: residential greenness during pregnancy, where the exposure window is short, dated per subject, and overlaps a seasonal cycle that is itself correlated with the outcome.

Problem Context & Constraints

The straightforward approach takes one annual greenness composite, samples it in a buffer around the maternal address, and calls it exposure. It fails on a birth cohort for a reason that is structural rather than technical: a birth in March and a birth in September have gestational windows covering different halves of the vegetation cycle, and birth month is associated with a long list of outcomes through pathways that have nothing to do with vegetation. Aligning the exposure to the gestational window therefore introduces seasonality into the exposure variable, and any association it produces is confounded by season unless season is handled explicitly.

The opposite error is equally common. Using a fixed summer composite for everyone removes the seasonal confound and also removes the within-subject timing that the trimester hypothesis depends on. If the study claims a second-trimester effect, the exposure must be measured in the second trimester; if it cannot be, the claim has to be about residential greenness generally rather than about timing.

Three constraints follow. The gestational window must be dated per subject, from an estimated conception date, not from birth year. Season must be a covariate, at minimum as birth month, and preferably as a term that allows a smooth seasonal effect. And leaf-off periods must be handled deliberately — a January window in a temperate deciduous region measures bare canopy, which is a real measurement of a different quantity, not a low-greenness neighbourhood.

Why Birth Month Becomes an Exposure Variable A vegetation index cycle over two calendar years, peaking each summer and falling each winter. Two gestational windows are drawn beneath it. The first, for a September birth, covers a full summer and receives a high mean greenness. The second, for a March birth, covers a full winter and receives a low mean greenness. Both mothers live in identical neighbourhoods, so the entire difference in assigned exposure comes from birth timing. A note states that unless season enters the model, the greenness coefficient is partly a birth-month coefficient. Same neighbourhood, two birth dates, two very different exposures summer peak winter trough gestation, September birth mean index 0.71 gestation, March birth mean index 0.34 Without a seasonal term, half the greenness coefficient is measuring birth month

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4, rasterio 1.3.10
  • A dated stack of cloud-masked greenness composites at sub-annual resolution — monthly or 16-day
  • Birth records with birth date, gestational age at birth, and a geocoded maternal address with match type
  • A land-cover or canopy layer if the leaf-off correction is to be applied

Step-by-Step Solution

Compute the gestational window per subject, average the dated composites that fall inside it within the chosen buffer, and return both the window-aligned exposure and a season-neutral annual exposure so the analysis can compare them.

# Gestational-window greenness exposure from a dated composite stack.
# Pinned: numpy==1.26.4, pandas==2.2.2, geopandas==0.14.4, rasterio==1.3.10
import logging

import numpy as np
import pandas as pd
import geopandas as gpd

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger("exposure.greenness")

DEFAULT_BUFFER_M = 300


def gestational_windows(births: pd.DataFrame) -> pd.DataFrame:
    """Conception-to-birth window and trimester boundaries per subject."""
    ga_days = births["gestational_age_weeks"].astype(float) * 7.0
    end = pd.to_datetime(births["birth_date"])
    start = end - pd.to_timedelta(ga_days, unit="D")
    out = pd.DataFrame({"study_id": births["study_id"].to_numpy(),
                        "win_start": start.to_numpy(), "win_end": end.to_numpy()})
    out["t1_end"] = out["win_start"] + (out["win_end"] - out["win_start"]) / 3.0
    out["t2_end"] = out["win_start"] + 2.0 * (out["win_end"] - out["win_start"]) / 3.0
    out["birth_month"] = end.dt.month.to_numpy()
    return out


def buffer_means_by_date(addresses: gpd.GeoDataFrame, composites: pd.DataFrame,
                         sample_fn, buffer_m: int = DEFAULT_BUFFER_M) -> pd.DataFrame:
    """Mean index inside the buffer for every composite date, per address.

    composites: date, raster_path (one row per dated composite)
    sample_fn:  callable(raster_path, GeoDataFrame_of_buffers) -> array of means
    """
    buf = addresses.copy()
    buf["geometry"] = buf.geometry.buffer(buffer_m)
    frames = []
    for rec in composites.sort_values("date").itertuples():
        vals = sample_fn(rec.raster_path, buf)
        frames.append(pd.DataFrame({"study_id": buf["study_id"].to_numpy(),
                                    "date": rec.date, "index": vals}))
    out = pd.concat(frames, ignore_index=True)
    log.info("sampled %d composite dates for %d addresses at %d m",
             composites.shape[0], len(buf), buffer_m)
    return out


def window_exposure(series: pd.DataFrame, windows: pd.DataFrame) -> pd.DataFrame:
    """Mean index inside the gestational window, per trimester, plus an annual mean."""
    m = series.merge(windows, on="study_id", how="inner")
    m["date"] = pd.to_datetime(m["date"])
    inside = m[(m["date"] >= m["win_start"]) & (m["date"] <= m["win_end"])]

    def agg(d):
        t1 = d[d["date"] <= d["t1_end"].iloc[0]]["index"]
        t2 = d[(d["date"] > d["t1_end"].iloc[0]) & (d["date"] <= d["t2_end"].iloc[0])]["index"]
        t3 = d[d["date"] > d["t2_end"].iloc[0]]["index"]
        return pd.Series({
            "green_gestation": float(d["index"].mean()),
            "green_t1": float(t1.mean()) if len(t1) else np.nan,
            "green_t2": float(t2.mean()) if len(t2) else np.nan,
            "green_t3": float(t3.mean()) if len(t3) else np.nan,
            "n_composites": int(d["index"].notna().sum()),
        })

    win = inside.groupby("study_id").apply(agg, include_groups=False).reset_index()

    # Season-neutral comparator: the same buffer averaged over ALL dates, which
    # carries no birth-timing information and isolates the neighbourhood effect.
    annual = (series.groupby("study_id")["index"].mean()
              .rename("green_annual").reset_index())
    out = win.merge(annual, on="study_id", how="outer").merge(
        windows[["study_id", "birth_month"]], on="study_id", how="left")
    thin = out["n_composites"] < 6
    out.loc[thin, ["green_gestation", "green_t1", "green_t2", "green_t3"]] = np.nan
    log.info("window exposure: n=%d thin_windows=%d median_composites=%.0f",
             len(out), int(thin.sum()), out["n_composites"].median())
    return out

Returning both the window-aligned and the season-neutral exposure is the design point. The annual mean over all dates at the same address contains no birth-timing information at all, so it estimates the neighbourhood greenness effect cleanly. The gestational mean adds timing and, inseparably, season. Running the analysis with both and comparing tells you how much of any association depends on the timing alignment, which is the question a reviewer will ask.

The trimester split uses the actual gestational age rather than a fixed 40 weeks, because preterm birth is frequently the outcome and assuming term length would make the exposure window a function of the outcome — a straightforward way to manufacture an association.

Three Exposure Definitions From One Address Three exposure variables derived from the same maternal address. The annual mean over all composite dates carries no birth-timing information and supports a claim about neighbourhood greenness. The gestational mean aligns to each subject's window and supports a claim about exposure during pregnancy, but is confounded with season and requires a seasonal term. The trimester-specific means support a timing claim and are the most confounded, because a trimester is a season for any given birth date. Each is annotated with the covariate adjustment it requires. More timing detail means more seasonal confounding, always Annual mean all dates, same buffer supports: neighbourhood greenness effect no seasonal adjustment needed Gestational mean dates inside the window supports: exposure during pregnancy requires a birth-month or smooth seasonal term Trimester means window split in three supports: a timing hypothesis a trimester is a season; adjust and expect little left Report all three; the pattern across them is more informative than any one estimate An effect present in the annual mean and absent in the trimester means is a neighbourhood effect

Validation & Edge Cases

Preterm births and truncated windows. A birth at 28 weeks has a shorter window with fewer composites in it. The composite-count threshold nulls the extreme cases, but the remainder still have systematically fewer observations, and that correlates with the outcome. Report the composite count as a covariate.

Maternal address is the address at delivery. Residential moves during pregnancy are common. Where a prenatal address is recorded, use it; where it is not, the exposure describes the delivery address and should be described as such.

Leaf-off windows in deciduous regions. A winter gestational window returns a genuinely low index, which is a correct measurement of bare canopy and a poor measure of the neighbourhood’s vegetation. Where the hypothesis is about vegetation presence rather than active vegetation, a canopy-cover product is a better exposure and carries no seasonality.

Composite cadence coarser than a trimester. A 16-day composite gives roughly six observations per trimester, which is adequate; a monthly composite gives three, and a quarterly one gives a single value that is not a mean of anything. Check the cadence against the shortest window the study will estimate before committing to a product, because no amount of care downstream recovers timing resolution the stack never had.

Address buffers that overlap water. A coastal or lakeside address buffer includes open water, which returns a strongly negative index and drags the mean down in a way unrelated to vegetation. Mask permanent water from the buffer and report the remaining land area, so two addresses with the same land greenness are not separated by how much sea their circle happened to include.

Evergreen and arid regions. In an evergreen forest the seasonal amplitude is small and the confounding is mild; in an arid city the amplitude is driven by irrigation rather than climate. The severity of the seasonal problem is regional, and it should be measured — plot the index cycle for the study area before choosing the design.

Buffer radius and the address itself. A 300 m buffer around a rural address includes farmland whose index dwarfs any residential vegetation. Consider excluding agricultural land cover where the hypothesis is about residential greenery.

INFO exposure.greenness sampled 46 composite dates for 22418 addresses at 300 m
INFO exposure.greenness window exposure: n=22418 thin_windows=412 median_composites=19

Measuring the seasonal amplitude of the study area before choosing a design takes minutes and tells you how severe the confounding will be:

The Seasonal Problem Is Regional The annual cycle of a greenness index in three regions. In a temperate deciduous region the index swings widely between summer and winter, so gestational windows falling in different seasons receive very different exposures and seasonal confounding is severe. In an evergreen region the swing is small and the confounding mild. In an arid irrigated city the cycle follows watering schedules rather than climate and may be out of phase with the other two. The amplitude should be measured before the analysis design is chosen. Plot your own study area's cycle before choosing the design temperate deciduous: large swing evergreen arid, irrigated: out of phase Jan Jul Dec greenness index through one year in three regions

Compliance Notes

Record the composite stack identifier and dates, the buffer radius and its justification, the gestational window derivation including the use of observed rather than assumed gestational age, the composite-count threshold, and both the window-aligned and season-neutral exposures with their distributions by birth month.

Birth records are protected health data and the maternal address is the identifying element, so the buffer sampling runs inside the protected environment and the outputs are keyed by study identifier. One study-design item belongs in the record as well: whether birth month or a smooth seasonal term entered the model. An analysis that assigns a season-dependent exposure without a seasonal adjustment has a known confounder left uncontrolled, and a reader needs to be able to check that rather than assume it.