⚙️ How-To Hydrology & watersheds Terrain & elevation Joins & relationships Hydrology & Watersheds · step 11 of 20 Updated

How to Summarise Catchment Attributes: Area, Slope and Land Cover

Problem statement

A delineated catchment is rarely the end product. Hydrological models, regional regressions and reports need its attributes: area, elevation range, mean slope, land cover shares, drainage density, shape. Each is one line of code, and each has a way of being quietly wrong: an area computed in the wrong projection, a slope averaged in the wrong units, a land cover share counted on a grid whose cells are not equal in area.

Measured on the Esopus Creek catchment above the Coldbrook gauge, New York (USGS 01362500, published 497.3 km²), delineated from the 10 m 3DEP DEM:

  • Area was 492.42 km² in UTM zone 18N, 492.77 km² in Albers equal-area and 492.77 km² geodesic — and 895.49 km² in Web Mercator, 82% too large. In EPSG:4326, .area returned 0.05362, in square degrees.
  • Mean elevation was 598.37 m, relief 1,084.8 m; exactextract and rasterstats agreed to the centimetre, in 0.26 s and 0.22 s.
  • Mean slope was 17.03°. Averaging slope in percent and converting afterwards gave 17.58°; NumPy's gradient on the geographic DEM, with degrees treated as metres, gave 89.91°.
  • Tree cover was 98.24% of ESA WorldCover, and three methods agreed to 0.001 percentage points.

Quick answer 

from exactextract import exact_extract

area_km2 = basin.to_crs(5070).area.iloc[0] / 1e6                       # an equal-area CRS, never EPSG:4326 or 3857
elevation = exact_extract("dem_utm.tif", basin, ["mean", "min", "max"], output="pandas")
slope = exact_extract("slope_degrees.tif", basin, ["mean", "median"], output="pandas")
cover = exact_extract("worldcover.tif", basin.to_crs(4326), ["unique", "frac"], output="pandas")

Compute area in an equal-area or local projected CRS, slope from a projected DEM in degrees, and categorical shares as area fractions of the polygon.

Bar chart of the Esopus catchment area computed in UTM 18N, Albers equal-area, geodesically, and in Web Mercator, against the published 497.3 km².
Three sensible methods agree within 0.1%; Web Mercator adds 403 km² of catchment that does not exist.

Step-by-step solution

1. Start from one polygon in a projected CRS

Delineate the catchment, polygonise it and dissolve it to a single geometry in the DEM's projected CRS; see delineating a watershed. Here: 4,924,175 cells of 10 m, from a filled DEM and a gauge snapped within 150 m.

2. Compute area in an equal-area CRS or geodesically

UTM 18N gave 492.42 km², within 0.07% of the geodesic 492.77 km², because the basin sits close to its zone's central meridian. Albers EPSG:5070 matched the geodesic value exactly. Web Mercator gave 895.49 km², since its scale at 42° N is 1.35 in each direction. The published 497.3 km² differs by a further 0.9% because it came from a different delineation. See choosing a projected CRS for area.

3. Summarise elevation with exact cell fractions

exact_extract weights each cell by the fraction of it inside the polygon; zonal_stats uses cells whose centres fall inside. For a basin of 4.9 million cells the two agreed exactly: mean 598.37 m, minimum 190.31 m, maximum 1,275.10 m. For small catchments a few hundred cells across, the edge cells matter much more.

4. Derive relief and the hypsometric integral

Relief is maximum minus minimum: 1,084.8 m. The hypsometric integral, (mean − min) / (max − min), was 0.376: most of the catchment lies in the lower half of its elevation range, typical of an eroded mountain valley.

5. Calculate slope on a projected DEM, in degrees

WhiteboxTools' slope on the 10 m UTM DEM gave a mean of 17.03° and median 16.40°. On the 30 m GLO-30 DEM, the mean fell to 16.71°, because coarser cells smooth steep ground. On the geographic 3DEP DEM, WhiteboxTools converted degrees to metres itself and gave 17.04°. A NumPy gradient on the same geographic grid, without that conversion, gave 89.91°: every slope near vertical. See why slope values come out wrong.

6. Average slope in the units you report

Slope in percent is not linear in angle. The mean of percent slope was 31.68%, which converts to 17.58° — 0.55° more than the mean of the degree values. Choose the unit first, then average.

7. Take land cover shares as area fractions

ESA WorldCover is on a geographic grid, where cells shrink towards the poles. Across a basin 0.2° tall the difference was negligible: count shares and cosine-weighted shares agreed to 0.001 percentage points. Tree cover 98.24%, grassland 1.28%, built-up 0.23%, water 0.09%. For basins spanning several degrees of latitude, weight by cell area.

8. Add network and shape attributes

