#! 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 math
import zipfile
import urllib.request
import urllib.parse
import Rhino
import Rhino.Geometry as rg
import System

DEFAULT_CACHE_DIR = os.path.expanduser("~/Library/Caches/Dossier/swisstopo")
CACHE_DIR = DEFAULT_CACHE_DIR
STAC_BASE = "https://data.geo.admin.ch/api/stac/v1"
SEARCH_API = "https://api3.geo.admin.ch/rest/services/api/SearchServer"


def get_cache_dir_for_doc(doc):
    """Cache-Pfad fuer ein Doc. Wenn das Doc auf Disk liegt: Subfolder neben
    der .3dm-Datei (`<dir>/<basename>_swisstopo/`). Damit reisen die Files
    mit dem Projekt — kann via SMB von anderen Maschinen geoeffnet werden
    solange der Mount-Pfad identisch ist. Falls Doc nicht gespeichert:
    globaler Fallback-Cache."""
    try:
        p = doc.Path if doc else None
        if p and os.path.isfile(p):
            doc_dir  = os.path.dirname(p)
            doc_base = os.path.splitext(os.path.basename(p))[0]
            return os.path.join(doc_dir, doc_base + "_swisstopo")
    except Exception: pass
    return DEFAULT_CACHE_DIR


def set_cache_dir(path):
    """Stellt das aktive Cache-Verzeichnis. Alle nachfolgenden Downloads
    landen dort. Aufrufer-Verantwortung: vor jedem Import den richtigen
    Cache setzen (per-Doc oder global)."""
    global CACHE_DIR
    CACHE_DIR = path
    _ensure_cache()


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",
                          version="v2"):
    """Holt swissBUILDINGS3D Tile-CAD-Files.

    version='v2': stabile 2.0-Variante (1km-Tiles, keine Solid/Separated-
                  Aufteilung — alle Kategorien auf eigenen DXF-Layern).
    version='v3': Beta 3.0-Variante mit Solid/Separated-Wahl. In Staedten
                  oft >700 MB pro Tile → auto-fallback auf v2 wenn v3
                  nichts brauchbares liefert."""
    bbox_wgs = lv95_bbox_to_wgs84_bbox(*bbox_lv95)
    if version == "v3":
        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...")
            paths = _fetch_buildings_from_collection(
                _BUILDINGS_COLLECTION_V2, bbox_wgs, variant, progress=progress)
    else:
        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


# --- TLM3D Vektor (Strassen / Gewaesser / Bahn / Vegetation) ----------------
# swisstopo bietet TLM3D unter mehreren Collection-IDs an (genaue Namen
# variieren). Wir probieren defensiv mehrere Kandidaten und nehmen DXF/DWG
# wenn verfuegbar (alles andere — GPKG/SHP — koennen wir nicht parsen).

# Echte swisstopo TLM-Collections (verifiziert via STAC API):
#   ch.swisstopo.swisstlm3d         — voller TLM3D Layer (~ganze CH)
#   ch.swisstopo.swisstlmregio      — kleinere Auflösung 1:200000
#   ch.swisstopo.swissboundaries3d  — Verwaltungsgrenzen
#   ch.swisstopo.swiss-map-vector25 — 1:25000 Vektor
# Achtung: ALLE liefern nur GDB/SHP/GPKG/XTF — KEIN DXF/DWG. Direkter Rhino-
# Import funktioniert nicht ohne Shapefile-/GPKG-Parser.
_TLM_COLLECTIONS = [
    "ch.swisstopo.swisstlm3d",
    "ch.swisstopo.swisstlmregio",
    "ch.swisstopo.swissboundaries3d",
]


