Downloading Coverage Subsets with WCS GetCoverage in Python

TL;DR: A Web Coverage Service GetCoverage request returns raw raster values rather than a rendered picture. In WCS 2.0.1 you carve out a spatial window by attaching one SUBSET=Lat(low,high) trim and one SUBSET=Long(low,high) trim to a plain HTTP GET, ask for FORMAT=image/tiff, and save the bytes as a GeoTIFF. Discover the correct axis labels with DescribeCoverage first, then open the result with rasterio to confirm the window and read a NumPy array.

The Core Challenge

A Web Coverage Service is the OGC protocol for retrieving coverages — the multidimensional grids behind elevation models, satellite imagery, and climate rasters. Unlike a Web Map Service, which paints pixels server-side and hands back a flattened image, WCS returns the underlying data values so you can run analysis on them. That distinction is covered in the WCS Coverage Service Fundamentals overview; if you only need a rendered picture, the Understanding OGC Web Map Service Specifications reference describes the lighter-weight alternative.

The GetCoverage operation is where subsetting happens, and it holds two traps for anyone arriving from a WMS background. First, WCS 2.0.1 does not accept a BBOX parameter at all. Extent is expressed with one SUBSET trim per axis, and each trim names the axis it constrains — Lat, Long, E, N, time, or whatever the coverage advertises. Second, those axis labels are not something you may assume: they are defined by the coverage’s own CRS and must be read from the service before you can build a valid request. Guessing Lat/Long against a projected coverage that actually uses E/N yields an InvalidAxisLabel exception, not a helpful default.

The named-axis design does have one clear payoff. Because a trim is keyed by axis name rather than by position, WCS sidesteps the latitude/longitude ordering ambiguity that makes WMS 1.3.0 BBOX construction so error-prone. The ordering question does not disappear entirely, though — it re-emerges in which CRS your numeric bounds are interpreted against, controlled by the optional subsettingCrs. That interaction is the same axis-order material examined in the SRS and Coordinate Reference System Handling guide.

WCS GetCoverage Subset Data Flow A five-stage pipeline: DescribeCoverage returns axis labels; those labels feed SUBSET trims in a GetCoverage KVP request; requests.get retrieves the coverage; the bytes are written as a GeoTIFF; rasterio opens the file to confirm the window and read a NumPy array. Subsetting a coverage: DescribeCoverage → GetCoverage → GeoTIFF → rasterio DescribeCoverage read axisLabels from gml:Envelope Lat, Long, time GetCoverage KVP COVERAGEID=srtm SUBSET=Lat(51.0,51.6) SUBSET=Long(-0.5,0.3) FORMAT=image/tiff SUBSETTINGCRS=…/4326 requests.get() HTTP GET → 200 image/tiff bytes subset.tif write_bytes() rasterio.open() confirm bounds · read NumPy GDAL / QGIS downstream

Production-Ready Code

The script below performs the whole round trip: it discovers axis labels with DescribeCoverage, builds a GetCoverage request with two spatial trims, guards against an exception report arriving in place of raster bytes, writes the GeoTIFF, and reopens it with rasterio to confirm the window. The only third-party dependencies are requests, lxml, rasterio, and numpy.

"""
wcs_getcoverage_subset.py — download a spatial subset of a WCS 2.0.1 coverage.

Dependencies (pip install):
    requests>=2.31, lxml>=5.0, rasterio>=1.3, numpy>=1.26
"""
from __future__ import annotations

from pathlib import Path

import numpy as np
import rasterio
import requests
from lxml import etree

# GML namespaces used by WCS 2.0.x DescribeCoverage responses.
WCS_NS = {
    "wcs": "http://www.opengis.net/wcs/2.0",
    "gml": "http://www.opengis.net/gml/3.2",
    "gmlcov": "http://www.opengis.net/gmlcov/1.0",
}


def discover_axis_labels(base_url: str, coverage_id: str, timeout: int = 30) -> list[str]:
    """
    Issue a DescribeCoverage request and read the envelope axis labels.

    The gml:Envelope carries an 'axisLabels' attribute whose order follows the
    coverage CRS — e.g. 'Lat Long' for EPSG:4326, 'E N' for a projected grid.
    Those exact strings are what SUBSET parameters must reference.
    """
    params = {
        "SERVICE": "WCS",
        "VERSION": "2.0.1",
        "REQUEST": "DescribeCoverage",
        "COVERAGEID": coverage_id,
    }
    resp = requests.get(base_url, params=params, timeout=timeout)
    resp.raise_for_status()

    root = etree.fromstring(resp.content)
    envelope = root.find(".//gml:Envelope", WCS_NS)
    if envelope is None:
        raise RuntimeError("No gml:Envelope in DescribeCoverage response")

    labels = envelope.get("axisLabels", "").split()
    print(f"CRS         : {envelope.get('srsName', '')}")
    print(f"Axis labels : {labels}")
    return labels


