xrspatial.contour.contours#

xrspatial.contour.contours(agg: DataArray, levels: Sequence[float] | ndarray | None = None, n_levels: int = 10, return_type: str = 'numpy') → List[Tuple[float, ndarray]] | gpd.GeoDataFrame[source]#

Extract contour lines (isolines) from a raster surface.

Uses the marching squares algorithm to trace isolines through the raster at specified elevation values. Each 2x2 cell quad is classified independently, so the algorithm parallelizes across Dask chunks.

Parameters:

  • agg (xr.DataArray) – 2D input raster (e.g. a DEM).

  • levels (sequence of float, optional) – Explicit contour levels to extract. If not provided, n_levels evenly spaced levels are chosen between the raster min and max.

  • n_levels (int, default 10) – Number of contour levels to generate when levels is not provided. Must be an integer >= 1; otherwise a TypeError (non-integer) or ValueError (< 1) is raised. Ignored when levels is given, in which case it is not validated.

  • return_type (str, default "numpy") – Output format. "numpy" returns a list of (level, coords) tuples where coords is an Nx2 array of (y, x) coordinates in the DataArray’s coordinate space. "geopandas" returns a GeoDataFrame with level and geometry columns (requires geopandas/shapely).

Returns:

Contour lines grouped by level.

Return type:

list of (float, ndarray) or GeoDataFrame

Notes

CuPy and Dask+CuPy arrays are accepted as input. Data is transferred to CPU for the tracing step because segment stitching is an inherently sequential graph traversal. For Dask inputs, chunking bounds the per-chunk scan buffers, but the global merge step materializes all contour segments at once to stitch polylines across chunk boundaries, so peak memory scales with total contour complexity rather than chunk size.

Corners are classified with >= level, so a corner exactly equal to the level is treated as above it. This can make a traced segment collapse onto a single point. Such degenerate (zero-length or single-distinct-point) segments and polylines are dropped before output, so no zero-length or invalid geometry reaches the numpy or geopandas result. The rule is applied identically across the numpy, cupy, dask+numpy, and dask+cupy backends.

Examples

>>> from xrspatial import contours
>>> lines = contours(dem, levels=[100, 500, 1000])
>>> # Each entry is (level_value, Nx2_coordinate_array)
>>> level, coords = lines[0]