Reading Cloud-Optimized GeoTIFFs for Impervious Surface Checks

Impervious cover limits are enforced as a percentage of lot area, which makes them a raster measurement against a vector boundary — and the naive implementation downloads a multi-gigabyte land cover raster to measure a quarter-acre lot. A cloud-optimized GeoTIFF is internally tiled and carries overviews, so a windowed read fetches only the tiles a parcel touches. This guide measures impervious cover per parcel from a remote COG, weights the boundary cells by their actual overlap rather than counting them in or out, and explains why the overviews that make the map draw quickly must never produce the number in a determination. It is the raster half of cloud-native geospatial formats for compliance pipelines.

Prerequisites

Step-by-step

Step 1: Confirm the file is genuinely cloud-optimized

A GeoTIFF served over HTTP is not a COG. Without internal tiling and overviews, every windowed read still transfers the whole file, and the pipeline appears to work while being orders of magnitude slower than it should be.

import rasterio

def cog_profile(url: str) -> dict:
    with rasterio.open(url) as src:
        return {
            "tiled": src.profile.get("tiled", False),
            "block": (src.profile.get("blockxsize"), src.profile.get("blockysize")),
            "overviews": src.overviews(1),
            "compression": str(src.compression),
            "crs": src.crs.to_string(),
            "res": src.res,
            "nodata": src.nodata,
            "dtype": src.dtypes[0],
        }

tiled: False means stripped layout, and a strip spans the full raster width — so reading a parcel-sized window pulls the entire width of the image for every row it touches. That single flag is the difference between a fast pipeline and a slow one.

Why the layout decides whether a windowed read is a readInternal tiling and overviews compared against a plain GeoTIFF served over HTTP, for a parcel-sized window.Internally tiledStripped layoutBytes fetched for aparcelA few blocksFull image width, every row touchedBatch over a countyNeighbouring parcels share blocksRe-reads the same stripsOverviewsPresent, for display onlyUsually absentSymptom when wrong—Correct answers, unreasonable runtime
A GeoTIFF on HTTP is not a COG. Without internal tiling every windowed read transfers far more than it needs and the pipeline is quietly slow.

Step 2: Reproject the parcels to the raster, not the raster to the parcels

Resampling a raster to match a vector layer’s CRS resamples every cell, changes the values at every boundary, and produces a measurement against data that no longer matches its source. Transforming the parcel geometry is exact and free.

def parcels_in_raster_crs(parcels, url):
    with rasterio.open(url) as src:
        return parcels.to_crs(src.crs)

There is a real trade-off here that has to be stated rather than ignored: areas measured in the raster’s CRS may differ slightly from areas in the project’s working CRS. The resolution is to take the fraction from the raster and apply it to the lot area computed in the working CRS, which is what step 4 does.

Step 3: Read only the window each parcel needs

The windowed read is what makes remote access viable. Ask for the parcel’s bounding box, masked to the parcel itself.

from rasterio.mask import mask

def parcel_window(url, geom, all_touched=False):
    """Cells within the parcel, read as a single windowed request."""
    with rasterio.open(url) as src:
        arr, transform = mask(src, [geom], crop=True, all_touched=all_touched,
                              filled=True, nodata=src.nodata)
        return arr[0], transform, src.nodata

For a batch, do not open the file once per parcel — group parcels by tile so that neighbouring parcels share the fetched blocks. GDAL’s block cache does this automatically within a session, so the practical rule is to sort parcels spatially before iterating, which the Hilbert ordering from converting parcel shapefiles to GeoParquet already provides.

Step 4: Weight the boundary cells by their actual overlap

A 1-metre cell straddling a parcel boundary is neither wholly in nor wholly out. Counting it in (all_touched=True) overstates the parcel; counting it out understates it — and for a quarter-acre lot the boundary cells are a substantial share of the total.

from exactextract import exact_extract

IMPERVIOUS_CLASSES = {21, 22, 23, 24}     # NLCD developed classes, as an example

def impervious_fraction(url, parcels, classes=IMPERVIOUS_CLASSES):
    """Area-weighted class fractions: boundary cells count for the part inside."""
    res = exact_extract(url, parcels, ["unique", "frac"], output="pandas")
    out = []
    for _, row in res.iterrows():
        vals, fracs = row["unique"], row["frac"]
        imp = sum(f for v, f in zip(vals, fracs) if int(v) in classes)
        out.append(imp)
    return out

The magnitude of the difference is worth internalising. On a 1-metre raster and a typical 20 by 40 metre suburban lot, about 120 of roughly 800 cells touch the boundary — so the gap between counting them all in and all out is around 15 percentage points of impervious cover, against a limit typically stated to the nearest whole percent. Exact area weighting is not a refinement here; it is the difference between a defensible number and an arbitrary one.

Boundary cells are not a rounding detailFor a 20 by 40 metre lot on a 1-metre raster, the share of cells that straddle the parcel boundary — and therefore the spread between counting them in and out.Cells wholly inside the parcel680 cellsCells straddling the boundary120 cellsAbout 15% of the parcel by cell count. all_touched=True and all_touched=False bracket the true value by roughly that margin.
Against a limit stated to the nearest whole percent, exact area weighting is the difference between a defensible number and an arbitrary one.

Step 5: Never measure from an overview

Overviews are decimated copies used for display. Reading class values from one does not average the classes — it picks a representative cell — so the impervious fraction it yields is a sample, not a measurement, and it changes with the zoom level requested.

