Roads

satprint.roads

OpenStreetMap roads -> one thin solid draped over the terrain.

Roads come from the transportation layer of the OpenFreeMap vector tiles the buildings and bridges use. Each line is drawn as a strip of a width per class (:data:satprint.bridges.BRIDGE_WIDTH_M), never narrower than a printable min_width_mm. At print scale a street is a fraction of a nozzle wide, so most roads are drawn wider than they are, and minor classes are left out of large areas, where they would cover the land (:data:ROAD_MIN_SCALE).

The strips are joined into one footprint, with water, building footprints and bridges cut out, and the footprint is cut along the triangles of a lattice of the terrain grid's nodes (:func:_lattice), every node on a coarse grid. Each piece lies on one flat lattice triangle, so its top sits a fixed height above the lattice surface; that height is height_mm plus however far the terrain rises above the lattice under the roads, so the top clears the terrain by at least height_mm and the floor sinks at least sink_mm into it. Pieces share their vertices, and walls stand only on the footprint's outer edges, as in :func:satprint.mesh.heightmap_split_solids, so each connected road network is a single closed solid. Tunnels are left out; bridges over water are :mod:satprint.bridges' job.

road_classes(detail, mm_per_m)

The classes drawn at detail for a model at mm_per_m.

Source code in satprint/roads.py

def road_classes(detail: Detail, mm_per_m: float) -> set[str]:
    """The classes drawn at ``detail`` for a model at ``mm_per_m``."""
    if detail == "major":
        return set(MAJOR_ROADS)
    if detail == "all":
        return set(ROAD_MIN_SCALE)
    return {cls for cls, s in ROAD_MIN_SCALE.items() if mm_per_m >= s}

road_mesh(lines, water, footprints, bbox, relief_mm, width_mm, depth_mm, base_mm, mm_per_m, detail='auto', height_mm=0.4, sink_mm=0.3, min_width_mm=0.8, simplify_mm=0.05, min_area_mm2=0.2, lattice_mm=0.6, lattice_tol_mm=0.3, bridges=None)

Roads as a thin solid on the terrain block built from relief_mm.

Parameters:

Name	Type	Description	Default
lines	list[tuple[LineString, str]]	(road line in lon/lat, OpenMapTiles class) pairs.	required
water	list[Polygon]	water polygons in lon/lat; roads stop at the water.	required
footprints	list[Polygon]	building footprints in lon/lat; roads stop at them.	required
bbox	BBox	area the terrain 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
base_mm	float	base thickness under the lowest terrain.	required
mm_per_m	float	horizontal model scale.	required
detail	Detail	which classes to draw (:func:road_classes).	'auto'
height_mm	float	height of the road top above the terrain.	0.4
sink_mm	float	depth of the road floor below the terrain.	0.3
min_width_mm	float	narrowest a road is drawn.	0.8
simplify_mm	float	footprint simplification tolerance.	0.05
min_area_mm2	float	connected pieces smaller than this are dropped.	0.2
lattice_mm	float	roads are cut on a lattice of terrain grid nodes about this far apart, not on every terrain cell (:func:_lattice).	0.6
lattice_tol_mm	float	the lattice is made finer until the terrain under the roads strays no more than this from it. Roads are raised and sunk by what strays, so they never sink into the terrain or float.	0.3
bridges	BuildingMesh | None	bridge solids from :func:satprint.bridges.bridge_mesh; roads stop at their outline, so a road meets a bridge's landing instead of running over it.	None

Returns:

Type	Description
BuildingMesh	the roads as one :class:BuildingMesh, the tops as roofs; count is the number of connected road networks.

Source code in satprint/roads.py

