DescribeCoverage publishes four things a GetCoverage request needs and cannot guess: the axis labels, the axis extents, the pixel size per axis, and the band names. Read all four before building any subset — an axis label differing only in case, or an interval written high-then-low, is rejected as InvalidSubsetting with no hint about which of the two it was.
There is no fixed vocabulary of WCS axis names. One service publishes Lat and Long, another lat and lon, a third E and N for a projected grid, and a fourth adds time and ansi for the temporal axis. The subset parameter of a GetCoverage request must use the coverage’s own labels exactly, including case, and a mismatch produces an exception naming the parameter rather than the label.
The document that carries them is well structured but deeply nested, and the values are spread across four separate blocks rather than gathered in one place. Reading them once into a typed record — as the subset download guide assumes has already happened — is what turns coverage requests from trial and error into arithmetic.
from __future__ import annotations
from dataclasses import dataclass
from xml.etree import ElementTree as ET
import requests
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",
"swe": "http://www.opengis.net/swe/2.0",
}
@dataclass(frozen=True)
class Axis:
label: str
uom: str
low: float
high: float
pixels: int
pixel_size: float # signed: negative means the axis descends
def clamp(self, lo: float, hi: float) -> tuple[float, float]:
"""Trim a requested interval into the declared extent, low first."""
lo, hi = min(lo, hi), max(lo, hi)
return max(lo, self.low), min(hi, self.high)
def pixels_for(self, lo: float, hi: float) -> int:
lo, hi = self.clamp(lo, hi)
return max(0, int(round((hi - lo) / abs(self.pixel_size))))
@dataclass(frozen=True)
class Band:
name: str
definition: str
uom: str
nil_values: tuple[float, ...]
@dataclass(frozen=True)
class CoverageProfile:
coverage_id: str
srs: str
axes: tuple[Axis, ...]
bands: tuple[Band, ...]
def axis(self, label: str) -> Axis:
for candidate in self.axes:
if candidate.label == label:
return candidate
raise KeyError(
f"{label!r} is not an axis of {self.coverage_id}; "
f"declared axes are {[a.label for a in self.axes]}")
def estimate_bytes(self, subsets: dict[str, tuple[float, float]],
bands: int | None = None,
bytes_per_sample: int = 4) -> int:
"""Payload estimate before committing to the download."""
total = 1
for axis in self.axes:
if axis.label in subsets:
total *= axis.pixels_for(*subsets[axis.label])
else:
total *= axis.pixels
return total * (bands or len(self.bands)) * bytes_per_sample
def describe(url: str, coverage_id: str, timeout: int = 60) -> CoverageProfile:
resp = requests.get(url, timeout=timeout, params={
"SERVICE": "WCS", "VERSION": "2.0.1",
"REQUEST": "DescribeCoverage", "COVERAGEID": coverage_id,
})
resp.raise_for_status()
return parse_describe(resp.content, coverage_id)
def parse_describe(xml: bytes, coverage_id: str) -> CoverageProfile:
root = ET.fromstring(xml)
envelope = root.find(".//gml:boundedBy/gml:Envelope", NS)
if envelope is None:
raise ValueError("no gml:Envelope — is this a DescribeCoverage response?")
labels = (envelope.get("axisLabels") or "").split()
uoms = (envelope.get("uomLabels") or "").split()
lower = [float(v) for v in envelope.findtext("gml:lowerCorner", "", NS).split()]
upper = [float(v) for v in envelope.findtext("gml:upperCorner", "", NS).split()]
grid = root.find(".//gml:RectifiedGrid", NS) or root.find(".//gml:Grid", NS)
high = [int(v) for v in (grid.findtext(".//gml:high", "", NS) or "").split()] if grid is not None else []
# gml:high is inclusive, so the pixel count is high + 1 per axis.
counts = [v + 1 for v in high]
offsets: list[float] = []
for i, vector in enumerate(grid.findall("gml:offsetVector", NS) if grid is not None else []):
components = [float(v) for v in (vector.text or "").split()]
# The magnitude on this axis is the component on its own dimension.
offsets.append(components[i] if i < len(components) else components[0])
axes = tuple(
Axis(label=labels[i],
uom=uoms[i] if i < len(uoms) else "",
low=lower[i], high=upper[i],
pixels=counts[i] if i < len(counts) else 0,
pixel_size=offsets[i] if i < len(offsets) else
((upper[i] - lower[i]) / counts[i] if i < len(counts) and counts[i] else 1.0))
for i in range(len(labels)))
bands = tuple(
Band(name=field.get("name") or "",
definition=field.findtext(".//swe:Quantity/swe:definition", "", NS),
uom=(field.find(".//swe:uom", NS).get("code")
if field.find(".//swe:uom", NS) is not None else ""),
nil_values=tuple(
float(n.text) for n in field.findall(".//swe:nilValue", NS)
if n.text and n.text.strip().lstrip("-").replace(".", "").isdigit()))
for field in root.findall(".//gmlcov:rangeType//swe:field", NS))
return CoverageProfile(coverage_id, envelope.get("srsName") or "", axes, bands)
def build_subset(profile: CoverageProfile,
requested: dict[str, tuple[float, float]]) -> list[str]:
"""Turn requested intervals into clamped, correctly ordered subset values."""
out = []
for label, (lo, hi) in requested.items():
axis = profile.axis(label) # raises with the real labels listed
low, high = axis.clamp(lo, hi)
out.append(f"{axis.label}({low:g},{high:g})")
return out
gml:high is inclusive. The grid envelope reports the index of the last pixel, not the count, so a coverage 4096 pixels wide publishes high of 4095. Adding one is a single character and getting it wrong shifts every size estimate by exactly one pixel per axis — invisible on a large subset and a real error on a small one.
Offset vectors are vectors, not scalars. A rectified grid publishes one offset vector per grid axis, each with a component on every coordinate axis. For an axis-aligned raster only the diagonal component is non-zero, which is what the code reads; for a rotated grid the off-diagonal components matter and a scalar pixel size is not meaningful at all. Reading the diagonal is correct for the overwhelming majority of published coverages and wrong in a way worth knowing about for the rest.
The sign of the offset carries orientation. A north-up raster has a negative offset on its northing axis, because row zero is the top. Axis.clamp normalises the interval to low-then-high regardless, which is what the subset parameter requires — the sign matters for pixel arithmetic, not for the request.
profile.axis() raises with the real labels. An unknown label is the single most common WCS mistake, and the exception naming the coverage’s actual axes turns a round trip through the server’s exception report into an immediate local answer.
Print the profile and estimate a subset before requesting it:
profile = describe(url, "elevation:dem_25m")
for axis in profile.axes:
print(f"{axis.label:<6} {axis.low:>12.4f} .. {axis.high:<12.4f} "
f"{axis.pixels:>7} px step {axis.pixel_size:+.5f} {axis.uom}")
print("bands:", [b.name for b in profile.bands])
wanted = {"Lat": (46.0, 47.0), "Long": (7.0, 9.0)}
print("estimate:", profile.estimate_bytes(wanted) / 1e6, "MB")
print("subset=", build_subset(profile, wanted))
Lat 45.8180 .. 47.8085 8000 px step -0.00025 deg
Long 5.9559 .. 10.4921 17000 px step +0.00025 deg
bands: ['elevation']
estimate: 128.0 MB
subset= ['Lat(46,47)', 'Long(7,9)']
The estimate is the number worth looking at before pressing send — see the streaming guide for what to do when it is large.
Axis labels are case-sensitive and per-coverage. Two coverages on the same server can use different labels, so caching a profile per coverage rather than per service is the correct granularity. Reusing one coverage’s labels against another is a rejected subset that looks like a server problem.
The temporal axis is not always named time. ansi, date and unix all appear in production, and its values are ISO 8601 strings that must be quoted in the subset parameter rather than written bare. Treating it as just another axis works as long as the value formatting branches on the unit.
nilValue matters more than it looks. A coverage’s nil value is frequently a large negative number such as -9999, and a consumer that averages the array without masking it produces an answer that is confidently wrong. Reading the declared nil values from the range type — rather than assuming the raster’s internal nodata tag will survive the transfer — is what makes the mask reliable.
Some servers publish an incomplete range type. Where band names are missing, rangesubset cannot be used at all and the whole coverage must be fetched. Detecting that at profile time is better than discovering it from an exception mid-batch.
Because the exception does not tell you which of several assumptions failed. InvalidSubsetting is returned for a wrong axis label, a reversed interval and an out-of-range bound alike, so trial and error costs several round trips per mistake. One DescribeCoverage call answers all three questions at once and gives you the payload estimate as well.
Good enough to decide. It assumes an uncompressed sample count times bytes per sample, which is an upper bound for a compressed GeoTIFF and roughly right for an uncompressed one. Since the decision it informs is whether a download is megabytes or gigabytes, an estimate within a factor of two is entirely sufficient.
Yes, per coverage and for the lifetime of a process or a short TTL. The document changes only when the coverage itself is republished, and it is large enough that fetching it per request is wasteful. Cache the parsed profile rather than the XML — the parse is the expensive part.
A referenceable grid has no constant pixel size, so no offset vector to read and no simple pixel arithmetic. The envelope and range type still parse, which is enough to build a subset; what you lose is the ability to estimate the response size in advance, and the practical substitute is requesting a small probe subset first.
Back to WCS Coverage Service Fundamentals
Related