def download_subset(
    base_url: str,
    coverage_id: str,
    out_path: Path,
    lat_range: tuple[float, float],
    long_range: tuple[float, float],
    lat_label: str = "Lat",
    long_label: str = "Long",
    subsetting_crs: str | None = None,
    output_crs: str | None = None,
    range_subset: str | None = None,
    timeout: int = 120,
) -> Path:
    """
    Request a trimmed subset of a coverage and save the returned GeoTIFF.

    Each SUBSET parameter is a trim: Label(low,high). Because axes are named,
    WCS avoids the lon/lat positional ambiguity that WMS BBOX suffers from.
    """
    # A list of tuples lets the SUBSET key repeat once per axis.
    params: list[tuple[str, str]] = [
        ("SERVICE", "WCS"),
        ("VERSION", "2.0.1"),
        ("REQUEST", "GetCoverage"),
        ("COVERAGEID", coverage_id),
        ("SUBSET", f"{lat_label}({lat_range[0]},{lat_range[1]})"),
        ("SUBSET", f"{long_label}({long_range[0]},{long_range[1]})"),
        ("FORMAT", "image/tiff"),
    ]
    if subsetting_crs:  # CRS the low/high numbers are expressed in
        params.append(("SUBSETTINGCRS", subsetting_crs))
    if output_crs:      # reproject the returned grid to this CRS
        params.append(("OUTPUTCRS", output_crs))
    if range_subset:    # Range Subsetting extension: pick bands, e.g. "red,nir"
        params.append(("RANGESUBSET", range_subset))

    resp = requests.get(base_url, params=params, timeout=timeout)
    resp.raise_for_status()

    # An error comes back as an XML ServiceExceptionReport, not a GeoTIFF.
    if "xml" in resp.headers.get("Content-Type", "").lower():
        raise RuntimeError(f"WCS returned an exception:\n{resp.text[:2000]}")

    out_path.write_bytes(resp.content)
    return out_path


def inspect_geotiff(path: Path) -> None:
    """Open the saved GeoTIFF and confirm the window, then read a band."""
    with rasterio.open(path) as ds:
        print(f"Driver      : {ds.driver}")
        print(f"CRS         : {ds.crs}")
        print(f"Size        : {ds.width} x {ds.height} ({ds.count} band(s))")
        print(f"Bounds      : {ds.bounds}")
        band1 = ds.read(1)  # NumPy ndarray, shape (height, width)
        print(f"Band 1      : dtype={band1.dtype}, "
              f"min={np.nanmin(band1)}, max={np.nanmax(band1)}")


if __name__ == "__main__":
    BASE_URL = "https://demo.example.org/geoserver/wcs"
    COVERAGE = "elevation__srtm"

    labels = discover_axis_labels(BASE_URL, COVERAGE)
    # Use the two horizontal labels the server actually advertises.
    lat_label, long_label = (labels + ["Lat", "Long"])[:2]

    out = download_subset(
        BASE_URL,
        COVERAGE,
        out_path=Path("subset.tif"),
        lat_range=(51.0, 51.6),
        long_range=(-0.5, 0.3),
        lat_label=lat_label,
        long_label=long_label,
        subsetting_crs="http://www.opengis.net/def/crs/EPSG/0/4326",
    )
    print(f"Saved {out.resolve()}")
    inspect_geotiff(out)

Step-by-Step Walkthrough

Axis-label discovery (discover_axis_labels). The function fires a DescribeCoverage request and parses the response with lxml. The value you need lives in the axisLabels attribute of the gml:Envelope element — a space-separated string such as Lat Long or E N. Reading it programmatically means the same script works against a geographic coverage and a projected one without hard-coding label names. The companion srsName attribute tells you which CRS those labels belong to, which is exactly the information you need when deciding whether to set subsettingCrs.

Repeating the SUBSET key (download_subset). WCS expresses a two-dimensional window as two separate SUBSET parameters, so the same query-string key appears twice. Passing a dict to requests cannot represent that — a dict collapses duplicate keys — so params is built as a list of (key, value) tuples. Each trim uses the Label(low,high) syntax; the low and high values are the minimum and maximum along that named axis.

Optional CRS controls. SUBSETTINGCRS declares the reference system your low/high numbers are expressed in. When omitted, the server interprets them in the coverage’s native CRS, so setting it explicitly to http://www.opengis.net/def/crs/EPSG/0/4326 removes any doubt about how your bounds are read. OUTPUTCRS is the mirror control: it asks the server to reproject the returned grid into a different CRS before encoding it. Keeping the two concerns separate is deliberate — mixing them up is the single most common source of a subset that lands in the wrong place, a failure mode explored in the SRS and Coordinate Reference System Handling guide.

