Computing Composite Curve Numbers from Land-Use Rasters in Python

A curve number is only as good as the spatial data behind it, and in practice that data arrives as two rasters that were never designed to be used together: a land-cover classification and a soil survey. This guide builds a single composite curve number per subbasin by joining a land-cover raster to a hydrologic soil group raster, mapping each pairing to a tabulated value, and area-weighting the result. It is the spatial-data companion to SCS curve number runoff estimation, the parent guide that covers the runoff equation itself, and it sits within the broader Rainfall-Runoff Modeling & Hydrologic Simulation domain. The worked example uses the National Land Cover Database (NLCD) for land cover and gridded gSSURGO for soils, but the pattern generalizes to any categorical land-cover and soil pair.

The end product is a number a modeler can defend: an area-weighted curve number for every hydrologic unit, traceable back to the specific land-cover class and soil group that produced it. Getting there means solving three concrete problems — aligning the two rasters onto one grid, encoding a two-dimensional lookup that survives missing and dual soil groups, and aggregating without letting nodata contaminate the average.

Inputs and Reference Values

Two categorical rasters and one polygon layer drive the workflow:

• Land cover — NLCD, a 30 m classification where each cell holds an integer class code (for example 82 for cultivated crops, 22 for low-intensity developed).
• Hydrologic soil group — gSSURGO, rasterized so each cell holds a group code A through D (or a dual group such as A/D).
• Subbasins — a polygon layer of the hydrologic units over which curve numbers are averaged, obtained from a delineation workflow.

The hydrologic soil group controls how much the same land cover runs off. The table below summarizes the groups and gives representative AMC II curve numbers for a single land cover — row crops in good hydrologic condition — to show the spread across soils.

Open water and fully impervious surfaces sit near 98 regardless of the underlying soil, and forest in good condition on group A soil can fall to the low 30s. Those extremes are what make the area-weighted average sensitive to how faithfully the lookup and the grid resolution preserve small, high-runoff features.

Step 1 — Align the Rasters onto One Grid

The land-cover and soil rasters almost never share a grid: NLCD ships in an Albers equal-area projection at 30 m, while a rasterized gSSURGO tile may be in a UTM zone at a different resolution and origin. Because the curve number lookup is applied cell by cell, the two rasters must be reprojected and resampled onto one identical grid — same CRS, resolution, transform, and shape — before anything else. Both are categorical, so resampling must use nearest-neighbour; any averaging or bilinear method would invent class codes that do not exist. A mismatch here is the most common source of a wrong composite curve number, and the general procedure for reconciling projections is covered in fixing CRS mismatches in watershed shapefiles.

import logging
import numpy as np
import rasterio
from rasterio.warp import reproject, Resampling

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


def align_to_reference(src_path: str, ref_path: str) -> np.ndarray:
    """Reproject a categorical raster onto the reference grid using nearest-neighbour."""
    with rasterio.open(ref_path) as ref:
        ref_profile = ref.profile
        ref_transform = ref.transform
        ref_crs = ref.crs
        dst = np.zeros((ref.height, ref.width), dtype="int32")

    with rasterio.open(src_path) as src:
        logger.info("Reprojecting %s (%s) onto reference grid (%s)",
                    src_path, src.crs, ref_crs)
        reproject(
            source=rasterio.band(src, 1),
            destination=dst,
            src_transform=src.transform,
            src_crs=src.crs,
            dst_transform=ref_transform,
            dst_crs=ref_crs,
            resampling=Resampling.nearest,  # categorical data: never interpolate
        )
    logger.info("Aligned raster classes present: %s", np.unique(dst)[:15])
    return dst

Use the land-cover raster as the reference grid so the composite curve number inherits the finest categorical resolution available for land use, which is where most of the runoff variability lives.

Step 2 — Encode the (Land Cover, HSG) Lookup

The curve number is a function of two categorical inputs, so the natural structure is a dictionary keyed by the pair. Keep the lookup explicit and version-controlled; it is the single most scrutinized artifact in a curve number study.