def fetch_tlm3d_vector(bbox_lv95, kinds, progress=None):
    """Versucht swissTLM3D-Daten als DXF/DWG zu holen. swisstopo liefert
    aktuell NUR GDB/SHP/GPKG-Formate — kein DXF. Diese Funktion findet
    daher in den meisten Faellen keine importierbaren Files; sie loggt
    aber sauber, was verfuegbar waere, falls wir spaeter einen
    Shapefile-Parser einbauen."""
    bbox_wgs = lv95_bbox_to_wgs84_bbox(*bbox_lv95)
    out = {}
    if progress:
        progress("TLM3D-Import: swisstopo bietet aktuell KEINE DXF-Assets")
        progress("  (nur GDB/SHP/GPKG — Rhino kann diese nicht nativ lesen)")
        progress("  Verfuegbare Collections (zur Info):")
    for coll in _TLM_COLLECTIONS:
        try:
            items = stac_query(coll, bbox_wgs,
                                asset_extensions=None)  # alle Assets
        except Exception as ex:
            if progress: progress("    {}: HTTP-fail ({})".format(coll, ex))
            continue
        if not items:
            if progress: progress("    {}: keine Items in der Region".format(coll))
            continue
        sample = items[0]
        formats = set()
        for k, a in (sample.get("assets") or {}).items():
            href = (a.get("href") or "").lower()
            for ext in (".gdb.zip", ".shp.zip", ".gpkg.zip", ".gpkg",
                         ".xtf.zip", ".dxf", ".dwg"):
                if href.endswith(ext): formats.add(ext.lstrip("."))
        if progress: progress("    {}: {} Items, Formate: {}".format(
            coll, len(items), ", ".join(sorted(formats)) or "?"))
    if progress:
        progress("→ TLM3D-Direct-Import nicht moeglich. Nutze OSM-Importer "
                 "fuer Vector-Daten (Strassen/Wasser/Gebaeude).")
    return out


# --- 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 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
    # 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
    # Phase relativ zum globalen LV95-Raster: swissALTI3D 0.5m liegt z.B.
    # auf .25-Marken (cell-center). Phase aus erstem Sample-Punkt ermitteln —
    # gilt global fuer alle Tiles, da das Raster LV95-aligned ist.
    e_phase = sample[0][0] - math.floor(sample[0][0] / raw_e_step) * raw_e_step
    n_phase = sample[0][1] - math.floor(sample[0][1] / raw_n_step) * raw_n_step
    if progress:
        progress("XYZ raw {:.2f}m → target {:.2f}m → sub-sample {}x{} ({:.2f}m actual, phase {:.2f}/{:.2f})".format(
            raw_e_step, target_step, factor_e, factor_n, actual_step_e, e_phase, n_phase))

    # --- 2. Pass: alle Punkte auf dem Sub-Raster behalten (+ optional clip)
    # WICHTIG: Sub-Sampling-Filter benutzt den GLOBALEN LV95-Raster-Index
    # (mit detected phase) statt eines tile-lokalen origin. Sonst waehlt
    # jedes Tile seine eigene Phase und am Tile-Boundary fehlen Faces /
    # das Mesh ist nicht durchgehend.
    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
        # Globaler Index im swissALTI3D-Raster: alle Tiles teilen Phase
        gi = int(round((e - e_phase) / raw_e_step))
        gj = int(round((n - n_phase) / raw_n_step))
        if gi % factor_e != 0 or gj % factor_n != 0: continue
        # Originale Koordinaten als Key (Tile A und Tile B teilen Phase,
        # also matchen ihre Keys am Boundary direkt). round(3) gegen
        # Float-Drift.
        e_snap = round(e, 3)
        n_snap = round(n, 3)
        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


