Geoid Models and Vertical Datum Automation

import numpy as np
from pyproj import CRS, Transformer


def to_orthometric(lon: np.ndarray, lat: np.ndarray, ellipsoidal_h: np.ndarray,
                   *, horizontal_epsg: int, vertical_epsg: int) -> dict:
    """Convert ellipsoidal to orthometric height using a geoid model.

    A compound CRS lets PROJ apply the geoid grid at each point rather than a
    constant, which is what removes the tilt a single subtraction introduces
    across a large site.
    """
    source = CRS.from_epsg(4979)                               # WGS84 3D, ellipsoidal
    target = CRS.from_user_input(f"EPSG:{horizontal_epsg}+{vertical_epsg}")
    transformer = Transformer.from_crs(source, target, always_xy=True)

    x, y, z = transformer.transform(lon, lat, ellipsoidal_h)
    undulation = np.asarray(ellipsoidal_h) - np.asarray(z)

    return {"x": x, "y": y, "orthometric_h": z,
            "undulation_m": undulation,
            "undulation_range_m": float(np.ptp(undulation)),
            "target_crs": target.to_string()}

from pyproj import CRS
from pyproj.transformer import TransformerGroup


def check_transformation_quality(horizontal_epsg: int, vertical_epsg: int) -> dict:
    """Is the best transformation available, and is it grid-based?

    TransformerGroup lists every candidate transformation and flags those
    whose grids are missing. A group whose best option is unavailable will
    silently use a worse one.
    """
    source = CRS.from_epsg(4979)
    target = CRS.from_user_input(f"EPSG:{horizontal_epsg}+{vertical_epsg}")
    group = TransformerGroup(source, target, always_xy=True)

    missing = [{"grid": g.short_name, "url": g.url, "available": g.available}
               for op in group.unavailable_operations for g in op.grids]

    best = group.transformers[0].description if group.transformers else None
    return {"best_available": best,
            "unavailable_operations": len(group.unavailable_operations),
            "missing_grids": [m for m in missing if not m["available"]],
            "best_is_available": group.best_available,
            "note": ("the best transformation is available" if group.best_available else
                     "the best transformation needs grids that are not installed — "
                     "results will use a lower-accuracy fallback")}

import rasterio
from pyproj import CRS


def write_with_compound_crs(src_path: str, dst_path: str,
                            horizontal_epsg: int, vertical_epsg: int) -> dict:
    """Re-write a raster with a compound CRS declaring both datums."""
    compound = CRS.from_user_input(f"EPSG:{horizontal_epsg}+{vertical_epsg}")
    with rasterio.open(src_path) as src:
        profile = src.profile
        profile.update(crs=compound)
        data = src.read()
        tags = src.tags()

    with rasterio.open(dst_path, "w", **profile) as dst:
        dst.write(data)
        dst.update_tags(**tags, VERTICAL_DATUM=f"EPSG:{vertical_epsg}",
                        HEIGHT_TYPE="orthometric")
    return {"crs": compound.to_string(), "path": dst_path}


def assert_vertical_datum(path: str, expected_vertical_epsg: int) -> None:
    """Refuse a product whose vertical datum is absent or unexpected."""
    with rasterio.open(path) as src:
        crs = src.crs
        tags = src.tags()
    if crs is None:
        raise ValueError(f"{path} has no CRS at all")
    parsed = CRS.from_user_input(crs)
    if not parsed.is_compound:
        raise ValueError(f"{path} declares a horizontal CRS only — its heights are "
                         "undefined")
    vertical = [sub for sub in parsed.sub_crs_list if sub.is_vertical]
    if not vertical or vertical[0].to_epsg() != expected_vertical_epsg:
        raise ValueError(f"{path} declares vertical {vertical[0].to_epsg() if vertical else None}, "
                         f"expected {expected_vertical_epsg}")

import numpy as np


