Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension


Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0

## [Unreleased]

### Added
- `ivert validate -V` accepts `mtl` (mean tide level) and `dtl` (diurnal tide level), passed to transformez as `vdatum:mtl` and `vdatum:dtl` (#235). transformez 1.0.1 added both references.

### Changed
- IVERT now needs `transformez>=1.0.1` (#235), for the `vdatum:mtl` and `vdatum:dtl` references. transformez 1.0.1 also extracts the NOAA VDatum model grids (such as xGEOID20B, about 200 MB) into its cache, where 1.0.0 extracted them into whatever folder `ivert` was run from.

### Fixed
- Vertical datum shifts are read from the right place in their shift grid (#231). transformez stores each grid value at its pixel centre, but IVERT read the grid as if its outer values sat on its edges, stretching it by half a pixel at each side. Heights moved by up to about 2 mm on IVERT's 3-arcsecond grids; points near the grid's edges were off the most.
- File sizes shown by `ivert database size` and `ivert database delete` no longer repeat a digit when they round up (#232). A 2047-byte file was listed as "11.0 KB" instead of "2.0 KB", and one just under 1 MB as "10231.0 KB" instead of "1.0 MB".
Expand Down
2 changes: 1 addition & 1 deletion docs/source/user_guide/validate.md
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@ Which CRS wins: the horizontal CRS comes from `-p`, else the DEM file. The verti

ICESat-2 photons are stored either in WGS84 'ellipsoid' (EPSG:4979) or EGM2008 'geoid' (EPSG:3855) elevations, depending on the setting used (see 'ivert options' for more information). IVERT automatically transforms them to the DEM's vertical datum before comparing. If the DEM's CRS already encodes the vertical datum, `-V` is not needed.

The vertical datum shift is done by [transformez](https://github.com/continuous-dems/transformez). Tidal datums (`mllw`, `mlw`, `msl`, `mhw`, `mhhw`) are passed to it as its own NOAA VDatum references (`vdatum:mllw`, …) rather than as EPSG codes, and transformez decides how each surface is applied. IVERT checks that transformez can transform the datum you give before it starts, so an unsupported one stops the run straight away.
The vertical datum shift is done by [transformez](https://github.com/continuous-dems/transformez). Tidal datums (`mllw`, `mlw`, `mtl`, `dtl`, `msl`, `mhw`, `mhhw`) are passed to it as its own NOAA VDatum references (`vdatum:mllw`, …) rather than as EPSG codes, and transformez decides how each surface is applied. IVERT checks that transformez can transform the datum you give before it starts, so an unsupported one stops the run straight away.

### Photon filtering

Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ dependencies = [
"xarray>=2023.10.0",
"fetchez>=0.9.1",
"globato>=0.3.9",
"transformez>=1.0.0",
"transformez>=1.0.1",
]

[dependency-groups]
Expand Down
35 changes: 9 additions & 26 deletions src/ivert/vdatum_lookup.py
Original file line number Diff line number Diff line change
Expand Up @@ -43,13 +43,11 @@
# Tidal datums go to transformez as its own IDs rather than EPSG codes: it builds
# them from NOAA VDatum grids, and it is the authority on whether each surface is a
# height or a depth.
#
# transformez's VDatum support also lists the mean tide level (mtl) and diurnal
# tide level (dtl) surfaces, but it has no vdatum:mtl or vdatum:dtl reference to
# reach them through yet. Once it does, add "mtl" and "dtl" here.
_TIDAL_NAMES: dict[str, str] = {
"mllw": "vdatum:mllw",
"mlw": "vdatum:mlw",
"mtl": "vdatum:mtl",
"dtl": "vdatum:dtl",
"msl": "vdatum:msl",
"mhw": "vdatum:mhw",
"mhhw": "vdatum:mhhw",
Expand Down Expand Up @@ -96,18 +94,6 @@ def resolve_vdatum(name: str | int | None) -> str | None:
return None if epsg is None else f"EPSG:{epsg}"


def _reference_errors() -> tuple[type[Exception], ...]:
"""Return the exceptions transformez raises for a reference it cannot use."""
# transformez derives these from the built-in ReferenceError, not ValueError.
from transformez.reference.parser import ( # noqa: PLC0415 - slow import
InvalidReferenceError,
ReferenceInputError,
UnsupportedReferenceError,
)

return InvalidReferenceError, ReferenceInputError, UnsupportedReferenceError


def check_vdatum(reference: str) -> None:
"""Check that transformez can transform heights to or from a vertical reference.

Expand All @@ -125,14 +111,12 @@ def check_vdatum(reference: str) -> None:
resolve_reference,
)

try:
parsed = parse_reference(reference)
if parsed.vertical is None:
msg = f"{reference!r} has no vertical component."
raise ValueError(msg)
resolve_reference(parsed)
except _reference_errors() as exc:
raise ValueError(str(exc)) from exc
# transformez raises a ValueError for a reference it cannot parse or transform.
parsed = parse_reference(reference)
if parsed.vertical is None:
msg = f"{reference!r} has no vertical component."
raise ValueError(msg)
resolve_reference(parsed)


def describe_vdatum(reference: str) -> str:
Expand All @@ -141,10 +125,9 @@ def describe_vdatum(reference: str) -> str:
parse_reference,
)

errors: tuple[type[Exception], ...] = (*_reference_errors(), ValueError)
try:
vertical = parse_reference(reference).vertical
except errors:
except ValueError:
return ""
return vertical.name if vertical is not None else ""

Expand Down
124 changes: 124 additions & 0 deletions tests/test_vdatum_lookup.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
"""Tests for translating vertical datum names into references transformez can use."""

import pytest

from ivert import vdatum_lookup

# --- resolve_vdatum -----------------------------------------------------------------


@pytest.mark.parametrize(
("name", "reference"),
[
("navd88", "EPSG:5703"),
("NAVD88", "EPSG:5703"),
(" egm2008 ", "EPSG:3855"),
("EGM96", "EPSG:5773"),
("prvd02", "EPSG:6641"),
("cgvd2013", "EPSG:6647"),
("navd88 height (ft)", "EPSG:8228"),
],
)
def test_geoid_datum_names_resolve_to_epsg_codes(name, reference):
"""Names are case-insensitive and ignore surrounding spaces."""
assert vdatum_lookup.resolve_vdatum(name) == reference


@pytest.mark.parametrize("name", ["ellipsoid", "wgs84", "WGS84"])
def test_ellipsoid_names_resolve_to_wgs84_3d(name):
"""Ellipsoidal heights mean WGS84, the photon database's frame (EPSG:4979)."""
assert vdatum_lookup.resolve_vdatum(name) == "EPSG:4979"


@pytest.mark.parametrize(
("name", "reference"),
[
("mllw", "vdatum:mllw"),
("MLW", "vdatum:mlw"),
("mtl", "vdatum:mtl"),
("dtl", "vdatum:dtl"),
("msl", "vdatum:msl"),
("mhw", "vdatum:mhw"),
("mhhw", "vdatum:mhhw"),
],
)
def test_tidal_datum_names_resolve_to_vdatum_references(name, reference):
"""Tidal datums go to transformez as its own VDatum references, not EPSG codes."""
assert vdatum_lookup.resolve_vdatum(name) == reference


@pytest.mark.parametrize(
("code", "reference"),
[
(5703, "EPSG:5703"),
("5703", "EPSG:5703"),
("EPSG:5703", "EPSG:5703"),
("vdatum:mllw", "vdatum:mllw"),
("epsg:3855", "epsg:3855"),
],
)
def test_codes_and_qualified_references_pass_through(code, reference):
"""A bare code gains its 'EPSG:'; a qualified reference is left as given."""
assert vdatum_lookup.resolve_vdatum(code) == reference


@pytest.mark.parametrize("name", [None, "", "not a datum", "navd89"])
def test_unknown_names_resolve_to_none(name):
"""A name the table doesn't know is None, not a guess."""
assert vdatum_lookup.resolve_vdatum(name) is None


# --- check_vdatum -------------------------------------------------------------------


@pytest.mark.parametrize("reference", sorted(vdatum_lookup.list_vdatums()))
def test_every_listed_datum_is_one_transformez_can_use(reference):
"""Each reference a name resolves to passes the check, so -V accepts it."""
vdatum_lookup.check_vdatum(reference)


@pytest.mark.parametrize(
("reference", "match"),
[
("EPSG:4326", "no vertical component"),
("EPSG:5868", "no operation binding"),
("vdatum:xyz", "Unknown custom reference"),
("not a datum", "Unsupported coordinate reference"),
],
)
def test_references_transformez_cannot_use_are_a_value_error(reference, match):
"""Horizontal-only, untransformable and unknown references are all ValueErrors."""
with pytest.raises(ValueError, match=match):
vdatum_lookup.check_vdatum(reference)


# --- describe_vdatum and list_vdatums -----------------------------------------------


@pytest.mark.parametrize(
("reference", "description"),
[
("EPSG:5703", "NAVD88 height"),
("vdatum:mtl", "NOAA VDatum Mean Tide Level"),
],
)
def test_references_are_described_by_transformez(reference, description):
"""The description is transformez's name for the vertical reference."""
assert vdatum_lookup.describe_vdatum(reference) == description


@pytest.mark.parametrize("reference", ["EPSG:4326", "vdatum:xyz", "not a datum"])
def test_a_reference_with_no_vertical_name_is_described_as_empty(reference):
"""No vertical part, or one transformez doesn't know, gives an empty string."""
assert vdatum_lookup.describe_vdatum(reference) == ""


def test_the_list_maps_each_reference_to_the_names_for_it():
"""Every listed name resolves to the reference it's listed under."""
listed = vdatum_lookup.list_vdatums()
assert "navd88" in listed["EPSG:5703"]
assert listed["vdatum:mtl"] == ["mtl"]
for reference, names in listed.items():
assert names
for name in names:
assert vdatum_lookup.resolve_vdatum(name) == reference
Loading