def generate_mesh_from_contours(doc, contour_curves, sample_step_m=2.0,
                                  m_to_unit=1.0, progress=None):
    """Baut ein TIN-Mesh aus Hoehenlinien-Curves. Jede Curve hat ihre echte
    Z-Hoehe — wir sampeln Vertices entlang der Curves und triangulieren
    sie via Rhinos _-MeshPatch / _-Delaunay Command. Resultat: Topographie-
    Mesh basierend auf den diskreten Hoehenlinien-Stufen.

    Liefert RhinoObject (Mesh) oder None."""
    import System
    if not contour_curves: return None
    pts = []
    for c in contour_curves:
        if c is None: continue
        # Polyline-Vertices wenn moeglich (exakt), sonst entlang Curve sampeln
        ok, poly = c.TryGetPolyline()
        if ok and poly is not None:
            for pt in poly: pts.append(rg.Point3d(pt))
        else:
            try:
                L = c.GetLength()
                n = max(2, int(L / (sample_step_m * m_to_unit)))
                params = c.DivideByCount(n, True)
                if params:
                    for t in params: pts.append(c.PointAt(t))
            except Exception: pass
    if len(pts) < 3:
        if progress: progress("Contour-Mesh: zu wenig Vertices ({})".format(len(pts)))
        return None
    if progress: progress("Contour-Mesh: trianguliere {} Vertices...".format(len(pts)))
    # Temp-Points erzeugen + selektieren
    temp_pids = []
    try:
        for p in pts:
            pid = doc.Objects.AddPoint(p)
            if pid and pid != System.Guid.Empty:
                temp_pids.append(pid)
        if not temp_pids:
            if progress: progress("Contour-Mesh: keine Temp-Points")
            return None
        doc.Objects.UnselectAll()
        for pid in temp_pids: doc.Objects.Select(pid)
        before = set(o.Id for o in doc.Objects
                      if o and not o.IsDeleted
                      and isinstance(o.Geometry, rg.Mesh))
        # Mehrere Commands probieren (Mac Rhino 8 vs neuere Versionen)
        cmd_tried = None
        for cmd in ['_-MeshPatch _Enter _Enter',
                     '_-Delaunay _Enter',
                     '_-DelaunayMesh _Enter',
                     '_-MeshFromPoints _Enter']:
            try:
                Rhino.RhinoApp.RunScript(cmd, False)
            except Exception: continue
            cmd_tried = cmd
            new_mesh = next((o for o in doc.Objects
                              if o and not o.IsDeleted
                              and isinstance(o.Geometry, rg.Mesh)
                              and o.Id not in before), None)
            if new_mesh:
                if progress: progress("→ Contour-Mesh via '{}'".format(cmd.split()[0]))
                return new_mesh
        if progress:
            progress("Contour-Mesh: kein Command lieferte ein Mesh "
                     "(zuletzt: {})".format(cmd_tried))
        return None
    finally:
        # Temp-Points wieder weg
        doc.Objects.UnselectAll()
        for pid in temp_pids:
            try: doc.Objects.Delete(pid, True)
            except Exception: pass


def generate_schichtenmodell(doc, contour_curves, progress=None):
    """Schichtenmodell: jede geschlossene Hoehenlinie wird zu einer planaren
    Flaeche auf ihrer Z-Hoehe. Stacked Discs — der architektonische
    'Pappmodell'-Look. Offene Konturen (typ. am bbox-Rand) werden
    uebersprungen.

    Liefert Liste von erzeugten RhinoObjects."""
    import System
    if not contour_curves: return []
    created = []
    tol = doc.ModelAbsoluteTolerance
    n_open = 0
    for c in contour_curves:
        if c is None: continue
        try:
            if not c.IsClosed:
                n_open += 1
                continue
            breps = rg.Brep.CreatePlanarBreps(c, tol)
        except Exception:
            continue
        if not breps: continue
        for brep in breps:
            gid = doc.Objects.AddBrep(brep)
            if gid and gid != System.Guid.Empty:
                obj = doc.Objects.Find(gid)
                if obj: created.append(obj)
    if progress:
        progress("→ {} Schichten-Flaechen ({} offene Konturen skipped)".format(
            len(created), n_open))
    return created


