Spatial Epi Public Health GIS

Combining Land-Use Regression With Satellite Aerosol Estimates

This guide, part of Land-Use Regression for Air Pollution Exposure, addresses the practical problem of a state-wide study where the land-use model was fitted in one metropolitan area and the rest of the state has neither monitors nor a campaign — but does have complete satellite coverage.

Problem Context & Constraints

Satellite aerosol optical depth is a column measurement: it describes how much light the whole atmospheric column scatters, from the ground to the top of the atmosphere. Surface particulate concentration, which is what people breathe, is a ground-level quantity. The two are related, and the relationship depends on the vertical profile of the aerosol, the humidity, and the season — all of which vary. A satellite product therefore carries genuine regional information about particulate loading and cannot be read directly as a surface concentration.

Land-use regression is the mirror image. It describes surface concentration well at street scale within the area sampled by its campaign, and it has nothing to say about a county 200 km away where its predictors take values outside the sampled range.

Combining them is attractive precisely because the errors are complementary: the satellite supplies the regional field, the land-use model supplies the local contrast. The combination fails in two ways if done carelessly. Double counting happens when a predictor already responds to regional background — population density within 5 km is effectively an urbanity index — and the satellite term adds the same signal again, splitting one effect across two coefficients whose sum is unstable. And resolution laundering happens when a 1 km satellite product is combined with 30 m land-use predictors and the result is presented at 30 m, implying a spatial precision the regional component does not have.

The constraint that resolves both is fitting the two components against measurements together, in one model, rather than fitting each separately and adding the surfaces.

What Each Component Contributes Along a Transect Three stacked profiles along the same transect running from a rural area through a city centre and out again. The satellite component is a broad smooth rise over the urban area with no street-scale detail, representing regional loading. The land-use component is a set of sharp peaks at road corridors with no regional trend, because the regional level is absorbed into the intercept. The hybrid is their sum, showing both the broad rise and the sharp peaks. A warning notes that a large-buffer population predictor would reproduce the broad rise on its own, so fitting the two components separately and adding them counts the urban background twice. Complementary components — but only if each keeps to its own scale Satellite: regional loading Land use: street-scale contrast Hybrid: the sum, fitted jointly Risk: a 5 km population predictor reproduces the top profile by itself fit both components in one regression so the shared signal is attributed once

Prerequisites

  • python 3.11
  • numpy 1.26.4, pandas 2.2.2, geopandas 0.14.4, rasterio 1.3.10, scikit-learn 1.5.0
  • A gap-filled satellite aerosol product on a regular grid, with its own quality flags retained
  • Buffer predictors built as in the sibling guide, restricted to radii below about 1 km so the families stay separable
  • Ground measurements at every site, which is what both components are fitted against

Step-by-Step Solution

Sample the satellite field at each site, add it as one more predictor alongside the local land-use candidates, and let the joint fit apportion the shared variance. Constrain the land-use candidate set to small buffers so it cannot compete with the satellite term for the regional signal.

# Hybrid exposure model: satellite regional term + local land-use predictors, fitted jointly.
# 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
import rasterio
from rasterio.sample import sample_gen
from sklearn.linear_model import LinearRegression

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

LOCAL_MAX_BUFFER_M = 1000     # keep land-use candidates below the satellite cell size


def sample_satellite(points: gpd.GeoDataFrame, raster_path: str,
                     band: int = 1, qa_band: int | None = 2) -> pd.DataFrame:
    """Sample the satellite field (and its QA band) at each point."""
    with rasterio.open(raster_path) as src:
        if points.crs != src.crs:
            points = points.to_crs(src.crs)
        coords = [(g.x, g.y) for g in points.geometry]
        vals = np.array(list(sample_gen(src, coords)), dtype=float)
        nodata = src.nodata
    out = pd.DataFrame({"sat": vals[:, band - 1]})
    if qa_band is not None and vals.shape[1] >= qa_band:
        out["sat_qa"] = vals[:, qa_band - 1]
    if nodata is not None:
        out.loc[out["sat"] == nodata, "sat"] = np.nan
    log.info("sampled satellite at %d points, %d missing", len(out), int(out["sat"].isna().sum()))
    return out


