🔧 Fix Hydrology & watersheds Terrain & elevation CRS & projections Hydrology & Watersheds · step 20 of 20 Updated

Fixing Catchment Areas That Are Wrong on a Latitude–Longitude DEM

Problem statement

Many DEMs are distributed in geographic coordinates: SRTM, Copernicus GLO-30, ASTER and the USGS 3DEP seamless tiles all come in degrees. Flow routing works on them — water still runs downhill from cell to cell — but every tool that turns cells into an area has to decide how big a cell is, and a degree is not a fixed distance. The result is a catchment area in square degrees, or one computed with the wrong cell size, typically without any warning.

Measured on the Esopus Creek catchment at the Coldbrook gauge, New York (published 497.3 km²), routed on the 1/3 arc-second 3DEP DEM in NAD83 geographic coordinates with WhiteboxTools:

  • The watershed had 6,253,412 cells, and WhiteboxTools' "catchment area" output at the outlet was 0.0536 — square degrees, 99.99% below the real area.
  • Summing geodesic cell areas row by row gave 492.67 km², 0.93% below the published figure; each cell was 78.65–78.90 m², varying by 0.31% from the south to the north of the basin.
  • Using 111.32 km per degree on both axes gave 664.37 km², 35% too large; using the east–west cell width at the gauge's latitude on both axes gave 366.76 km², 26% too small.
  • A single cell area taken at the basin's middle row gave 492.62 km², 0.01% from the row-by-row answer — adequate for a catchment 0.2° tall, and increasingly wrong for taller ones.

Quick answer 

import numpy as np
from pyproj import Geod


def row_cell_areas_m2(transform, height, ellps="GRS80"):
    geod, dx = Geod(ellps=ellps), abs(transform.a)
    tops = transform.f + np.arange(height) * transform.e
    return np.array([abs(geod.polygon_area_perimeter([0, dx, dx, 0], [t + transform.e] * 2 + [t] * 2)[0]) for t in tops])


area_km2 = (watershed_mask.sum(axis=1) * row_cell_areas_m2(transform, watershed_mask.shape[0])).sum() / 1e6

Count the watershed's cells in each row, multiply by that row's geodesic cell area, and sum. Or reproject the DEM to an equal-area or local projected CRS before routing.

Bar chart of the Esopus catchment area from the same 6.25 million geographic cells using square degrees, equator metres, the gauge-latitude cell width on both axes, a constant mid-row cell area and row-by-row geodesic cell areas.
Same cells, same watershed: the area depends entirely on how a degree is turned into metres.

Step-by-step solution

1. Check whether the DEM is geographic

rasterio.open(dem).crs.is_geographic answers it; so does a cell size like 9.26 × 10⁻⁵ or 0.000277. The 3DEP tile here had cells of 9.259 × 10⁻⁵ degrees, 1/3 arc-second.

2. Route in cells, not in areas

Flow direction and accumulation in cells are unaffected by the cell's size in metres: D8 compares elevation drops to the eight neighbours. The distance used for diagonal slopes differs slightly on a geographic grid, which can change a few directions on flat ground, but the counts remain valid. Filling, pointer and accumulation on the geographic DEM took 2.6 s.

3. Do not trust "catchment area" outputs on geographic grids

WhiteboxTools' d8_flow_accumulation with out_type="catchment area" multiplied cells by the cell size in the raster's units: 6,253,412 × (9.259 × 10⁻⁵)² = 0.0536. Other tools do the same, or silently use a fixed number of metres per degree. Convert cells to area yourself.

4. Compute cell area per row

A cell on a geographic grid is a quadrilateral on the ellipsoid whose area depends only on latitude. pyproj.Geod.polygon_area_perimeter gives its exact area. For the rows covering the Esopus basin, the area ran from 78.654 m² at the northern edge to 78.899 m² at the southern edge.

5. Multiply per-row cell counts by per-row areas

Counting the watershed's cells in each row and multiplying by that row's cell area gave 492.67 km². That is the geodesic area of the delineated watershed, equivalent to projecting the cells to an equal-area CRS.

6. Avoid the common shortcuts

  • 111,319.49 m per degree on both axes ignores the narrowing of meridians: at 42° N it overstated the area by 34.85%.
  • The east–west width at the gauge's latitude on both axes shrinks the north–south side too: 25.56% low.
  • One geodesic cell area for the whole catchment is fine for short catchments: 0.01% error at the middle row and 0.12% at the gauge row here.