Guarding against exception reports. A WCS answers a malformed request with HTTP 200 and an XML ServiceExceptionReport body rather than a non-2xx status in many server configurations, so raise_for_status() alone will not catch it. The Content-Type check intercepts an xml response before it is written to disk and surfaces the server’s error text instead of leaving you with a .tif file that is really an error document.

Range/band subsetting. The range_subset argument feeds the optional RANGESUBSET parameter from the WCS Range Subsetting extension. It selects range components — bands — by the names published in the gmlcov:rangeType block of DescribeCoverage, either as a comma list (red,green,blue) or a colon interval (band1:band3). It is orthogonal to spatial trimming and can be combined with the SUBSET trims in one request.

Reopening with rasterio (inspect_geotiff). Once the bytes are on disk, rasterio.open() reads the GeoTIFF’s embedded georeferencing. Printing ds.bounds confirms the returned window matches the trim you asked for, and ds.read(1) hands back the first band as a NumPy array ready for analysis, reprojection, or ingestion into a downstream pipeline.

Verification

Run the script against a coverage endpoint. The console output confirms both that the discovered labels were used and that the saved file carries the expected window:

python wcs_getcoverage_subset.py

Expected output (values depend on the coverage’s resolution):

CRS         : http://www.opengis.net/def/crs/EPSG/0/4326
Axis labels : ['Lat', 'Long']
Saved /home/user/subset.tif
Driver      : GTiff
CRS         : EPSG:4326
Size        : 960 x 720 (1 band(s))
Bounds      : BoundingBox(left=-0.5, bottom=51.0, right=0.3, top=51.6)
Band 1      : dtype=int16, min=-4, max=245

Cross-check the file independently with the GDAL command-line tools if rasterio is unavailable:

gdalinfo subset.tif | grep -E "Size is|Lower Left|Upper Right"

The reported corner coordinates should bracket the Lat/Long ranges you supplied. If the bounds are transposed or land in the wrong hemisphere, the numeric values were interpreted against a different CRS than you intended — set SUBSETTINGCRS explicitly and retry.

Gotchas & Edge Cases

My SUBSET request returns an InvalidAxisLabel exception. What is wrong?

The label inside SUBSET must exactly match an axis label advertised in the coverage’s DescribeCoverage envelope, and the match is case-sensitive. A coverage in EPSG:4326 usually exposes Lat and Long, but a projected coverage may expose E and N, or X and Y. Never assume — run DescribeCoverage first, read the axisLabels attribute on gml:Envelope, and use those strings verbatim, which is exactly what discover_axis_labels automates.

Do I have to swap latitude and longitude like I do for a WMS 1.3.0 BBOX?

No. SUBSET parameters are keyed by named axis rather than by position, so SUBSET=Lat(...) and SUBSET=Long(...) are unambiguous regardless of the CRS’s declared axis order. The ordering concern instead moves to which CRS your numbers are read in: they are interpreted in the subsettingCrs, which defaults to the coverage’s native CRS. Set SUBSETTINGCRS explicitly whenever your bounds are expressed in a different system, and consult the SRS and Coordinate Reference System Handling guide when the returned window lands somewhere unexpected.

What is the difference between a trim and a slice in a WCS subset?

A trim supplies two bounds — SUBSET=Lat(51.0,51.6) — and returns a coverage that keeps that axis but with a reduced extent. A slice supplies a single value — SUBSET=Lat(51.3) — and collapses the axis entirely, lowering the coverage’s dimensionality by one. Use trims to cut a spatial window; use slices to pick a single timestamp, elevation band, or layer from a coverage that has more than two dimensions.

How do I download only specific bands from a multi-band coverage?

Use the Range Subsetting extension parameter RANGESUBSET. Pass a comma-separated list of range component (band) names, such as RANGESUBSET=red,green,blue, or a contiguous interval with a colon, such as RANGESUBSET=band1:band3. The component names are published in the gmlcov:rangeType section of DescribeCoverage. Band selection is independent of spatial SUBSET trimming, and both can appear in the same GetCoverage request to fetch a smaller window of fewer bands.

The server ignored my FORMAT=image/tiff and returned XML. Why?

Two causes are common. First, the request may have failed validation, in which case the XML is a ServiceExceptionReport describing the fault — the Content-Type guard in download_subset catches this and prints the message. Second, the coverage may not advertise image/tiff as a supported output format; check the wcs:ServiceMetadata/wcs:formatSupported entries in the GetCapabilities document and request one of the listed MIME types instead.

Back to WCS Coverage Service Fundamentals

Related