Water and 3MF parts

satprint.water

Water map for multi-color prints: which terrain grid cells are water.

Water shapes come from the water layer of the same OpenFreeMap vector tiles the buildings use (sea, rivers, lakes, ponds). Where the sea was flattened to 0 m, the flattened cells count as water too, so a coast still colors when the tiles are unavailable.

level_water(hm, tiles)

hm with each lake, pond, reservoir and sea surface flattened to its shore.

Elevation data over water is unreliable: radar returns and void filling leave noise, and a lake in a steep valley can read tens of meters above its shore, which prints as a raised plateau. Each water body is set to the 10th percentile of the ground along its shore, the lower shore, and only ever lowered, so a body the data already shows flat is unchanged.

Parameters:

Name	Type	Description	Default
hm	Heightmap	a heightmap with a bounding box.	required
tiles	list[tuple[int, int, int, bytes]]	vector tiles covering it, as :meth:VectorTileClient.tiles returns.	required

Returns:

Type	Description
Heightmap	a new heightmap, or hm itself when there is no still water.

Source code in satprint/water.py

def level_water(hm: Heightmap, tiles: list[tuple[int, int, int, bytes]]) -> Heightmap:
    """``hm`` with each lake, pond, reservoir and sea surface flattened to its shore.

    Elevation data over water is unreliable: radar returns and void filling
    leave noise, and a lake in a steep valley can read tens of meters above
    its shore, which prints as a raised plateau. Each water body is set to
    the 10th percentile of the ground along its shore, the lower shore, and
    only ever lowered, so a body the data already shows flat is unchanged.

    :param hm: a heightmap with a bounding box.
    :param tiles: vector tiles covering it, as :meth:`VectorTileClient.tiles` returns.
    :return: a new heightmap, or ``hm`` itself when there is no still water.
    """
    if hm.bbox is None:
        return hm
    rows, cols = hm.data.shape
    shapes = [
        poly
        for poly, props in vector_tile_features(tiles, "water")
        if props.get("class", "lake") in STILL_WATER
    ]
    if not shapes:
        return hm
    # one pixel per unit, x east and y north, so a node (r, c) is at (c, rows-1-r);
    # pieces cut at a tile edge are joined first so each body has one shore
    to_px = model_projection(hm.bbox, cols - 1, rows - 1)
    bodies = _polygons(unary_union([shapely.transform(p, to_px) for p in shapes]))
    # The tiles outline a coast far finer than the grid; buffering that is
    # slow and uses gigabytes, so keep the detail to a fraction of a node.
    bodies = [b for body in bodies for b in _polygons(body.simplify(0.25))]
    data = hm.data.copy()
    changed = False
    for body in bodies:
        x0, y0, x1, y1 = body.buffer(SHORE_PX[1], quad_segs=2).bounds
        c0, c1 = max(int(np.floor(x0)), 0), min(int(np.ceil(x1)), cols - 1)
        r0, r1 = (
            max(int(np.floor(rows - 1 - y1)), 0),
            min(int(np.ceil(rows - 1 - y0)), rows - 1),
        )
        if c1 < c0 or r1 < r0:
            continue
        gx, gy = np.meshgrid(np.arange(c0, c1 + 1), rows - 1 - np.arange(r0, r1 + 1))
        inside = shapely.contains_xy(body, gx, gy)
        shore = shapely.contains_xy(
            body.buffer(SHORE_PX[1], quad_segs=2), gx, gy
        ) & ~shapely.contains_xy(body.buffer(SHORE_PX[0], quad_segs=2), gx, gy)
        if inside.sum() < 4 or not shore.any():
            continue
        window = data[r0 : r1 + 1, c0 : c1 + 1]
        level = np.percentile(window[shore], 10)
        lowered = inside & (window > level)
        window[lowered] = level
        changed |= bool(lowered.any())
    return replace(hm, data=data) if changed else hm

multicolor_parts(polygons, bbox, relief_mm, width_mm, depth_mm, base_mm, sea_level_flat=False, buildings=None, frame=None, roads=None, bridges=None)

The 3MF parts, land, water, buildings, roads and border, and the water mask.

Bridges go in the roads part, so they print in the road color. Parts without geometry are left out; :func:write_3mf numbers the filaments by color in this order: land 1, water 2, then buildings and the gray of roads and border as present.

Returns:

Type	Description
tuple[list[tuple[str, Mesh, str]], ndarray]	([(name, mesh, color), ...] for :func:write_3mf, mask).

Source code in satprint/water.py