7. Watch latitude span, not just latitude

For a basin 0.2° tall, the difference between the cell areas at its top and bottom was 0.00% at the equator, 0.20% at 30°, 0.31% at 42°, 0.60% at 60° and 0.96% at 70°. A catchment spanning 10° of latitude — a large river basin — has cells differing by tens of per cent, and a single cell area no longer works.

8. Or reproject before routing

Reprojecting to a projected CRS makes every cell the same size and every downstream area calculation simple. On the Esopus basin the UTM 10 m DEM gave 492.45 km² at the gauge, 0.04% from the geodesic row-by-row answer. Use nearest-neighbour or bilinear resampling at a similar resolution, and route after reprojecting, not before.

Bar chart of how much the equator-metres shortcut overstates catchment area at latitudes 0, 30, 42, 60 and 70 degrees.
The degree-to-metre shortcut is harmless only near the equator; at 60° it doubles the area.

Code examples

Example 1 — route a geographic DEM and read what the tool reports 

import os

import numpy as np
import rasterio
import whitebox

wbt = whitebox.WhiteboxTools()
wbt.set_verbose_mode(False)
wbt.set_working_dir(os.path.abspath("."))
wbt.fill_depressions("esopus_3dep13_geo.tif", "geo_filled.tif", fix_flats=True)
wbt.d8_pointer("geo_filled.tif", "geo_d8.tif")
wbt.d8_flow_accumulation("geo_d8.tif", "geo_cells.tif", out_type="cells", pntr=True)
wbt.d8_flow_accumulation("geo_d8.tif", "geo_area.tif", out_type="catchment area", pntr=True)

with rasterio.open("geo_cells.tif") as src:
    cells, transform, crs = src.read(1), src.transform, src.crs
with rasterio.open("geo_area.tif") as src:
    reported = src.read(1)
col, row = ~transform * (-74.2701944, 42.0144722)
r, c, k = int(row), int(col), 20
window = cells[r - k:r + k + 1, c - k:c + k + 1]
i = np.unravel_index(window.argmax(), window.shape)
R, C = r - k + i[0], c - k + i[1]
print(f"{crs.to_string()}, geographic {crs.is_geographic}, cell {transform.a:.4e} degrees")
print(f"outlet ({R}, {C}): {cells[R, C]:,.0f} cells; tool's catchment area {reported[R, C]:.6f}")
EPSG:4269, geographic True, cell 9.2593e-05 degrees
outlet (2455, 3123): 6,253,412 cells; tool's catchment area 0.053613

Example 2 — geodesic area from per-row cell areas 

import math

import geopandas as gpd
from pyproj import Geod
from shapely.geometry import Point

lon, lat = transform * (C + 0.5, R + 0.5)
gpd.GeoDataFrame(geometry=[Point(lon, lat)], crs=crs).to_file("geo_outlet.shp")
wbt.watershed("geo_d8.tif", "geo_outlet.shp", "geo_ws.tif")
with rasterio.open("geo_ws.tif") as src:
    ws = src.read(1)
    inside = (ws != src.nodata) & (ws > 0)


def row_cell_areas_m2(transform, height, ellps="GRS80"):
    geod, dx = Geod(ellps=ellps), abs(transform.a)
    tops = transform.f + np.arange(height) * transform.e
    return np.array([abs(geod.polygon_area_perimeter([0, dx, dx, 0], [t + transform.e] * 2 + [t] * 2)[0]) for t in tops])


