Terrain and elevation

satprint.terrain

Elevation data sources and heightmap processing.

Sources

• fetch_terrarium — satellite-derived global elevation (SRTM/ASTER/GMTED merged) from the AWS Open Data "Terrain Tiles" set, in Mapzen terrarium PNG encoding. No API key is required.
• fetch_imagery -- Esri World Imagery tiles for the same bbox, used as the texture on the GLB export.
• synthetic_heightmap — procedural mountains for offline use and tests.
• load_heightmap_file — user-supplied PNG / TIFF / GeoTIFF heightmaps.

BBox(south, west, north, east) dataclass

ground_size_m()

(width, height) in meters along the bbox's central lines.

Source code in satprint/terrain.py

def ground_size_m(self) -> tuple[float, float]:
    """(width, height) in meters along the bbox's central lines."""
    w = haversine_m(self.mid_lat, self.west, self.mid_lat, self.east)
    h = haversine_m(self.south, self.mid_lon, self.north, self.mid_lon)
    return w, h

ImageryFetcher(cache_dir=None, session=None, url_template=IMAGERY_URL, timeout=30.0)

Bases: TileFetcher

Fetch + disk-cache satellite imagery tiles as (256, 256, 3) uint8 RGB.

Source code in satprint/terrain.py

def __init__(
    self,
    cache_dir: str | None = None,
    session: requests.Session | None = None,
    url_template: str = IMAGERY_URL,
    timeout: float = 30.0,
):
    super().__init__(
        cache_dir
        or os.path.join(os.path.expanduser("~"), ".cache", "satprint", "imagery"),
        session,
        url_template,
        timeout,
    )

TileFetcher(cache_dir=None, session=None, url_template=TERRARIUM_URL, timeout=30.0)

Fetch + disk-cache terrarium tiles.

Source code in satprint/terrain.py

def __init__(
    self,
    cache_dir: str | None = None,
    session: requests.Session | None = None,
    url_template: str = TERRARIUM_URL,
    timeout: float = 30.0,
):
    self.cache_dir = cache_dir or os.path.join(
        os.path.expanduser("~"), ".cache", "satprint", "tiles"
    )
    self.session = session or requests.Session()
    # Assign, not setdefault: a Session already carries python-requests'
    # own User-Agent, which some tile servers refuse.
    self.session.headers["User-Agent"] = "satprint/0.1 (+terrain relief models)"
    self.url_template = url_template
    self.timeout = timeout

choose_zoom(bbox, target_cols, max_zoom=MAX_ZOOM)

Smallest zoom at which the bbox spans at least target_cols pixels.

Source code in satprint/terrain.py

def choose_zoom(bbox: BBox, target_cols: int, max_zoom: int = MAX_ZOOM) -> int:
    """Smallest zoom at which the bbox spans at least ``target_cols`` pixels."""
    for z in range(0, max_zoom + 1):
        x0, _ = lonlat_to_global_px(bbox.west, bbox.north, z)
        x1, _ = lonlat_to_global_px(bbox.east, bbox.south, z)
        if x1 - x0 >= target_cols:
            return z
    return max_zoom

decode_terrarium(png_bytes)

Decode a terrarium PNG into float32 meters.

Source code in satprint/terrain.py

def decode_terrarium(png_bytes: bytes) -> np.ndarray:
    """Decode a terrarium PNG into float32 meters."""
    img = Image.open(io.BytesIO(png_bytes)).convert("RGB")
    rgb = np.asarray(img, dtype=np.float32)
    return rgb[..., 0] * 256.0 + rgb[..., 1] + rgb[..., 2] / 256.0 - 32768.0

encode_terrarium(elev)

Inverse of :func:decode_terrarium (used by tests and fixtures).

Source code in satprint/terrain.py

def encode_terrarium(elev: np.ndarray) -> bytes:
    """Inverse of :func:`decode_terrarium` (used by tests and fixtures)."""
    v = np.asarray(elev, dtype=np.float64) + 32768.0
    r = np.floor(v / 256.0)
    g = np.floor(v - r * 256.0)
    b = np.floor((v - r * 256.0 - g) * 256.0)
    rgb = np.stack([r, g, b], axis=-1).clip(0, 255).astype(np.uint8)
    buf = io.BytesIO()
    Image.fromarray(rgb, "RGB").save(buf, format="PNG")
    return buf.getvalue()

