feat(extract): Küste, Straßenflächen, Brückenerkennung, Mittellinien

Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>

tools/extract/leonida/vectors.py

@@ -0,0 +1,97 @@
+import numpy as np
+from scipy import ndimage as ndi
+from shapely.geometry import LineString, Polygon
+from skimage.measure import find_contours
+from skimage.morphology import closing, disk, skeletonize
+
+from .classify import CLASSES
+from .coords import SLICE
+
+ROAD_CLASSES = ("strasse", "weg", "spekulativ")
+
+
+def _to_world(contour: np.ndarray) -> list[tuple[float, float]]:
+    """Konturpunkte (Zeile, Spalte) → Weltkoordinaten in Metern."""
+    return [(SLICE["x0"] + c, SLICE["y1"] - r) for r, c in contour]
+
+
+def mask_polygons(mask: np.ndarray, min_area: float, simplify: float) -> list[Polygon]:
+    padded = np.pad(mask, 1)
+    outers, holes = [], []
+    for contour in find_contours(padded.astype(float), 0.5):
+        local = contour - 1
+        ring = _to_world(local)
+        if len(ring) < 4:
+            continue
+        poly = Polygon(ring).buffer(0)
+        if poly.is_empty or poly.area < min_area:
+            continue
+        (outers if _signed_area(local) > 0 else holes).append(poly)
+    result = []
+    for outer in outers:
+        inner = [h.exterior.coords for h in holes if outer.contains(h)]
+        result.append(Polygon(outer.exterior.coords, inner).buffer(0).simplify(simplify))
+    return [p for p in result if p.geom_type == "Polygon" and p.area >= min_area]
+
+
+def _signed_area(contour: np.ndarray) -> float:
+    """Shoelace-Summe in Pixelkoordinaten (Zeile, Spalte).
+
+    find_contours liefert Außenkonturen und Löcher mit entgegengesetztem Umlaufsinn;
+    das Vorzeichen trennt sie zuverlässig. Eine Stichprobe der Maske am Mittelpunkt der
+    Kontur (naheliegende Alternative) versagt bei konzentrischen Formen wie einer Insel
+    mit Lagune: der Mittelpunkt der äußeren Kontur liegt dann selbst in der Lagune.
+    """
+    row, col = contour[:, 0], contour[:, 1]
+    return float(0.5 * np.sum(col * np.roll(row, -1) - np.roll(col, -1) * row))
+
+
+def land_polygons(labels: np.ndarray) -> list[Polygon]:
+    land = labels != CLASSES["wasser"]
+    land = ndi.binary_closing(ndi.binary_opening(land, disk(2)), disk(2))
+    return mask_polygons(land, min_area=50, simplify=1.0)
+
+
+def road_mask(labels: np.ndarray) -> np.ndarray:
+    mask = np.zeros(labels.shape, dtype=bool)
+    for name in ROAD_CLASSES:
+        mask |= labels == CLASSES[name]
+    # Schließen überbrückt die eingezeichneten Straßennamen.
+    return closing(mask, disk(3))
+
+
+BRIDGE_MARGIN_PX = 3  # Rand um die Straßenfläche; die Fläche selbst ist nie "wasser".
+
+
+def road_polygons(labels: np.ndarray) -> list[tuple[Polygon, bool]]:
+    water = labels == CLASSES["wasser"]
+    polys = mask_polygons(road_mask(labels), min_area=20, simplify=0.8)
+    out = []
+    for poly in polys:
+        r0, c0, r1, c1 = _bbox_px(poly)
+        r0 = max(r0 - BRIDGE_MARGIN_PX, 0)
+        c0 = max(c0 - BRIDGE_MARGIN_PX, 0)
+        r1 = min(r1 + BRIDGE_MARGIN_PX, water.shape[0] - 1)
+        c1 = min(c1 + BRIDGE_MARGIN_PX, water.shape[1] - 1)
+        window = water[r0:r1 + 1, c0:c1 + 1]
+        out.append((poly, bool(window.mean() > 0.3)))
+    return out
+
+
+def _bbox_px(poly: Polygon) -> tuple[int, int, int, int]:
+    minx, miny, maxx, maxy = poly.bounds
+    return (int(SLICE["y1"] - maxy), int(minx - SLICE["x0"]),
+            int(SLICE["y1"] - miny), int(maxx - SLICE["x0"]))
+
+
+def centerlines(labels: np.ndarray) -> list[LineString]:
+    """Skelett der Straßenflächen, nur zum Platzieren von Laternen und Palmen."""
+    skeleton = skeletonize(road_mask(labels))
+    lines = []
+    for contour in find_contours(skeleton.astype(float), 0.5):
+        pts = _to_world(contour)
+        if len(pts) >= 2:
+            line = LineString(pts).simplify(1.5)
+            if line.length >= 15:
+                lines.append(line)
+    return lines

tools/extract/tests/test_vectors.py

@@ -0,0 +1,33 @@
+import numpy as np
+from leonida.classify import CLASSES
+from leonida.vectors import land_polygons, mask_polygons, road_polygons
+
+def test_square_becomes_one_polygon_with_right_area():
+    mask = np.zeros((100, 100), dtype=bool)
+    mask[20:60, 30:80] = True                       # 40 x 50 Pixel = 2000 m²
+    polys = mask_polygons(mask, min_area=10, simplify=1.0)
+    assert len(polys) == 1
+    assert abs(polys[0].area - 2000) < 60           # Vereinfachung darf etwas abweichen
+
+def test_island_keeps_its_lagoon():
+    labels = np.full((80, 80), CLASSES["wasser"], dtype=np.uint8)
+    labels[10:70, 10:70] = CLASSES["land"]
+    labels[30:50, 30:50] = CLASSES["wasser"]        # Lagune
+    polys = land_polygons(labels)
+    assert len(polys) == 1
+    assert len(polys[0].interiors) == 1
+
+def test_road_over_water_is_flagged_as_bridge():
+    labels = np.full((60, 60), CLASSES["wasser"], dtype=np.uint8)
+    labels[0:10, :] = CLASSES["land"]
+    labels[25:31, :] = CLASSES["strasse"]           # Damm quer durchs Wasser
+    roads = road_polygons(labels)
+    assert len(roads) == 1
+    assert roads[0][1] is True
+
+def test_crossing_roads_form_a_single_polygon():
+    labels = np.full((60, 60), CLASSES["land"], dtype=np.uint8)
+    labels[28:32, :] = CLASSES["strasse"]
+    labels[:, 28:32] = CLASSES["strasse"]
+    roads = road_polygons(labels)
+    assert len(roads) == 1                          # Kreuzung, kein Zerfall in vier Arme