isobands

An easy way to make filled contour maps with Python.

isobands converts a regular two-dimensional xarray DataArray into contoured polygons in a GeoPandas GeoDataFrame.

Installation

By default, isobands requires an installed, matching osgeo.gdal Python binding. In a Conda or system-managed environment that provides such a binding, you can install isobands directly:

pip install isobands

Lacking that, you can let pip build the binding for you by specifying the GDAL version. The PyPI GDAL package is source-only: it compiles against the native GDAL library, which must already be installed with its development headers at the matching version. The latest tested version is GDAL 3.13.2.

pip install "isobands[gdal313]"

Installers for gdal310, gdal311 and gdal312 are also available.

If the install fails with gdal-config: No such file or directory, the native library is missing. Install it first — for example brew install gdal on macOS or apt install libgdal-dev on Debian and Ubuntu — then choose the extra matching gdal-config --version, or use conda-forge bindings instead.

Verify the installed bindings, their native-library version, and in-memory contour support after installation:

python -m isobands check

The command exits nonzero when a check fails and includes installation guidance. Add --json for machine-readable output, or call isobands.check() to receive the same structured report in Python.

Quick start

Pass an in-memory xarray.DataArray and receive a geopandas.GeoDataFrame:

import numpy as np
import xarray as xr

import isobands

data = xr.DataArray(
    np.array([[0.0, 1.0, 2.0], [1.0, 2.0, 3.0], [2.0, 3.0, 4.0]]),
    dims=("y", "x"),
    coords={"x": [0.0, 1.0, 2.0], "y": [2.0, 1.0, 0.0]},
)
bands = isobands.from_raster(data, levels=[1.5, 2.5], crs="EPSG:4326")
print(bands[["min_value", "max_value", "geometry"]])

The levels input supplies the interior breakpoints, dividing the data into bands at 1.5 and 2.5.

The result always has the stable min_value, max_value, and geometry columns and a GeoPandas CRS.

The returned GeoDataFrame of the example above would look like so:

min_value	max_value	geometry
0.0	1.5	MULTIPOLYGON (...)
1.5	2.5	MULTIPOLYGON (...)
2.5	4.0	MULTIPOLYGON (...)

Equal interval breaks

You can use the intervals input to specify uniform threshold breaks. This example uses interval=5 to create five-degree temperature bands.

import xarray as xr

import isobands

temperature = xr.open_dataarray("west-coast-daily-highs.nc")
bands = isobands.from_raster(temperature, interval=5, crs="EPSG:4326")

For this field, the returned GeoDataFrame begins like this:

min_value	max_value	geometry
12.3	15.0	MULTIPOLYGON (...)
15.0	20.0	MULTIPOLYGON (...)
20.0	25.0	MULTIPOLYGON (...)

Set offset when equal-width bands need a nonzero starting point. For example, interval=5, offset=2.5 creates thresholds at 2.5, 7.5, 12.5, and so on.

The map shows high temperatures on August 16, 2020, when Death Valley recorded a record high of 54.4°C in Furnace Creek.

Threshold breaks

The levels input allows you to specify whatever threshold you like. This example uses the EPA’s air quality risk categories for particulate matter.

import xarray as xr

import isobands

pm25 = xr.open_dataarray("nyc-smoke-pm25.nc")
bands = isobands.from_raster(
    pm25,
    levels=[12.0, 35.4, 55.4, 125.4, 225.4],
    crs="EPSG:4326",
)

For this field, the returned GeoDataFrame begins like this:

min_value	max_value	geometry
1.4	12.0	MULTIPOLYGON (...)
12.0	35.4	MULTIPOLYGON (...)
35.4	55.4	MULTIPOLYGON (...)

The map shows EPA AirData risk scores on June 7, 2023, when Canadian wildfire smoke blanketed the East Coast.

Data-derived breaks

Pass a callable to levels to determine breaks based on the data. The function receives a one-dimensional array of valid values and should return the interior thresholds. This can be useful for quintiles, natural breaks, Jenks breaks and other data-driven strategies.

import numpy as np
import xarray as xr

import isobands