def generate_patch_from_contours(doc, contour_curves, progress=None):
    """Erzeugt eine NURBS-Patch-Flaeche durch alle Hoehenlinien — der
    kanonische Rhino-Workflow fuer Terrain aus Konturen (Brep.CreatePatch
    fittet eine Surface durch Curves + Points). Loft funktioniert hier
    nicht, weil benachbarte Konturen unterschiedliche Topologie haben
    (Inseln, Halbinseln, geschlossen vs. offen am Rand).

    UV-Spans werden aus dem Bounding-Box-Aspect-Ratio abgeleitet, so
    dass ein rechteckiger Site mehr Spans in der laengeren Richtung
    bekommt.

    Liefert das erstellte RhinoObject oder None."""
    import System
    if not contour_curves: return None
    valid = [c for c in contour_curves if c is not None]
    if len(valid) < 2:
        if progress: progress("Patch braucht mindestens 2 Hoehenlinien")
        return None
    bb = rg.BoundingBox.Unset
    for c in valid:
        bb_c = c.GetBoundingBox(True)
        if not bb_c.IsValid: continue
        if not bb.IsValid:
            bb = bb_c
        else:
            bb.Union(bb_c)
    if not bb.IsValid:
        if progress: progress("Patch: ungueltige Bounding-Box")
        return None
    dx = bb.Max.X - bb.Min.X
    dy = bb.Max.Y - bb.Min.Y
    base_span = 40
    if dx >= dy and dy > 1e-6:
        aspect = dx / dy
        u_spans = int(round(base_span * math.sqrt(aspect)))
        v_spans = base_span
    elif dx > 1e-6:
        aspect = dy / dx
        u_spans = base_span
        v_spans = int(round(base_span * math.sqrt(aspect)))
    else:
        u_spans = v_spans = base_span
    u_spans = max(8, min(200, u_spans))
    v_spans = max(8, min(200, v_spans))
    if progress:
        progress("Patch aus {} Hoehenlinien ({}x{} UV-Spans)...".format(
            len(valid), u_spans, v_spans))
    geom_list = System.Collections.Generic.List[rg.GeometryBase]()
    for c in valid: geom_list.Add(c)
    try:
        brep = rg.Brep.CreatePatch(
            geom_list, u_spans, v_spans, doc.ModelAbsoluteTolerance)
        if brep is None:
            if progress: progress("Patch fehlgeschlagen (None zurueck)")
            return None
        gid = doc.Objects.AddBrep(brep)
        if gid and gid != System.Guid.Empty:
            obj = doc.Objects.Find(gid)
            if progress: progress("→ Patch-Terrain erzeugt")
            return obj
    except Exception as ex:
        if progress: progress("Patch-Fehler: {}".format(ex))
    return None