fetch_imagery(bbox, max_px=2048, fetcher=None, workers=8, progress=None)

Satellite image of bbox, cropped exactly like :func:fetch_terrarium.

The zoom is the smallest that reaches max_px columns, lowered until the area fits in MAX_TILES. The image is downscaled so its longer side is at most max_px. Row 0 is north, so grid node (r, c) of the heightmap sits at texture coordinates (c / (cols-1), r / (rows-1)).

Source code in satprint/terrain.py

def fetch_imagery(
    bbox: BBox,
    max_px: int = 2048,
    fetcher: ImageryFetcher | None = None,
    workers: int = 8,
    progress: Progress | None = None,
) -> tuple[Image.Image, dict]:
    """Satellite image of ``bbox``, cropped exactly like :func:`fetch_terrarium`.

    The zoom is the smallest that reaches ``max_px`` columns, lowered until
    the area fits in ``MAX_TILES``. The image is downscaled so its longer
    side is at most ``max_px``. Row 0 is north, so grid node (r, c) of the
    heightmap sits at texture coordinates (c / (cols-1), r / (rows-1)).
    """
    fetcher = fetcher or ImageryFetcher()
    zoom = choose_zoom(bbox, max_px, max_zoom=MAX_IMAGERY_ZOOM)
    while zoom > 0 and _tile_count(bbox, zoom) > MAX_TILES:
        zoom -= 1
    crop, n_tiles = _mosaic(bbox, zoom, fetcher.fetch, workers, progress, "imagery")
    img = Image.fromarray(np.ascontiguousarray(crop), "RGB")
    if max(img.size) > max_px:
        scale = max_px / max(img.size)
        size = (max(1, round(img.width * scale)), max(1, round(img.height * scale)))
        img = img.resize(size, Image.Resampling.LANCZOS)
    return img, {"zoom": zoom, "tiles": n_tiles, "px": [img.height, img.width]}

fetch_terrarium(bbox, target_cols=256, fetcher=None, workers=8, progress=None)

Download, mosaic and crop terrain tiles covering bbox.

The result is resampled so that its pixel aspect ratio matches the ground aspect ratio (width/height in meters), with target_cols columns.

Source code in satprint/terrain.py

def fetch_terrarium(
    bbox: BBox,
    target_cols: int = 256,
    fetcher: TileFetcher | None = None,
    workers: int = 8,
    progress: Progress | None = None,
) -> Heightmap:
    """Download, mosaic and crop terrain tiles covering ``bbox``.

    The result is resampled so that its pixel aspect ratio matches the
    ground aspect ratio (width/height in meters), with ``target_cols`` columns.
    """
    fetcher = fetcher or TileFetcher()
    zoom = choose_zoom(bbox, target_cols)
    crop, n_tiles = _mosaic(bbox, zoom, fetcher.fetch, workers, progress, "elevation")

    gw, gh = bbox.ground_size_m()
    rows = max(2, int(round(target_cols * gh / gw)))
    data = resample(crop, rows, target_cols)
    return Heightmap(
        data=data,
        ground_width_m=gw,
        ground_height_m=gh,
        source="terrarium",
        bbox=bbox,
        meta={"zoom": zoom, "tiles": n_tiles, "native_px": list(crop.shape)},
    )

gaussian_smooth(arr, sigma)

Separable Gaussian blur with edge replication (no SciPy needed).

Source code in satprint/terrain.py

def gaussian_smooth(arr: np.ndarray, sigma: float) -> np.ndarray:
    """Separable Gaussian blur with edge replication (no SciPy needed)."""
    if sigma <= 0:
        return arr
    radius = max(1, int(3 * sigma))
    k = np.exp(-0.5 * (np.arange(-radius, radius + 1) / sigma) ** 2)
    k /= k.sum()
    padded = np.pad(arr.astype(np.float64), radius, mode="edge")
    tmp = np.apply_along_axis(lambda m: np.convolve(m, k, mode="valid"), 1, padded)
    out = np.apply_along_axis(lambda m: np.convolve(m, k, mode="valid"), 0, tmp)
    return out.astype(np.float32)

heightmap_png(hm, clamp_sea_level=True)