def road_mesh(
    lines: list[tuple[LineString, str]],
    water: list[Polygon],
    footprints: list[Polygon],
    bbox: BBox,
    relief_mm: np.ndarray,
    width_mm: float,
    depth_mm: float,
    base_mm: float,
    mm_per_m: float,
    detail: Detail = "auto",
    height_mm: float = 0.4,
    sink_mm: float = 0.3,
    min_width_mm: float = 0.8,
    simplify_mm: float = 0.05,
    min_area_mm2: float = 0.2,
    lattice_mm: float = 0.6,
    lattice_tol_mm: float = 0.3,
    bridges: BuildingMesh | None = None,
) -> BuildingMesh:
    """Roads as a thin solid on the terrain block built from ``relief_mm``.

    :param lines: (road line in lon/lat, OpenMapTiles class) pairs.
    :param water: water polygons in lon/lat; roads stop at the water.
    :param footprints: building footprints in lon/lat; roads stop at them.
    :param bbox: area the terrain block covers.
    :param relief_mm: the relief passed to :func:`heightmap_to_mesh`.
    :param width_mm: block width.
    :param depth_mm: block depth.
    :param base_mm: base thickness under the lowest terrain.
    :param mm_per_m: horizontal model scale.
    :param detail: which classes to draw (:func:`road_classes`).
    :param height_mm: height of the road top above the terrain.
    :param sink_mm: depth of the road floor below the terrain.
    :param min_width_mm: narrowest a road is drawn.
    :param simplify_mm: footprint simplification tolerance.
    :param min_area_mm2: connected pieces smaller than this are dropped.
    :param lattice_mm: roads are cut on a lattice of terrain grid nodes about
        this far apart, not on every terrain cell (:func:`_lattice`).
    :param lattice_tol_mm: the lattice is made finer until the terrain under
        the roads strays no more than this from it. Roads are raised and
        sunk by what strays, so they never sink into the terrain or float.
    :param bridges: bridge solids from :func:`satprint.bridges.bridge_mesh`;
        roads stop at their outline, so a road meets a bridge's landing
        instead of running over it.
    :return: the roads as one :class:`BuildingMesh`, the tops as roofs;
        ``count`` is the number of connected road networks.
    """
    relief = np.asarray(relief_mm, dtype=np.float64)
    if np.min(relief) < 0:  # match heightmap_to_mesh
        relief = relief - np.min(relief)
    rows, cols = relief.shape
    to_model = model_projection(bbox, width_mm, depth_mm)
    inset = 0.05  # keep walls off the block's own walls
    frame = box(inset, inset, width_mm - inset, depth_mm - inset)

    keep = road_classes(detail, mm_per_m)
    kept = [(line, cls) for line, cls in lines if cls in keep]
    if not kept:
        return BuildingMesh.empty()
    strips = shapely.buffer(
        np.array([shapely.transform(line, to_model) for line, _ in kept]),
        np.array(
            [
                max(BRIDGE_WIDTH_M.get(cls, DEFAULT_WIDTH_M) * mm_per_m, min_width_mm)
                / 2
                for _, cls in kept
            ]
        ),
        quad_segs=2,
    )
    footprint = shapely.union_all(strips).intersection(frame)
    cut = [shapely.transform(p, to_model) for p in [*water, *footprints]]
    if bridges is not None and bridges.count:
        # The bridges' outline in plan: their tops, seen from above.
        tris = shapely.polygons(bridges.vertices[bridges.roof_faces][:, :, :2])
        cut.append(shapely.union_all(tris[shapely.area(tris) > 0]))
    if cut:
        footprint = footprint.difference(shapely.union_all(cut))
    footprint = shapely.make_valid(footprint.simplify(simplify_mm))
    # Shrinking by a hair parts outlines that touch at a point, where four
    # walls would share one edge.
    footprint = footprint.buffer(-1e-4, join_style="mitre")
    parts = [p for p in _polygons(footprint) if p.area >= min_area_mm2]
    if not parts:
        return BuildingMesh.empty()
    footprint = shapely.MultiPolygon(parts)

    # Grid node (r, c) is at (xs[c], ys[r]); cell (r, c) is split along its
    # NW-SE diagonal into the two triangles heightmap_to_mesh draws. Roads
    # are cut on a coarser lattice of those nodes, cut the same way.
    xs = np.linspace(0.0, width_mm, cols)
    ys = np.linspace(depth_mm, 0.0, rows)
    ir, ic, above, below = _lattice(
        relief, xs, ys, footprint, lattice_mm, lattice_tol_mm
    )
    xs, ys, relief = xs[ic], ys[ir], relief[np.ix_(ir, ic)]
    pieces, planes = _pieces(footprint, xs, ys)
    if not pieces:
        return BuildingMesh.empty()

    index: dict[tuple[float, float], int] = {}
    xyz: list[tuple[float, float, float]] = []
    tops: list[list[int]] = []
    for piece, (r, c, upper) in zip(pieces, planes, strict=True):
        # Triangle corners: NW, SE and SW (lower) or NE (upper).
        corner = (r, c + 1) if upper else (r + 1, c)
        p = np.array(
            [
                [xs[c], ys[r], relief[r, c]],
                [xs[c + 1], ys[r + 1], relief[r + 1, c + 1]],
                [xs[corner[1]], ys[corner[0]], relief[corner]],
            ]
        )
        normal = np.cross(p[1] - p[0], p[2] - p[0])

        def vid(x: float, y: float, p=p, normal=normal) -> int:
            key = (round(x, KEY_DECIMALS), round(y, KEY_DECIMALS))
            i = index.get(key)
            if i is None:
                z = (
                    p[0, 2]
                    - (normal[0] * (x - p[0, 0]) + normal[1] * (y - p[0, 1]))
                    / normal[2]
                )
                i = index[key] = len(xyz)
                xyz.append((x, y, z))
            return i

        for tri in shapely.constrained_delaunay_triangles(piece).geoms:
            pts = np.asarray(tri.exterior.coords)[:3]
            a, b, c3 = pts
            if (b[0] - a[0]) * (c3[1] - a[1]) - (b[1] - a[1]) * (c3[0] - a[0]) <= 0:
                pts = pts[[0, 2, 1]]
            ids = [vid(float(x), float(y)) for x, y in pts]
            if len(set(ids)) == 3:
                tops.append(ids)
    if not tops:
        return BuildingMesh.empty()

    # The top clears the terrain by height_mm where it rises most above the
    # lattice, and the floor sinks sink_mm where it dips most below it.
    top = np.asarray(xyz)
    top[:, 2] += base_mm + height_mm + above
    bottom = top.copy()
    bottom[:, 2] -= height_mm + sink_mm + above + below
    n = top.shape[0]
    roof = np.asarray(tops, dtype=np.int64)
    # Boundary edges are the directed top edges whose reverse is absent;
    # each gets a wall facing out, as in heightmap_split_solids.
    u = roof.ravel()
    v = np.roll(roof, -1, axis=1).ravel()
    fwd = u * (2 * n) + v
    edge = ~np.isin(fwd, v * (2 * n) + u)
    u, v = u[edge], v[edge]
    walls = np.concatenate(
        [np.column_stack([u + n, v + n, v]), np.column_stack([u + n, v, u])]
    )
    floor = roof[:, ::-1] + n
    return BuildingMesh(
        vertices=np.vstack([top, bottom]).astype(np.float32),
        roof_faces=roof,
        wall_faces=np.vstack([walls, floor]),
        count=len(parts),
    )

roads_from_vector_tiles(tiles)

Road lines, in lon/lat, from the tiles' transportation layer.

Tunnels are left out, and so are classes not in :data:ROAD_MIN_SCALE, such as rail, ferries and roads under construction.

Returns:

Type	Description
list[tuple[LineString, str]]	(line, OpenMapTiles class) pairs.

Source code in satprint/roads.py

def roads_from_vector_tiles(
    tiles: list[tuple[int, int, int, bytes]],
) -> list[tuple[LineString, str]]:
    """Road lines, in lon/lat, from the tiles' ``transportation`` layer.

    Tunnels are left out, and so are classes not in :data:`ROAD_MIN_SCALE`,
    such as rail, ferries and roads under construction.

    :return: (line, OpenMapTiles class) pairs.
    """
    return [
        (line, props["class"])
        for line, props in vector_tile_features(tiles, "transportation", lines=True)
        if props.get("class") in ROAD_MIN_SCALE and props.get("brunnel") != "tunnel"
    ]