API Reference

The public API is the two accessors plus the function behind them. import xcorduroy registers .dem on DataArray and Dataset; in most cases that is all you need, and compute_terrain is the escape hatch for calling the machinery directly.

Terrain modes (Slope, Aspect, Hillshade) are frozen dataclasses carrying the output's name, units and long_name, and in the hillshade case the light position. They are what compute_terrain dispatches on.

The kernel

All three methods come from one function, _terrain_kernel in src/xcorduroy/DEM.py. It takes a 1-pixel-padded 2D elevation array and returns an array of the same shape as the unpadded input.

Gradients use the Horn (1981) 3x3 kernel — a weighted finite difference over the eight neighbours, with the centre-row and centre-column neighbours weighted 2x:

dz/dcol = ((z[-1,+1] + 2*z[0,+1] + z[+1,+1])
         - (z[-1,-1] + 2*z[0,-1] + z[+1,-1])) / (8 * res_x)

and the transpose for dz/drow. This is what gdaldem uses, which is why the outputs match it.

Array orientation

The kernel works in array index space, but slope and aspect are defined in world space. compute_terrain takes the signed median spacing of each coordinate axis and passes x_sign / y_sign into the kernel, which flips the gradients so they are positive eastward and northward regardless of how the raster is stored. A north-up raster (descending y) and a south-up one give the same answer.

Longitude is differenced with a dateline guard: a step larger than 180 degrees is wrapped, so an array crossing the antimeridian does not read as one enormous cell.

Conventions

Output	Range	Notes
`slope`	0–90 degrees	`atan(hypot(dz_dx, dz_dy))`
`aspect`	0–360 degrees	Downslope direction, clockwise from north. Flat cells are `NaN`.
`hillshade`	0–1	Clipped, `float32`

Hillshade is the standard Lambertian term:

cos(zenith) * cos(slope) + sin(zenith) * sin(slope) * cos(azimuth - aspect)

with zenith = 90 - altitude, and azimuth converted from compass degrees to the counter-clockwise-from-east convention that arctan2 returns.

Why flat cells are NaN

A perfectly flat cell has no downslope direction. arctan2(0, 0) returns 0, which would report every flat cell as facing due north — visible as a false wall of colour across hydro-flattened lakes and plateaus. xcorduroy returns NaN there instead, matching gdaldem. Slope and hillshade are unaffected: a flat cell has zero slope and takes the illumination of a horizontal surface.

Degrees to metres

Elevation is in metres, but a geographic CRS gives coordinate spacing in degrees. Dividing one by the other is meaningless, so for a geographic CRS the spacing is scaled to metres before the gradient is taken:

res_y = dy * 111320
res_x = dx * 111320 * cos(mean_latitude)

111320 is the WGS84 mean metres per degree of latitude. Using the mean latitude of the array is the approximation: the true longitude scale varies across the array, and the error grows with its north–south extent. For a tile this is negligible; for a continental extent it is not, and the right answer is to reproject to an equal-area or local projected CRS first.

Near a pole cos(lat) collapses toward zero, which would make every cell read as a cliff. The scaling is floored at its 89.9-degree value and a RuntimeWarning is raised.

An explicit resolution= bypasses all of this and is taken in elevation units.

Chunked arrays

A 3x3 kernel needs one cell of context beyond each chunk boundary. For dask-backed input the kernel is applied through map_overlap with depth=1, boundary="nearest" and trim=False: dask grows each chunk by one cell on every side, the kernel consumes that halo, and the output chunk comes back the same shape as the input chunk. Seam values are therefore identical to the in-memory result — tests/test_DEM.py::test_dask_chunk_seams asserts exactly that.

The in-memory path does the same thing with np.pad(..., mode="edge"), so the two code paths agree at array edges as well as at chunk seams.

NaN handling

NaN is not treated as nodata; it propagates. A single NaN cell in the input produces a 3x3 block of NaN in the output: every neighbour of that cell had it in its stencil, and the centre cell is masked explicitly, since the Horn kernel gives it weight 0 and would otherwise return a finite, meaningless value. Mask or fill nodata before calling if that is not what you want.

Testing

Beyond unit tests against hand-derived analytic expectations, tests/test_hypothesis.py asserts properties that must hold for any input: a flat surface has zero slope, a linear ramp has constant slope, slope is invariant under translation and rotation of the elevation field, aspect stays in 0–360, hillshade stays in 0–1, and halving the resolution doubles the gradient.

Signatures

xcorduroy.DEMDataArrayAccessor

DEMDataArrayAccessor(xarray_obj: DataArray)

Xarray accessor for terrain analysis on DataArrays.

Provides methods to compute slope, aspect, and hillshade from elevation data. Access via da.dem.slope(), da.dem.aspect(), or da.dem.hillshade().

crs `property`

crs: Any

Get coordinate reference system from DataArray.

Returns:

Any –

CRS object from xproj accessor

Raises:

ValueError –

If no CRS is found

aspect

aspect(
    x: str | None = None,
    y: str | None = None,
    resolution: Resolution = None,
    **kwargs: Any,
) -> DataArray

Compute aspect: downslope direction in degrees clockwise from north (0-360). Flat cells are NaN.

Parameters:

x (str | None, default: None ) –

X dimension name (auto-detected if None)
y (str | None, default: None ) –

Y dimension name (auto-detected if None)
resolution (Resolution, default: None ) –

Cell size or (y_res, x_res) in elevation units. Derived from coordinates if None (degrees converted to metres for geographic CRSs)
**kwargs (Any, default: {} ) –

Additional arguments (e.g., z_factor vertical exaggeration)