def reconcile_heights(gcp_heights: np.ndarray, gcp_datum: str,
                      camera_heights: np.ndarray, camera_datum: str,
                      undulation_m: np.ndarray) -> dict:
    """Bring control and camera heights into one vertical frame.

    The undulation is per point, not a constant, because the control points
    and the camera positions are at different places and the model differs
    between them.
    """
    if gcp_datum == camera_datum:
        return {"converted": False, "note": "both are already in the same datum"}

    if gcp_datum == "orthometric" and camera_datum == "ellipsoidal":
        cameras = camera_heights - undulation_m
        return {"converted": True, "camera_heights": cameras,
                "note": "camera heights converted to orthometric"}
    if gcp_datum == "ellipsoidal" and camera_datum == "orthometric":
        gcps = gcp_heights - undulation_m
        return {"converted": True, "gcp_heights": gcps,
                "note": "control heights converted to orthometric"}
    raise ValueError(f"unhandled datum pair: {gcp_datum} and {camera_datum}")

from pyproj import CRS, Transformer


def transform_with_epoch(lon, lat, h, *, source_epsg: int, source_epoch: float,
                         target_epsg: int, target_epoch: float) -> dict:
    """Transform between dynamic datum realisations at stated epochs.

    Both epochs must be supplied. A transformation between dynamic frames
    without them uses a default that is unlikely to match either the survey
    or the deliverable, and the resulting shift is a systematic error nobody
    can attribute later.
    """
    source = CRS.from_epsg(source_epsg)
    target = CRS.from_epsg(target_epsg)
    transformer = Transformer.from_crs(source, target, always_xy=True)
    x, y, z = transformer.transform(lon, lat, h, source_epoch, target_epoch)
    return {"x": x, "y": y, "z": z,
            "source_epoch": source_epoch, "target_epoch": target_epoch,
            "note": "epochs recorded; a transformation without them uses a default"}


def epoch_warning(survey_epoch: float, control_epoch: float,
                  *, max_gap_years: float = 2.0) -> str | None:
    """Flag a survey whose control was measured at a materially different epoch."""
    gap = abs(survey_epoch - control_epoch)
    if gap > max_gap_years:
        return (f"control was measured {gap:.1f} years from the survey epoch — "
                "in a dynamic frame this is a real coordinate difference, "
                "not a measurement error")
    return None

def vertical_statement(vertical_epsg: int, geoid_model: str,
                       epoch: float | None, benchmark_check: dict | None) -> str:
    """The sentence that belongs in every survey report."""
    from pyproj import CRS
    name = CRS.from_epsg(vertical_epsg).name
    parts = [f"Heights are orthometric, referenced to {name} (EPSG:{vertical_epsg}), "
             f"converted from ellipsoidal heights using the {geoid_model} geoid model."]
    if epoch is not None:
        parts.append(f"Horizontal coordinates are at epoch {epoch:.1f}.")
    if benchmark_check and benchmark_check.get("benchmarks"):
        parts.append(f"Verified against {benchmark_check['benchmarks']} published "
                     f"benchmarks with a mean difference of "
                     f"{benchmark_check['bias_m']:+.3f} m.")
    return " ".join(parts)

import numpy as np


def verify_against_benchmarks(measured_h: np.ndarray, published_h: np.ndarray,
                              *, tolerance_m: float = 0.05) -> dict:
    """Compare converted heights against published benchmark values.

    The bias is the number that matters: a systematic offset means the wrong
    datum or a missing grid, while scatter means ordinary measurement error.
    """
    d = np.asarray(measured_h) - np.asarray(published_h)
    d = d[np.isfinite(d)]
    if d.size < 2:
        return {"note": "need at least two benchmarks"}

    bias = float(np.mean(d))
    scatter = float(np.std(d, ddof=1))
    return {"benchmarks": int(d.size), "bias_m": bias, "scatter_m": scatter,
            "within_tolerance": abs(bias) < tolerance_m,
            "diagnosis": ("heights agree with the published datum" if abs(bias) < tolerance_m
                          else f"systematic {bias:+.3f} m — check the vertical datum, "
                               "the geoid grid and the antenna height")}