Hillshaded 8-bit preview PNG of the heightmap.

Source code in satprint/terrain.py

def heightmap_png(hm: Heightmap, clamp_sea_level: bool = True) -> bytes:
    """Hillshaded 8-bit preview PNG of the heightmap."""
    elev = hm.data.astype(np.float64)
    if clamp_sea_level:
        elev = np.maximum(elev, 0)
    lo, hi = np.min(elev), np.max(elev)
    norm = (elev - lo) / max(hi - lo, 1e-6)
    gy, gx = np.gradient(elev)
    # light from the north-west, 45 deg elevation
    scale = 2.0 * (elev.shape[1] / hm.ground_width_m)
    nx, ny, nz = -gx * scale, gy * scale, np.ones_like(elev)
    n = np.sqrt(nx * nx + ny * ny + nz * nz)
    light = np.array([-0.5, 0.5, 0.7071])
    shade = np.clip((nx * light[0] + ny * light[1] + nz * light[2]) / n, 0, 1)
    img = np.clip(255 * (0.25 + 0.45 * norm + 0.3 * shade), 0, 255).astype(np.uint8)
    buf = io.BytesIO()
    Image.fromarray(img, "L").save(buf, format="PNG")
    return buf.getvalue()

load_heightmap_file(data, filename='', ground_width_m=None, max_cols=1024)

Load a heightmap image. GeoTIFFs use rasterio when available, otherwise the file is treated as a plain raster whose values are meters (16-bit PNGs are interpreted as meters as well, 8-bit as 0-255 relative units).

Source code in satprint/terrain.py

def load_heightmap_file(
    data: bytes,
    filename: str = "",
    ground_width_m: float | None = None,
    max_cols: int = 1024,
) -> Heightmap:
    """Load a heightmap image. GeoTIFFs use rasterio when available, otherwise
    the file is treated as a plain raster whose values are meters (16-bit PNGs
    are interpreted as meters as well, 8-bit as 0-255 relative units)."""
    arr: np.ndarray | None = None
    meta: dict = {"filename": filename}
    gw = ground_width_m
    gh = None

    if filename.lower().endswith((".tif", ".tiff")):
        try:
            # rasterio is the optional `geotiff` extra.
            from rasterio.io import MemoryFile  # ty: ignore[unresolved-import]

            with MemoryFile(data) as mem, mem.open() as ds:
                arr = ds.read(1).astype(np.float32)
                if ds.nodata is not None:
                    arr[arr == ds.nodata] = np.nan
                if ds.crs and ds.bounds:
                    b = ds.bounds
                    if ds.crs.is_geographic:
                        bbox = BBox(b.bottom, b.left, b.top, b.right)
                        gw, gh = bbox.ground_size_m()
                        meta["bbox"] = bbox.__dict__
                    else:
                        gw, gh = float(b.right - b.left), float(b.top - b.bottom)
                meta["georeferenced"] = gw is not None
        except ImportError:
            arr = None

    if arr is None:
        img = Image.open(io.BytesIO(data))
        if img.mode in ("I;16", "I;16B", "I", "F"):
            arr = np.asarray(img.convert("F"), dtype=np.float32)
        else:
            arr = np.asarray(img.convert("L"), dtype=np.float32)

    if arr.ndim != 2:
        raise ValueError("heightmap must be single-band")
    # Fill nodata holes with the minimum valid elevation.
    if np.isnan(arr).any():
        arr = np.where(np.isnan(arr), np.nanmin(arr), arr)
    rows, cols = arr.shape
    if cols > max_cols:
        rows = max(2, int(round(rows * max_cols / cols)))
        arr = resample(arr, rows, max_cols)
        cols = max_cols
    if gw is None:
        gw = 30.0 * cols  # assume ~30 m pixels (SRTM-like) if nothing better
        meta["assumed_pixel_m"] = 30.0
    if gh is None:
        gh = gw * rows / cols
    return Heightmap(
        data=arr.astype(np.float32),
        ground_width_m=float(gw),
        ground_height_m=float(gh),
        source="upload",
        meta=meta,
    )

lonlat_to_global_px(lon, lat, zoom)

Web-Mercator pixel coordinates (x right, y down) at zoom.

Source code in satprint/terrain.py