Returns:

DataArray –

DataArray with aspect values in degrees

hillshade

hillshade(
    x: str | None = None,
    y: str | None = None,
    resolution: Resolution = None,
    azimuth: float = 315.0,
    altitude: float = 45.0,
    **kwargs: Any,
) -> DataArray

Compute hillshade (0-1 shaded relief).

Parameters:

x (str | None, default: None ) –

X dimension name (auto-detected if None)
y (str | None, default: None ) –

Y dimension name (auto-detected if None)
resolution (Resolution, default: None ) –

Cell size or (y_res, x_res) in elevation units. Derived from coordinates if None (degrees converted to metres for geographic CRSs)
azimuth (float, default: 315.0 ) –

Light source direction in degrees (0-360)
altitude (float, default: 45.0 ) –

Light source angle above horizon in degrees (0-90)
**kwargs (Any, default: {} ) –

Additional arguments (e.g., z_factor vertical exaggeration)

Returns:

DataArray –

DataArray with hillshade values (0-1)

slope

slope(
    x: str | None = None,
    y: str | None = None,
    resolution: Resolution = None,
    **kwargs: Any,
) -> DataArray

Compute slope in degrees.

Parameters:

x (str | None, default: None ) –

X dimension name (auto-detected if None)
y (str | None, default: None ) –

Y dimension name (auto-detected if None)
resolution (Resolution, default: None ) –

Cell size or (y_res, x_res) in elevation units. Derived from coordinates if None (degrees converted to metres for geographic CRSs)
**kwargs (Any, default: {} ) –

Additional arguments (e.g., z_factor vertical exaggeration)

Returns:

DataArray –

DataArray with slope values in degrees

xcorduroy.DEMDatasetAccessor

DEMDatasetAccessor(xarray_obj: Dataset)

Xarray accessor for terrain analysis on Datasets.

Auto-detects elevation variable or access specific variables via ds.dem('var_name'). Provides methods to compute slope, aspect, and hillshade.

call

__call__(name: str | None = None) -> DEMDataArrayAccessor

Get DEM accessor for a specific variable or auto-detect elevation variable.

Parameters:

name (str | None, default: None ) –

Variable name or None to auto-detect

Returns:

DEMDataArrayAccessor –

DEMDataArrayAccessor for the selected variable

Raises:

ValueError –

If multiple candidate variables found
AttributeError –

If no suitable elevation variable found

aspect

aspect(**kwargs: Any) -> DataArray

Compute aspect. See DEMDataArrayAccessor.aspect for details.

hillshade

hillshade(**kwargs: Any) -> DataArray

Compute hillshade. See DEMDataArrayAccessor.hillshade for details.

slope

slope(**kwargs: Any) -> DataArray

Compute slope. See DEMDataArrayAccessor.slope for details.

xcorduroy.compute_terrain

compute_terrain(
    da: DataArray,
    mode: ModeType,
    resolution: float | int | tuple[float, float] | None = None,
    crs: Any = None,
    x_dim: str = "x",
    y_dim: str = "y",
    **kwargs: Any,
) -> DataArray

Compute terrain analysis from a 2D elevation DataArray.

Parameters:

da (DataArray) –

Input elevation DataArray with dims y_dim and x_dim
mode (ModeType) –

Terrain mode (Slope, Aspect, or Hillshade instance)
resolution (float | int | tuple[float, float] | None, default: None ) –

Cell size as scalar or (y_res, x_res) in elevation units. Derived from coordinates if None. For geographic CRSs the degree spacing is converted to metres using the mean latitude; within 0.1 degrees of a pole that scaling is clamped and a RuntimeWarning is raised.
crs (Any, default: None ) –

pyproj-like CRS with is_geographic. Used for degree-to-metre scaling
x_dim (str, default: 'x' ) –

Name of x dimension
y_dim (str, default: 'y' ) –

Name of y dimension
**kwargs (Any, default: {} ) –

Passed to the kernel (e.g. z_factor vertical exaggeration)

Returns:

DataArray –

DataArray with computed terrain values

xcorduroy.types.TerrainMode `dataclass`

TerrainMode()

Bases: ABC

Base class for all terrain analysis methods.

xcorduroy.types.Slope `dataclass`

Slope(
    name: str = "slope",
    units: str = "degrees",
    long_name: str = "terrain_slope",
)

Bases: TerrainMode

xcorduroy.types.Aspect `dataclass`

Aspect(
    name: str = "aspect",
    units: str = "degrees",
    long_name: str = "terrain_aspect",
)

Bases: TerrainMode

xcorduroy.types.Hillshade `dataclass`

Hillshade(
    azimuth: float = 315.0,
    altitude: float = 45.0,
    name: str = "hillshade",
    units: str = "0-1",
    long_name: str = "hillshade",
)

Bases: TerrainMode

API Reference

The kernel

Array orientation

Conventions

Why flat cells are NaN

Degrees to metres

Chunked arrays

NaN handling

Testing

Signatures

xcorduroy.DEMDataArrayAccessor

crs property

aspect

hillshade

slope

xcorduroy.DEMDatasetAccessor

__call__

aspect

hillshade

slope

xcorduroy.compute_terrain

xcorduroy.types.TerrainMode dataclass

xcorduroy.types.Slope dataclass

xcorduroy.types.Aspect dataclass

xcorduroy.types.Hillshade dataclass

crs `property`

call

xcorduroy.types.TerrainMode `dataclass`

xcorduroy.types.Slope `dataclass`

xcorduroy.types.Aspect `dataclass`

xcorduroy.types.Hillshade `dataclass`