Converting GML Responses to GeoDataFrames

pip install "geopandas>=1.0" "pyogrio>=0.9" "pyproj>=3.6" "lxml>=5.0"

from __future__ import annotations

import logging
import re
import tempfile
from pathlib import Path

import geopandas as gpd
from pyproj import CRS

logger = logging.getLogger(__name__)

SRS_PATTERN = re.compile(rb'srsName="([^"]+)"')


def _srs_from_payload(payload: bytes) -> str | None:
    """Read the first srsName declaration out of the raw GML."""
    match = SRS_PATTERN.search(payload)
    return match.group(1).decode() if match else None


def _needs_axis_swap(srs: str) -> bool:
    """URN-form EPSG references follow authority axis order (lat, lon for 4326)."""
    if not srs:
        return False
    if "CRS84" in srs.upper():
        return False
    try:
        crs = CRS.from_user_input(srs)
    except Exception:  # noqa: BLE001 - unknown authority strings are common in the wild
        logger.warning("unrecognised srsName %r; assuming lon/lat order", srs)
        return False
    order = [axis.abbrev.lower() for axis in crs.axis_info][:2]
    return order == ["lat", "lon"]


def gml_to_geodataframe(
    payload: bytes,
    target_crs: int = 4326,
    flatten_sep: str = "_",
) -> gpd.GeoDataFrame:
    """Parse a WFS GML payload into a flat, CRS-correct GeoDataFrame.

    The payload is written to a temporary file because GDAL's GML driver
    resolves the application schema from the document and caches it beside it.
    """
    srs = _srs_from_payload(payload)
    swap = _needs_axis_swap(srs or "")

    with tempfile.TemporaryDirectory() as tmpdir:
        gml_path = Path(tmpdir) / "response.gml"
        gml_path.write_bytes(payload)

        # SWAP_COORDINATES=YES makes GDAL apply the authority order itself,
        # which is more reliable than transforming the geometries afterwards.
        gdf = gpd.read_file(
            gml_path,
            driver="GML",
            SWAP_COORDINATES="YES" if swap else "NO",
        )

    if gdf.empty:
        logger.warning("GML payload parsed to zero features (srsName=%s)", srs)
        return gpd.GeoDataFrame(geometry=[], crs=target_crs)

    if gdf.crs is None and srs:
        gdf = gdf.set_crs(CRS.from_user_input(srs), allow_override=True)
    if gdf.crs is None:
        raise ValueError("GML response declared no CRS and none could be inferred")

    # Flatten any nested property columns the driver surfaced as dicts.
    for column in list(gdf.columns):
        if column == gdf.geometry.name:
            continue
        sample = gdf[column].dropna().head(1)
        if len(sample) and isinstance(sample.iloc[0], dict):
            expanded = gdf[column].apply(lambda value: value or {}).apply(pd.Series)
            expanded.columns = [f"{column}{flatten_sep}{sub}" for sub in expanded.columns]
            gdf = gdf.drop(columns=[column]).join(expanded)

    result = gdf.to_crs(target_crs)
    logger.info("parsed %d features from GML (srsName=%s, swap=%s)", len(result), srs, swap)
    return result

import pandas as pd  # imported alongside the recipe for the flattening step