def volumize_terrain_object(doc, top_obj, depth_doc, progress=None):
    """Wandelt ein offenes Terrain (Mesh ODER Brep) in ein geschlossenes
    Mesh-Volumen um: Skirt um den Boundary + planarer Boden bei
    (min_z - depth_doc). Resultat hat eine Section beim Schneiden mit
    einer Clipping-Plane.

    Strategie:
      1. Mesh-Source ermitteln (Brep → Mesh.CreateFromBrep, Mesh → direkt)
      2. GetNakedEdges() liefert die Boundary-Loop(s) als Polylines
      3. Pro Loop: Skirt-Quads zwischen Top-Edge und Bottom-Vertices
      4. Pro Loop: Bottom-Cap via Mesh.CreateFromClosedPolyline (Rhino
         triangliert auch nicht-konvexe Boundaries sauber)
      5. CombineIdentical schweisst Top + Skirt-Top zusammen

    Ersetzt das Original im Doc (Delete+Add mit gleichen Attributes).
    Liefert das neue RhinoObject oder None bei Fehler."""
    import System
    if top_obj is None or top_obj.IsDeleted: return None
    geom = top_obj.Geometry
    if geom is None: return None
    # 1) Top-Mesh ermitteln (Brep meshen wenn noetig)
    top_mesh = None
    if isinstance(geom, rg.Mesh):
        top_mesh = geom.Duplicate()
    elif isinstance(geom, rg.Brep):
        try:
            mp = rg.MeshingParameters.Default
            meshes = rg.Mesh.CreateFromBrep(geom, mp)
            if meshes and len(meshes) > 0:
                joined = rg.Mesh()
                for m in meshes: joined.Append(m)
                top_mesh = joined
        except Exception as ex:
            if progress: progress("Volumize: Brep-Meshing-Fehler: {}".format(ex))
            return None
    elif isinstance(geom, rg.Extrusion):
        try:
            brep = geom.ToBrep(False)
            mp = rg.MeshingParameters.Default
            meshes = rg.Mesh.CreateFromBrep(brep, mp)
            if meshes and len(meshes) > 0:
                joined = rg.Mesh()
                for m in meshes: joined.Append(m)
                top_mesh = joined
        except Exception: pass
    if top_mesh is None or top_mesh.Vertices.Count < 3:
        if progress: progress("Volumize: kein Mesh-Top")
        return None
    # 2) Boundary-Loops
    naked = top_mesh.GetNakedEdges()
    if naked is None or len(naked) == 0:
        if progress: progress("Volumize: keine Boundary — Terrain schon geschlossen")
        return None
    # 3) Bottom-Z = min_z des Top - depth
    bb = top_mesh.GetBoundingBox(True)
    if not bb.IsValid:
        if progress: progress("Volumize: ungueltige BoundingBox")
        return None
    bottom_z = bb.Min.Z - float(depth_doc)
    if progress:
        progress("Volumize: {} Boundary-Loop(s), Boden bei Z={:.3f}".format(
            len(naked), bottom_z))
    # 4) Volumen-Mesh aufbauen: top + Skirt + Bottom-Cap
    vol = top_mesh.Duplicate()
    for loop in naked:
        try:
            if loop is None or loop.Count < 3: continue
            # Polyline-Punkte (offene Form — closing point ggf. entfernen)
            pts = [rg.Point3d(p) for p in loop]
            if len(pts) > 1 and pts[0].DistanceTo(pts[-1]) < 1e-6:
                pts = pts[:-1]
            n = len(pts)
            if n < 3: continue
            # Top + Bottom Vertices anfuegen
            top_idx = []
            bot_idx = []
            for p in pts:
                top_idx.append(vol.Vertices.Add(p.X, p.Y, p.Z))
                bot_idx.append(vol.Vertices.Add(p.X, p.Y, bottom_z))
            # Skirt: Quads zwischen aufeinanderfolgenden Top/Bottom-Paaren.
            # Faces sind "innen orientiert" — bei Bedarf normals
            # umdrehen via ComputeNormals + RebuildNormals.
            for i in range(n):
                j = (i + 1) % n
                vol.Faces.AddFace(top_idx[i], top_idx[j],
                                   bot_idx[j], bot_idx[i])
            # Bottom-Cap via planar Polyline → Mesh
            bot_pts = [rg.Point3d(p.X, p.Y, bottom_z) for p in pts]
            bot_pts.append(bot_pts[0])  # schliessen
            bot_poly = rg.Polyline(bot_pts)
            cap = rg.Mesh.CreateFromClosedPolyline(bot_poly)
            if cap is not None and cap.Vertices.Count >= 3:
                vol.Append(cap)
        except Exception as ex:
            if progress: progress("Volumize: Loop-Fehler: {}".format(ex))
    # 5) Cleanup
    try: vol.Vertices.CombineIdentical(True, True)
    except Exception: pass
    try: vol.Compact()
    except Exception: pass
    try:
        vol.Normals.ComputeNormals()
        vol.FaceNormals.ComputeFaceNormals()
        # Topologie pruefen + Naked-Edges-Anzahl loggen
        post_naked = vol.GetNakedEdges()
        if progress:
            n_naked = len(post_naked) if post_naked else 0
            progress("Volumize: Resultat {} naked-edge-loops (0 = closed)".format(
                n_naked))
    except Exception: pass
    # 6) Original ersetzen — Attributes + LayerIndex behalten
    try:
        attrs = top_obj.Attributes.Duplicate()
        old_id = top_obj.Id
        new_gid = doc.Objects.AddMesh(vol, attrs)
        if new_gid is None or new_gid == System.Guid.Empty:
            if progress: progress("Volumize: AddMesh fehlgeschlagen")
            return None
        doc.Objects.Delete(old_id, True)
        new_obj = doc.Objects.Find(new_gid)
        if progress: progress("→ Terrain-Volumen erzeugt")
        return new_obj
    except Exception as ex:
        if progress: progress("Volumize: Replace-Fehler: {}".format(ex))
        return None


