Plugin_SN_Verfahrensgebiet/functions/grenzpunkt_alkis.py

"""sn_verfahrensgebiet/functions/grenzpunkt_alkis.py – ALKIS-Grenzpunktabfrage.

Lädt Grenzpunkte aus dem ALKIS-NAS-WFS in zwei Schritten:

1. ``adv:AX_PunktortTA`` per BBox laden.
   Diese Objektart enthält die eigentliche Koordinate (Geometrie) sowie das
   Attribut ``gehoertZu``, das per xlink-Referenz auf die GML-ID des
   zugehörigen ``AX_Grenzpunkt``-Objekts verweist.

2. GML-IDs aus ``gehoertZu`` extrahieren → ``AX_Grenzpunkt`` per
   ``gml_id IN (...)`` abfragen (batching, ~25 IDs/Batch).

3. Beide Layer werden über ``gehoertZu ↔ gml_id`` verknüpft (Join) und als
   einzelner angereicherter Layer zurückgegeben.

**Einschränkungen des WFS-Dienstes:**

- ``AX_Grenzpunkt`` hat **keine eigene Geometrie** → kein BBOX-Filter möglich.
- Das Attribut ``gehoertZu`` in ``AX_PunktortTA`` ist serverseitig **nicht
  filterbar** (bestätigt durch WFS-Fehlermeldung ``ErrorInPropertyName``).
- Der serverseitige GML-ID-Filter (``ResourceId``) ist aktiviert
  (``ImplementsResourceId=TRUE`` in den Capabilities).

**Strategie:** BBox → PunktortTA → GML-IDs → Grenzpunkt-Batches →
client-seitiger Join.
"""

from __future__ import annotations

import math
from typing import Callable, Dict, List, Optional, Tuple

from sn_basis.functions import QgsVectorLayer, QgsProject
from sn_basis.functions import QMessageBox
from sn_basis.functions.message_wrapper import info, warning
from sn_basis.functions.dialog_wrapper import create_progress_dialog

ALKIS_SF_URL = "https://geodienste.sachsen.de/aaa/public_alkis/sf/wfs"
PUNKTORT_TA_TYPENAME = "adv:AX_PunktortTA"
GRENZPUNKT_TYPENAME = "adv:AX_Grenzpunkt"

PUNKTORT_LAYER_NAME = "PunktortTA"
GRENZPUNKT_LAYER_NAME = "Grenzpunkte"
JOINED_LAYER_NAME = "Grenzpunkte (ALKIS)"

# Attributname in AX_PunktortTA, der auf AX_Grenzpunkt verweist
_ATTR_GEHOERT_ZU = "gehoertZu"
# Attributname für die GML-ID in AX_Grenzpunkt
_ATTR_GML_ID = "gml_id"

# Batch-Größe für ResourceId-Abfragen per GET
# WAF auf geodienste.sachsen.de blockiert POST und GML-BBOX-XML in URLs.
# GET+ResourceId-Filter funktioniert, aber URL-Länge ist begrenzt (~8192 Bytes).
# Bei ~54 Bytes/ID (URL-kodiert) passen sicher 100 IDs pro Request.
_BATCH_GROESSE = 100
# Maximale Features pro WFS-Abruf (Serverlimit)
_MAX_FEATURES = 20000

# Präfix für Grenzpunkt-Attribute im angereicherten Join-Layer
_GP_PREFIX = "gp_"

# Bekannte Feldnamen für die xlink-Referenz auf AX_Grenzpunkt (ALKIS 6 + 7)
# AX_PunktortTA NAS 6.x: "istTeilVon" verweist auf den übergeordneten Grenzpunkt
_GEHOERT_ZU_KANDIDATEN = ["istTeilVon", "gehoertZu", "gehoertZuAX_Grenzpunkt"]

# Vereinfachter ALKIS-WFS für Gemarkungs-BBox-Abfragen
_ALKIS_VEREINF_URL = "https://geodienste.sachsen.de/aaa/public_alkis/vereinf/wfs"
_KATASTERBEZIRK_TYPENAME = "ave:KatasterBezirk"
_KB_ATTR_SCHLUESSEL = "schluessel"


# ---------------------------------------------------------------------------
# Interne Hilfsfunktionen
# ---------------------------------------------------------------------------

def _erkenne_gehoert_zu_feld(layer) -> Optional[str]:
    """Erkennt den Feldnamen im ``AX_PunktortTA``-Layer, der auf
    ``AX_Grenzpunkt`` verweist – ALKIS-6/7-kompatibel.

    Prüft Kandidaten in fester Reihenfolge:

    1. Bekannte Feldnamen aus :data:`_GEHOERT_ZU_KANDIDATEN`.
    2. Fallback: erstes Feld, dessen Name ``gehoert`` (case-insensitive)
       enthält.

    Parameters
    ----------
    layer :
        Geladener ``AX_PunktortTA``-Layer.

    Returns
    -------
    str | None
        Gefundener Feldname oder ``None`` wenn kein Kandidat vorhanden.
    """
    feld_namen = [f.name() for f in layer.fields()]
    for kandidat in _GEHOERT_ZU_KANDIDATEN:
        if kandidat in feld_namen:
            return kandidat
    # Fallback: erstes Feld, das 'gehoert' oder 'istteil' oder 'besteht' enthält
    for name in feld_namen:
        low = name.lower()
        if "gehoert" in low or "istteil" in low or "bestehtaus" in low:
            return name
    warning(
        "Grenzpunkte (Diagnose)",
        "Kein Referenzfeld in AX_PunktortTA gefunden. "
        "Vorhandene Felder: " + ", ".join(feld_namen),
        duration=15,
    )
    return None


def _punktort_ta_bbox_param(bbox) -> str:
    """Gibt den ``BBOX``-KVP-Wert für eine PunktortTA-Abfrage zurück.

    **Warum BBOX-Parameter statt FILTER-XML?**
    Die WAF auf geodienste.sachsen.de blockiert URLs, die GML-Geometrie-XML
    im ``FILTER``-Parameter enthalten (``xmlns:gml``, ``gml:Envelope`` usw.).
    Der native WFS-2.0-KVP-Parameter ``BBOX=minx,miny,maxx,maxy,crs`` wird
    dagegen akzeptiert.

    Parameters
    ----------
    bbox :
        ``QgsRectangle`` in EPSG:25833.

    Returns
    -------
    str
        Wert für den ``BBOX``-URL-Parameter.
    """
    return (
        f"{bbox.xMinimum()},{bbox.yMinimum()}"
        f",{bbox.xMaximum()},{bbox.yMaximum()}"
        ",urn:ogc:def:crs:EPSG::25833"
    )