# AMC II curve numbers keyed by (nlcd_class, hsg_code). Abbreviated for illustration;
# a production table covers every NLCD class present in the study area.
CN_LOOKUP = {
    (82, "A"): 67, (82, "B"): 78, (82, "C"): 85, (82, "D"): 89,  # cultivated crops
    (81, "A"): 49, (81, "B"): 69, (81, "C"): 79, (81, "D"): 84,  # pasture/hay
    (41, "A"): 30, (41, "B"): 55, (41, "C"): 70, (41, "D"): 77,  # deciduous forest
    (71, "A"): 30, (71, "B"): 58, (71, "C"): 71, (71, "D"): 78,  # grassland
    (21, "A"): 49, (21, "B"): 69, (21, "C"): 79, (21, "D"): 84,  # developed, open
    (22, "A"): 61, (22, "B"): 75, (22, "C"): 83, (22, "D"): 87,  # developed, low
    (23, "A"): 77, (23, "B"): 85, (23, "C"): 90, (23, "D"): 92,  # developed, medium
    (24, "A"): 89, (24, "B"): 92, (24, "C"): 94, (24, "D"): 95,  # developed, high
    (11, "A"): 98, (11, "B"): 98, (11, "C"): 98, (11, "D"): 98,  # open water
}

HSG_CODES = {1: "A", 2: "B", 3: "C", 4: "D"}  # integer soil raster -> group letter
CN_NODATA = -9999

Two design choices matter here. First, water and high-intensity development carry curve numbers near the ceiling for every soil group, so they must never be diluted by an area-weighted average that treats them as ordinary land. Second, dual hydrologic soil groups (A/D, B/D, C/D) reported for high-water-table soils force a decision: use the drained member where artificial drainage exists, otherwise the wetter undrained member. Encode that decision as its own key rather than collapsing it silently.

Step 3 — Map to a Curve Number Raster

With aligned rasters and a lookup, build the curve number grid. The efficient pattern combines the two categorical rasters into a single integer key per cell, then applies the lookup over the unique keys rather than looping over every pixel.

def build_cn_raster(
    landcover: np.ndarray,
    hsg: np.ndarray,
    lookup: dict,
    hsg_codes: dict,
    nodata: int = CN_NODATA,
) -> np.ndarray:
    """Combine land-cover and HSG rasters into a curve number grid via a paired lookup."""
    if landcover.shape != hsg.shape:
        raise ValueError(
            f"Shape mismatch: landcover {landcover.shape} vs hsg {hsg.shape}. "
            "Align both rasters onto one grid before mapping."
        )

    cn = np.full(landcover.shape, nodata, dtype="int32")
    # Iterate over the distinct (landcover, hsg) combinations actually present
    lc_codes = np.unique(landcover)
    for lc in lc_codes:
        for hsg_int, hsg_letter in hsg_codes.items():
            key = (int(lc), hsg_letter)
            cn_value = lookup.get(key)
            if cn_value is None:
                continue  # unmapped pairing stays nodata for later auditing
            mask = (landcover == lc) & (hsg == hsg_int)
            cn[mask] = cn_value

    unmapped = int(np.sum((cn == nodata) & (landcover > 0)))
    if unmapped > 0:
        logger.warning(
            "%d cells have land cover but no CN (unmapped pairing or missing HSG). "
            "Assign a default group or mask before area-weighting.", unmapped
        )
    logger.info("CN raster built: min=%d max=%d over %d valid cells",
                int(cn[cn != nodata].min()), int(cn[cn != nodata].max()),
                int(np.sum(cn != nodata)))
    return cn

The unmapped count is the workflow’s early-warning signal. A large value means either the lookup is missing a land-cover class present in the study area or the soil raster has holes where land cover does not — both must be resolved before the average is trustworthy.

Step 4 — Area-Weight per Subbasin

A curve number raster is rarely the deliverable; a single representative value per subbasin is. Because every cell in a projected equal-area grid covers the same area, the area-weighted mean reduces to a simple masked mean of the CN raster within each polygon, which rasterstats computes directly.

import rasterio
from rasterstats import zonal_stats


def composite_cn_per_subbasin(
    cn_raster_path: str,
    subbasins_path: str,
    nodata: int = CN_NODATA,
) -> list:
    """Area-weighted mean curve number per subbasin polygon."""
    logger.info("Computing zonal composite CN for %s", subbasins_path)
    stats = zonal_stats(
        subbasins_path,
        cn_raster_path,
        stats=["mean", "count"],
        nodata=nodata,          # nodata cells are excluded from the mean
        geojson_out=True,
    )
    for feature in stats:
        props = feature["properties"]
        props["composite_cn"] = round(props["mean"], 1) if props["mean"] else None
        logger.info("Subbasin %s: composite CN = %s over %s cells",
                    props.get("id", "?"), props["composite_cn"], props["count"])
    return stats