def multicolor_parts(
    polygons: list[Polygon],
    bbox: BBox,
    relief_mm: np.ndarray,
    width_mm: float,
    depth_mm: float,
    base_mm: float,
    sea_level_flat: bool = False,
    buildings: BuildingMesh | None = None,
    frame: Mesh | None = None,
    roads: BuildingMesh | None = None,
    bridges: BuildingMesh | None = None,
) -> tuple[list[tuple[str, Mesh, str]], np.ndarray]:
    """The 3MF parts, land, water, buildings, roads and border, and the water
    mask.

    Bridges go in the roads part, so they print in the road color. Parts
    without geometry are left out; :func:`write_3mf` numbers the filaments by
    color in this order: land 1, water 2, then buildings and the gray of roads
    and border as present.

    :return: ([(name, mesh, color), ...] for :func:`write_3mf`, mask).
    """
    mask = water_mask(polygons, bbox, relief_mm, width_mm, depth_mm, sea_level_flat)
    water, land = heightmap_split_solids(
        relief_mm, width_mm, depth_mm, base_mm, mask, skin_mm=WATER_SKIN_MM
    )
    parts = [
        ("land", land, PART_COLORS["land"]),
        ("water", water, PART_COLORS["water"]),
    ]
    if buildings is not None and buildings.count:
        parts.append(("buildings", buildings.as_mesh(), PART_COLORS["buildings"]))
    deck = merge_building_meshes(*(m for m in (roads, bridges) if m is not None))
    if deck.count:
        parts.append(("roads", deck.as_mesh(), PART_COLORS["roads"]))
    if frame is not None:
        parts.append(("border", frame, PART_COLORS["border"]))
    return parts, mask

water_from_vector_tiles(tiles)

Water polygons, in lon/lat, from the tiles' water layer.

Source code in satprint/water.py

def water_from_vector_tiles(tiles: list[tuple[int, int, int, bytes]]) -> list[Polygon]:
    """Water polygons, in lon/lat, from the tiles' ``water`` layer."""
    return [
        poly
        for poly, props in vector_tile_features(tiles, "water")
        if props.get("class", "lake") in WATER_CLASSES
    ]

water_mask(polygons, bbox, relief_mm, width_mm, depth_mm, sea_level_flat=False)

(rows-1, cols-1) boolean grid: True where a terrain cell is water.

Parameters:

Name	Type	Description	Default
polygons	list[Polygon]	water shapes in lon/lat.	required
bbox	BBox	area the block covers.	required
relief_mm	ndarray	the relief passed to :func:heightmap_to_mesh.	required
width_mm	float	block width.	required
depth_mm	float	block depth.	required
sea_level_flat	bool	the sea was flattened to the lowest relief, so cells lying flat there count as water.	False

Returns:

Type	Description
ndarray	the mask, row 0 north, as the heightmap.

Source code in satprint/water.py

def water_mask(
    polygons: list[Polygon],
    bbox: BBox,
    relief_mm: np.ndarray,
    width_mm: float,
    depth_mm: float,
    sea_level_flat: bool = False,
) -> np.ndarray:
    """(rows-1, cols-1) boolean grid: True where a terrain cell is water.

    :param polygons: water shapes in lon/lat.
    :param bbox: area the block covers.
    :param relief_mm: the relief passed to :func:`heightmap_to_mesh`.
    :param width_mm: block width.
    :param depth_mm: block depth.
    :param sea_level_flat: the sea was flattened to the lowest relief, so cells
        lying flat there count as water.
    :return: the mask, row 0 north, as the heightmap.
    """
    relief = np.asarray(relief_mm, dtype=np.float64)
    rows, cols = relief.shape
    xs = (np.arange(cols - 1) + 0.5) / (cols - 1) * width_mm
    ys = depth_mm - (np.arange(rows - 1) + 0.5) / (rows - 1) * depth_mm
    gx, gy = np.meshgrid(xs, ys)
    mask = np.zeros((rows - 1, cols - 1), dtype=bool)
    if polygons:
        to_model = model_projection(bbox, width_mm, depth_mm)
        water = shapely.union_all([shapely.transform(p, to_model) for p in polygons])
        shapely.prepare(water)
        mask = shapely.contains_xy(water, gx, gy)
    if sea_level_flat:
        lo = np.min(relief)
        corners = np.max(
            np.stack(
                [relief[:-1, :-1], relief[:-1, 1:], relief[1:, :-1], relief[1:, 1:]]
            ),
            axis=0,
        )
        mask |= corners <= lo + 1e-6
    return _fix_diagonals(mask)

water_zoom(bbox)

Highest vector-tile zoom (8-14) covering bbox in MAX_WATER_TILES.

A city reuses the zoom-14 tiles its buildings came from.

Source code in satprint/water.py

def water_zoom(bbox: BBox) -> int:
    """Highest vector-tile zoom (8-14) covering ``bbox`` in MAX_WATER_TILES.

    A city reuses the zoom-14 tiles its buildings came from.
    """
    zoom = 14
    while zoom > 8 and _tile_count(bbox, zoom) > MAX_WATER_TILES:
        zoom -= 1
    return zoom