Estimate the response size from DescribeCoverage before requesting anything, split the extent into bounded tiles, stream each response to disk with iter_content, and read the files back through rasterio windows. Peak memory then tracks one tile, the run is resumable, and a failure costs one request rather than an hour.
Feature services page. OGC API - Features has limit and a rel="next" link; WFS has count and startIndex. A coverage service has neither. GetCoverage returns exactly the subset you asked for as one response, so the only mechanism for bounding a response is to ask for less — which means the client has to compute what “less” means before it asks.
That computation needs the axis extents and pixel sizes from DescribeCoverage. Without it, the usual sequence is: request a plausible-looking bounding box, wait several minutes, watch the process get killed, and repeat with a smaller box. With it, the size is arithmetic and the tiling falls out of a target byte budget.
from __future__ import annotations
import os
from dataclasses import dataclass
from typing import Iterator
import numpy as np
import rasterio
import requests
from rasterio.windows import Window
@dataclass(frozen=True)
class Tile:
index: tuple[int, int]
subsets: dict[str, tuple[float, float]]
@property
def name(self) -> str:
return f"tile_{self.index[0]:03d}_{self.index[1]:03d}.tif"
def plan_tiles(profile, extent: dict[str, tuple[float, float]],
budget_bytes: int = 256 * 1024 * 1024,
bytes_per_sample: int = 4) -> list[Tile]:
"""Split a requested extent into tiles no larger than the byte budget.
`profile` is the CoverageProfile from DescribeCoverage: it supplies the
pixel size per axis, which is what turns a byte budget into a degree or
metre span.
"""
x_label, y_label = list(extent)[:2]
x_axis, y_axis = profile.axis(x_label), profile.axis(y_label)
x_lo, x_hi = x_axis.clamp(*extent[x_label])
y_lo, y_hi = y_axis.clamp(*extent[y_label])
bands = max(1, len(profile.bands))
samples = budget_bytes // (bytes_per_sample * bands)
side_px = max(256, int(samples ** 0.5))
x_span = side_px * abs(x_axis.pixel_size)
y_span = side_px * abs(y_axis.pixel_size)
tiles: list[Tile] = []
row = 0
y = y_lo
while y < y_hi:
col = 0
x = x_lo
while x < x_hi:
tiles.append(Tile((row, col), {
x_label: (x, min(x + x_span, x_hi)),
y_label: (y, min(y + y_span, y_hi)),
}))
x += x_span
col += 1
y += y_span
row += 1
return tiles
def fetch_tile(url: str, coverage_id: str, tile: Tile, directory: str,
fmt: str = "image/tiff", timeout: int = 600) -> str:
"""Stream one tile to disk. Never touches resp.content."""
path = os.path.join(directory, tile.name)
if os.path.exists(path) and os.path.getsize(path) > 0:
return path # already fetched — the run resumes
params = [("SERVICE", "WCS"), ("VERSION", "2.0.1"),
("REQUEST", "GetCoverage"), ("COVERAGEID", coverage_id),
("FORMAT", fmt)]
for label, (lo, hi) in tile.subsets.items():
params.append(("SUBSET", f"{label}({lo:g},{hi:g})"))
partial = path + ".part"
with requests.get(url, params=params, stream=True, timeout=timeout) as resp:
resp.raise_for_status()
ctype = resp.headers.get("content-type", "")
if "xml" in ctype:
raise RuntimeError(f"service exception: {resp.text[:400]}")
with open(partial, "wb") as fh:
for chunk in resp.iter_content(chunk_size=1 << 20):
fh.write(chunk)
os.replace(partial, path) # atomic: no half file is ever seen
return path
def windows(path: str, size: int = 1024) -> Iterator[tuple[Window, np.ndarray]]:
"""Yield (window, array) pairs — never the whole band."""
with rasterio.open(path) as src:
for row in range(0, src.height, size):
for col in range(0, src.width, size):
window = Window(col, row,
min(size, src.width - col),
min(size, src.height - row))
yield window, src.read(1, window=window)
@dataclass
class Running:
"""A fold, not a collection — this is what keeps memory flat."""
count: int = 0
total: float = 0.0
minimum: float = float("inf")
maximum: float = float("-inf")
def update(self, values: np.ndarray, nodata: float | None) -> None:
data = values if nodata is None else values[values != nodata]
if data.size == 0:
return
self.count += int(data.size)
self.total += float(data.sum(dtype="float64"))
self.minimum = min(self.minimum, float(data.min()))
self.maximum = max(self.maximum, float(data.max()))
@property
def mean(self) -> float:
return self.total / self.count if self.count else float("nan")
def summarise(paths: list[str], nodata: float | None = -9999.0) -> Running:
stats = Running()
for path in paths:
for _, block in windows(path):
stats.update(block, nodata)
return stats
The byte budget drives the tile size, not the other way round. plan_tiles converts a target response size into a pixel side length, then into a coordinate span using the pixel size from the coverage profile. Choosing a tile size in degrees directly is the common alternative and it produces wildly different response sizes across coverages of different resolutions.
iter_content and an atomic rename. Streaming to a .part file and renaming on completion means a killed process never leaves a truncated GeoTIFF that a later run would treat as done. Combined with the existence check at the top of fetch_tile, the whole job becomes resumable: re-running skips what completed and retries only what did not.
Check the content type before writing. A WCS exception arrives as XML with an HTTP 200, so a client that streams straight to disk saves an error document with a .tif extension and discovers it later when rasterio refuses to open it. One header check converts that into an immediate, readable error.
Running folds rather than collects. The whole point of windowed reads is defeated by appending each window to a list and concatenating at the end. Accumulating a count, a sum and the extremes gives mean, min and max at constant memory; percentiles need more care, and a t-digest or a histogram is the right structure when they are required.
Watch peak memory while the pipeline runs over a coverage far larger than RAM:
import resource
profile = describe(url, "elevation:dem_25m")
tiles = plan_tiles(profile, {"Long": (5.9, 10.5), "Lat": (45.8, 47.8)})
print(f"{len(tiles)} tiles at ~256 MB each")
paths = [fetch_tile(url, "elevation:dem_25m", t, "work/") for t in tiles]
stats = summarise(paths)
peak_mb = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
print(f"mean={stats.mean:.2f} min={stats.minimum} max={stats.maximum} "
f"peak_rss={peak_mb:.0f} MB")
54 tiles at ~256 MB each
mean=1043.71 min=193.0 max=4478.0 peak_rss=214 MB
A resident set of a couple of hundred megabytes while summarising tens of gigabytes is the property being tested. If it tracks the data volume instead, one of the three lines in the diagnosis below is present.
src.read() without a window is the whole band. This is the single most common way a windowed pipeline regresses, because the call looks almost identical to the correct one. A band of a large coverage is gigabytes; the window argument is what makes the read bounded.
Tiles overlap at their edges or they do not. The planner above produces abutting, non-overlapping tiles, which is correct for statistics and wrong for any focal operation — a slope or a filter needs a halo of pixels beyond the tile edge, or the result has seams. Add the halo deliberately when the operation needs it, and remember it means the tiles no longer sum to the whole.
Server-side limits are usually lower than your budget. Many deployments cap the pixel count of a single GetCoverage, returning an exception rather than a large response. Where that cap is lower than the tile size you computed, the run fails on the first tile — so reading the advertised limits from capabilities, or simply probing with one tile before planning the rest, saves a wasted batch.
Scale rather than tile when full resolution is not needed. scalefactor asks the server to resample before responding, which for a visualisation or an overview reduces the payload by the square of the factor. Tiling is for when the full resolution genuinely matters.
Because a single response has no restart point. Streaming to disk keeps memory flat but a connection dropped at ninety per cent of a twenty-gigabyte download costs the whole transfer. Tiling turns that into losing one tile, and it lets the fetch run concurrently with a bounded pool.
Large enough that per-request overhead is negligible and small enough that a failure is cheap. A couple of hundred megabytes is a reasonable default: at that size the request overhead is a rounding error and a retry costs seconds. Server-side pixel caps often decide it for you.
Not with a simple fold, since percentiles need the distribution rather than a running total. Accumulate a histogram with fixed bins, or use a streaming quantile structure such as a t-digest. Both keep memory bounded; what you cannot do is sort the full array, which is the operation windowing exists to avoid.
The windowed reading half does, and more so — a COG’s internal tiling means rasterio can fetch only the byte ranges it needs over HTTP, so the download step often disappears entirely. The tiling half is specific to WCS, where the server decides what a response contains and the client’s only lever is asking for less.
Back to WCS Coverage Service Fundamentals
Related