def _wfs_get_gml(url: str, typename: str, extra_params: dict) -> Optional[bytes]:
    """Sendet einen WFS 2.0 GetFeature GET-Request und gibt die rohe GML-Antwort zurück.

    **Warum GET?**
    Die WAF auf geodienste.sachsen.de blockiert POST-Requests und GET-Requests
    mit GML-Geometrie-XML im ``FILTER``-Parameter (``xmlns:gml``, ``gml:Envelope``).
    GET-Requests mit ``BBOX=``-KVP-Parameter oder ResourceId-``FILTER`` (nur
    ``xmlns:fes``, kein GML) werden dagegen akzeptiert.

    Parameters
    ----------
    url :
        WFS-Endpoint-URL.
    typename :
        Vollqualifizierter Feature-Typ, z.B. ``adv:AX_PunktortTA``.
    extra_params :
        Zusätzliche KVP-Parameter, z.B. ``{'BBOX': 'minx,miny,maxx,maxy,crs'}``
        oder ``{'FILTER': '<fes:Filter...>'}``.

    Returns
    -------
    bytes | None
        Rohe GML-Antwort oder ``None`` bei Fehler.
    """
    import urllib.request
    import urllib.parse
    import urllib.error
    from qgis.core import QgsMessageLog, Qgis

    params = {
        'SERVICE': 'WFS',
        'VERSION': '2.0.0',
        'REQUEST': 'GetFeature',
        'TYPENAMES': typename,
        'COUNT': str(_MAX_FEATURES),
    }
    params.update(extra_params)

    full_url = url + '?' + urllib.parse.urlencode(params)

    try:
        with urllib.request.urlopen(full_url, timeout=120) as resp:
            gml_bytes = resp.read()
    except urllib.error.HTTPError as exc:
        body_preview = exc.read(200).decode('utf-8', errors='replace')
        QgsMessageLog.logMessage(
            f"WFS GET HTTP {exc.code} ({typename}): {exc.reason}\n"
            f"URL: {full_url[:300]}\nAntwort: {body_preview}",
            "sn_verfahrensgebiet",
            Qgis.MessageLevel.Critical,
        )
        warning("WFS GET", f"HTTP {exc.code}: {exc.reason}", duration=10)
        return None
    except Exception as exc:
        QgsMessageLog.logMessage(
            f"WFS GET fehlgeschlagen ({typename}): {exc}\nURL: {full_url[:300]}",
            "sn_verfahrensgebiet",
            Qgis.MessageLevel.Critical,
        )
        warning("WFS GET", f"Netzwerkfehler: {exc}", duration=10)
        return None

    if not gml_bytes:
        return None

    # WFS-Fehlermeldung erkennen
    if b"ExceptionReport" in gml_bytes[:800] or b"ExceptionText" in gml_bytes[:800]:
        import xml.etree.ElementTree as _ET
        try:
            _root = _ET.fromstring(gml_bytes)
            _exc = _root.find(".//{http://www.opengis.net/ows/1.1}ExceptionText")
            msg = _exc.text.strip() if _exc is not None and _exc.text else "Unbekannter WFS-Fehler"
        except Exception:
            msg = gml_bytes[:500].decode("utf-8", errors="replace")
        QgsMessageLog.logMessage(
            f"WFS ExceptionReport ({typename}):\n{msg}",
            "sn_verfahrensgebiet",
            Qgis.MessageLevel.Critical,
        )
        warning("WFS GET", f"WFS-Fehler: {msg}", duration=10)
        return None

    return gml_bytes


def _parse_alkis_gml_memory_layer(
    gml_bytes: bytes,
    feature_local_name: str,
    layer_name: str,
    geom_prop: Optional[str] = "position",
) -> Optional[object]:
    """Parst ALKIS-GML-Bytes und erstellt einen QGIS Memory-Layer.

    Verarbeitet die GML-Antwort eines WFS-POST-Requests und extrahiert
    Features des angegebenen Typs als Memory-Layer.

    Parameters
    ----------
    gml_bytes :
        Rohe GML-Bytes der WFS-Antwort.
    feature_local_name :
        Lokaler Name des Feature-Typs, z.B. ``AX_PunktortTA``.
    layer_name :
        Anzeigename des resultierenden Memory-Layers.
    geom_prop :
        Lokaler Name des Geometrie-Attributs (``None`` = kein Geometrie-Layer).

    Returns
    -------
    QgsVectorLayer | None
        Memory-Layer oder ``None`` wenn keine Features gefunden.
    """
    import xml.etree.ElementTree as ET
    from sn_basis.functions import QgsFeature, QgsField
    from sn_basis.functions.qt_wrapper import QVariant
    from qgis.core import QgsGeometry, QgsPointXY

    _ADV = "http://www.adv-online.de/namespaces/adv/gid/6.0"
    _GML = "http://www.opengis.net/gml/3.2"
    _XLINK = "http://www.w3.org/1999/xlink"

    try:
        root = ET.fromstring(gml_bytes)
    except ET.ParseError as exc:
        warning("GML-Parser", f"GML konnte nicht geparst werden: {exc}", duration=8)
        return None

    feat_tag = f"{{{_ADV}}}{feature_local_name}"
    feat_elems = root.findall(f".//{feat_tag}")
    if not feat_elems:
        return None

    # Felder aus ALLEN Features ermitteln (2-Pass), damit optionale Felder
    # wie "istTeilVon" nicht verloren gehen, wenn sie im ersten Feature fehlen.
    feld_namen: List[str] = ["gml_id"]
    for feat_el in feat_elems:
        for child in feat_el:
            local = child.tag.split("}")[-1] if "}" in child.tag else child.tag
            if local not in feld_namen and local != geom_prop:
                feld_namen.append(local)

    # Memory-Layer anlegen
    mem_uri = f"Point?crs=EPSG:25833" if geom_prop else "none"
    mem_layer = QgsVectorLayer(mem_uri, layer_name, "memory")
    pr = mem_layer.dataProvider()
    pr.addAttributes([QgsField(n, QVariant.String) for n in feld_namen if n])
    mem_layer.updateFields()

    neue_features: List = []
    for feat_elem in feat_elems:
        gml_id = feat_elem.get(f"{{{_GML}}}id", "")

        # Geometrie extrahieren
        geom = None
        if geom_prop:
            pos_el = feat_elem.find(f"{{{_ADV}}}{geom_prop}")
            if pos_el is not None:
                pe = pos_el.find(f".//{{{_GML}}}pos")
                if pe is not None and pe.text:
                    coords = pe.text.strip().split()
                    if len(coords) >= 2:
                        try:
                            geom = QgsGeometry.fromPointXY(
                                QgsPointXY(float(coords[0]), float(coords[1]))
                            )
                        except (ValueError, Exception):
                            pass

        qgs_feat = QgsFeature(mem_layer.fields())
        if geom:
            qgs_feat.setGeometry(geom)
        qgs_feat["gml_id"] = gml_id

        for child in feat_elem:
            local = child.tag.split("}")[-1] if "}" in child.tag else child.tag
            if local == geom_prop or local not in feld_namen:
                continue
            href = child.get(f"{{{_XLINK}}}href", "")
            val = href if href else (child.text.strip() if child.text else "")
            try:
                qgs_feat[local] = val
            except Exception:
                pass

        neue_features.append(qgs_feat)

    pr.addFeatures(neue_features)
    mem_layer.updateExtents()
    return mem_layer


def _extrahiere_feature_gml_id(feat, feld_namen: List[str]) -> str:
    """Extrahiert die GML-ID aus einem WFS-Feature robust.

    Versucht mehrere Strategien in dieser Reihenfolge:

    1. Bekannte Feldnamen (``gml_id``, ``identifier``).
    2. Scan aller String-Felder auf GML-ID-Muster (``DE`` + 14-20 alphanumerische Zeichen).

    Parameters
    ----------
    feat :
        WFS-Feature.
    feld_namen :
        Liste der Feldnamen des Layers.

    Returns
    -------
    str
        GML-ID oder leerer String.
    """
    import re as _re
    _GML_ID_PATTERN = _re.compile(r'^DE[0-9A-Za-z_]{10,24}$')

    for candidate in ("gml_id", "identifier", "gmlId", "featureId", "fid"):
        if candidate in feld_namen:
            try:
                val = feat[candidate]
                if val:
                    s = str(val).strip()
                    if s:
                        return s
            except Exception:
                pass

    # Scan aller Felder auf GML-ID-Muster
    for name in feld_namen:
        try:
            val = feat[name]
            if val:
                s = str(val).strip()
                if _GML_ID_PATTERN.match(s):
                    return s
        except Exception:
            pass

    return ""