def fit_hybrid(sites: gpd.GeoDataFrame, y: np.ndarray, local_predictors: pd.DataFrame,
               sat: pd.Series, select_fn) -> dict:
    """Joint fit. The satellite term is forced in; local terms are selected around it."""
    cand = local_predictors[[c for c in local_predictors.columns
                             if _buffer_of(c) is not None
                             and _buffer_of(c) <= LOCAL_MAX_BUFFER_M]].copy()
    cand.insert(0, "sat", sat.to_numpy(float))
    ok = np.isfinite(cand["sat"].to_numpy()) & np.isfinite(y)
    if (~ok).any():
        log.warning("dropping %d sites with no satellite value", int((~ok).sum()))
    cand, yy = cand.loc[ok], y[ok]

    # Fit the regional term alone first: this is the variance the satellite explains.
    sat_only = LinearRegression().fit(cand[["sat"]], yy)
    r2_sat = sat_only.score(cand[["sat"]], yy)

    # Then select local terms on the RESIDUAL of the regional fit, so a large-buffer
    # urbanity proxy cannot claim signal the satellite has already accounted for.
    resid = yy - sat_only.predict(cand[["sat"]])
    local_terms = select_fn(cand.drop(columns=["sat"]), resid)

    terms = ["sat"] + list(local_terms)
    joint = LinearRegression().fit(cand[terms], yy)
    out = {"terms": terms,
           "coefficients": dict(zip(terms, map(float, joint.coef_))),
           "intercept": float(joint.intercept_),
           "r2_satellite_only": float(r2_sat),
           "r2_joint": float(joint.score(cand[terms], yy)),
           "n": int(len(yy))}
    log.info("hybrid n=%d r2_sat=%.3f r2_joint=%.3f terms=%s",
             out["n"], out["r2_satellite_only"], out["r2_joint"], terms)
    return out


def _buffer_of(name: str) -> int | None:
    tail = name.rsplit("_", 1)[-1]
    return int(tail) if tail.isdigit() else None

Selecting the local terms on the residual of the satellite-only fit is what prevents double counting. A 5 km population-density predictor would otherwise be chosen first, explaining the same urban–rural contrast the satellite already carries, and the joint model would then hold two large coefficients of opposite sign whose difference happens to fit the campaign. Restricting the local candidates to buffers no larger than the satellite cell size is the second, cruder guard on the same problem.

The reported pair of fits is the useful summary. If the satellite alone reaches 0.55 and the joint model reaches 0.58, the local component is adding almost nothing and the hybrid is not worth its complexity. If the satellite reaches 0.30 and the joint model 0.68, the two components are doing genuinely different work, which is the case the method exists for.

Where Each Model Wins, and Where the Hybrid Matters Held-out prediction error plotted against distance from the metropolitan campaign area. The land-use model has the lowest error inside the campaign area and degrades sharply beyond it, because its predictors move outside the sampled range. The satellite-only model has a higher but almost flat error at all distances, because the regional field is measured everywhere. The hybrid tracks the land-use model inside the campaign area and converges to the satellite model far from it, so it is never worse than the better of the two at any distance. A shaded band marks the campaign extent. The hybrid is never worse than the better component at any distance campaign extent land use alone satellite alone hybrid 0 100 km 250 km distance from the campaign area held-out error

Validation & Edge Cases

Satellite gaps that are not random. Retrievals fail under cloud, over bright surfaces and at high viewing angles, and cloud is correlated with season and terrain. A gap-filled product hides that with an imputation whose uncertainty is rarely propagated. Keep the quality band, and report what fraction of each area’s value came from filled cells.

Resolution claims. The hybrid inherits the coarser resolution for its regional component. Publish it at the land-use resolution if you like, but state in the metadata that variation beyond the satellite cell size comes entirely from the land-use term.

Validation must span the distance axis. A hybrid validated only at campaign sites reports the land-use model’s strength and tells you nothing about the case the hybrid was built for. Hold out distant sites explicitly, using the blocked design from Cross-Validating an Exposure Surface With Held-Out Monitors.

Temporal mismatch between components. A satellite composite for a calendar year combined with a campaign run over three summer months describes two different exposures. Align the windows or restrict the hybrid to the overlap, and record which was done.

Sites with no satellite value. Dropping them silently biases the fit toward cloud-free regions. Count them, report them, and check whether they cluster.

INFO lur.hybrid sampled satellite at 58 points, 3 missing
WARNING lur.hybrid dropping 3 sites with no satellite value
INFO lur.hybrid hybrid n=55 r2_sat=0.312 r2_joint=0.681 terms=['sat', 'road_major_100', 'industrial_300']

The residual-selection order is the whole difference between a stable hybrid and an unstable one, and it is visible in the coefficients each ordering produces:

Offsetting Coefficients Are the Signature of Double Counting Two sets of standardised coefficients. In the naive fit, where a five-kilometre population density term competes with the satellite term for the same regional signal, the satellite coefficient is large and positive while the population term is large and negative, two big numbers cancelling to fit the campaign. In the residual-selection fit, the satellite term carries the regional signal alone, and the local terms selected on its residual are small, positive and physically interpretable. Two large coefficients of opposite sign mean one signal counted twice naive joint fit residual selection satellite popdens 5 km road 100 m satellite road 100 m industrial cancelling terms, unstable small, positive, interpretable standardised coefficients; bars to the left of each axis are negative

Compliance Notes

A hybrid model’s provenance must name both parents and the rule that joined them. Record the satellite product identifier, version, compositing window and gap-filling method with the quality-band summary; the land-use candidate set and the buffer ceiling applied; the residual-selection order that prevented double counting; the joint coefficients; the satellite-only and joint fits; and held-out error reported separately inside and outside the campaign area.

Two statements belong in any methods text. That the local component was selected on the regional residual, which is what makes the coefficients interpretable rather than offsetting. And the effective resolution of each component, so a reader does not attribute street-scale meaning to a term that varies over a kilometre.