def lonlat_to_global_px(lon: float, lat: float, zoom: int) -> tuple[float, float]:
    """Web-Mercator pixel coordinates (x right, y down) at ``zoom``."""
    n = TILE_SIZE * (2**zoom)
    x = (lon + 180.0) / 360.0 * n
    lat_r = math.radians(lat)
    y = (1.0 - math.log(math.tan(lat_r) + 1.0 / math.cos(lat_r)) / math.pi) / 2.0 * n
    return x, y

prepare_relief(hm, p)

Convert elevation (m) to relief heights (mm above the base) and report the scale actually used.

Source code in satprint/terrain.py

def prepare_relief(hm: Heightmap, p: PrintParams) -> tuple[np.ndarray, dict]:
    """Convert elevation (m) to relief heights (mm above the base) and report
    the scale actually used."""
    elev = hm.data.astype(np.float32)
    if p.clamp_sea_level:
        elev = np.maximum(elev, 0.0)
    if p.smoothing > 0:
        elev = gaussian_smooth(elev, p.smoothing)
    lo, hi = float(np.min(elev)), float(np.max(elev))
    span = max(hi - lo, 1e-6)

    mm_per_m_plan = p.width_mm / hm.ground_width_m  # horizontal scale
    if p.relief_mm is not None:
        z_scale = p.relief_mm / span
        exaggeration = z_scale / mm_per_m_plan
    else:
        z_scale = mm_per_m_plan * p.exaggeration
        exaggeration = p.exaggeration
    relief = (elev - lo) * z_scale
    depth_mm = p.width_mm * hm.ground_height_m / hm.ground_width_m
    info = {
        "min_elev_m": lo,
        "max_elev_m": hi,
        "relief_m": span,
        "width_mm": p.width_mm,
        "depth_mm": depth_mm,
        "base_mm": p.base_mm,
        "relief_mm": float(np.max(relief)),
        "height_mm": float(np.max(relief)) + p.base_mm,
        "plan_scale": f"1:{int(round(1 / mm_per_m_plan * 1000)):,}",
        "mm_per_m_plan": mm_per_m_plan,
        "mm_per_m_vertical": z_scale,
        "exaggeration": exaggeration,
        "rows": int(elev.shape[0]),
        "cols": int(elev.shape[1]),
    }
    return relief.astype(np.float32), info

synthetic_heightmap(rows=200, cols=256, seed=0, ground_width_m=20000.0)

Procedural alpine-looking terrain with ridges, valleys and a lake.

Source code in satprint/terrain.py

def synthetic_heightmap(
    rows: int = 200, cols: int = 256, seed: int = 0, ground_width_m: float = 20_000.0
) -> Heightmap:
    """Procedural alpine-looking terrain with ridges, valleys and a lake."""
    rng = np.random.default_rng(seed)
    y, x = np.mgrid[0:rows, 0:cols].astype(np.float64)
    x /= cols
    y /= rows

    elev = np.zeros((rows, cols))
    # Fractal value noise: sum of upsampled random grids.
    amp, cells = 1.0, 4
    for _ in range(6):
        grid = rng.random((cells + 1, cells + 1)).astype(np.float32)
        layer = np.asarray(
            Image.fromarray(grid, "F").resize((cols, rows), Image.Resampling.BICUBIC)
        )
        elev += amp * (layer - 0.5)
        amp *= 0.5
        cells *= 2
    # A couple of big peaks.
    for _ in range(3):
        px, py, s = rng.random(), rng.random(), 0.12 + 0.1 * rng.random()
        elev += 1.2 * np.exp(-((x - px) ** 2 + (y - py) ** 2) / (2 * s * s))
    # Ridged transform for sharper crests.
    elev = 1.0 - np.abs(elev - elev.mean())
    elev = (elev - elev.min()) / (elev.max() - elev.min())
    elev = 400.0 + 2400.0 * elev**1.6
    # Flat lake in the lowest basin.
    lake = np.percentile(elev, 8)
    elev = np.where(elev < lake, lake, elev)
    gh = ground_width_m * rows / cols
    return Heightmap(
        data=elev.astype(np.float32),
        ground_width_m=ground_width_m,
        ground_height_m=gh,
        source="synthetic",
        meta={"seed": seed},
    )