def _extrahiere_gml_ids(
    punktort_layer,
    gehoert_zu_feld: Optional[str] = None,
) -> List[str]:
    """Extrahiert GML-IDs aus dem Referenzfeld des PunktortTA-Layers.

    Der WFS gibt das Referenzfeld üblicherweise als ``xlink:href``-Referenz im
    Format ``urn:adv:oid:DEXXXXXXXXXXXXXXXXXX`` oder als flachen Wert.
    Es wird der letzte Bestandteil nach dem letzten ``:`` verwendet.

    Parameters
    ----------
    punktort_layer :
        Geladener ``AX_PunktortTA``-WFS-Layer.
    gehoert_zu_feld :
        Feldname der Referenz auf ``AX_Grenzpunkt``. Wenn ``None``, wird
        :data:`_ATTR_GEHOERT_ZU` als Fallback verwendet
        (Rückwärtskompatibilität).

    Returns
    -------
    List[str]
        Deduplizierte, sortierte Liste von GML-IDs.
    """
    feld = gehoert_zu_feld if gehoert_zu_feld is not None else _ATTR_GEHOERT_ZU
    ids: set[str] = set()
    for feat in punktort_layer.getFeatures():
        try:
            val = feat[feld]
            if val:
                raw = str(val).strip()
                # xlink:href hat Format "urn:adv:oid:DE..." oder reiner GML-ID-Wert
                gml_id = raw.rsplit(":", 1)[-1] if ":" in raw else raw
                if gml_id:
                    ids.add(gml_id)
        except (KeyError, Exception):
            continue
    return sorted(ids)


def _lade_grenzpunkte_als_dict(
    id_liste: List[str],
    progress_callback: Optional[Callable[[int], None]] = None,
) -> Tuple[Dict[str, object], List[str]]:
    """Lädt ``AX_Grenzpunkt``-Features als Dict für O(1)-Lookup.

    Parameters
    ----------
    id_liste :
        GML-IDs der gesuchten ``AX_Grenzpunkt``-Objekte.
    progress_callback :
        Optionaler Callback ``fn(batch_index)`` der nach jedem Batch aufgerufen
        wird (0-basierter Batch-Index). Dient zur Fortschrittsanzeige.

    Returns
    -------
    Tuple[Dict[str, QgsFeature], List[str]]
        ``(gp_dict, feld_namen)`` – Dict ``{gml_id: QgsFeature}`` und
        Liste der GP-Feldnamen (leer wenn keine Features geladen wurden).
    """
    gp_dict: Dict[str, object] = {}
    feld_namen: List[str] = []

    anzahl_batches = math.ceil(len(id_liste) / _BATCH_GROESSE)
    for i in range(anzahl_batches):
        batch = id_liste[i * _BATCH_GROESSE: (i + 1) * _BATCH_GROESSE]
        info(
            "Grenzpunkte (Diagnose)",
            f"AX_Grenzpunkt Batch {i + 1}/{anzahl_batches}: {len(batch)} IDs per ResourceId-Filter …",
            duration=2,
        )

        rid_elemente = "".join(f'<fes:ResourceId rid="{gid}"/>' for gid in batch)
        filter_xml = (
            '<fes:Filter xmlns:fes="http://www.opengis.net/fes/2.0">'
            + rid_elemente
            + "</fes:Filter>"
        )
        _gml_bytes = _wfs_get_gml(ALKIS_SF_URL, GRENZPUNKT_TYPENAME, {'FILTER': filter_xml})
        if not _gml_bytes:
            warning(
                "Grenzpunkte (Diagnose)",
                f"Batch {i + 1}/{anzahl_batches}: GET-Anfrage fehlgeschlagen – Batch wird übersprungen.",
                duration=4,
            )
        else:
            tmp = _parse_alkis_gml_memory_layer(_gml_bytes, "AX_Grenzpunkt", "_TmpGP", geom_prop=None)
            if not tmp:
                warning(
                    "Grenzpunkte (Diagnose)",
                    f"Batch {i + 1}/{anzahl_batches}: Keine Features in GML-Antwort – Batch wird übersprungen.",
                    duration=4,
                )
            else:
                batch_feld_namen = [f.name() for f in tmp.fields()]
                if not feld_namen and batch_feld_namen:
                    feld_namen = batch_feld_namen
                    info(
                        "Grenzpunkte (Diagnose)",
                        f"AX_Grenzpunkt-Felder: {', '.join(feld_namen)}",
                        duration=8,
                    )

                batch_count = 0
                for feat in tmp.getFeatures():
                    gml_id = _extrahiere_feature_gml_id(feat, batch_feld_namen or feld_namen)
                    if gml_id:
                        gp_dict[gml_id] = feat
                        batch_count += 1

                info(
                    "Grenzpunkte (Diagnose)",
                    f"Batch {i + 1}/{anzahl_batches}: {batch_count}/{len(batch)} Grenzpunkte empfangen "
                    f"(gesamt bisher: {len(gp_dict)}).",
                    duration=2,
                )

        if progress_callback is not None:
            progress_callback(i)

    return gp_dict, feld_namen


def _lade_grenzpunkte_per_ids(id_liste: List[str]) -> Optional[object]:
    """Lädt ``AX_Grenzpunkt``-Features für die angegebenen GML-IDs.

    Die ID-Liste wird in Batches der Größe :data:`_BATCH_GROESSE` aufgeteilt,
    um URL-Längenlimits zu umgehen. Alle Batch-Ergebnisse werden zu einem
    Memory-Layer zusammengeführt.

    Delegiert intern an :func:`_lade_grenzpunkte_als_dict`.

    Parameters
    ----------
    id_liste :
        GML-IDs der gesuchten ``AX_Grenzpunkt``-Objekte.

    Returns
    -------
    QgsVectorLayer | None
        Memory-Layer mit allen gefundenen Grenzpunkten oder ``None`` bei Fehler.
    """
    from sn_basis.functions import QgsField
    from sn_basis.functions.qt_wrapper import QVariant

    gp_dict, feld_namen = _lade_grenzpunkte_als_dict(id_liste)
    if not gp_dict:
        return None

    alle_features = list(gp_dict.values())

    # Ergebnis-Memory-Layer aufbauen (ohne Geometrie, da AX_Grenzpunkt keine hat)
    mem_uri = "Point?crs=EPSG:25833" if _hat_geometrie(alle_features[0]) else "none"
    result_layer = QgsVectorLayer(mem_uri, GRENZPUNKT_LAYER_NAME, "memory")
    pr = result_layer.dataProvider()

    felder = [QgsField(n, QVariant.String) for n in feld_namen if n]
    pr.addAttributes(felder)
    result_layer.updateFields()
    pr.addFeatures(alle_features)
    result_layer.updateExtents()

    return result_layer


def _hat_geometrie(feature) -> bool:
    """Gibt ``True`` zurück, wenn das Feature eine nicht-leere Geometrie hat."""
    try:
        g = feature.geometry()
        return g is not None and not g.isEmpty()
    except Exception:
        return False


def _lade_gemarkung_extents(schluessel_liste: List[str]) -> Dict[str, object]:
    """Lädt BBox-Rechtecke der Gemarkungen aus dem vereinfachten ALKIS-WFS.

    Da der vereinfachte WFS keine serverseitigen Filter unterstützt, werden
    alle Katasterbezirke geladen und client-seitig nach den angegebenen
    Gemarkungsnummern gefiltert. Der vollständige Katasterbezirksschlüssel hat
    ≥ 4 Stellen; die letzten 4 Stellen entsprechen der Gemarkungsnummer.

    Parameters
    ----------
    schluessel_liste :
        Liste der 4-stelligen Gemarkungsnummern.

    Returns
    -------
    Dict[str, QgsRectangle]
        ``{gemarkungsnummer: QgsRectangle}`` – BBox je Gemarkung.
        Leer bei Fehler oder keinen Treffern.
    """
    from qgis.core import QgsRectangle

    uri = (
        f"url='{_ALKIS_VEREINF_URL}' typename='{_KATASTERBEZIRK_TYPENAME}' "
        f"srsname='EPSG:25833'"
    )
    layer = QgsVectorLayer(uri, "_TmpKB", "WFS")
    if not layer or not layer.isValid():
        return {}

    schluessel_set = set(schluessel_liste)
    extents: Dict[str, object] = {}

    for feat in layer.getFeatures():
        try:
            val = feat[_KB_ATTR_SCHLUESSEL]
            schluessel = str(val).strip() if val else ""
        except Exception:
            continue

        nummer = schluessel[-4:] if len(schluessel) >= 4 else schluessel
        if nummer not in schluessel_set:
            continue

        if not feat.hasGeometry() or feat.geometry().isEmpty():
            continue

        bbox = feat.geometry().boundingBox()
        if nummer in extents:
            extents[nummer].combineExtentWith(bbox)
        else:
            extents[nummer] = QgsRectangle(bbox)

    return extents