row_area = row_cell_areas_m2(transform, inside.shape[0])
per_row = inside.sum(axis=1)
geodesic = (per_row * row_area).sum() / 1e6
used = np.nonzero(per_row)[0]
n = int(inside.sum())
dx = abs(transform.a)
shortcuts = {
    "square degrees": n * dx * dx,
    "111,319.49 m per degree": n * (111_319.49 * dx) ** 2 / 1e6,
    "width at outlet latitude, squared": n * (111_319.49 * dx * math.cos(math.radians(lat))) ** 2 / 1e6,
    "one cell area, middle row": n * row_area[(used.min() + used.max()) // 2] / 1e6,
    "row-by-row geodesic": geodesic,
}
print(f"{n:,} cells; cell area {row_area[used].min():.3f}-{row_area[used].max():.3f} m2")
for name, value in shortcuts.items():
    print(f"{name:34} {value:10.4f}  {value / geodesic - 1:+8.2%}")
6,253,412 cells; cell area 78.654-78.899 m2
square degrees                         0.0536   -99.99%
111,319.49 m per degree              664.3730   +34.85%
width at outlet latitude, squared    366.7637   -25.56%
one cell area, middle row            492.6232    -0.01%
row-by-row geodesic                  492.6694    +0.00%

Example 3 — accumulate area instead of cells 

from pysheds.grid import Grid
from pysheds.sview import Raster

grid = Grid.from_raster("geo_d8.tif")
pointer = grid.read_raster("geo_d8.tif")
fdir = Raster(np.asarray(pointer).astype("int64"), pointer.viewfinder)
weights = Raster(np.repeat(row_area[:, None] / 1e6, pointer.shape[1], axis=1), pointer.viewfinder)
area_km2 = grid.accumulation(fdir, weights=weights, dirmap=(128, 1, 2, 4, 8, 16, 32, 64))
print(f"area-weighted accumulation at the outlet: {area_km2[R, C]:.2f} km2")
area-weighted accumulation at the outlet: 492.67 km2

Weighting accumulation by cell area gives drainage area in km² at every cell at once, which is what snapping and stream thresholds on a geographic grid need.

Explanation

Why a degree is not a fixed distance

A degree of latitude is about 111 km everywhere, varying by 1% between equator and pole because the Earth is flattened. A degree of longitude is 111 km at the equator and shrinks with the cosine of latitude: 82.7 km at 42° N, 55.8 km at 60° N. A geographic cell is therefore a tall, narrow rectangle away from the equator — 7.7 m by 10.3 m for a 1/3 arc-second cell at the Esopus basin.

Why each shortcut fails the way it does

Squaring the equatorial length treats a 7.7 m × 10.3 m cell as 10.3 m × 10.3 m, overstating by 1 / cos(42°) = 1.35. Squaring the east–west width at the gauge treats it as 7.7 m × 7.7 m, understating by cos(42°) = 0.74. The right answer needs both sides, and the north–south side barely changes while the east–west side changes with latitude.

Why a single cell area often works

The error from using one cell area is set by how much cos(latitude) changes across the catchment. Across 0.2°, near 42° N, it changed by 0.31%, and using the middle row splits that error in half. For national-scale basins, the change is large, and per-row areas are needed.

Why routing in degrees is still acceptable

D8 picks the steepest descent among eight neighbours, dividing elevation drop by distance. On a geographic grid, tools either use true distances or treat cells as square; the choice changes directions only where drops are nearly equal. The areas above come from counting cells that drain to the outlet, which is unaffected by how a cell's size is expressed.

Two panels contrasting a geographic grid cell at 42 degrees north, 7.7 metres wide and 10.3 metres tall, with the square cells assumed by the two common shortcuts.
A 1/3 arc-second cell at 42° N is not square, and neither shortcut uses its real shape.

Edge cases or notes

  • Snap distances and stream thresholds in geographic units need converting too; a threshold of 1,000 cells means different areas at different latitudes.
  • Slope on a geographic DEM needs a z-factor or unit conversion; see why slope values come out wrong.
  • Polygonised watersheds in EPSG:4326 need a geodesic or projected area, not .area.
  • Very large basins cross many latitudes; use per-row areas or an equal-area projection.
  • Datums: GRS80 and WGS84 areas differ by far less than 0.01%.
  • Polar regions make geographic grids strongly anisotropic; reproject before routing.
  • Mixed-resolution mosaics need cell areas computed per source, not per row alone.

FAQ

Why is my catchment area a tiny decimal like 0.05?

It is in square degrees. WhiteboxTools' catchment area output on the geographic Esopus DEM was 0.0536; the real area was 492.67 km².

How do I calculate watershed area from a geographic DEM?

Count the watershed's cells in each row, multiply by that row's geodesic cell area from pyproj, and sum. That gave 492.67 km² for Esopus Creek.

Can I use 111 km per degree to convert cell area?

Not for longitude away from the equator. At 42° N it overstated the Esopus catchment by 35%.

Should I reproject a DEM before delineating a watershed?

Usually, yes: routing on a projected DEM makes every cell the same area. The UTM 10 m DEM gave 492.45 km², within 0.04% of the geodesic answer.

Is flow accumulation wrong on a geographic DEM?

Accumulation in cells is valid; converting cells to area is where errors come in. Weight accumulation by per-row cell area to get km² directly.

When is a single cell area good enough?

When the catchment spans a small range of latitude. For the 0.2°-tall Esopus basin, a mid-row cell area was within 0.01%.