The NHDPlus HR network inside the basin was 800.1 km long, a drainage density of 1.62 km/km². The perimeter of the raw polygon, following every cell edge, was 152.2 km; simplified at 30 m it was 117.7 km. Shape indices built on perimeter depend on that choice: the compactness coefficient was 1.93 on the staircase and 1.50 on the simplified outline.

Bar chart of the Esopus catchment's mean slope computed five ways: 10 m UTM degrees, 30 m GLO-30 degrees, geographic 3DEP with WhiteboxTools, the mean of percent slope converted to degrees, and a NumPy gradient on a degrees grid.
Resolution moved the mean by 0.3°, averaging percent by 0.55°, and ignoring the degree grid by 73°.

Code examples

Example 1 — delineate and measure area four ways 

import os

import geopandas as gpd
import rasterio
import whitebox
from pyproj import Geod, Transformer
from rasterio.features import shapes
from shapely.geometry import shape

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))
wbt.fill_depressions("esopus_3dep13_utm.tif", "filled.tif", fix_flats=True)
wbt.d8_pointer("filled.tif", "d8.tif")
wbt.d8_flow_accumulation("d8.tif", "acc.tif", pntr=True)
x, y = Transformer.from_crs("EPSG:4269", "EPSG:32618", always_xy=True).transform(-74.2701944, 42.0144722)
gpd.GeoDataFrame(geometry=gpd.points_from_xy([x], [y]), crs=32618).to_file("gauge.shp")
wbt.snap_pour_points("gauge.shp", "acc.tif", "gauge_snapped.shp", snap_dist=150)
wbt.watershed("d8.tif", "gauge_snapped.shp", "watershed.tif")

with rasterio.open("watershed.tif") as src:
    ws = src.read(1)
    inside = (ws != src.nodata) & (ws > 0)
    parts = [shape(g) for g, v in shapes(inside.astype("uint8"), mask=inside, transform=src.transform)]
basin = gpd.GeoDataFrame(geometry=[gpd.GeoSeries(parts).union_all()], crs=32618)

geodesic = abs(Geod(ellps="WGS84").geometry_area_perimeter(basin.to_crs(4326).geometry.iloc[0])[0]) / 1e6
areas = {"UTM 18N": basin.area.iloc[0] / 1e6, "Albers EPSG:5070": basin.to_crs(5070).area.iloc[0] / 1e6,
         "geodesic": geodesic, "Web Mercator": basin.to_crs(3857).area.iloc[0] / 1e6}
print(f"{int(inside.sum()):,} cells")
for name, km2 in areas.items():
    print(f"{name:17} {km2:7.2f} km2 ({km2 / geodesic - 1:+.2%})")
print(f"EPSG:4326 .area   {basin.to_crs(4326).area.iloc[0]:.5f} (square degrees)")
4,924,175 cells
UTM 18N            492.42 km2 (-0.07%)
Albers EPSG:5070   492.77 km2 (+0.00%)
geodesic           492.77 km2 (+0.00%)
Web Mercator       895.49 km2 (+81.72%)
EPSG:4326 .area   0.05362 (square degrees)

Example 2 — elevation, relief and slope 

import numpy as np
from exactextract import exact_extract

elevation = exact_extract("esopus_3dep13_utm.tif", basin, ["mean", "min", "max"], output="pandas").iloc[0]
relief = elevation["max"] - elevation["min"]
print(f"elevation mean {elevation['mean']:.2f} m, range {elevation['min']:.2f}-{elevation['max']:.2f} m, "
      f"relief {relief:.1f} m, hypsometric integral {(elevation['mean'] - elevation['min']) / relief:.3f}")

for dem, crs_basin in (("esopus_3dep13_utm.tif", basin), ("esopus_glo30_utm.tif", basin)):
    out = dem.replace(".tif", "_slope.tif")
    wbt.slope(dem, out, units="degrees")
    slope = exact_extract(out, crs_basin, ["mean", "median"], output="pandas").iloc[0]
    print(f"{dem}: slope mean {slope['mean']:.2f} deg, median {slope['median']:.2f} deg")

wbt.slope("esopus_3dep13_utm.tif", "slope_percent.tif", units="percent")
percent = exact_extract("slope_percent.tif", basin, ["mean"], output="pandas")["mean"].iloc[0]
print(f"mean of percent slope {percent:.2f}% = {np.degrees(np.arctan(percent / 100)):.2f} deg")
elevation mean 598.37 m, range 190.31-1275.10 m, relief 1084.8 m, hypsometric integral 0.376
esopus_3dep13_utm.tif: slope mean 17.03 deg, median 16.40 deg
esopus_glo30_utm.tif: slope mean 16.71 deg, median 16.28 deg
mean of percent slope 31.68% = 17.58 deg

Example 3 — land cover, drainage density and shape 

import pandas as pd

WORLDCOVER = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up",
              60: "bare or sparse", 80: "water", 90: "herbaceous wetland"}