def _lade_gemarkung_extents_mit_puffer(
    schluessel_liste: List[str],
    vg_geom,
    puffer_m: float = 200.0,
) -> Dict[str, object]:
    """Lädt Gemarkungsgeometrien, verschneidet sie mit einem Puffer um das
    Verfahrensgebiet und liefert BBoxen der Schnittflächen.

    Statt der vollen Gemarkungs-BBox wird nur die BBox des Anteils verwendet,
    der innerhalb des gepufferten Verfahrensgebiets liegt. Das reduziert die
    Anzahl der vom WFS gelieferten PunktortTA-Objekte erheblich.

    Parameters
    ----------
    schluessel_liste :
        Liste der 4-stelligen Gemarkungsnummern.
    vg_geom :
        ``QgsGeometry`` des Verfahrensgebiets (EPSG:25833).
    puffer_m :
        Pufferdistanz in Metern (Standard: 200 m).

    Returns
    -------
    Dict[str, QgsRectangle]
        ``{gemarkungsnummer: QgsRectangle}`` – BBox der Schnittfläche.
        Gemarkungen ohne Schnitt werden ausgelassen.
    """
    from qgis.core import QgsRectangle

    uri = (
        f"url='{_ALKIS_VEREINF_URL}' typename='{_KATASTERBEZIRK_TYPENAME}' "
        f"srsname='EPSG:25833'"
    )
    layer = QgsVectorLayer(uri, "_TmpKB", "WFS")
    if not layer or not layer.isValid():
        return {}

    puffer_geom = vg_geom.buffer(puffer_m, 16) if vg_geom and not vg_geom.isEmpty() else None

    schluessel_set = set(schluessel_liste)
    extents: Dict[str, object] = {}

    for feat in layer.getFeatures():
        try:
            val = feat[_KB_ATTR_SCHLUESSEL]
            schluessel = str(val).strip() if val else ""
        except Exception:
            continue

        nummer = schluessel[-4:] if len(schluessel) >= 4 else schluessel
        if nummer not in schluessel_set:
            continue

        if not feat.hasGeometry() or feat.geometry().isEmpty():
            continue

        geom = feat.geometry()

        # Mit gepuffertem VG verschneiden
        if puffer_geom and not puffer_geom.isEmpty():
            schnitt = geom.intersection(puffer_geom)
            if schnitt is None or schnitt.isEmpty():
                continue
            bbox = schnitt.boundingBox()
        else:
            bbox = geom.boundingBox()

        if bbox.isEmpty():
            continue

        if nummer in extents:
            extents[nummer].combineExtentWith(bbox)
        else:
            extents[nummer] = QgsRectangle(bbox)

    return extents


def _wfs_hits_count(url: str, typename: str, extra_params: dict) -> int:
    """Fragt die Anzahl der Features per ``resultType=hits`` ab.

    Parameters
    ----------
    url :
        WFS-Endpoint-URL.
    typename :
        Vollqualifizierter Feature-Typ.
    extra_params :
        Zus\u00e4tzliche KVP-Parameter (z.B. ``{'BBOX': '...'}``).

    Returns
    -------
    int
        Anzahl der Features oder ``-1`` wenn nicht ermittelbar.
    """
    import urllib.request
    import urllib.parse
    import urllib.error
    import xml.etree.ElementTree as ET

    params = {
        'SERVICE': 'WFS',
        'VERSION': '2.0.0',
        'REQUEST': 'GetFeature',
        'TYPENAMES': typename,
        'RESULTTYPE': 'hits',
    }
    params.update(extra_params)

    full_url = url + '?' + urllib.parse.urlencode(params)
    try:
        with urllib.request.urlopen(full_url, timeout=30) as resp:
            gml_bytes = resp.read()
    except Exception:
        return -1

    try:
        root = ET.fromstring(gml_bytes)
        matched = root.get('numberMatched', '')
        if matched and matched != 'unknown':
            return int(matched)
        returned = root.get('numberReturned', '')
        if returned and returned != 'unknown':
            return int(returned)
    except Exception:
        pass
    return -1


def _lade_punktort_ta_batches_per_gemarkung(
    schluessel_liste: List[str],
) -> Optional[object]:
    """Lädt ``AX_PunktortTA`` für jede Gemarkung separat per BBox.

    Pro Gemarkung wird eine eigene WFS-Abfrage durchgeführt, was präzisere
    Ergebnisse als eine kombinierte BBox liefert. Doppelte Features an
    Gemarkungsgrenzen werden anhand der ``gml_id`` dedupliziert.

    Parameters
    ----------
    schluessel_liste :
        Liste der 4-stelligen Gemarkungsnummern.

    Returns
    -------
    QgsVectorLayer | None
        Memory-Layer mit allen PunktortTA-Features oder ``None`` bei Fehler.
    """
    from sn_basis.functions import QgsField
    from sn_basis.functions.qt_wrapper import QVariant

    extents = _lade_gemarkung_extents(schluessel_liste)
    if not extents:
        return None

    alle_features: List = []
    gesehene_ids: set[str] = set()
    feld_namen: List[str] = []

    n_gemarkungen_gesamt = len(extents)
    for g_idx, (_gemarkungsnummer, bbox) in enumerate(extents.items(), 1):
        info(
            "Grenzpunkte (Diagnose)",
            f"PunktortTA Gemarkung {_gemarkungsnummer} ({g_idx}/{n_gemarkungen_gesamt}): "
            f"BBox {bbox.xMinimum():.0f},{bbox.yMinimum():.0f} – "
            f"{bbox.xMaximum():.0f},{bbox.yMaximum():.0f} anfragen …",
            duration=2,
        )
        _gml_bytes = _wfs_get_gml(ALKIS_SF_URL, PUNKTORT_TA_TYPENAME, {'BBOX': _punktort_ta_bbox_param(bbox)})
        if not _gml_bytes:
            warning(
                "Grenzpunkte (Diagnose)",
                f"Gemarkung {_gemarkungsnummer}: GET-Anfrage fehlgeschlagen – wird übersprungen.",
                duration=4,
            )
            continue
        tmp = _parse_alkis_gml_memory_layer(_gml_bytes, "AX_PunktortTA", "_TmpPTA", geom_prop="position")
        if not tmp:
            warning(
                "Grenzpunkte (Diagnose)",
                f"Gemarkung {_gemarkungsnummer}: Keine Features in GML-Antwort – wird übersprungen.",
                duration=4,
            )
            continue

        if not feld_namen:
            feld_namen = [f.name() for f in tmp.fields()]

        cur_feld_namen = [f.name() for f in tmp.fields()]

        neu = 0
        duplikat = 0
        for feat in tmp.getFeatures():
            gml_id = _extrahiere_feature_gml_id(feat, cur_feld_namen)
            if gml_id and gml_id in gesehene_ids:
                duplikat += 1
                continue
            if gml_id:
                gesehene_ids.add(gml_id)
            alle_features.append(feat)
            neu += 1

        info(
            "Grenzpunkte (Diagnose)",
            f"Gemarkung {_gemarkungsnummer}: {neu} neue PunktortTA-Objekte "
            f"({duplikat} Duplikate übersprungen, gesamt bisher: {len(alle_features)}).",
            duration=2,
        )

    if not alle_features:
        return None

    result_layer = QgsVectorLayer("Point?crs=EPSG:25833", PUNKTORT_LAYER_NAME, "memory")
    pr = result_layer.dataProvider()

    felder = [QgsField(n, QVariant.String) for n in feld_namen if n]
    pr.addAttributes(felder)
    result_layer.updateFields()
    pr.addFeatures(alle_features)
    result_layer.updateExtents()

    return result_layer