def generate_contour_curves(grid, shift_lv95, m_to_unit, interval=2.0,
                              progress=None):
    """Generiert Hoehenlinien (Contour-Curves) aus dem Terrain-Grid via
    Mesh.CreateContourCurves.

    interval: Hoehenabstand in REALEN METERN (1.0/2.0/5.0 typisch).
    Liefert Liste von rg.Curve-Objekten in Doc-Units. Caller macht
    doc.Objects.AddCurve + Layer-Move."""
    if not grid or not grid.get("points"): return []
    # Temp-Mesh aus Grid (gleicher Pipeline wie mesh_from_grid)
    mesh = mesh_from_grid(grid, origin_shift=shift_lv95, unit_scale=m_to_unit)
    if mesh.Vertices.Count < 3: return []
    bb = mesh.GetBoundingBox(True)
    z_min_doc = bb.Min.Z
    z_max_doc = bb.Max.Z
    interval_doc = interval * m_to_unit
    if interval_doc <= 0: return []
    if progress:
        z_min_m = z_min_doc / m_to_unit + shift_lv95[2]
        z_max_m = z_max_doc / m_to_unit + shift_lv95[2]
        progress("Hoehenlinien: Z {:.1f}–{:.1f} m.ü.M, Abstand {} m".format(
            z_min_m, z_max_m, interval))
    try:
        curves = rg.Mesh.CreateContourCurves(
            mesh,
            rg.Point3d(0, 0, z_min_doc),
            rg.Point3d(0, 0, z_max_doc),
            interval_doc)
    except Exception as ex:
        if progress: progress("Contour fail: {}".format(ex))
        return []
    if not curves:
        if progress: progress("Keine Hoehenlinien erzeugt")
        return []
    out = list(curves)
    if progress: progress("→ {} Hoehenlinien-Kurven".format(len(out)))
    return out