def assert_full_resolution(url, requested_res):
    with rasterio.open(url) as src:
        if abs(requested_res[0] - src.res[0]) > 1e-9:
            raise ValueError(
                f"measurement requested at {requested_res} but the source is {src.res}; "
                "overviews are for display, not for determination")

This is the mistake most likely to survive review, because an overview-derived number looks entirely reasonable — it is in the right range, it varies sensibly between parcels, and it is simply wrong by a few points in an unpredictable direction. Assert the resolution at the point of measurement and the class of error disappears.

Step 6: Map the raster’s classes to the ordinance’s definition, explicitly

The raster’s idea of impervious and the code’s are different, and the mapping is a regulatory decision. A gravel driveway, a deck, a swimming pool and a green roof are each treated differently by different ordinances, and none of those distinctions exists in a land cover product.

CLASS_MAP = {
    21: {"name": "developed, open space", "impervious_share": 0.10},
    22: {"name": "developed, low intensity", "impervious_share": 0.35},
    23: {"name": "developed, medium intensity", "impervious_share": 0.65},
    24: {"name": "developed, high intensity", "impervious_share": 0.90},
}

def modelled_impervious(class_fractions: dict) -> float:
    """The product's own class-to-impervious model, stated rather than implied."""
    return sum(frac * CLASS_MAP[c]["impervious_share"]
               for c, frac in class_fractions.items() if c in CLASS_MAP)
From a land cover class to the ordinance’s definitionThe raster reports classes; the code defines impervious. The mapping between them is a stated regulatory model, not a property of the data.Read at full resolution onlyoverviews are for display, not determinationArea-weighted class fractionsboundary cells count for the part insideApply the stated class-to-impervious modela decision, recorded with the resultReport, or flag for site measurementwith the cell count that supports it
Where the threshold is close, this modelling assumption dominates the uncertainty and the correct output is a flag for site measurement.

Where the ordinance’s threshold is close to the measurement, this modelling assumption dominates the uncertainty, and the correct output is a flag for site-specific measurement rather than a verdict — the same three-valued treatment used in computing building height from lidar-derived surfaces.

Verification

Check that the read was selective, that the fractions are well-formed, and that the parcels with too few cells are excluded rather than reported.

assert cog_profile(URL)["tiled"], "source is not internally tiled — reads are not windowed"

f = results["impervious_fraction"]
assert ((f >= 0) & (f <= 1)).all(), "fraction outside [0, 1] — check nodata handling"

# Parcels too small for the raster to say anything about.
cells = results["cell_count"]
too_small = cells < 30
print(f"{too_small.sum()} parcels have fewer than 30 cells — report, do not measure")

# Reconcile against a known total for a sample area.
sample_imp = (results["impervious_fraction"] * results["lot_area"]).sum()
print(f"modelled impervious area {sample_imp:,.0f} vs published {PUBLISHED_TOTAL:,.0f}")

The cell-count check matters more than it looks. On a 30-metre land cover product, a suburban lot is about one cell, and a per-parcel impervious percentage derived from one cell is not a measurement of that parcel at all — it is the value of the cell the parcel happens to sit in.

Common Pitfalls

  • Assuming a GeoTIFF on HTTP is a COG. Without internal tiling, every windowed read transfers far more than it needs and the pipeline is quietly slow.
  • Resampling the raster to the vector CRS. It changes every cell value; transform the parcels instead and apply the fraction to the working-CRS lot area.
  • Counting boundary cells all-in or all-out. On typical lots this moves the answer by more than the precision the limit is stated to.
  • Measuring from an overview. The number looks plausible and is wrong by an amount that depends on the zoom level requested.
  • Ignoring nodata. Nodata cells inside a parcel reduce the denominator; treating them as pervious understates impervious cover systematically over water and cloud gaps.
  • Using a 30-metre product for parcel-level determinations. The resolution is not adequate to the question, whatever the arithmetic says.

Frequently Asked Questions

What resolution is adequate?

One metre or finer for parcel-level work; a 30-metre national product is suitable for screening and for aggregate reporting only. The practical test is the cell count in the verification step — a few hundred cells per parcel gives a stable fraction, and a handful does not.

Should the impervious fraction be cached?

Yes, keyed by the parcel identifier and the raster’s content hash. Impervious cover changes when the imagery is refreshed, not when the pipeline runs, so recomputing it every run is pure cost — and caching against the raster hash means a new imagery vintage invalidates the cache automatically rather than by anyone remembering.

How do I handle parcels that span two raster tiles or two rasters?

Tiles are handled transparently by the windowed read. Two separate rasters — adjacent tiles of a tiled product — need either a VRT or a mosaic read, and the thing to check is that they share a resolution, a CRS and a classification vintage. Mixing vintages across a parcel is a real hazard on multi-year products and is invisible in the output.

Is exact extraction much slower than zonal statistics?

Marginally, and it is the wrong thing to economise on. Exact area weighting costs perhaps 20–30% more than cell-centre counting and removes the largest source of error in the whole measurement.

Can the impervious fraction feed a density or coverage rule directly?

It can feed a lot coverage check directly, since that is exactly what it measures. Density rules use dwelling units rather than surface, so the connection is indirect — though impervious cover is often a useful cross-check on a density figure, since the two disagree in ways that usually indicate a stale parcel attribute. Computing dwelling unit density per acre covers the density side.

How should the measurement be recorded?

With the raster’s identity, its resolution, the class-to-impervious model applied, the cell count and the extraction method. The class model is the part that a reviewer will question, and it is the part that cannot be reconstructed from the number alone — which is exactly what capturing CRS provenance in validation logs argues for on the vector side.

Part of: Cloud-native geospatial formats for compliance pipelines