def _erstelle_joined_layer(
    punktort_layer,
    gp_dict: Dict[str, object],
    gp_feld_namen: List[str],
    gehoert_zu_feld: str,
) -> object:
    """Erstellt einen angereicherten Memory-Layer aus PunktortTA + GP-Attributen.

    Geometrie und Attribute von ``AX_PunktortTA`` werden mit den Attributen
    des zugehörigen ``AX_Grenzpunkt`` zusammengeführt. GP-Felder erhalten den
    Präfix :data:`_GP_PREFIX` (``gp_``), um Namenskollisionen zu vermeiden.
    Features ohne GP-Match werden mit leeren ``gp_``-Feldern aufgenommen.

    Parameters
    ----------
    punktort_layer :
        Memory-Layer mit ``AX_PunktortTA``-Features.
    gp_dict :
        ``{gml_id: QgsFeature}`` der geladenen ``AX_Grenzpunkt``-Features.
    gp_feld_namen :
        Feldnamen aus dem ``AX_Grenzpunkt``-WFS.
    gehoert_zu_feld :
        Feldname in ``AX_PunktortTA``, der die xlink-Referenz enthält.

    Returns
    -------
    QgsVectorLayer
        Angereicherter Memory-Layer ``Point?crs=EPSG:25833``.
    """
    from sn_basis.functions import QgsFeature, QgsField
    from sn_basis.functions.qt_wrapper import QVariant

    pta_feld_namen = [f.name() for f in punktort_layer.fields()]

    result_layer = QgsVectorLayer("Point?crs=EPSG:25833", JOINED_LAYER_NAME, "memory")
    pr = result_layer.dataProvider()

    felder: List = []
    for n in pta_feld_namen:
        if n:
            felder.append(QgsField(n, QVariant.String))
    for n in gp_feld_namen:
        if n:
            felder.append(QgsField(f"{_GP_PREFIX}{n}", QVariant.String))

    pr.addAttributes(felder)
    result_layer.updateFields()

    neue_features: List = []
    for pta_feat in punktort_layer.getFeatures():
        try:
            raw = pta_feat[gehoert_zu_feld]
            raw = str(raw).strip() if raw else ""
            gml_id = raw.rsplit(":", 1)[-1] if ":" in raw else raw
        except Exception:
            gml_id = ""

        gp_feat = gp_dict.get(gml_id)

        new_feat = QgsFeature(result_layer.fields())
        new_feat.setGeometry(pta_feat.geometry())

        # PunktortTA-Attribute
        for n in pta_feld_namen:
            if n:
                try:
                    new_feat[n] = pta_feat[n]
                except Exception:
                    pass

        # Grenzpunkt-Attribute mit gp_-Präfix
        if gp_feat is not None:
            for n in gp_feld_namen:
                if n:
                    try:
                        new_feat[f"{_GP_PREFIX}{n}"] = gp_feat[n]
                    except Exception:
                        pass

        neue_features.append(new_feat)

    pr.addFeatures(neue_features)
    result_layer.updateExtents()
    return result_layer


def lade_punktort_ta(
    schluessel_liste: List[str],
    vg_geom,
    progress=None,
) -> Optional[object]:
    """Lädt ``AX_PunktortTA`` gemarkungsweise mit STARTINDEX-Paginierung.

    Verwendet dieselbe Strategie wie :func:`grenzpunkte_gemarkungsweise`:

    1. Gemarkungsextents laden und mit einem 200-m-VG-Puffer verschneiden.
    2. Hits-Abfrage pro Gemarkung zur Batch-Planung.
    3. PunktortTA je Gemarkung in 5000-Feature-Batches laden (STARTINDEX-Paging).
    4. gml_id-Deduplizierung über Gemarkungsgrenzen hinweg.

    Es werden keine Grenzpunkt-Batches geladen – die Funktion gibt ausschließlich
    den PunktortTA-Layer zurück.

    Parameters
    ----------
    schluessel_liste :
        4-stellige Gemarkungsnummern (aus dem aktuellen Workflow bereits bekannt).
    vg_geom :
        ``QgsGeometry`` des BROR-Objekts (EPSG:25833). Wird für den 200-m-Puffer
        zur Einschränkung der Gemarkungsextents verwendet.
    progress :
        Optionales Fortschrittsdialog-Objekt mit ``set_label()``, ``set_value()``
        und ``is_canceled()``.

    Returns
    -------
    QgsVectorLayer | None
        Memory-Layer mit allen PunktortTA-Features oder ``None`` bei Fehler.
    """
    from sn_basis.functions import QgsField
    from sn_basis.functions import QgsFeature as _QgsFeature
    from sn_basis.functions.qt_wrapper import QVariant

    _BBOX_BATCH = 5000

    if not schluessel_liste:
        return None

    if vg_geom and not vg_geom.isEmpty():
        extents = _lade_gemarkung_extents_mit_puffer(schluessel_liste, vg_geom, 200.0)
    else:
        extents = _lade_gemarkung_extents(schluessel_liste)

    if not extents:
        return None

    n_extents = len(extents)
    alle_records: List[Dict] = []
    gesehene_ids: set = set()
    feld_namen: List[str] = []

    for g_idx, (gemarkungsnummer, bbox) in enumerate(extents.items(), 1):
        if progress is not None and progress.is_canceled():
            return None

        bbox_param = _punktort_ta_bbox_param(bbox)
        n_hits = _wfs_hits_count(ALKIS_SF_URL, PUNKTORT_TA_TYPENAME, {'BBOX': bbox_param})
        n_batches = max(1, math.ceil(n_hits / _BBOX_BATCH)) if n_hits > 0 else 1

        info(
            "Grenzpunkte (Diagnose)",
            f"Gemarkung {gemarkungsnummer} ({g_idx}/{n_extents}): "
            f"~{n_hits} PunktortTA → {n_batches} Batch(es).",
            duration=2,
        )

        for b_idx in range(n_batches):
            if progress is not None:
                progress.set_label(
                    f"PunktortTA: Gemarkung {gemarkungsnummer} "
                    f"({g_idx}/{n_extents}), Batch {b_idx + 1}/{n_batches} …"
                )
                if progress.is_canceled():
                    return None

            extra: dict = {'BBOX': bbox_param}
            if n_batches > 1:
                extra['COUNT'] = str(_BBOX_BATCH)
                extra['STARTINDEX'] = str(b_idx * _BBOX_BATCH)

            _gml_bytes = _wfs_get_gml(ALKIS_SF_URL, PUNKTORT_TA_TYPENAME, extra)
            if not _gml_bytes:
                warning(
                    "Grenzpunkte (Diagnose)",
                    f"Gemarkung {gemarkungsnummer} Batch {b_idx + 1}: "
                    f"GET-Anfrage fehlgeschlagen – übersprungen.",
                    duration=4,
                )
                continue

            tmp = _parse_alkis_gml_memory_layer(
                _gml_bytes, "AX_PunktortTA", "_TmpPTA", geom_prop="position"
            )
            if not tmp:
                continue

            cur_felder = [f.name() for f in tmp.fields()]
            for n in cur_felder:
                if n not in feld_namen:
                    feld_namen.append(n)

            for feat in tmp.getFeatures():
                gml_id = _extrahiere_feature_gml_id(feat, cur_felder)
                if gml_id and gml_id in gesehene_ids:
                    continue
                if gml_id:
                    gesehene_ids.add(gml_id)
                attrs: Dict = {}
                for n in cur_felder:
                    if n:
                        try:
                            attrs[n] = feat[n]
                        except Exception:
                            pass
                alle_records.append({'geom': feat.geometry(), 'attrs': attrs})

        info(
            "Grenzpunkte (Diagnose)",
            f"Gemarkung {gemarkungsnummer}: {len(gesehene_ids)} PunktortTA bisher.",
            duration=2,
        )

    if not alle_records:
        return None

    result = QgsVectorLayer("Point?crs=EPSG:25833", PUNKTORT_LAYER_NAME, "memory")
    pr = result.dataProvider()
    pr.addAttributes([QgsField(n, QVariant.String) for n in feld_namen if n])
    result.updateFields()

    neue_features = []
    for record in alle_records:
        f = _QgsFeature(result.fields())
        f.setGeometry(record['geom'])
        for n, val in record['attrs'].items():
            if n in feld_namen:
                try:
                    f[n] = val
                except Exception:
                    pass
        neue_features.append(f)

    pr.addFeatures(neue_features)
    result.updateExtents()
    return result