cover = exact_extract("esopus_worldcover2021.tif", basin.to_crs(4326), ["unique", "frac"], output="pandas").iloc[0]
shares = pd.Series({WORLDCOVER[int(k)]: 100 * f for k, f in zip(cover["unique"], cover["frac"])}).sort_values(ascending=False)
print(shares.round(3).to_string())

flowlines = gpd.read_file("nhdplushr_flowlines_esopus.gpkg").to_crs(32618)
network_km = gpd.clip(flowlines[flowlines.ftype.isin([460, 558, 334])], basin).length.sum() / 1000
area_km2 = areas["Albers EPSG:5070"]
attributes = basin.assign(
    area_km2=area_km2, elev_mean_m=elevation["mean"], relief_m=relief,
    tree_cover_pct=shares["tree cover"], network_km=network_km, drainage_density=network_km / area_km2,
    perimeter_km=basin.simplify(30).length / 1000,
)
attributes["compactness"] = attributes.perimeter_km / (2 * np.sqrt(np.pi * area_km2))
attributes.to_file("esopus_attributes.gpkg")
print(attributes.drop(columns="geometry").round(3).T.to_string(header=False))
tree cover            98.241
grassland              1.278
built-up               0.235
water                  0.091
herbaceous wetland     0.079
cropland               0.062
bare or sparse         0.014
area_km2           492.774
elev_mean_m        598.375
relief_m          1084.794
tree_cover_pct      98.241
network_km         800.111
drainage_density     1.624
perimeter_km       117.721
compactness          1.496

Explanation

Why projections change area and not much else

Every map projection distorts either area, shape or both. UTM is conformal but its scale error within a zone is below 0.1% near the central meridian; Albers is designed for equal area; a geodesic calculation works on the ellipsoid directly. Web Mercator preserves shape at the cost of area, and at 42° N its area scale is 1 / cos²(42°) = 1.81. The 0.05362 from EPSG:4326 is not an area in any useful unit.

Why exact fractions matter less on big catchments

Only boundary cells differ between centre-in and fractional coverage. The Esopus boundary crosses roughly 15,000 of its 4.9 million cells — 0.3% — and the boundary cells are no different from the interior on average. For a 5 km² sub-catchment of 50,000 cells, the share of boundary cells is ten times higher.

Why averaging order matters

The mean of a non-linear transform is not the transform of the mean. Percent slope is tan(angle) × 100, which grows faster than the angle as slopes steepen, so the steep cells pull the mean of percent up more than they pull the mean of degrees. Averaging percent and converting gave a steeper catchment than averaging degrees directly.

Why the geographic gradient failed

The DEM's cells are 1/3 arc-second, about 7.6 m by 10.3 m at this latitude, but a gradient in "cell units" of degrees treats 9.26 × 10⁻⁵ as the horizontal distance. A 1 m rise over that distance is a gradient of about 10,800, and the arctangent is 89.99°. WhiteboxTools detects geographic coordinates and converts them.

Table of the Esopus Creek catchment attributes: area, elevation, relief, hypsometric integral, slope, tree cover, drainage density and compactness, with the method used for each.
Each value with the method that produced it: the method is part of the attribute.

Edge cases or notes

  • Catchments crossing UTM zones should use an equal-area CRS or geodesic area.
  • NoData inside the polygon lowers counts silently; compare the cell count with the area.
  • Land cover on projected grids needs no latitude weighting; on geographic grids spanning many degrees it does.
  • Categorical rasters resampled bilinearly produce codes that do not exist; reproject them with nearest neighbour.
  • Main-stem length and elongation need a longest flow path, not the perimeter.
  • Very small catchments should use exact fractions; centre-in counts can miss thin polygons entirely.
  • Rasters in a different CRS from the polygon need the polygon reprojected, not the raster.

FAQ

How do I calculate catchment area in Python?

Reproject the catchment polygon to an equal-area CRS such as EPSG:5070, or compute a geodesic area with pyproj. The Esopus catchment was 492.77 km² both ways; Web Mercator gave 895.49 km².

Should I use exactextract or rasterstats for catchment statistics?

exactextract weights boundary cells by coverage, which matters for small catchments. For the 4.9-million-cell Esopus basin, both gave a mean elevation of 598.37 m.

How do I calculate mean catchment slope?

Calculate slope in degrees on a projected DEM, then average it over the polygon. The Esopus mean was 17.03° at 10 m and 16.71° at 30 m.

Why is my mean slope close to 90 degrees?

The DEM is in degrees and the gradient treated them as metres. A NumPy gradient on the geographic DEM gave 89.91°; use a projected DEM or a tool that converts units.

How do I calculate land cover percentages for a catchment?

Use exact_extract with the unique and frac operations on the categorical raster. The Esopus basin was 98.24% tree cover in ESA WorldCover 2021.

What is the hypsometric integral?

It is (mean elevation − minimum) / (maximum − minimum), a measure of how much of a catchment lies high in its range. The Esopus value was 0.376.