if __name__ == "__main__":
    lc = align_to_reference("nlcd_2021.tif", "nlcd_2021.tif")  # reference = itself
    soils = align_to_reference("gssurgo_hsg.tif", "nlcd_2021.tif")
    cn_grid = build_cn_raster(lc, soils, CN_LOOKUP, HSG_CODES)

    with rasterio.open("nlcd_2021.tif") as ref:
        profile = ref.profile
    profile.update(dtype="int32", nodata=CN_NODATA, count=1)
    with rasterio.open("curve_number.tif", "w", **profile) as dst:
        dst.write(cn_grid.astype("int32"), 1)

    composite_cn_per_subbasin("curve_number.tif", "subbasins.gpkg")

Passing nodata=CN_NODATA to zonal_stats is what keeps unmapped and out-of-basin cells out of the average. Omit it and every nodata cell is read as its raw sentinel value, which either crashes the mean or, worse, drags it toward an implausible number. The count field reports how many valid cells backed each subbasin’s average — a low count relative to the polygon area is a flag that too much of the basin fell to nodata.

Validating the Composite Curve Number

A composite curve number is a single scalar per subbasin, which makes it easy to accept without scrutiny and easy to get quietly wrong. Three checks catch the large majority of errors before the value propagates into a runoff estimate.

Range and distribution. Every subbasin composite should land between the minimum and maximum curve number present in its land cover, typically somewhere in the 55 to 95 band for mixed watersheds. A composite below every value in the lookup is the classic symptom of nodata leaking into the mean; a composite pinned near 98 means water or high-intensity development dominates the basin, which should be corroborated by the land-cover fractions. Histogram the CN raster and confirm the spike at 98 corresponds to genuinely mapped water and pavement rather than an unmapped class that fell back to a ceiling default.

Areal fractions. Tabulate the fraction of each subbasin occupied by each land-cover class and each soil group, and confirm those fractions match the source data independently. If cultivated crops cover 60 percent of a basin in the NLCD but only 20 percent of the cells feeding the curve number average, the alignment or the lookup dropped a class. This cross-check is the most reliable way to detect a silently missing lookup entry, because the affected cells simply vanish from the count rather than raising an error.

Sensitivity to the aggregation order. Recall that averaging the curve numbers is only exact for a linear response. For a subbasin that mixes very high and very low curve numbers, compute a trial runoff depth two ways for a representative design storm — once from the area-weighted composite curve number, and once by computing runoff per cell and then averaging the depths. If the two disagree by more than a few percent, the basin is heterogeneous enough that per-cell runoff aggregation is the more defensible path, and the composite value should be flagged as an approximation. The runoff arithmetic behind both approaches is covered in the parent guide on the curve number method.

Document the resulting composite curve numbers alongside the land-cover and soil vintages that produced them. A curve number derived from 2021 land cover and a 2020 soil survey is a dated snapshot, and recording the provenance is what lets a later analyst reproduce or revise the value rather than re-deriving it from scratch.

Gotchas and Edge Cases

• Nodata masquerading as data. If the CN sentinel falls inside the plausible 30–100 range, or if zonal_stats is not told the nodata value, unmapped cells silently enter the average and pull the composite curve number toward a wrong value. Use a sentinel well outside the valid range, such as -9999, and pass it to every downstream call.

• Missing soil groups. gSSURGO leaves the hydrologic soil group blank for open water and some mapped urban land. Decide a policy up front — default to group D for a conservative estimate, borrow the neighbourhood mode, or mask the cell — and log how many cells the policy touched. Silent fallback to a default is the failure that most often survives review undetected.

• Dual groups. A/D, B/D, and C/D soils behave as the drier member only when artificially drained. For undrained conditions use the wetter D member. Splitting them requires a drainage layer or a documented assumption, not a coin flip.

• Categorical resampling. Reprojecting either raster with bilinear or cubic resampling invents class codes between the real ones (a cell halfway between class 41 and 42 becomes a nonexistent 41.5 that rounds unpredictably). Always resample categorical rasters with nearest-neighbour.

• Impervious dilution. A 30 m grid smears small connected impervious surfaces across mostly pervious cells, lowering the composite curve number below the true runoff response. Where connected impervious area matters, model it as a separate class or move to a finer land-cover product before averaging.

• CRS drift between the raster and the subbasins. The polygon layer must share the CRS of the CN raster, or rasterstats will read cells from the wrong location. Reproject the subbasins to match the raster before the zonal call.

Related Topics

• SCS Curve Number Runoff Estimation — the parent guide covering the runoff equation, retention, and antecedent moisture that consume the composite curve number produced here
• Estimating Time of Concentration with the NRCS Velocity Method — the sibling guide deriving the basin timing parameter that pairs with the curve number in a hydrograph model
• Rainfall-Runoff Modeling & Hydrologic Simulation — the domain overview linking loss methods, transforms, and calibration