def grenzpunkt_alkis(tab_widget, extent=None) -> Tuple[Optional[object], Optional[object]]:
    """Lädt Grenzpunkte aus dem ALKIS-WFS für den angegebenen Begrenzungsrahmen.

    Ablauf:

    1. ``AX_PunktortTA`` per BBox laden (enthält Koordinaten + Referenz auf
       ``AX_Grenzpunkt`` via ``gehoertZu``).
    2. GML-IDs aus ``gehoertZu`` extrahieren.
    3. ``AX_Grenzpunkt`` per ID-Liste in Batches laden.
    4. Beide Layer ins Projekt einfügen. Der Caller ist verantwortlich dafür,
       den Join in der QGIS-Projektkonfiguration einzurichten oder die
       Verknüpfung client-seitig vorzunehmen.

    Wenn kein ``extent`` übergeben wird, fragt die Funktion nach dem aktuellen
    Verfahrensgebiet-Layer und nutzt dessen BBox.

    Parameters
    ----------
    tab_widget :
        Eltern-Widget für Dialoge.
    extent :
        Optional. ``QgsRectangle`` in EPSG:25833. Wenn ``None``, wird der
        aktive Verfahrensgebiet-Layer gesucht und dessen Extent genutzt.

    Returns
    -------
    Tuple[Optional[punktort_layer], Optional[grenzpunkt_layer]]
        Beide geladenen Layer (bereits im Projekt) oder ``(None, None)``
        bei Fehler.
    """
    project = QgsProject.instance()

    # BBox bestimmen
    if extent is None:
        vg_layers = project.mapLayersByName("Verfahrensgebiet")
        if not vg_layers:
            QMessageBox.warning(
                tab_widget,
                "Kein Verfahrensgebiet",
                "Bitte zuerst ein Verfahrensgebiet laden, um die BBox abzuleiten.",
            )
            return None, None
        extent = vg_layers[0].extent()

    xmin = extent.xMinimum()
    ymin = extent.yMinimum()
    xmax = extent.xMaximum()
    ymax = extent.yMaximum()

    # -----------------------------------------------------------------------
    # Schritt 1: AX_PunktortTA per BBox laden
    # -----------------------------------------------------------------------
    info("Grenzpunkte", "Schritt 1/3: PunktortTA per BBox laden …", duration=3)
    uri_pta = (
        f"url='{ALKIS_SF_URL}' typename='{PUNKTORT_TA_TYPENAME}' "
        f"srsname='EPSG:25833' "
        f"bbox={xmin},{ymin},{xmax},{ymax},EPSG:25833"
    )
    punktort_layer = QgsVectorLayer(uri_pta, PUNKTORT_LAYER_NAME, "WFS")

    if not punktort_layer or not punktort_layer.isValid():
        QMessageBox.critical(
            tab_widget,
            "Fehler",
            "PunktortTA-Layer konnte nicht geladen werden.\n"
            "Bitte Netzwerkverbindung und Ausdehnung prüfen.",
        )
        return None, None

    anzahl_pta = punktort_layer.featureCount()
    if anzahl_pta == 0:
        QMessageBox.warning(
            tab_widget,
            "Keine Grenzpunkte",
            "Im Bereich des Verfahrensgebiets wurden keine PunktortTA-Objekte gefunden.",
        )
        return None, None

    # -----------------------------------------------------------------------
    # Schritt 2: Referenzfeld erkennen + GML-IDs extrahieren (ALKIS-6/7)
    # -----------------------------------------------------------------------
    info("Grenzpunkte", "Schritt 2/3: Referenzfeld erkennen und GML-IDs extrahieren …", duration=3)
    gehoert_zu_feld = _erkenne_gehoert_zu_feld(punktort_layer)
    if not gehoert_zu_feld:
        vorhandene_felder = ", ".join(f.name() for f in punktort_layer.fields())
        QMessageBox.warning(
            tab_widget,
            "Kein Referenzfeld",
            f"{anzahl_pta} PunktortTA-Objekte geladen, aber kein Referenzfeld auf "
            f"AX_Grenzpunkt gefunden.\n\n"
            f"Gesuchte Feldnamen: {', '.join(_GEHOERT_ZU_KANDIDATEN)}\n"
            f"Vorhandene Felder: {vorhandene_felder}",
        )
        project.addMapLayer(punktort_layer)
        return punktort_layer, None

    gml_ids = _extrahiere_gml_ids(punktort_layer, gehoert_zu_feld)

    if not gml_ids:
        QMessageBox.warning(
            tab_widget,
            "Keine Verknüpfungen",
            f"{anzahl_pta} PunktortTA-Objekte gefunden, aber Feld '{gehoert_zu_feld}' "
            "enthält keine gültigen IDs.\nBitte Attributwerte im WFS prüfen.",
        )
        # PunktortTA trotzdem zurückgeben
        project.addMapLayer(punktort_layer)
        return punktort_layer, None

    # -----------------------------------------------------------------------
    # Schritt 3: AX_Grenzpunkt per ID-Listen-Batches laden
    # -----------------------------------------------------------------------
    info(
        "Grenzpunkte",
        f"Schritt 3/3: {len(gml_ids)} Grenzpunkte in "
        f"{math.ceil(len(gml_ids) / _BATCH_GROESSE)} Batches laden …",
        duration=3,
    )
    grenzpunkt_layer = _lade_grenzpunkte_per_ids(gml_ids)

    if not grenzpunkt_layer:
        warning(
            "Grenzpunkte",
            "Grenzpunkt-Daten konnten nicht geladen werden (IDs gefunden, aber WFS lieferte kein Ergebnis).",
            duration=6,
        )
        project.addMapLayer(punktort_layer)
        return punktort_layer, None

    # Beide Layer ins Projekt
    project.addMapLayer(punktort_layer)
    project.addMapLayer(grenzpunkt_layer)

    return punktort_layer, grenzpunkt_layer


# ---------------------------------------------------------------------------
# Gemarkungsweise Abfrage
# ---------------------------------------------------------------------------

