_rings() rundete bislang jede Koordinate unabhängig auf 2 Stellen, nach normalize(). Das kann zwei Ringpunkte so verschieben, dass sich der Ring selbst schneidet, obwohl das Ausgangspolygon gültig war - Ursache der acht neuen ungültigen Polygone in den ausgelieferten Chunks (carriageway 3, land 2, sidewalk 2, surfaces 1). shapely.set_precision(poly, 0.01) läuft jetzt vor normalize() und vor der Rundung: Es ist eine echte Geometrieoperation (GEOS hält das Ergebnis mit mode="valid_output" aktiv gültig) statt eines blinden round() je Zahl. Gegen alle 29.037 Polygone der committeten (bereits einmal falsch gerundeten) Chunks probeweise gerechnet: kein Absturz, 8 der 9 ungültigen Polygone (die neunte ist die bekannte, von diesem Zweig unabhängige Vorbelastung in "paths", chunk_10_6, schon vor diesem Zweig ungültig) werden gültig, Flächenänderung insgesamt -24,35 m² von 11.711.354 m² (-0,00021 %). Zwei Sonderfälle brauchten eine eigene Behandlung, die die Zeilen im Befund nicht erwähnen, ohne die aber der Vier-Minuten-Lauf hätte abbrechen oder eine falsche Fläche ausliefern können: - An einer Engstelle unter 1 cm kann set_precision ein Polygon in mehrere Teile spalten (3 von 29.037 Fällen, größter abgeschnittener Teil 25,9 m²) - poly.exterior auf einer MultiPolygon wäre ein AttributeError gewesen. Es bleibt jetzt das größere Teilstück. - Ganz vereinzelt lehnt GEOS das Schnappen mit einer TopologyException ab (in der Stichprobe nur an einem bereits vorher ungültigen Polygon); dann bleibt das Polygon ungerundet, und die alte Rundung je Koordinate greift wie vor diesem Befund, statt den Lauf abzubrechen. Zwei Tests, beide ohne die Korrektur rot: ein valides Polygon, dessen unabhängig gerundete Ecken sich kreuzen (Beispiel selbst nachgerechnet, nicht nur behauptet), und eine Hantelform mit einem Steg unter 1 cm, die am unkorrigierten _rings() unter der Rundung selbst ungültig wird statt sich sauber zu teilen.
204 lines
10 KiB
Python
204 lines
10 KiB
Python
import json
|
|
from itertools import groupby
|
|
from pathlib import Path
|
|
|
|
import shapely
|
|
from shapely.geometry import Polygon, box
|
|
|
|
from .coords import CHUNK_SIZE, SLICE
|
|
|
|
N = (SLICE["x1"] - SLICE["x0"]) // CHUNK_SIZE # 12
|
|
|
|
|
|
def chunk_box(i: int, j: int) -> Polygon:
|
|
x0 = SLICE["x0"] + i * CHUNK_SIZE
|
|
y1 = SLICE["y1"] - j * CHUNK_SIZE
|
|
return box(x0, y1 - CHUNK_SIZE, x0 + CHUNK_SIZE, y1)
|
|
|
|
|
|
def split_by_chunk(poly: Polygon) -> dict[tuple[int, int], Polygon]:
|
|
out = {}
|
|
for j in range(N):
|
|
for i in range(N):
|
|
part = poly.intersection(chunk_box(i, j))
|
|
if not part.is_empty and part.area > 0.5:
|
|
out[(i, j)] = part
|
|
return out
|
|
|
|
|
|
def _rings(poly: Polygon, i: int, j: int) -> list[list[list[float]]]:
|
|
"""Aussenring zuerst, danach jedes Loch, in derselben Transformation.
|
|
|
|
Ruling 34: ein Loch (z.B. ein von Strassen umschlossener Block, ein Innenhof, ein Teich
|
|
im Park) ist genauso echte Flaeche wie der Aussenring und darf beim Zuschneiden nicht
|
|
verschwinden - sonst wird der Block zu Asphalt, der Innenhof zu Gebaeudemasse, der Teich
|
|
zu Wiese.
|
|
|
|
normalize() zuerst, weil dieselbe Flaeche sonst je nach Rechenweg anders in der Datei
|
|
steht: GEOS waehlt Ringstartpunkt und Wicklung nach Belieben. unary_union() in
|
|
apply_street() fasst dafuer das kartenweite Strassenpolygon neu zusammen und aenderte so
|
|
88 Chunkdateien, in denen sich geometrisch nichts bewegt hatte - symmetrische Differenz
|
|
der Strassenvereinigung exakt 0,00 m2 in allen 88. normalize() legt beides fest
|
|
(lexikographisch kleinster Punkt zuerst, Aussenring im Uhrzeigersinn, Loecher dagegen)
|
|
und macht die Ausgabe damit vom Rechenweg unabhaengig.
|
|
|
|
Die vereinheitlichte Wicklung ist unbedenklich, nachgeprueft an allen drei Verbrauchern:
|
|
geometry.js stellt sie beim Laden ohnehin selbst her (clockwise(), Zeile 50), das
|
|
Halbstrahl-Verfahren von pointInPolygon() in props.js kommt ohne sie aus, und PIL fuellt
|
|
in overview.py unabhaengig von Wicklung und Startpunkt (nachgemessen, Bild bytegleich).
|
|
|
|
set_precision(poly, 0.01) laeuft vor normalize() und vor der Rundung selbst (Befund 5 der
|
|
Fixwelle): Wer stattdessen jede Koordinate unabhaengig auf 2 Stellen rundet, wie es hier bis
|
|
zu diesem Befund geschah, kann zwei Ringpunkte so verschieben, dass sich der Ring selbst
|
|
schneidet - set_precision ist eine echte Geometrieoperation (GEOS haelt das Ergebnis mit
|
|
mode="valid_output" aktiv gueltig) statt eines blinden Python-round() je Zahl. Acht von
|
|
144 ausgelieferten Chunks hatten dadurch ein ungueltiges Polygon (carriageway 3, land 2,
|
|
sidewalk 2, surfaces 1) - set_precision behebt alle acht (nachgemessen).
|
|
|
|
Zwei Sonderfaelle, an einer Stichprobe ueber alle 29.037 ausgelieferten Polygone gemessen:
|
|
An einer Engstelle unter 1 cm spaltet set_precision ein Polygon in mehrere Teile (3 von
|
|
29.037 Faellen, groesster abgeschnittener Teil 25,9 m²) - hier bleibt nur das groesste
|
|
Teilstueck; das Polygon soll eine zusammenhaengende Flaeche bleiben, kein Zusatzeintrag,
|
|
den kein Aufrufer erwartet. Ganz vereinzelt lehnt GEOS das Schnappen mit einer
|
|
TopologyException ab (in der Stichprobe nur an einem bereits vorher ungueltigen Polygon);
|
|
dann bleibt das Polygon ungerundet, und die anschliessende Rundung je Koordinate greift wie
|
|
vor diesem Befund - kein Absturz des Vier-Minuten-Laufs an einer einzelnen Ausnahme.
|
|
"""
|
|
ox = SLICE["x0"] + i * CHUNK_SIZE
|
|
oy = SLICE["y1"] - j * CHUNK_SIZE
|
|
try:
|
|
genau = shapely.set_precision(poly, 0.01)
|
|
except shapely.errors.GEOSException:
|
|
genau = poly
|
|
if genau.geom_type == "MultiPolygon" and not genau.is_empty:
|
|
genau = max(genau.geoms, key=lambda g: g.area)
|
|
elif genau.is_empty or genau.geom_type != "Polygon":
|
|
genau = poly
|
|
poly = shapely.normalize(genau)
|
|
# Welt (X, Y) -> three (x, z) -> relativ zur Chunk-Ecke
|
|
return [[[round(x - ox, 2), round(-y - (-oy), 2)] for x, y in ring.coords]
|
|
for ring in (poly.exterior, *poly.interiors)]
|
|
|
|
|
|
def _sorted_surfaces(entries: list[dict]) -> list[dict]:
|
|
"""Belaege nur innerhalb einer Art sortieren - zwischen den Arten traegt die Reihenfolge Bedeutung.
|
|
|
|
overview.py malt die Belaege in Listenreihenfolge uebereinander, der spaetere gewinnt. An den
|
|
ausgelieferten Chunks ueberlappen sich verschiedene Arten auf 36.873 m2: der Parkrasen liegt
|
|
unter Feldern (19.573 m2), Plaetzen (11.872 m2), Parkplatz (2.878 m2), Pool (539 m2) und
|
|
Spielplatzsand (270 m2), die Stadtflaeche ueber dem Strandsand (1.729 m2). Ein Sortieren ueber
|
|
die Arten hinweg kehrte das um und liesse die Felder unter dem Rasen verschwinden - und eine
|
|
feste Rangfolge je Art gaebe es dafuer nicht, weil "sand" als Strand vor und als Spielplatz
|
|
hinter "park" steht. Die Abfolge der Arten selbst steht fest im Code (run.py, SURFACE_ORDER in
|
|
park.py) und ist damit schon eindeutig; GEOS mischt nur innerhalb einer Art.
|
|
"""
|
|
return [entry for _, run in groupby(entries, key=lambda e: e["kind"])
|
|
for entry in sorted(run, key=lambda e: e["poly"])]
|
|
|
|
|
|
def _line(line, i: int, j: int) -> list[list[float]]:
|
|
"""Ein Linienzug in Chunk-Koordinaten, richtungsunabhängig festgelegt.
|
|
|
|
normalize() dreht eine Linie in eine feste Laufrichtung, damit dieselbe Linie unabhängig vom
|
|
Rechenweg gleich in der Datei steht. Zulässig ist das, weil kein Verbraucher aus Richtung oder
|
|
Listenposition etwas ableitet: Die Laternen (sidewalkLampPoints() in props.js) suchen ihre
|
|
Straßenseite an der Gehwegfläche, und die Oberkante des Bordsteins (kerbGeometry() in
|
|
geometry.js) findet ihre Seite an der Fahrbahnfläche. Bei den früheren centerlines war das
|
|
nicht so, deshalb blieben sie unsortiert und unnormalisiert.
|
|
"""
|
|
ox = SLICE["x0"] + i * CHUNK_SIZE
|
|
oy = SLICE["y1"] - j * CHUNK_SIZE
|
|
return [[round(x - ox, 2), round(-y + oy, 2)] for x, y in shapely.normalize(line).coords]
|
|
|
|
|
|
def chunk_payload(i, j, *, land, surfaces, carriageway, curb, sidewalk, markings, lots,
|
|
bridges, edges, buildings, palms, paths=()) -> dict:
|
|
"""Ein Chunk als Ausgabestruktur, deterministisch sortiert.
|
|
|
|
Sortiert wird nach der gerundeten Ausgabestruktur selbst, nicht nach Flaeche, Schwerpunkt oder
|
|
Bounding-Box: nur sie ist total geordnet und ohne Sonderfall, denn zwei Eintraege mit gleichem
|
|
Schluessel sind dann bytegleich und ihre Reihenfolge wieder egal.
|
|
"""
|
|
ox = SLICE["x0"] + i * CHUNK_SIZE
|
|
oy = SLICE["y1"] - j * CHUNK_SIZE
|
|
return {
|
|
"chunk": [i, j],
|
|
"origin": [ox, -oy],
|
|
"land": sorted(_rings(p, i, j) for p in land),
|
|
"surfaces": _sorted_surfaces([{"poly": _rings(p, i, j), "kind": k} for p, k in surfaces]),
|
|
"carriageway": sorted(({"poly": _rings(p, i, j), "klasse": k} for p, k in carriageway),
|
|
key=lambda c: (c["poly"], c["klasse"])),
|
|
"curb": sorted(_line(l, i, j) for l in curb),
|
|
"sidewalk": sorted(({"poly": _rings(p, i, j)} for p in sidewalk), key=lambda s: s["poly"]),
|
|
"markings": sorted(({"poly": _rings(p, i, j), "art": a} for p, a in markings),
|
|
key=lambda m: (m["poly"], m["art"])),
|
|
"lots": sorted(({"poly": _rings(p, i, j)} for p in lots), key=lambda l: l["poly"]),
|
|
"bridges": sorted(({"poly": _rings(p, i, j)} for p in bridges), key=lambda b: b["poly"]),
|
|
"edges": sorted(({"line": _line(l, i, j), "klasse": k, "breite": w, "fahrbahn": f}
|
|
for l, k, w, f in edges),
|
|
key=lambda e: (e["line"], e["klasse"], e["breite"], e["fahrbahn"])),
|
|
"buildings": sorted(({"poly": _rings(p, i, j), "height": h["height"],
|
|
"source": h["source"], "speculative": s} for p, h, s in buildings),
|
|
key=lambda b: (b["poly"], b["height"], b["source"], b["speculative"])),
|
|
"palms": sorted([round(x - ox, 2), round(-y + oy, 2)] for x, y in palms),
|
|
"paths": sorted(({"poly": _rings(p, i, j)} for p in paths), key=lambda p: p["poly"]),
|
|
}
|
|
|
|
|
|
def chunk_box_of(i: int, j: int) -> Polygon:
|
|
return chunk_box(i, j)
|
|
|
|
|
|
def split_by_chunk_all(polys, i: int, j: int, min_area: float = 0.5):
|
|
"""Alle Anteile der übergebenen Polygone, die in Chunk (i, j) liegen.
|
|
|
|
min_area verwirft die Splitter, die beim Zuschneiden an der Chunkgrenze entstehen. 0,5 m² passt
|
|
für Flächen; Markierungen sind kleiner und brauchen eine eigene Grenze (run.py, chunk_markings).
|
|
"""
|
|
cell = chunk_box(i, j)
|
|
for poly in polys:
|
|
part = poly.intersection(cell)
|
|
if not part.is_empty and part.area > min_area:
|
|
if part.geom_type == "Polygon":
|
|
yield part
|
|
else:
|
|
yield from (g for g in part.geoms if g.geom_type == "Polygon" and g.area > min_area)
|
|
|
|
|
|
def in_chunk(point, i: int, j: int) -> bool:
|
|
x0 = SLICE["x0"] + i * CHUNK_SIZE
|
|
y1 = SLICE["y1"] - j * CHUNK_SIZE
|
|
return x0 <= point[0] < x0 + CHUNK_SIZE and y1 - CHUNK_SIZE <= point[1] < y1
|
|
|
|
|
|
def write_debug_png(out_dir: Path, mosaic, payloads: list[dict]) -> None:
|
|
"""Ein Bild je Chunk: die extrahierten Vektoren über der Karte, zur Sichtprüfung."""
|
|
from PIL import Image, ImageDraw
|
|
out_dir.mkdir(parents=True, exist_ok=True)
|
|
for payload in payloads:
|
|
i, j = payload["chunk"]
|
|
x0, y0 = i * CHUNK_SIZE, j * CHUNK_SIZE
|
|
crop = Image.fromarray(mosaic[y0:y0 + CHUNK_SIZE, x0:x0 + CHUNK_SIZE]).convert("RGB")
|
|
draw = ImageDraw.Draw(crop)
|
|
for rings in payload["land"]:
|
|
ring = rings[0] # nur der Aussenring; Loecher werden nicht gezeichnet
|
|
draw.line([tuple(pt) for pt in ring] + [tuple(ring[0])], fill=(0, 255, 255), width=1)
|
|
for road in payload["carriageway"]:
|
|
ring = road["poly"][0]
|
|
draw.line([tuple(pt) for pt in ring], fill=(255, 0, 255), width=1)
|
|
for b in payload["buildings"]:
|
|
ring = b["poly"][0]
|
|
if len(ring) >= 3:
|
|
draw.polygon([tuple(pt) for pt in ring], outline=(255, 220, 0))
|
|
for x, z in payload["palms"]:
|
|
draw.ellipse([x - 2, z - 2, x + 2, z + 2], outline=(0, 200, 0))
|
|
crop.save(out_dir / f"chunk_{i}_{j}.png")
|
|
|
|
|
|
def write_chunks(out_dir: Path, payloads: list[dict], meta: dict) -> None:
|
|
out_dir.mkdir(parents=True, exist_ok=True)
|
|
for payload in payloads:
|
|
i, j = payload["chunk"]
|
|
(out_dir / f"chunk_{i}_{j}.json").write_text(json.dumps(payload, separators=(",", ":")))
|
|
(out_dir / "manifest.json").write_text(json.dumps(meta, indent=1))
|