def quintiles(values):
    return np.quantile(values, [0.2, 0.4, 0.6, 0.8])


rainfall = xr.open_dataarray("harvey-daily-rainfall.nc")
bands = isobands.from_raster(rainfall, levels=quintiles, crs="EPSG:4326")

For this field, the returned GeoDataFrame begins like this:

min_value	max_value	geometry
0.0	0.1	MULTIPOLYGON (...)
0.1	0.9	MULTIPOLYGON (...)
0.9	4.3	MULTIPOLYGON (...)

The map shows rainfall on August 27, 2017, as Hurricane Harvey stalled over Texas.

No-data cells

Use nodata to leave unavailable cells out of the contours instead of treating them as measured values.

import xarray as xr

import isobands

temperature = xr.open_dataarray("iowa-land-surface-temperature.nc")
bands = isobands.from_raster(
    temperature,
    levels=[20, 25, 30, 35],
    nodata=-9999,
    crs="EPSG:4326",
)

For this field, the returned GeoDataFrame begins like this:

min_value	max_value	geometry
12.1	20.0	MULTIPOLYGON (...)
20.0	25.0	MULTIPOLYGON (...)
25.0	30.0	MULTIPOLYGON (...)

The map shows land-surface temperatures captured by NASA satellites on August 12, 2020, two days after a historic derecho swept across Iowa. The no-data cells indicate pixels without valid data, typically due to cloud cover obscuring the satellite’s sensors.

API reference

isobands.check()

Check whether the active environment can generate isobands.

Ordinary installation failures are returned in the report rather than raised.

Return type:

CheckReport

class isobands.CheckReport(ok, checks)

Bases: object

Structured GDAL diagnostic report returned by check().

Parameters:

• ok (bool)

• checks (tuple[CheckResult, ...])

ok: bool

checks: tuple[CheckResult, ...]

to_dict()

Return a JSON-serializable representation.

Return type:

dict[str, object]

__init__(ok, checks)

Parameters:

• ok (bool)

• checks (tuple[CheckResult, ...])

Return type:

None

class isobands.CheckResult(name, ok, observed, message, guidance)

Bases: object

One result returned by check().

Parameters:

• name (Literal['python_bindings', 'gdal_versions', 'supported_gdal_version', 'contour_smoke'])

• ok (bool)

• observed (dict[str, str])

• message (str)

• guidance (str)

name: Literal['python_bindings', 'gdal_versions', 'supported_gdal_version', 'contour_smoke']

ok: bool

observed: dict[str, str]

message: str

guidance: str

__init__(name, ok, observed, message, guidance)

Parameters:

• name (Literal['python_bindings', 'gdal_versions', 'supported_gdal_version', 'contour_smoke'])

• ok (bool)

• observed (dict[str, str])

• message (str)

• guidance (str)

Return type:

None

isobands.from_raster(data, *, levels=None, interval=None, offset=0.0, crs=None, nodata=None)

Create finite, filled contour polygons from a two-dimensional raster.

Exactly one of levels or interval is required. Explicit levels are interior thresholds; a callable receives a one-dimensional array of valid raster values and returns those thresholds. The returned outer bands begin and end at the valid raster extrema. Interval thresholds are integral multiples of the supplied interval, optionally shifted by offset. A constant raster returns one covering band whose min_value and max_value are necessarily equal. Interval requests are limited to 100,000 interior thresholds; use a larger interval for wider value ranges. Integer samples must be within Float64’s exact consecutive-integer range because GDAL contours using Float64 values.

Parameters:

• data (DataArray)

• levels (Sequence[float] | ndarray | Callable[[ndarray], Sequence[float] | ndarray] | None)

• interval (float | None)

• offset (float)

• crs (Any | None)

• nodata (float | None)

Return type:

GeoDataFrame

Links

• Source code

• Issue tracker

• Changelog

• PyPI package

About

Ben Welsh first released this module in August 2026 as a spinoff of the Reuters Climate Monitor. GitHub’s Copilot, an AI-powered text generator, helped draft this documentation. Map examples use MapLibre, OpenFreeMap, and © OpenStreetMap contributors.