karim
afb59b6626
Swisstopo
- swissBUILDINGS3D 3.0 + Variant-Toggle (separated/solid) im Dialog
- Auto-Fallback auf 2.0 wenn 3.0-Tiles ueber 200 MB sind (Stadt-Fall)
- Defensiver Variant-Filter auf 3 Ebenen (Item, Asset, ZIP-Extract) — keine
Doppelimporte mehr
- Auto-Skala korrigiert jetzt die importierten Objekte (×1000) statt die
User-bbox zu schrumpfen — Buildings bleiben in m-Doc-Skala
- merge_grids: XYZ-Tiles werden vor dem Mesh-Bau vereint, kein 1m-Streifen
zwischen Tiles mehr
- Layer-Konsolidierung: Build_*/Roof_*/Wall_*/Floor_* DWG-Source-Layer
werden auf Sub-Sub-Layer unter 81_Swissbuildings/{Build,Roof,Wall,Floor}
gemappt; solid-Variante landet flach direkt auf dem Parent
- 0-Kote m.ü.M (Projekt-Nullpunkt) wird beim Import als Z-Offset angewandt
Hierarchische Ebenen
- dossier_ebenen unterstuetzt jetzt 'children'-Array (rekursiv)
- layer_builder.build_layers rekursiv (Parent + Children unter jedem Geschoss)
- apply_visibility/update_layer_style/set_ebene_visible/set_ebene_locked
walken den Tree (Sub-Sub-Layer mit gleichem Code-Prefix werden mit-gepflegt)
- EbenenManager mit Chevron-Toggle + Indent pro Level + Context-Menue-Item
'Sub-Ebene hinzufuegen'
- rhinoBridge.applyVisibility schickt Children-Tree (nicht nur Top-Level) —
sonst kommen Sub-Toggles nicht beim Backend an
- Visibility-Key in App.jsx rekursiv durch Children — useEffect feuert jetzt
auch bei Sub-Eye-Toggles
0-Kote m.ü.M
- Eingabefeld im Geschoss-Settings-Dialog (projektweit)
- Speicherung als dossier_project_zero_mum in doc.Strings
- Wird im Swisstopo-Import als Z-Offset (m + doc-units) angewandt
Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
801 lines
33 KiB
Python
801 lines
33 KiB
Python
#! python 3
|
||
# -*- coding: utf-8 -*-
|
||
"""
|
||
swisstopo.py
|
||
STAC-API-Client + GeoTIFF/XYZ-Parser + Mesh-Builder fuer swisstopo-Daten.
|
||
Alle APIs sind offen + ohne Auth, ohne Key.
|
||
|
||
Collections die wir nutzen:
|
||
- ch.swisstopo.swissbuildings3d_3 : 3D-Gebaeudemodelle (DWG/OBJ/DAE/IFC)
|
||
- ch.swisstopo.swissalti3d : Hoehenmodell DTM (GeoTIFF/XYZ)
|
||
- ch.swisstopo.swissimage-dop10 : Orthofoto 10 cm (GeoTIFF, RGB)
|
||
"""
|
||
import os
|
||
import re
|
||
import json
|
||
import zipfile
|
||
import urllib.request
|
||
import urllib.parse
|
||
import Rhino
|
||
import Rhino.Geometry as rg
|
||
|
||
CACHE_DIR = os.path.expanduser("~/Library/Caches/Dossier/swisstopo")
|
||
STAC_BASE = "https://data.geo.admin.ch/api/stac/v1"
|
||
SEARCH_API = "https://api3.geo.admin.ch/rest/services/api/SearchServer"
|
||
|
||
|
||
def _ensure_cache():
|
||
if not os.path.isdir(CACHE_DIR):
|
||
try: os.makedirs(CACHE_DIR)
|
||
except Exception as ex: print("[SWISSTOPO] cache mkdir:", ex)
|
||
|
||
|
||
def _http_get_json(url, timeout=30):
|
||
"""HTTP GET + JSON-decode. Wirft bei Fehler."""
|
||
req = urllib.request.Request(url, headers={"Accept": "application/json"})
|
||
with urllib.request.urlopen(req, timeout=timeout) as resp:
|
||
return json.loads(resp.read().decode("utf-8"))
|
||
|
||
|
||
_MAX_SAFE_DOWNLOAD_MB = 200 # ueber diesem Wert: Download skippen + Warn
|
||
|
||
|
||
def _yield_ui():
|
||
"""Gibt der UI eine Chance zu repainten waehrend einer langen Operation.
|
||
Ohne diesen Aufruf friert die WebView ein bis _run_import fertig ist."""
|
||
try: Rhino.RhinoApp.Wait()
|
||
except Exception: pass
|
||
|
||
|
||
def _http_download(url, dest_path, timeout=120, progress=None, status=None):
|
||
"""Download URL → dest_path. Skippt wenn Datei schon im Cache existiert
|
||
(gleicher Name, > 0 Bytes). Liefert True bei Erfolg. Bricht ab wenn
|
||
Content-Length > _MAX_SAFE_DOWNLOAD_MB (Schutz vor versehentlichen
|
||
Mega-Downloads wie der gesamt-CH-GDB).
|
||
progress: optional callable(bytes_done, bytes_total).
|
||
status: optional callable(msg) fuer Log-Updates."""
|
||
if os.path.isfile(dest_path) and os.path.getsize(dest_path) > 0:
|
||
if status: status("Cache: {} ({:.1f} MB)".format(
|
||
os.path.basename(dest_path), os.path.getsize(dest_path) / 1e6))
|
||
return True
|
||
try:
|
||
req = urllib.request.Request(url, headers={"User-Agent": "Dossier/0.6"})
|
||
with urllib.request.urlopen(req, timeout=timeout) as resp:
|
||
total = int(resp.headers.get("Content-Length") or -1)
|
||
if total > 0:
|
||
size_mb = total / 1e6
|
||
if status: status("Download {:.1f} MB...".format(size_mb))
|
||
if size_mb > _MAX_SAFE_DOWNLOAD_MB:
|
||
if status: status("→ Abgebrochen: > {} MB (vermutlich kein "
|
||
"Per-Tile-Asset)".format(_MAX_SAFE_DOWNLOAD_MB))
|
||
return False
|
||
done = 0
|
||
last_log_mb = 0
|
||
with open(dest_path + ".part", "wb") as f:
|
||
while True:
|
||
chunk = resp.read(65536)
|
||
if not chunk: break
|
||
f.write(chunk)
|
||
done += len(chunk)
|
||
if progress:
|
||
try: progress(done, total)
|
||
except Exception: pass
|
||
# Alle 2 MB einen Status + UI-Yield
|
||
mb_done = done // (2 * 1024 * 1024)
|
||
if status and mb_done > last_log_mb:
|
||
last_log_mb = mb_done
|
||
if total > 0:
|
||
pct = 100.0 * done / total
|
||
status(" {:.0f}% ({:.1f}/{:.1f} MB)".format(
|
||
pct, done/1e6, total/1e6))
|
||
else:
|
||
status(" {:.1f} MB...".format(done/1e6))
|
||
_yield_ui()
|
||
os.rename(dest_path + ".part", dest_path)
|
||
return True
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] download {}: {}".format(url, ex))
|
||
try:
|
||
if os.path.isfile(dest_path + ".part"): os.remove(dest_path + ".part")
|
||
except Exception: pass
|
||
return False
|
||
|
||
|
||
# --- Koordinaten-Transformationen ------------------------------------------
|
||
# Schweizer Koordinaten (LV95 = EPSG:2056) zu/von WGS84.
|
||
# Approximations-Formeln aus dem swisstopo-Dokument
|
||
# "Naeherungsloesungen fuer die direkte Transformation CH1903<->WGS84".
|
||
# Genau auf ca. 1m im gesamten Schweizer Gebiet — fuer unsere Zwecke ueppig.
|
||
|
||
def lv95_to_wgs84(e, n):
|
||
"""LV95 (E, N) -> WGS84 (lon, lat)."""
|
||
e_n = (e - 2600000) / 1_000_000.0
|
||
n_n = (n - 1200000) / 1_000_000.0
|
||
lon_p = (2.6779094
|
||
+ 4.728982 * e_n
|
||
+ 0.791484 * e_n * n_n
|
||
+ 0.1306 * e_n * n_n * n_n
|
||
- 0.0436 * e_n * e_n * e_n)
|
||
lat_p = (16.9023892
|
||
+ 3.238272 * n_n
|
||
- 0.270978 * e_n * e_n
|
||
- 0.002528 * n_n * n_n
|
||
- 0.0447 * e_n * e_n * n_n
|
||
- 0.0140 * n_n * n_n * n_n)
|
||
lon = lon_p * 100 / 36
|
||
lat = lat_p * 100 / 36
|
||
return (lon, lat)
|
||
|
||
|
||
def wgs84_to_lv95(lon, lat):
|
||
"""WGS84 (lon, lat) -> LV95 (E, N)."""
|
||
lon_p = (lon * 3600 - 26782.5) / 10000.0
|
||
lat_p = (lat * 3600 - 169028.66) / 10000.0
|
||
e = (2600072.37
|
||
+ 211455.93 * lon_p
|
||
- 10938.51 * lon_p * lat_p
|
||
- 0.36 * lon_p * lat_p * lat_p
|
||
- 44.54 * lon_p * lon_p * lon_p)
|
||
n = (1200147.07
|
||
+ 308807.95 * lat_p
|
||
+ 3745.25 * lon_p * lon_p
|
||
+ 76.63 * lat_p * lat_p
|
||
- 194.56 * lon_p * lon_p * lat_p
|
||
+ 119.79 * lat_p * lat_p * lat_p)
|
||
return (e, n)
|
||
|
||
|
||
def lv95_bbox_to_wgs84_bbox(e_min, n_min, e_max, n_max):
|
||
"""4 Ecken transformieren + min/max in WGS84 zurueck."""
|
||
corners = [
|
||
lv95_to_wgs84(e_min, n_min),
|
||
lv95_to_wgs84(e_max, n_min),
|
||
lv95_to_wgs84(e_max, n_max),
|
||
lv95_to_wgs84(e_min, n_max),
|
||
]
|
||
lons = [c[0] for c in corners]
|
||
lats = [c[1] for c in corners]
|
||
return (min(lons), min(lats), max(lons), max(lats))
|
||
|
||
|
||
# --- Geocoding -------------------------------------------------------------
|
||
|
||
def geocode(text):
|
||
"""Sucht Adresse via swisstopo Search-API. Liefert dict mit:
|
||
{label, e, n, lon, lat, origin} oder None.
|
||
"""
|
||
if not text or not text.strip():
|
||
return None
|
||
try:
|
||
params = {
|
||
"searchText": text.strip(),
|
||
"type": "locations",
|
||
"origins": "address,gg25,gazetteer,parcel",
|
||
"sr": "2056", # LV95-Koords im Response
|
||
"limit": "1",
|
||
}
|
||
url = SEARCH_API + "?" + urllib.parse.urlencode(params)
|
||
data = _http_get_json(url, timeout=15)
|
||
results = data.get("results") or []
|
||
if not results: return None
|
||
attrs = results[0].get("attrs") or {}
|
||
# Im LV95-Modus liefert die API y=East, x=North (Geo-Admin-Konvention).
|
||
e = attrs.get("y")
|
||
n = attrs.get("x")
|
||
label = attrs.get("label") or attrs.get("detail") or text
|
||
if e is None or n is None: return None
|
||
e = float(e); n = float(n)
|
||
lon, lat = lv95_to_wgs84(e, n)
|
||
return {
|
||
"label": _strip_html(label),
|
||
"e": e,
|
||
"n": n,
|
||
"lon": lon,
|
||
"lat": lat,
|
||
"origin": attrs.get("origin"),
|
||
}
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] geocode '{}':".format(text), ex)
|
||
return None
|
||
|
||
|
||
def _strip_html(s):
|
||
"""Search-API liefert Labels mit <b>-Tags fuer Highlighting."""
|
||
try:
|
||
import re
|
||
return re.sub(r"<[^>]+>", "", s or "")
|
||
except Exception:
|
||
return s
|
||
|
||
|
||
# --- STAC-Queries ----------------------------------------------------------
|
||
|
||
def stac_query(collection_id, bbox_wgs84, asset_extensions=None, limit=100):
|
||
"""STAC-Items in einer WGS84-bbox. Liefert Liste von:
|
||
{id, bbox_wgs84, assets: {key: {href, type, size?}}}
|
||
asset_extensions: optional list von suffixen (z.B. ['.dwg', '.tif']) zum
|
||
filtern. Default: alle Assets behalten.
|
||
"""
|
||
if not bbox_wgs84 or len(bbox_wgs84) != 4: return []
|
||
try:
|
||
url = ("{base}/collections/{cid}/items"
|
||
"?bbox={mn_lon},{mn_lat},{mx_lon},{mx_lat}"
|
||
"&limit={lim}").format(
|
||
base=STAC_BASE, cid=collection_id,
|
||
mn_lon=bbox_wgs84[0], mn_lat=bbox_wgs84[1],
|
||
mx_lon=bbox_wgs84[2], mx_lat=bbox_wgs84[3],
|
||
lim=int(limit))
|
||
data = _http_get_json(url, timeout=30)
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] STAC '{}': {}".format(collection_id, ex))
|
||
return []
|
||
out = []
|
||
for feat in (data.get("features") or []):
|
||
assets_raw = feat.get("assets") or {}
|
||
assets = {}
|
||
for k, v in assets_raw.items():
|
||
href = v.get("href")
|
||
if not href: continue
|
||
if asset_extensions:
|
||
low = href.lower()
|
||
if not any(low.endswith(ext) for ext in asset_extensions): continue
|
||
assets[k] = {
|
||
"href": href,
|
||
"type": v.get("type"),
|
||
"size": v.get("file:size"),
|
||
}
|
||
if not assets: continue
|
||
out.append({
|
||
"id": feat.get("id"),
|
||
"bbox_wgs84": feat.get("bbox"),
|
||
"assets": assets,
|
||
})
|
||
return out
|
||
|
||
|
||
# --- Download helpers ------------------------------------------------------
|
||
|
||
def download_asset(href, subdir="misc", progress=None, status=None):
|
||
"""Laedt asset, liefert lokalen Pfad oder None bei Fehler.
|
||
status: optional callable(msg) fuer Live-Progress-Log."""
|
||
_ensure_cache()
|
||
sub = os.path.join(CACHE_DIR, subdir)
|
||
if not os.path.isdir(sub):
|
||
try: os.makedirs(sub)
|
||
except Exception: pass
|
||
fn = os.path.basename(urllib.parse.urlparse(href).path) or "asset"
|
||
dest = os.path.join(sub, fn)
|
||
if _http_download(href, dest, progress=progress, status=status):
|
||
return dest
|
||
return None
|
||
|
||
|
||
# --- Helpers: Tile-Dedupe + ZIP-Extract ------------------------------------
|
||
|
||
# Item-IDs haben verschiedene Year-Patterns:
|
||
# swissalti3d: "_YYYY_NNNN-NNNN" (z.B. _2022_2664-1212)
|
||
# swissbuildings3d_2: "_YYYY-MM_NNNN-NN" (z.B. _2022-12_1150-23)
|
||
# Beide gleichzeitig matchen — Year-Sort dann via String-Compare (lex == chrono).
|
||
_TILE_PAT = re.compile(r"_(\d{4}(?:-\d{2})?)_(\d{3,4}-\d{2,4})")
|
||
# Filter fuer per-Tile Assets — die URL/Filename MUSS eine Tile-Coord-
|
||
# Markierung haben (z.B. `_1150-23_`). Bewahrt vor versehentlichem Download
|
||
# der gesamt-CH-Geodatabase (>1 GB).
|
||
_TILE_FILE_PAT = re.compile(r"_\d{3,4}-\d{2,4}[_.]")
|
||
|
||
|
||
def _dedupe_latest(items):
|
||
"""Bei Items mit gleicher Tile-Coord aber unterschiedlicher Year-Markierung:
|
||
nur das neueste behalten (Year-String-Compare: '2022-12' > '2020-03' ✓)."""
|
||
keep = {}
|
||
for it in items:
|
||
m = _TILE_PAT.search(it.get("id") or "")
|
||
if not m:
|
||
keep[it["id"]] = it
|
||
continue
|
||
year = m.group(1); tile = m.group(2)
|
||
prev = keep.get(tile)
|
||
if prev is None or year > prev["_year"]:
|
||
it["_year"] = year
|
||
keep[tile] = it
|
||
return list(keep.values())
|
||
|
||
|
||
def _extract_zip_to_dir(zip_path, dest_dir):
|
||
"""Entpackt alle Files aus einem ZIP nach dest_dir. Liefert Liste der
|
||
extrahierten Pfade."""
|
||
if not os.path.isdir(dest_dir):
|
||
try: os.makedirs(dest_dir)
|
||
except Exception: pass
|
||
paths = []
|
||
try:
|
||
with zipfile.ZipFile(zip_path, "r") as zf:
|
||
for name in zf.namelist():
|
||
if name.endswith("/"): continue # Verzeichnis
|
||
out = os.path.join(dest_dir, os.path.basename(name))
|
||
if not os.path.isfile(out) or os.path.getsize(out) == 0:
|
||
with zf.open(name) as src, open(out, "wb") as dst:
|
||
dst.write(src.read())
|
||
paths.append(out)
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] zip extract {}: {}".format(zip_path, ex))
|
||
return paths
|
||
|
||
|
||
# --- Buildings: 3D-Gebaeude DWG --------------------------------------------
|
||
|
||
# swissBUILDINGS3D 3.0: liefert mehrere Formate (DXF/DWG/OBJ/IFC/3DTiles) und
|
||
# variant-Filter (solid/separated). In Staedten sind die 3.0-Tiles aber riesig
|
||
# (>700 MB), weil nicht 1km-strukturiert — dann fallen wir auf 2.0 zurueck
|
||
# (verlaesslich 1km-Tiles, ~50 MB).
|
||
_BUILDINGS_COLLECTION_V3 = "ch.swisstopo.swissbuildings3d_3_0"
|
||
_BUILDINGS_COLLECTION_V2 = "ch.swisstopo.swissbuildings3d_2"
|
||
|
||
|
||
def _fetch_buildings_from_collection(collection_id, bbox_wgs, variant,
|
||
progress=None):
|
||
"""Holt Tile-CAD-Files aus EINER STAC-Collection. Liefert Liste Pfade
|
||
oder [] wenn nichts brauchbar geladen werden konnte (z.B. alle ueber
|
||
Size-Limit)."""
|
||
if progress: progress("STAC-Query {} (variant={})...".format(
|
||
collection_id, variant))
|
||
items = stac_query(collection_id, bbox_wgs,
|
||
asset_extensions=[".dwg", ".dxf", ".obj", ".ifc",
|
||
".dwg.zip", ".dxf.zip",
|
||
".obj.zip", ".ifc.zip"])
|
||
items = _dedupe_latest(items)
|
||
if not items:
|
||
if progress: progress(" Keine Tiles in der Region")
|
||
return []
|
||
variant_marker = "_{}_".format(variant.lower())
|
||
items_v = [it for it in items
|
||
if variant_marker in (it.get("id") or "").lower()
|
||
or any(variant_marker in (a.get("href") or "").lower()
|
||
for a in it.get("assets", {}).values())]
|
||
if items_v and len(items_v) < len(items):
|
||
if progress: progress(" Item-Filter: {}/{} matchen variant '{}'".format(
|
||
len(items_v), len(items), variant))
|
||
items = items_v
|
||
paths = []
|
||
for i, item in enumerate(items):
|
||
if progress: progress("Lade Tile {}/{}: {}".format(
|
||
i+1, len(items), item["id"]))
|
||
# Asset waehlen: nur Per-Tile-Files (Pattern `_NNNN-NN_`), bevorzugt
|
||
# direkt-DXF/DWG, sonst ZIP. Filter eliminiert auch national-GDB.
|
||
per_tile = [(k, a) for k, a in item["assets"].items()
|
||
if _TILE_FILE_PAT.search(a["href"])]
|
||
if not per_tile:
|
||
if progress: progress("→ kein Per-Tile-Asset, skip")
|
||
continue
|
||
# Varianten-Filter: `_solid_` vs `_separated_` im Filename. Default
|
||
# ist 'separated'. Falls keine Asset mit dem Marker matcht (alte
|
||
# Collection-Version o.ae.), fallen wir auf alle per-tile zurueck.
|
||
variant_marker = "_{}_".format(variant.lower())
|
||
per_tile_v = [(k, a) for k, a in per_tile
|
||
if variant_marker in a["href"].lower()]
|
||
if per_tile_v:
|
||
per_tile = per_tile_v
|
||
if progress: progress(" → {} Asset(s) matchen variant '{}'".format(
|
||
len(per_tile), variant))
|
||
else:
|
||
if progress:
|
||
hrefs_short = ", ".join(os.path.basename(a["href"])
|
||
for _, a in per_tile[:3])
|
||
progress(" → kein '{}' Marker, nehme erstes Asset (verfuegbar: {})".format(
|
||
variant, hrefs_short))
|
||
# Priorisierung: DXF/DWG (am stabilsten in Rhino) > OBJ > IFC
|
||
chosen = None
|
||
for prio_ext in [(".dxf", ".dwg"), (".obj",), (".ifc",),
|
||
(".dxf.zip", ".dwg.zip"), (".obj.zip",), (".ifc.zip",)]:
|
||
for k, a in per_tile:
|
||
low = a["href"].lower()
|
||
if low.endswith(prio_ext):
|
||
chosen = a["href"]; break
|
||
if chosen is not None: break
|
||
if chosen is None:
|
||
chosen = per_tile[0][1]["href"]
|
||
|
||
p = download_asset(chosen, subdir="buildings3d_dwg", status=progress)
|
||
if not p: continue
|
||
|
||
# ZIP-Wrapper aufloesen + Variant-Filter (ZIP kann beide DWGs enthalten)
|
||
if p.lower().endswith(".zip"):
|
||
extracted = _extract_zip_to_dir(
|
||
p, os.path.join(CACHE_DIR, "buildings3d_dwg", "_unzipped"))
|
||
cads_all = [e for e in extracted
|
||
if e.lower().endswith((".dwg", ".dxf", ".obj", ".ifc"))]
|
||
cads_v = [e for e in cads_all
|
||
if variant_marker in os.path.basename(e).lower()]
|
||
cads = cads_v if cads_v else cads_all
|
||
if cads_v and len(cads_v) < len(cads_all):
|
||
if progress: progress(" ZIP-Filter: {}/{} Files matchen '{}'".format(
|
||
len(cads_v), len(cads_all), variant))
|
||
paths.extend(cads)
|
||
else:
|
||
paths.append(p)
|
||
return paths
|
||
|
||
|
||
def fetch_buildings_dwg(bbox_lv95, progress=None, variant="separated"):
|
||
"""Holt swissBUILDINGS3D Tile-CAD-Files. Versucht erst v3.0 (separated/
|
||
solid Varianten), faellt automatisch auf v2.0 zurueck wenn v3.0 in der
|
||
Region keine brauchbaren Files liefert (typisch in Staedten — die 3.0-
|
||
Tiles sind dort >700 MB pro Stueck und werden vom Size-Limit geblockt)."""
|
||
bbox_wgs = lv95_bbox_to_wgs84_bbox(*bbox_lv95)
|
||
paths = _fetch_buildings_from_collection(
|
||
_BUILDINGS_COLLECTION_V3, bbox_wgs, variant, progress=progress)
|
||
if not paths:
|
||
if progress: progress("v3.0 lieferte keine Tiles — fallback auf v2.0 (1km-Tiles)...")
|
||
# v2.0 hat keine variant-Marker im Filename, ist immer "separated"-
|
||
# artig (Kategorien auf eigenen DXF-Layern innerhalb einer DWG).
|
||
paths = _fetch_buildings_from_collection(
|
||
_BUILDINGS_COLLECTION_V2, bbox_wgs, variant, progress=progress)
|
||
if progress: progress("{} CAD-Datei(en) bereit".format(len(paths)))
|
||
return paths
|
||
|
||
|
||
# --- Terrain: swissALTI3D via XYZ ASCII -------------------------------------
|
||
|
||
def fetch_terrain_xyz(bbox_lv95, resolution="2.0", progress=None):
|
||
"""Holt swissALTI3D-XYZ-Punktwolke. resolution: '0.5' oder '2.0' (m).
|
||
Tiles kommen als .xyz.zip — wir packen lokal aus. Liefert Liste lokaler
|
||
.xyz-Pfade (ausgepackt)."""
|
||
bbox_wgs = lv95_bbox_to_wgs84_bbox(*bbox_lv95)
|
||
if progress: progress("STAC-Query swissALTI3D...")
|
||
items = stac_query("ch.swisstopo.swissalti3d",
|
||
bbox_wgs, asset_extensions=[".xyz.zip", ".xyz"])
|
||
items = _dedupe_latest(items)
|
||
if not items:
|
||
if progress: progress("Keine ALTI-Tiles in der Region gefunden")
|
||
return []
|
||
target_tag = "_{}_".format(resolution)
|
||
paths = []
|
||
for i, item in enumerate(items):
|
||
if progress: progress("Lade Terrain {}/{}: {}".format(
|
||
i+1, len(items), item["id"]))
|
||
# Asset mit gewuenschter Resolution suchen, sonst irgendeines
|
||
chosen_href = None
|
||
for k, a in item["assets"].items():
|
||
if target_tag in a["href"]:
|
||
chosen_href = a["href"]; break
|
||
if chosen_href is None:
|
||
chosen_href = next(iter(item["assets"].values()))["href"]
|
||
p = download_asset(chosen_href, subdir="alti3d_xyz", status=progress)
|
||
if not p: continue
|
||
# ZIP-Wrapper aufloesen
|
||
if p.lower().endswith(".zip"):
|
||
extracted = _extract_zip_to_dir(
|
||
p, os.path.join(CACHE_DIR, "alti3d_xyz", "_unzipped"))
|
||
xyzs = [e for e in extracted if e.lower().endswith(".xyz")]
|
||
paths.extend(xyzs)
|
||
else:
|
||
paths.append(p)
|
||
if progress: progress("{} XYZ-Datei(en) bereit".format(len(paths)))
|
||
return paths
|
||
|
||
|
||
def xyz_to_grid(path, target_step=2.0, clip_bbox=None, progress=None):
|
||
"""Parse swissALTI3D-XYZ-ASCII. Format: jede Zeile 'E N Z' (LV95).
|
||
Liest binaer ein (XYZ-Header/Footer kann Sonderzeichen enthalten).
|
||
|
||
target_step: gewuenschter Vertex-Abstand in Metern. Sub-Sampling laeuft
|
||
RAEUMLICH (jeder N-te Punkt im E×N-Raster), nicht ueber
|
||
Zeilen-Indizes — sonst entstehen Aliasing-Streifen weil
|
||
die Zeilenzahl in der Datei nicht zwingend ein ganzzahliges
|
||
Vielfaches der Reihenbreite ist.
|
||
clip_bbox: optionales (e_min, n_min, e_max, n_max) — Punkte ausserhalb
|
||
werden verworfen. Reduziert das Mesh auf das User-Gebiet
|
||
statt der ganzen 1km×1km Tile."""
|
||
# Binaer + ASCII-decode mit errors='ignore' damit BOM/Sonderzeichen am
|
||
# File-Start nicht crashen
|
||
with open(path, "rb") as fb:
|
||
raw = fb.read()
|
||
text = raw.decode("ascii", errors="ignore")
|
||
lines = text.split("\n")
|
||
# Header-Detection
|
||
start_idx = 0
|
||
if lines:
|
||
first = lines[0].split()
|
||
if len(first) < 3:
|
||
start_idx = 1
|
||
else:
|
||
try: float(first[0])
|
||
except Exception: start_idx = 1
|
||
|
||
# --- 1. Pass: raw Step + Origin aus ersten ~200 Punkten erkennen
|
||
sample = []
|
||
for ln in lines[start_idx:start_idx + 500]:
|
||
parts = ln.split()
|
||
if len(parts) < 3: continue
|
||
try:
|
||
e = float(parts[0]); n = float(parts[1])
|
||
except Exception: continue
|
||
sample.append((e, n))
|
||
if len(sample) >= 200: break
|
||
if len(sample) < 4:
|
||
if progress: progress("XYZ leer / nicht parsebar")
|
||
return None
|
||
# Smallest non-zero E-diff = raw_e_step
|
||
e_diffs = sorted({round(abs(sample[i+1][0] - sample[i][0]), 3)
|
||
for i in range(len(sample) - 1)
|
||
if abs(sample[i+1][0] - sample[i][0]) > 0.001})
|
||
raw_e_step = e_diffs[0] if e_diffs else 0.5
|
||
n_diffs = sorted({round(abs(sample[i+1][1] - sample[i][1]), 3)
|
||
for i in range(len(sample) - 1)
|
||
if abs(sample[i+1][1] - sample[i][1]) > 0.001})
|
||
raw_n_step = n_diffs[0] if n_diffs else 0.5
|
||
origin_e = min(p[0] for p in sample)
|
||
origin_n = min(p[1] for p in sample)
|
||
# Sub-Sampling-Faktoren — nur ganzzahlig damit das Raster regulaer bleibt
|
||
factor_e = max(1, int(round(target_step / raw_e_step)))
|
||
factor_n = max(1, int(round(target_step / raw_n_step)))
|
||
actual_step_e = raw_e_step * factor_e
|
||
actual_step_n = raw_n_step * factor_n
|
||
if progress:
|
||
progress("XYZ raw {:.2f}m → target {:.2f}m → sub-sample {}x{} ({:.2f}m actual)".format(
|
||
raw_e_step, target_step, factor_e, factor_n, actual_step_e))
|
||
|
||
# --- 2. Pass: alle Punkte auf dem Sub-Raster behalten (+ optional clip)
|
||
points = {}
|
||
es = set(); ns = set()
|
||
cb = clip_bbox
|
||
for ln in lines[start_idx:]:
|
||
parts = ln.split()
|
||
if len(parts) < 3: continue
|
||
try:
|
||
e = float(parts[0]); n = float(parts[1]); z = float(parts[2])
|
||
except Exception: continue
|
||
if cb is not None:
|
||
if e < cb[0] or e > cb[2] or n < cb[1] or n > cb[3]: continue
|
||
# Raster-Pruefung: nur jeden factor_e-ten E-Schritt + factor_n-ten N-Schritt
|
||
di = int(round((e - origin_e) / raw_e_step))
|
||
dj = int(round((n - origin_n) / raw_n_step))
|
||
if di % factor_e != 0 or dj % factor_n != 0: continue
|
||
# Auf snapped Koords runden um Float-Drift zu vermeiden
|
||
e_snap = origin_e + di * raw_e_step
|
||
n_snap = origin_n + dj * raw_n_step
|
||
points[(e_snap, n_snap)] = z
|
||
es.add(e_snap); ns.add(n_snap)
|
||
|
||
if not points:
|
||
if progress: progress("Keine Punkte im Clipping-Gebiet / nach Sub-Sample")
|
||
return None
|
||
es_sorted = sorted(es); ns_sorted = sorted(ns)
|
||
if progress: progress("XYZ → {} Vertices ({}×{} Raster)".format(
|
||
len(points), len(es_sorted), len(ns_sorted)))
|
||
return {
|
||
"bbox": (es_sorted[0], ns_sorted[0], es_sorted[-1], ns_sorted[-1]),
|
||
"step": (actual_step_e, actual_step_n),
|
||
"es": es_sorted,
|
||
"ns": ns_sorted,
|
||
"points": points,
|
||
}
|
||
|
||
|
||
def merge_grids(grids):
|
||
"""Vereint mehrere Grid-Dicts (eines pro XYZ-Tile) zu einem zusammen-
|
||
haengenden Grid. swissALTI3D-Tiles liefern jeweils 1km×1km Punkte —
|
||
benachbarte Tiles teilen KEINE Rand-Punkte (Tile A endet z.B. bei
|
||
e=2700999.5, Tile B startet bei e=2701000.0). Beim getrennten Meshen
|
||
entsteht dadurch ein 1m-Streifen ohne Faces. Hier mergen wir die
|
||
Punkte VORHER zu einem unified Grid, dann verbindet mesh_from_grid
|
||
die Tile-Grenze automatisch (benachbarte Spalten = ein step Abstand
|
||
nach Sub-Sampling).
|
||
|
||
Annahme: alle Grids teilen sich denselben step (gleicher target_step
|
||
+ Quell-Resolution). Origin-Alignment ist gegeben, weil swissALTI3D
|
||
auf einem globalen 0.5m-Raster liegt."""
|
||
if not grids: return None
|
||
grids = [g for g in grids if g is not None]
|
||
if not grids: return None
|
||
if len(grids) == 1: return grids[0]
|
||
step = grids[0]["step"]
|
||
all_points = {}
|
||
all_es = set(); all_ns = set()
|
||
for g in grids:
|
||
for (e, n), z in g["points"].items():
|
||
all_points[(e, n)] = z
|
||
all_es.add(e); all_ns.add(n)
|
||
if not all_points: return None
|
||
es_sorted = sorted(all_es); ns_sorted = sorted(all_ns)
|
||
return {
|
||
"bbox": (es_sorted[0], ns_sorted[0], es_sorted[-1], ns_sorted[-1]),
|
||
"step": step,
|
||
"es": es_sorted,
|
||
"ns": ns_sorted,
|
||
"points": all_points,
|
||
}
|
||
|
||
|
||
def mesh_from_grid(grid, origin_shift=(0, 0, 0), unit_scale=1.0):
|
||
"""Baut ein Rhino-Mesh aus dem XYZ-Grid. origin_shift wird auf jeden
|
||
Vertex angewendet (typisch: bbox-Center zu Welt-0/0/0 schieben).
|
||
unit_scale: Skalierung von Meter (Quelle XYZ) auf Doc-Units. Bei
|
||
mm-Doc = 1000, bei m-Doc = 1.0 ."""
|
||
es = grid["es"]; ns = grid["ns"]
|
||
pts = grid["points"]
|
||
sx, sy, sz = origin_shift
|
||
mesh = rg.Mesh()
|
||
idx_for = {}
|
||
for j, ny in enumerate(ns):
|
||
for i, ex in enumerate(es):
|
||
z = pts.get((ex, ny))
|
||
if z is None: continue
|
||
v_idx = mesh.Vertices.Add(
|
||
(ex - sx) * unit_scale,
|
||
(ny - sy) * unit_scale,
|
||
(z - sz) * unit_scale)
|
||
idx_for[(i, j)] = v_idx
|
||
# Quads
|
||
for j in range(len(ns) - 1):
|
||
for i in range(len(es) - 1):
|
||
a = idx_for.get((i, j))
|
||
b = idx_for.get((i+1, j))
|
||
c = idx_for.get((i+1, j+1))
|
||
d = idx_for.get((i, j+1))
|
||
if a is None or b is None or c is None or d is None: continue
|
||
mesh.Faces.AddFace(a, b, c, d)
|
||
mesh.Normals.ComputeNormals()
|
||
mesh.Compact()
|
||
return mesh
|
||
|
||
|
||
# --- Orthofoto: SWISSIMAGE 10cm via GeoTIFF --------------------------------
|
||
|
||
def fetch_orthophoto(bbox_lv95, resolution="2.0", progress=None):
|
||
"""Holt SWISSIMAGE-10cm-Tiles fuer LV95-bbox. resolution: '0.1' (10 cm,
|
||
sehr gross!), '0.5' (50 cm), '2.0' (2 m, Default — Material-Texture
|
||
braucht keine 10cm)."""
|
||
bbox_wgs = lv95_bbox_to_wgs84_bbox(*bbox_lv95)
|
||
if progress: progress("STAC-Query SWISSIMAGE...")
|
||
items = stac_query("ch.swisstopo.swissimage-dop10",
|
||
bbox_wgs, asset_extensions=[".tif"])
|
||
items = _dedupe_latest(items)
|
||
if not items:
|
||
if progress: progress("Kein Luftbild in der Region gefunden")
|
||
return []
|
||
target_tag = "_{}_".format(resolution)
|
||
paths = []
|
||
for i, item in enumerate(items):
|
||
if progress: progress("Lade Luftbild {}/{}: {}".format(
|
||
i+1, len(items), item["id"]))
|
||
chosen_href = None
|
||
for k, a in item["assets"].items():
|
||
if target_tag in a["href"]:
|
||
chosen_href = a["href"]; break
|
||
if chosen_href is None:
|
||
chosen_href = next(iter(item["assets"].values()))["href"]
|
||
p = download_asset(chosen_href, subdir="swissimage_tif", status=progress)
|
||
if p: paths.append(p)
|
||
if progress: progress("{} Luftbild-Tile(s) bereit".format(len(paths)))
|
||
return paths
|
||
|
||
|
||
def _geotiff_to_png(tif_path, max_dim=2048):
|
||
"""SWISSIMAGE kommt als GeoTIFF — Rhinos Material-Bitmap kann GeoTIFF nicht
|
||
direkt lesen. Konvertiere zu PNG. Zwei Wege:
|
||
1) Pillow (wenn in Rhinos CPython verfuegbar) — universell + downsample
|
||
2) Eto.Drawing.Bitmap (Mac: NSImage liest TIFF nativ) — Fallback
|
||
Liefert PNG-Pfad oder None bei Fehler."""
|
||
if not tif_path: return None
|
||
base, _ = os.path.splitext(tif_path)
|
||
png_path = base + "_2k.png"
|
||
if os.path.isfile(png_path) and os.path.getsize(png_path) > 0:
|
||
print("[SWISSTOPO] PNG-Cache:", os.path.basename(png_path))
|
||
return png_path
|
||
# --- Variante 1: Pillow
|
||
try:
|
||
from PIL import Image
|
||
img = Image.open(tif_path)
|
||
if max(img.width, img.height) > max_dim:
|
||
scale = max_dim / float(max(img.width, img.height))
|
||
new_w = max(1, int(img.width * scale))
|
||
new_h = max(1, int(img.height * scale))
|
||
img = img.resize((new_w, new_h), Image.LANCZOS)
|
||
if img.mode not in ("RGB", "RGBA"):
|
||
img = img.convert("RGB")
|
||
img.save(png_path, "PNG", optimize=False)
|
||
print("[SWISSTOPO] Pillow: {} → {} ({}x{}px)".format(
|
||
os.path.basename(tif_path), os.path.basename(png_path),
|
||
img.width, img.height))
|
||
return png_path
|
||
except ImportError:
|
||
print("[SWISSTOPO] Pillow nicht verfuegbar — versuche Eto.Drawing")
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] Pillow-convert fail:", ex)
|
||
# --- Variante 2: Eto.Drawing (Mac NSImage liest TIFF)
|
||
try:
|
||
import Eto.Drawing as _ed
|
||
bmp_src = _ed.Bitmap(tif_path)
|
||
if bmp_src is None:
|
||
print("[SWISSTOPO] Eto konnte TIFF nicht laden")
|
||
return None
|
||
# Downsample falls > max_dim
|
||
w, h = bmp_src.Width, bmp_src.Height
|
||
if max(w, h) > max_dim:
|
||
scale = max_dim / float(max(w, h))
|
||
new_w = max(1, int(w * scale))
|
||
new_h = max(1, int(h * scale))
|
||
target = _ed.Bitmap(new_w, new_h, _ed.PixelFormat.Format32bppRgba)
|
||
g = _ed.Graphics(target)
|
||
try:
|
||
try: g.AntiAlias = True
|
||
except Exception: pass
|
||
g.DrawImage(bmp_src, 0, 0, new_w, new_h)
|
||
finally: g.Dispose()
|
||
bmp_src = target
|
||
w, h = new_w, new_h
|
||
try: bmp_src.Save(png_path, _ed.ImageFormat.Png)
|
||
except Exception:
|
||
# Eto.ImageFormat-Variante kann je nach Eto-Version variieren
|
||
bmp_src.Save(png_path)
|
||
print("[SWISSTOPO] Eto: {} → {} ({}x{}px)".format(
|
||
os.path.basename(tif_path), os.path.basename(png_path), w, h))
|
||
return png_path
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] Eto-convert fail:", ex)
|
||
return None
|
||
|
||
|
||
def add_ortho_plane(doc, ortho_path, tile_bbox_lv95, shift_lv95, m_to_unit,
|
||
z_doc=0.0):
|
||
"""Erzeugt eine planare Brep-Flaeche mit dem SWISSIMAGE-Foto als Material,
|
||
direkt sichtbar in Top/Shaded/Rendered Display-Mode.
|
||
|
||
tile_bbox_lv95: (e_min, n_min, e_max, n_max) in LV95-Metern der Tile-Region
|
||
shift_lv95: (sx, sy, sz) — Origin-Shift in LV95-Metern (typisch eC,nC)
|
||
m_to_unit: Skalierung m → doc-units (z.B. 0.001 fuer km-Doc)
|
||
z_doc: Z-Hoehe der Plane in Doc-Units (typisch max-Terrain-Z + Epsilon)
|
||
|
||
Liefert den RhinoObject der erzeugten Plane (oder None)."""
|
||
if not (ortho_path and os.path.isfile(ortho_path)): return None
|
||
# GeoTIFF → PNG damit Rhino's Material-Bitmap es als Diffuse nehmen kann
|
||
if ortho_path.lower().endswith((".tif", ".tiff")):
|
||
png = _geotiff_to_png(ortho_path)
|
||
if not png: return None
|
||
ortho_path = png
|
||
# bbox in Doc-Units (nach Shift + Scale)
|
||
e_min, n_min, e_max, n_max = tile_bbox_lv95
|
||
sx, sy, sz = shift_lv95
|
||
x_min = (e_min - sx) * m_to_unit
|
||
x_max = (e_max - sx) * m_to_unit
|
||
y_min = (n_min - sy) * m_to_unit
|
||
y_max = (n_max - sy) * m_to_unit
|
||
# Mesh-Quad mit expliziten Per-Vertex-UV-Koordinaten — bombensicher
|
||
# fuer Cycles/Raytraced. Eine Brep-Plane braucht erst Render-Mesh-
|
||
# Erzeugung + TextureMapping, was diverse Fallstricke hat.
|
||
mesh = rg.Mesh()
|
||
mesh.Vertices.Add(x_min, y_min, z_doc) # 0 → UV (0,0)
|
||
mesh.Vertices.Add(x_max, y_min, z_doc) # 1 → UV (1,0)
|
||
mesh.Vertices.Add(x_max, y_max, z_doc) # 2 → UV (1,1)
|
||
mesh.Vertices.Add(x_min, y_max, z_doc) # 3 → UV (0,1)
|
||
mesh.Faces.AddFace(0, 1, 2, 3)
|
||
mesh.TextureCoordinates.Add(0.0, 0.0)
|
||
mesh.TextureCoordinates.Add(1.0, 0.0)
|
||
mesh.TextureCoordinates.Add(1.0, 1.0)
|
||
mesh.TextureCoordinates.Add(0.0, 1.0)
|
||
mesh.Normals.ComputeNormals()
|
||
mesh.Compact()
|
||
gid = doc.Objects.AddMesh(mesh)
|
||
obj = doc.Objects.Find(gid)
|
||
if obj is None: return None
|
||
# Material: Legacy + ToPhysicallyBased + PBR_BaseColor-Texture.
|
||
# Bekannt instabil unter Mac Rhino 8 für Raytraced (Cycles greift den
|
||
# Shim nicht zuverlaessig); zumindest Shaded zeigt die Textur.
|
||
try:
|
||
mat = Rhino.DocObjects.Material()
|
||
mat.Name = "swisstopo_ortho"
|
||
mat.SetBitmapTexture(ortho_path)
|
||
mat.ToPhysicallyBased()
|
||
tex = Rhino.DocObjects.Texture()
|
||
tex.FileName = ortho_path
|
||
tex.Enabled = True
|
||
mat.SetTexture(tex, Rhino.DocObjects.TextureType.PBR_BaseColor)
|
||
midx = doc.Materials.Add(mat)
|
||
attrs = obj.Attributes.Duplicate()
|
||
attrs.MaterialSource = Rhino.DocObjects.ObjectMaterialSource.MaterialFromObject
|
||
attrs.MaterialIndex = midx
|
||
doc.Objects.ModifyAttributes(obj, attrs, True)
|
||
except Exception as ex:
|
||
print("[SWISSTOPO] ortho-material:", ex)
|
||
return obj
|