# --- 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_draped_mesh(doc, ortho_path, tile_bbox_lv95, terrain_grid,
                            shift_lv95, m_to_unit, z_lift=0.05,
                            target_layer_idx=-1):
    """Erzeugt ein Mesh, das der Topographie folgt — textured mit dem Ortho-
    Foto. Statt einer flachen Plane: Per-Tile-Sub-Mesh aus dem Terrain-Grid
    mit Per-Vertex-UV (0..1 ueber die Tile-Breite). Material kommt von einem
    temporaeren PictureFrame (das ist der einzige Weg auf Mac Rhino 8 die
    embedded Bitmap in Cycles zur Anzeige zu bringen) — der PictureFrame
    wird hinterher geloescht, nur das Drape-Mesh bleibt.

    terrain_grid: dict aus merge_grids() — wir extrahieren daraus die Punkte
                  innerhalb der Tile-bbox.
    z_lift:       kleiner Z-Offset (in doc-units) gegen Z-Fighting mit dem
                  darunterliegenden Terrain-Mesh."""
    if not (ortho_path and os.path.isfile(ortho_path)): return None
    # TIF direkt verwenden — Rhino's _Picture liest GeoTIFF nativ ueber
    # NSImage (Mac) und behaelt 10cm-Aufloesung (10000×10000 px statt 2k PNG).
    e_min, n_min, e_max, n_max = tile_bbox_lv95
    sx, sy, sz = shift_lv95
    # Terrain-Punkte innerhalb des Tiles aus dem Merged-Grid extrahieren
    es = sorted(e for e in terrain_grid["es"]
                if e_min - 0.01 <= e <= e_max + 0.01)
    ns = sorted(n for n in terrain_grid["ns"]
                if n_min - 0.01 <= n <= n_max + 0.01)
    if len(es) < 2 or len(ns) < 2:
        print("[SWISSTOPO] drape: zu wenig Terrain-Punkte fuer Tile")
        return None
    pts = terrain_grid["points"]
    span_e = e_max - e_min
    span_n = n_max - n_min
    # Half-Pixel-Inset: bei 10000×10000 px Tiles wuerde ein Sample exakt an
    # u=0 oder u=1 auf der Pixel-Grenze landen; mit clamp-to-border kann das
    # weisse Linien an den Tile-Boundaries erzeugen. Wir verschieben UV
    # minimal nach innen.
    UV_INSET = 0.5 / 10000.0   # halbe Pixel-Breite im UV-Raum
    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) * m_to_unit,
                (ny - sy) * m_to_unit,
                (z - sz) * m_to_unit + z_lift)
            u = UV_INSET + (ex - e_min) / span_e * (1.0 - 2 * UV_INSET)
            v = UV_INSET + (ny - n_min) / span_n * (1.0 - 2 * UV_INSET)
            mesh.TextureCoordinates.Add(u, v)
            idx_for[(i, j)] = v_idx
    n_faces = 0
    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)
            n_faces += 1
    if n_faces == 0:
        print("[SWISSTOPO] drape: keine Faces erzeugt")
        return None
    mesh.Normals.ComputeNormals()
    mesh.Compact()
    # Temp-PictureFrame off-screen erzeugen — ergibt working RenderMaterial
    # mit Bitmap-Texture, das wir auf das Mesh uebertragen.
    # embedBitmap=False: Pfad-Referenz statt 70MB-TIF-Embedding ins .3dm.
    # Cache ist persistent (~/Library/Caches), Pfad bleibt gueltig.
    pf_plane = rg.Plane(rg.Point3d(-1e6, -1e6, -1e6),
                         rg.Vector3d.XAxis, rg.Vector3d.YAxis)
    try:
        pf_gid = doc.Objects.AddPictureFrame(
            pf_plane, ortho_path, False, 1.0, 1.0, True, False)
    except Exception as ex:
        print("[SWISSTOPO] drape: PictureFrame-create fail:", ex)
        return None
    if not pf_gid or pf_gid == System.Guid.Empty:
        print("[SWISSTOPO] drape: PictureFrame Empty-GUID")
        return None
    pf_obj = doc.Objects.Find(pf_gid)
    pf_mat_idx = pf_obj.Attributes.MaterialIndex
    # Mesh ins Doc + Material vom PictureFrame uebernehmen
    mesh_gid = doc.Objects.AddMesh(mesh)
    mesh_obj = doc.Objects.Find(mesh_gid)
    if mesh_obj is None:
        try: doc.Objects.Delete(pf_gid, True)
        except Exception: pass
        return None
    attrs = mesh_obj.Attributes.Duplicate()
    attrs.MaterialSource = Rhino.DocObjects.ObjectMaterialSource.MaterialFromObject
    attrs.MaterialIndex = pf_mat_idx
    if target_layer_idx >= 0:
        attrs.LayerIndex = target_layer_idx
    doc.Objects.ModifyAttributes(mesh_obj, attrs, True)
    # Temp-PictureFrame loeschen — das Mesh hat jetzt das Material
    try: doc.Objects.Delete(pf_gid, True)
    except Exception: pass
    print("[SWISSTOPO] drape mesh: {}x{} grid, {} faces, mat={}".format(
        len(es), len(ns), n_faces, pf_mat_idx))
    return mesh_obj