def grenzpunkte_gemarkungsweise(
    tab_widget,
    schluessel_liste: List[str],
    vg_geom=None,
    progress: Optional[Callable[[int, int], None]] = None,
) -> Optional[object]:
    """Lädt Grenzpunkte aus dem ALKIS-WFS gemarkungsweise mit per-Gemarkung-Tracking.

    Ablauf:

    1. Gemarkungsgeometrien laden, mit 200-m-Puffer um das Verfahrensgebiet
       verschneiden, BBoxen der Schnittflächen bestimmen.
    2. Pro Gemarkung: Hits-Abfrage für ``AX_PunktortTA``, dann ggf. in
       Batches à 5000 Features laden, Referenzfeld erkennen, GP-IDs
       extrahieren, PunktortTA deduplizieren.
    3. ``AX_Grenzpunkt`` für alle eindeutigen IDs per ResourceId-Batches laden.
    4. Per-Gemarkung-Statistik ausgeben.
    5. Angereicherten Join-Layer erzeugen und ins Projekt einfügen.

    Parameters
    ----------
    tab_widget :
        Eltern-Widget für Fehlerdialoge.
    schluessel_liste :
        Liste der 4-stelligen Gemarkungsnummern (z.B. ``['1234', '5678']``).
    vg_geom :
        ``QgsGeometry`` des Verfahrensgebiets (EPSG:25833). Wird für den
        200-m-Puffer zur Verschneidung der Gemarkungen verwendet.
        ``None`` = keine Einschränkung (volle Gemarkungs-BBox).
    progress :
        Reserviert (wird nicht ausgewertet; interner ``ProgressDialog`` wird
        verwendet).

    Returns
    -------
    QgsVectorLayer | None
        Angereicherter Memory-Layer mit Geometrie und allen GP-Attributen
        (Prefix ``gp_``) oder ``None`` bei Fehler.
    """
    import urllib.parse
    from sn_basis.functions import QgsField
    from sn_basis.functions.qt_wrapper import QVariant

    # Maximale Features pro einzelnem PunktortTA-GET-Aufruf
    _BBOX_BATCH = 5000

    project = QgsProject.instance()

    if not schluessel_liste:
        QMessageBox.warning(
            tab_widget,
            "Keine Gemarkungen",
            "Bitte mindestens einen Gemarkungsschlüssel angeben.",
        )
        return None

    n_gemarkungen = len(schluessel_liste)

    # Fortschrittsschritte werden nach den Hits-Abfragen präzisiert.
    # Vorab-Schätzung: 1 (Extents) + N*2 (Hits+Laden) + 1 (GP) + 1 (Join)
    prog_vorab = n_gemarkungen * 2 + 3
    prog = create_progress_dialog(
        prog_vorab,
        "Grenzpunkte laden",
        f"Lade Gemarkungs-Geometrien für {n_gemarkungen} Gemarkung(en) …",
    )
    abgeschlossen = False

    try:
        prog.set_value(0)
        if prog.is_canceled():
            return None

        # ------------------------------------------------------------------
        # Schritt 1: Gemarkungsgeometrien + Verschneidung mit VG-Puffer
        # ------------------------------------------------------------------
        if vg_geom and not vg_geom.isEmpty():
            extents = _lade_gemarkung_extents_mit_puffer(schluessel_liste, vg_geom, 200.0)
        else:
            extents = _lade_gemarkung_extents(schluessel_liste)

        if not extents:
            prog.close()
            QMessageBox.critical(
                tab_widget,
                "Fehler",
                "Gemarkungs-BBoxen konnten nicht geladen werden.\n"
                "Bitte Netzwerkverbindung und Gemarkungsschlüssel prüfen.\n"
                f"Gesuchte Schlüssel: {', '.join(schluessel_liste)}",
            )
            return None

        n_extents = len(extents)
        if n_extents < n_gemarkungen:
            warning(
                "Grenzpunkte",
                f"Nur {n_extents}/{n_gemarkungen} Gemarkungen im WFS gefunden. "
                f"Nicht gefunden: {set(schluessel_liste) - set(extents.keys())}",
                duration=8,
            )

        # ------------------------------------------------------------------
        # Schritt 2a: Hits-Abfragen für alle BBoxen
        # ------------------------------------------------------------------
        prog.set_label(f"Anzahl PunktortTA-Objekte abfragen ({n_extents} Gemarkungen) …")
        prog.set_value(1)
        if prog.is_canceled():
            return None

        # {gemarkungsnummer: anzahl_features} – -1 = unbekannt
        hits_pro_gemarkung: Dict[str, int] = {}
        n_pta_batches_gesamt = 0

        for gemarkungsnummer, bbox in extents.items():
            bbox_param = _punktort_ta_bbox_param(bbox)
            n_hits = _wfs_hits_count(ALKIS_SF_URL, PUNKTORT_TA_TYPENAME, {'BBOX': bbox_param})
            hits_pro_gemarkung[gemarkungsnummer] = n_hits
            batches = max(1, math.ceil(n_hits / _BBOX_BATCH)) if n_hits > 0 else 1
            n_pta_batches_gesamt += batches
            info(
                "Grenzpunkte (Diagnose)",
                f"Gemarkung {gemarkungsnummer}: ~{n_hits} PunktortTA-Objekte "
                f"→ {batches} Batch(es).",
                duration=2,
            )

        # Gesamtschritte jetzt präzise:
        # 1 (Extents/Hits) + n_pta_batches_gesamt (PunktortTA) + n_gp_batches (GP) + 1 (Join)
        # GP-Batches kennen wir noch nicht – vorläufig 0, wird unten korrigiert.
        prog_gesamt_pta = 1 + n_pta_batches_gesamt
        prog.set_total(prog_gesamt_pta + 10)  # Platzhalter für GP + Join

        # ------------------------------------------------------------------
        # Schritt 2b: PunktortTA pro Gemarkung + Batching
        # ------------------------------------------------------------------
        # Statt rohe QgsFeature-Objekte (positionsbasierte Attribute!) zu sammeln,
        # speichern wir (geom, attrs-dict) – vollständig namensbasiert.
        # Das verhindert Feldwert-Verschiebungen, wenn verschiedene Batches
        # die Felder in unterschiedlicher Reihenfolge liefern.
        alle_pta_records: List[Dict] = []  # {'geom': QgsGeometry, 'attrs': {name: val}}
        gesehene_pta_ids: set[str] = set()
        pta_feld_namen: List[str] = []
        gehoert_zu_feld: Optional[str] = None

        per_gemarkung_gp_ids: Dict[str, set] = {}
        prog_schritt = 1  # nach Hits-Phase

        for g_idx, (gemarkungsnummer, bbox) in enumerate(extents.items(), 1):
            n_hits = hits_pro_gemarkung.get(gemarkungsnummer, -1)
            n_batches_diese_gemarkung = max(1, math.ceil(n_hits / _BBOX_BATCH)) if n_hits > 0 else 1

            if prog.is_canceled():
                return None

            gp_ids_dieser_gemarkung: set[str] = set()
            neu_gesamt = 0
            duplikat_gesamt = 0
            bbox_param = _punktort_ta_bbox_param(bbox)

            for b_idx in range(n_batches_diese_gemarkung):
                start_index = b_idx * _BBOX_BATCH
                prog.set_label(
                    f"Gemarkung {gemarkungsnummer} ({g_idx}/{n_extents}), "
                    f"Batch {b_idx + 1}/{n_batches_diese_gemarkung} "
                    f"(~{n_hits if n_hits >= 0 else '?'} Objekte) …"
                )
                prog.set_value(prog_schritt)
                prog_schritt += 1
                if prog.is_canceled():
                    return None

                extra: dict = {'BBOX': bbox_param}
                if n_batches_diese_gemarkung > 1:
                    extra['COUNT'] = str(_BBOX_BATCH)
                    extra['STARTINDEX'] = str(start_index)

                info(
                    "Grenzpunkte (Diagnose)",
                    f"Gemarkung {gemarkungsnummer} ({g_idx}/{n_extents}), "
                    f"Batch {b_idx + 1}/{n_batches_diese_gemarkung}: "
                    f"STARTINDEX={start_index}, COUNT={_BBOX_BATCH} …",
                    duration=2,
                )

                _gml_bytes = _wfs_get_gml(ALKIS_SF_URL, PUNKTORT_TA_TYPENAME, extra)
                if not _gml_bytes:
                    warning(
                        "Grenzpunkte (Diagnose)",
                        f"Gemarkung {gemarkungsnummer} Batch {b_idx + 1}: "
                        f"GET-Anfrage fehlgeschlagen – Batch wird übersprungen.",
                        duration=6,
                    )
                    continue

                tmp = _parse_alkis_gml_memory_layer(
                    _gml_bytes, "AX_PunktortTA", "_TmpPTA", geom_prop="position"
                )
                if not tmp:
                    continue

                cur_feld_namen = [f.name() for f in tmp.fields()]

                # Globale Feldliste inkrementell um neue Felder erweitern
                # (Union über alle Batches/Gemarkungen).
                for n in cur_feld_namen:
                    if n not in pta_feld_namen:
                        pta_feld_namen.append(n)

                # Referenzfeld beim ersten Batch erkennen, der es enthält
                if gehoert_zu_feld is None:
                    gehoert_zu_feld = _erkenne_gehoert_zu_feld(tmp)
                    if gehoert_zu_feld:
                        info(
                            "Grenzpunkte",
                            f"Referenzfeld erkannt: '{gehoert_zu_feld}'. "
                            f"PunktortTA-Felder bisher: {', '.join(pta_feld_namen)}",
                            duration=5,
                        )
                    elif not pta_feld_namen:
                        warning(
                            "Grenzpunkte (Diagnose)",
                            "Kein Referenzfeld in AX_PunktortTA erkannt. "
                            "Vorhandene Felder: " + ", ".join(cur_feld_namen),
                            duration=15,
                        )

                neu = 0
                duplikat = 0
                for feat in tmp.getFeatures():
                    pta_id = _extrahiere_feature_gml_id(feat, cur_feld_namen)
                    if pta_id and pta_id in gesehene_pta_ids:
                        duplikat += 1
                        continue
                    if pta_id:
                        gesehene_pta_ids.add(pta_id)

                    # Attribute namensbasiert als Dict speichern – verhindert
                    # Positions-Verschiebung beim späteren Aufbau von punktort_layer.
                    _attrs: Dict[str, str] = {}
                    for n in cur_feld_namen:
                        if n:
                            try:
                                _attrs[n] = feat[n]
                            except Exception:
                                pass
                    alle_pta_records.append({'geom': feat.geometry(), 'attrs': _attrs})
                    neu += 1

                    if gehoert_zu_feld:
                        try:
                            raw = feat[gehoert_zu_feld]
                            raw = str(raw).strip() if raw else ""
                            gp_id = raw.rsplit(":", 1)[-1] if ":" in raw else raw
                            if gp_id:
                                gp_ids_dieser_gemarkung.add(gp_id)
                        except Exception:
                            pass

                neu_gesamt += neu
                duplikat_gesamt += duplikat

            per_gemarkung_gp_ids[gemarkungsnummer] = gp_ids_dieser_gemarkung
            info(
                "Grenzpunkte",
                f"Gemarkung {gemarkungsnummer}: {neu_gesamt} PunktortTA-Objekte "
                f"({duplikat_gesamt} Duplikate übersprungen), "
                f"{len(gp_ids_dieser_gemarkung)} GP-Referenzen.",
                duration=4,
            )

        # Abbruch wenn keine Features
        if not alle_pta_records:
            prog.close()
            QMessageBox.warning(
                tab_widget,
                "Keine Grenzpunkte",
                "Für die angegebenen Gemarkungen wurden keine PunktortTA-Objekte gefunden.",
            )
            return None

        if not gehoert_zu_feld:
            prog.close()
            QMessageBox.warning(
                tab_widget,
                "Kein Referenzfeld",
                f"{len(alle_pta_records)} PunktortTA-Objekte geladen, "
                f"aber kein Referenzfeld auf AX_Grenzpunkt gefunden.\n\n"
                f"Gesuchte Feldnamen: {', '.join(_GEHOERT_ZU_KANDIDATEN)}\n"
                f"Vorhandene Felder: {', '.join(pta_feld_namen)}",
            )
            return None

        punktort_layer = QgsVectorLayer("Point?crs=EPSG:25833", PUNKTORT_LAYER_NAME, "memory")
        pr_pta = punktort_layer.dataProvider()
        pr_pta.addAttributes([QgsField(n, QVariant.String) for n in pta_feld_namen if n])
        punktort_layer.updateFields()

        # Features mit vollständig namensbasierter Attributzuweisung aufbauen –
        # vermeidet Positions-Verschiebungen, die entstehen würden, wenn
        # QgsFeatures aus tmp-Layern mit unterschiedlicher Feldreihenfolge
        # positionsbasiert in punktort_layer übertragen werden.
        neue_pta_features = []
        from sn_basis.functions import QgsFeature as _QgsFeature
        for record in alle_pta_records:
            new_f = _QgsFeature(punktort_layer.fields())
            new_f.setGeometry(record['geom'])
            for n, val in record['attrs'].items():
                if n in pta_feld_namen:
                    try:
                        new_f[n] = val
                    except Exception:
                        pass
            neue_pta_features.append(new_f)
        pr_pta.addFeatures(neue_pta_features)
        punktort_layer.updateExtents()

        alle_gp_ids: set[str] = set()
        for ids in per_gemarkung_gp_ids.values():
            alle_gp_ids.update(ids)

        n_pta_gesamt = len(alle_pta_records)
        n_gp_ids = len(alle_gp_ids)
        info(
            "Grenzpunkte",
            f"Gesamt: {n_pta_gesamt} PunktortTA-Objekte, {n_gp_ids} eindeutige GP-Referenzen "
            f"aus {n_extents} Gemarkung(en).",
            duration=5,
        )

        if not alle_gp_ids:
            prog.close()
            QMessageBox.warning(
                tab_widget,
                "Keine Verknüpfungen",
                f"{n_pta_gesamt} PunktortTA-Objekte geladen, aber Referenzfeld "
                f"'{gehoert_zu_feld}' enthält keine gültigen IDs.",
            )
            return None

        # ------------------------------------------------------------------
        # Schritt 3: AX_Grenzpunkt per ResourceId-Batches
        # ------------------------------------------------------------------
        anzahl_gp_batches = math.ceil(n_gp_ids / _BATCH_GROESSE)
        # Gesamtschritte jetzt vollständig bekannt
        prog_gesamt = prog_schritt + anzahl_gp_batches + 1
        prog.set_total(prog_gesamt)
        prog.set_label(
            f"Grenzpunkte laden: {n_gp_ids} IDs in {anzahl_gp_batches} Batch(es) …"
        )
        prog.set_value(prog_schritt)
        if prog.is_canceled():
            return None

        def _gp_progress(batch_idx: int) -> None:
            prog.set_label(
                f"AX_Grenzpunkt Batch {batch_idx + 1}/{anzahl_gp_batches} …"
            )
            prog.set_value(prog_schritt + batch_idx + 1)

        gp_dict, gp_feld_namen = _lade_grenzpunkte_als_dict(
            sorted(alle_gp_ids), progress_callback=_gp_progress
        )

        if not gp_dict:
            prog.close()
            warning(
                "Grenzpunkte",
                "Grenzpunkt-Daten konnten nicht geladen werden "
                "(IDs vorhanden, WFS lieferte kein Ergebnis).",
                duration=6,
            )
            return None

        info(
            "Grenzpunkte",
            f"{len(gp_dict)}/{n_gp_ids} Grenzpunkte vom WFS empfangen "
            f"({n_gp_ids - len(gp_dict)} IDs nicht aufgelöst).",
            duration=4,
        )

        for gemarkungsnummer, gp_ids in per_gemarkung_gp_ids.items():
            if not gp_ids:
                continue
            n_gefunden = sum(1 for gid in gp_ids if gid in gp_dict)
            info(
                "Grenzpunkte",
                f"Gemarkung {gemarkungsnummer}: {n_gefunden}/{len(gp_ids)} "
                f"Grenzpunkte geladen.",
                duration=6,
            )

        # ------------------------------------------------------------------
        # Schritt 4: Angereicherten Join-Layer erstellen
        # ------------------------------------------------------------------
        prog.set_label("Angereicherten Layer erstellen …")
        prog.set_value(prog_gesamt - 1)
        if prog.is_canceled():
            return None

        joined_layer = _erstelle_joined_layer(
            punktort_layer, gp_dict, gp_feld_namen, gehoert_zu_feld
        )

        if not joined_layer or not joined_layer.isValid():
            prog.close()
            QMessageBox.critical(
                tab_widget,
                "Fehler",
                "Angereicherter Grenzpunkt-Layer konnte nicht erstellt werden.",
            )
            return None

        project.addMapLayer(joined_layer)

        prog.set_value(prog_gesamt)
        prog.close()
        abgeschlossen = True

        info(
            "Grenzpunkte",
            f"Fertig: {joined_layer.featureCount()} Grenzpunkte im Layer "
            f"({len(gp_dict)} mit GP-Attributen, "
            f"{joined_layer.featureCount() - len(gp_dict)} ohne GP-Match).",
            duration=10,
        )

        return joined_layer

    finally:
        if not abgeschlossen:
            prog.close()