def add_ortho_plane(doc, ortho_path, tile_bbox_lv95, shift_lv95, m_to_unit,
                      z_doc=0.0, target_layer_idx=-1):
    """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
    # TIF direkt — Rhino's Picture-Pfad liest GeoTIFF nativ (NSImage auf Mac).
    # Behaelt die volle 10cm-Aufloesung statt auf 2k PNG runter zu skalieren.
    # 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
    # PictureFrame mit embedded Bitmap. AddPictureFrame interpretiert
    # plane.Origin als BOTTOM-LEFT corner der Picture (nicht als Zentrum!).
    # Width geht in +X-Richtung, Height in +Y-Richtung des Planes.
    width  = abs(x_max - x_min)
    height = abs(y_max - y_min)
    plane = rg.Plane(rg.Point3d(x_min, y_min, z_doc),
                      rg.Vector3d.XAxis, rg.Vector3d.YAxis)
    try:
        size_mb = os.path.getsize(ortho_path) / 1e6
        print("[SWISSTOPO] PictureFrame src: {} ({:.1f} MB)".format(
            os.path.basename(ortho_path), size_mb))
    except Exception:
        print("[SWISSTOPO] file nicht lesbar:", ortho_path)
        return None
    try:
        gid = doc.Objects.AddPictureFrame(
            plane, ortho_path,
            False,          # asMesh=False (Brep) — Mac Rhino 8 ignoriert die
                            # Plane bei asMesh=True, alle Pictures landen
                            # uebereinander
            width, height,
            True,           # selfIllumination=True — Textur unabhaengig von
                            # Lighting sichtbar (sonst evtl. dunkel in modes
                            # ohne Lichtquellen)
            False)          # embedBitmap=False — Pfad-Referenz (Cache bleibt
                            # persistent, kein 70MB-Embedding pro Tile)
        if gid == System.Guid.Empty:
            print("[SWISSTOPO] AddPictureFrame: Empty-GUID")
            return None
    except Exception as ex:
        print("[SWISSTOPO] AddPictureFrame exception:", ex)
        return None
    obj = doc.Objects.Find(gid)
    if obj is None: return None
    # Auf Ziel-Layer schieben (nachträglich; Material bleibt auf Object).
    if target_layer_idx >= 0:
        try:
            at = obj.Attributes.Duplicate()
            at.LayerIndex = target_layer_idx
            doc.Objects.ModifyAttributes(obj, at, True)
        except Exception as ex:
            print("[SWISSTOPO] Layer-Move fail:", ex)
    # Diagnose: hat das Material tatsaechlich eine Bitmap-Textur drin?
    try:
        o2 = doc.Objects.Find(gid)
        a = o2.Attributes
        print("[SWISSTOPO] PictureFrame OK id={} layer='{}' MatSrc={} MatIdx={} hidden={}".format(
            gid, doc.Layers[a.LayerIndex].FullPath,
            a.MaterialSource, a.MaterialIndex, o2.IsHidden))
        # Material-Inspect
        mat = None
        try:
            if a.MaterialIndex >= 0 and a.MaterialIndex < doc.Materials.Count:
                mat = doc.Materials[a.MaterialIndex]
        except Exception: pass
        if mat is not None:
            try: bmp_fn = mat.GetBitmapTexture().FileName if mat.GetBitmapTexture() else None
            except Exception: bmp_fn = None
            try: tex = mat.GetTexture(Rhino.DocObjects.TextureType.Bitmap)
            except Exception: tex = None
            print("[SWISSTOPO] material[{}].Name='{}' bitmap='{}' tex={} textures={}".format(
                a.MaterialIndex, mat.Name, bmp_fn, tex,
                mat.GetTextures().Length if hasattr(mat, "GetTextures") else "?"))
        # RenderMaterial-Inspect
        try:
            rm = a.RenderMaterial
            print("[SWISSTOPO] RenderMaterial: {}".format(
                rm.Name if rm else "None"))
        except Exception as ex:
            print("[SWISSTOPO] RenderMaterial-check fail:", ex)
    except Exception as ex:
        print("[SWISSTOPO] diag fail:", ex)
    return obj