From 47c55e4231005dfbd6acaa45dca11df4135ccdf8 Mon Sep 17 00:00:00 2001 From: Romain Hugonnet Date: Fri, 28 Apr 2023 12:24:05 -0700 Subject: [PATCH] Rework `DEM.vcrs` to support any 3D transformation (and fix with PROJ update) (#350) * Rework vref into vcrs, clean behaviour with ccrs and fix with pyproj update * Incremental commit, waiting for pyproj to_2d solution * Incremental commit on tests * Move vcrs to its own module, write more tests * Linting * Add future import of annotations for python 3.8 * Add pyproj version condition to tests * Linting * Skip vertical transform is dest compound CRS is the same as source, and add tests for it * Remove proj-data dependency * Add automatic download of PROJ grids with grid user input, and through TransformationGroup during transform * Linting * Fix test return for pyproj < 3.5.1 * Add error to build_ccrs and tests * Make axis length checks robust to older pyproj versions * Linting * Make URL error more user-friendly, remove Windows test skipping * Add a couple useful functions and VCRS setting based on CRS during DEM instantiation * Linting * Remove vcrs_equals and override from_array * Linting * Comment other test based on speed * Linting * Incremental commit on doc * Linting * Finalize documentation page * Fix tests and refactor vcrs_from_user_input to not be public * Linting * Try to fix PROJ path on readthedocs * Is it a problem with the datadir then? * Comment the pre_job * Retry with unset PROJ * Try to unset during post-build instead * Revert to default proj data dir in transformer group * Fix re-initiation of trans_group object after grid downloading in transform_zz * Linting * Accout for eriks comments * Linting --- dev-environment.yml | 3 +- doc/source/api.md | 11 +- doc/source/vertical_ref.md | 237 ++++++++++++++++++++++- environment.yml | 3 +- tests/test_dem.py | 331 ++++++++++++++++++++------------ tests/test_spatialstats.py | 21 +-- tests/test_vcrs.py | 201 ++++++++++++++++++++ xdem/ddem.py | 14 +- xdem/dem.py | 376 ++++++++++++++++++++----------------- xdem/vcrs.py | 333 ++++++++++++++++++++++++++++++++ 10 files changed, 1213 insertions(+), 317 deletions(-) create mode 100644 tests/test_vcrs.py create mode 100644 xdem/vcrs.py diff --git a/dev-environment.yml b/dev-environment.yml index 379f7e24..a93563ab 100644 --- a/dev-environment.yml +++ b/dev-environment.yml @@ -12,12 +12,11 @@ dependencies: - matplotlib - opencv - openh264 - - pyproj + - pyproj>=3.4 - rasterio>=1.3 - scipy - tqdm - scikit-image - - proj-data - scikit-gstat>=0.6.8 - pytransform3d - geoutils==0.0.11 diff --git a/doc/source/api.md b/doc/source/api.md index 7b695692..d8597160 100644 --- a/doc/source/api.md +++ b/doc/source/api.md @@ -49,7 +49,7 @@ A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Rast .. autosummary:: :toctree: gen_modules/ - DEM.vref + DEM.vcrs ``` ### Derived attributes @@ -61,6 +61,15 @@ A {class}`~xdem.DEM` inherits four unique attributes from {class}`~geoutils.Rast DEM.ccrs ``` +### Vertical referencing + +```{eval-rst} +.. autosummary:: + :toctree: gen_modules/ + + DEM.set_vcrs + DEM.to_vcrs +``` ## dDEM diff --git a/doc/source/vertical_ref.md b/doc/source/vertical_ref.md index 71d12e7e..96f965f8 100644 --- a/doc/source/vertical_ref.md +++ b/doc/source/vertical_ref.md @@ -1,5 +1,240 @@ +--- +file_format: mystnb +jupytext: + formats: md:myst + text_representation: + extension: .md + format_name: myst +kernelspec: + display_name: xdem-env + language: python + name: xdem +--- (vertical-ref)= # Vertical referencing -TODO: In construction +xDEM supports the use of **vertical coordinate reference systems (vertical CRSs)** and vertical transformations for DEMs +by conveniently wrapping PROJ pipelines through [Pyproj](https://pyproj4.github.io/pyproj/stable/) in the {class}`~xdem.DEM` class. + +```{important} +**A {class}`~xdem.DEM` already possesses a {class}`~xdem.DEM.crs` attribute that defines its 2- or 3D CRS**, inherited from +{class}`~geoutils.Raster`. Unfortunately, most DEM products do not yet come with a 3D CRS in their raster metadata, and +vertical CRSs often have to be set by the user. See {ref}`vref-setting` below. +``` + +## What is a vertical CRS? + +A vertical CRS is a **1D, often gravity-related, coordinate reference system of surface elevation** (or height), used to expand a [2D CRS](https://en.wikipedia.org/wiki/Spatial_reference_system) to a 3D CRS. + +There are two types of 3D CRSs, related to two types of definition of the vertical axis: +- **Ellipsoidal heights** CRSs, that are simply 2D CRS promoted to 3D by explicitly adding an elevation axis to the ellipsoid used by the 2D CRS, +- **Geoid heights** CRSs, that are a compound of a 2D CRS and a vertical CRS (2D + 1D), where the vertical CRS of the geoid is added relative to the ellipsoid. + +In xDEM, we merge these into a single vertical CRS attribute {class}`DEM.vcrs` that takes two types of values: +- the string `"Ellipsoid"` for any ellipsoidal CRS promoted to 3D (e.g., the WGS84 ellipsoid), or +- a {class}`pyproj.CRS` with only a vertical axis for a CRS based on geoid heights (e.g., the EGM96 geoid). + +In practice, a {class}`pyproj.CRS` with only a vertical axis is either a {class}`~pyproj.crs.BoundCRS` (when created from a grid) or a {class}`~pyproj.crs.VerticalCRS` (when created in any other manner). + +## Methods to manipulate vertical CRSs + +The parsing, setting and transformation of vertical CRSs revolves around **three methods**, all described in details further below: +- The **instantiation** of {class}`~xdem.DEM` that implicitly tries to set the vertical CRS from the metadata (or explicitly through the `vcrs` argument), +- The **setting** method {func}`~xdem.DEM.set_vcrs` to explicitly set the vertical CRS of a {class}`~xdem.DEM`, +- The **transformation** method {func}`~xdem.DEM.to_vcrs` to explicitly transform the vertical CRS of a {class}`~xdem.DEM`. + +```{caution} +As of now, **[Rasterio](https://rasterio.readthedocs.io/en/stable/) does not support vertical transformations during CRS reprojection** (even if the CRS +provided contains a vertical axis). +We therefore advise to perform horizontal transformation and vertical transformation independently using {func}`DEM.reproject` and {func}`DEM.to_vcrs`, respectively. +``` + +(vref-setting)= +## Automated vertical CRS detection + +During instantiation of a {class}`~xdem.DEM`, the vertical CRS {attr}`~xdem.DEM.vcrs` is tentatively set with the following priority order: + +1. **From the {attr}`~xdem.DEM.crs` of the DEM**, if 3-dimensional, + +```{code-cell} ipython3 +:tags: [remove-cell] + +import xdem + +# Replace this with a new DEM in xdem-data +import numpy as np +import pyproj +import rasterio as rio +dem = xdem.DEM.from_array(data=np.ones((2,2)), + transform=rio.transform.from_bounds(0, 0, 1, 1, 2, 2), + crs=pyproj.CRS("EPSG:4326+5773"), + nodata=None) +dem.save("mydem_with3dcrs.tif") +``` + +```{code-cell} ipython3 +import xdem + +# Open a DEM with a 3D CRS +dem = xdem.DEM("mydem_with3dcrs.tif") +# First, let's look at what was the 3D CRS +pyproj.CRS(dem.crs) +``` + +```{code-cell} ipython3 +# The vertical CRS is extracted automatically +dem.vcrs +``` + +```{code-cell} ipython3 +:tags: [remove-cell] + +import os +os.remove("mydem_with3dcrs.tif") +``` + +2. **From the {attr}`~xdem.DEM.product` name of the DEM**, if parsed from the filename by {class}`geoutils.SatelliteImage`. + + +```{see-also} +The {class}`~geoutils.SatelliteImage` parent class that parses the product metadata is described in [a dedicated section of GeoUtils' documentation](https://geoutils.readthedocs.io/en/latest/satimg_class.html). +``` + +```{code-cell} ipython3 +:tags: [remove-cell] + +# Replace this with a new DEM in xdem-data +import rasterio as rio +dem = xdem.DEM.from_array(data=np.ones((2,2)), + transform=rio.transform.from_bounds(0, 0, 1, 1, 2, 2), + crs=pyproj.CRS("EPSG:4326"), + nodata=None) +# Save with the name of an ArcticDEM strip file +dem.save("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem.tif") +``` + +```{code-cell} ipython3 +# Open an ArcticDEM strip file, recognized as an ArcticDEM product by SatelliteImage +dem = xdem.DEM("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem.tif") +# The vertical CRS is set as "Ellipsoid" from knowing that is it an ArcticDEM product +dem.vcrs +``` + +```{code-cell} ipython3 +:tags: [remove-cell] + +os.remove("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem.tif") +``` + +**Currently recognized DEM products**: + +```{list-table} + :widths: 1 1 + :header-rows: 1 + + * - **DEM** + - **Vertical CRS** + * - [ArcticDEM](https://www.pgc.umn.edu/data/arcticdem/) + - Ellipsoid + * - [REMA](https://www.pgc.umn.edu/data/arcticdem/) + - Ellipsoid + * - [EarthDEM](https://www.pgc.umn.edu/data/earthdem/) + - Ellipsoid + * - [TanDEM-X global DEM](https://geoservice.dlr.de/web/dataguide/tdm90/) + - Ellipsoid + * - [NASADEM-HGTS](https://lpdaac.usgs.gov/documents/592/NASADEM_User_Guide_V1.pdf) + - Ellipsoid + * - [NASADEM-HGT](https://lpdaac.usgs.gov/documents/592/NASADEM_User_Guide_V1.pdf) + - EGM96 + * - [ALOS World 3D](https://www.eorc.jaxa.jp/ALOS/en/aw3d30/aw3d30v11_format_e.pdf) + - EGM96 + * - [SRTM v4.1](http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1) + - EGM96 + * - [ASTER GDEM v2-3](https://lpdaac.usgs.gov/documents/434/ASTGTM_User_Guide_V3.pdf) + - EGM96 + * - [Copernicus DEM](https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198) + - EGM08 +``` + +If your DEM does not have a `.vcrs` after instantiation, none of those steps worked. You can define the vertical CRS +explicitly during {class}`~xdem.DEM` instantiation with the `vcrs` argument or with {func}`~xdem.DEM.set_vcrs`, +with user inputs described below. + +## Setting a vertical CRS with convenient user inputs + +The vertical CRS of a {class}`~xdem.DEM` can be set or re-set manually at any point using {func}`~xdem.DEM.set_vcrs`. + +The `vcrs` argument, consistent across the three functions {class}`~xdem.DEM`, {func}`~xdem.DEM.set_vcrs` and {func}`~xdem.DEM.to_vcrs`, +accepts a **wide variety of user inputs**: + +- **Simple strings for the three most common references: `"Ellipsoid"`, `"EGM08"` or `"EGM96"`**, + +```{code-cell} ipython3 +# Set a geoid vertical CRS based on a string +dem.set_vcrs("EGM96") +dem.vcrs +``` + +```{code-cell} ipython3 +# Set a vertical CRS extended from the ellipsoid of the DEM's CRS +dem.set_vcrs("Ellipsoid") +dem.vcrs +``` + +- **Any PROJ grid name available at [https://cdn.proj.org/](https://cdn.proj.org/)**, + +```{tip} +**No need to download the grid!** This is done automatically during the setting operation, if the grid does not already exist locally. +``` + +```{code-cell} ipython3 +# Set a geoid vertical CRS based on a grid +dem.set_vcrs("us_noaa_geoid06_ak.tif") +dem.vcrs +``` + +- **Any EPSG code as {class}`int`**, + +```{code-cell} ipython3 +# Set a geoid vertical CRS based on an EPSG code +dem.set_vcrs(5773) +dem.vcrs +``` + +- **Any {class}`~pyproj.crs.CRS` that possesses a vertical axis**. + +```{code-cell} ipython3 +# Set a vertical CRS based on a pyproj.CRS +import pyproj +dem.set_vcrs(pyproj.CRS("EPSG:3855")) +dem.vcrs +``` + +## Transforming to another vertical CRS + +To transform a {class}`~xdem.DEM` to a different vertical CRS, {func}`~xdem.DEM.to_vcrs` is used. + +```{note} +If your transformation requires a grid that is not available locally, it will be **downloaded automatically**. +xDEM uses only the best available (i.e. best accuracy) transformation returned by {class}`pyproj.transformer.TransformerGroup`, considering the area-of-interest as the DEM extent {class}`~xdem.DEM.bounds`. +``` + +```{code-cell} ipython3 +# Open a DEM and set its CRS +filename_dem = xdem.examples.get_path("longyearbyen_ref_dem") +dem = xdem.DEM(filename_dem, vcrs="Ellipsoid") +dem.to_vcrs("EGM96") +dem.vcrs +``` + +The operation updates the DEM array **in-place**, shifting each pixel by the transformation at their coordinates: + +```{code-cell} ipython3 +import numpy as np + +# Mean difference after transformation (about 30 m in Svalbard) +dem_orig = xdem.DEM(filename_dem) +diff = dem_orig - dem +np.nanmean(diff) +``` diff --git a/environment.yml b/environment.yml index 798d2ed3..7f239c4c 100644 --- a/environment.yml +++ b/environment.yml @@ -12,12 +12,11 @@ dependencies: - matplotlib - opencv - openh264 - - pyproj + - pyproj>=3.4 - rasterio>=1.3 - scipy - tqdm - scikit-image - - proj-data - scikit-gstat>=0.6.8 - pytransform3d - geoutils==0.0.11 diff --git a/tests/test_dem.py b/tests/test_dem.py index aaf5b301..0179d06d 100644 --- a/tests/test_dem.py +++ b/tests/test_dem.py @@ -1,16 +1,22 @@ """ Functions to test the DEM tools.""" +from __future__ import annotations + import os +import tempfile import warnings +from typing import Any import geoutils.raster as gr import geoutils.raster.satimg as si import numpy as np -import pyproj import pytest import rasterio as rio from geoutils.raster.raster import _default_rio_attrs +from pyproj import CRS +from pyproj.transformer import Transformer import xdem +import xdem.vcrs from xdem.dem import DEM DO_PLOT = False @@ -75,6 +81,92 @@ def test_init(self) -> None: nodata=None, ) + def test_init__vcrs(self) -> None: + """Test that vcrs is set properly during instantiation.""" + + # Tests 1: instantiation with a file that has a 2D CRS + + # First, check a DEM that does not have any vertical CRS set + fn_img = xdem.examples.get_path("longyearbyen_ref_dem") + dem = DEM(fn_img) + assert dem.vcrs is None + + # Setting a vertical CRS during instantiation should work here + dem = DEM(fn_img, vcrs="EGM96") + assert dem.vcrs_name == "EGM96 height" + + # Tests 2: instantiation with a file that has a 3D CRS + # Create such a file + dem = DEM(fn_img) + dem_reproj = dem.reproject(dst_crs=4979) + + # Save to temporary folder + temp_dir = tempfile.TemporaryDirectory() + temp_file = os.path.join(temp_dir.name, "test.tif") + dem_reproj.save(temp_file) + + # Check opening a DEM with a 3D CRS sets the vcrs + dem_3d = DEM(temp_file) + assert dem_3d.vcrs == "Ellipsoid" + + # Check that a warning is raised when trying to override with user input + with pytest.warns( + UserWarning, + match="The CRS in the raster metadata already has a vertical component, " + "the user-input 'EGM08' will override it.", + ): + DEM(temp_file, vcrs="EGM08") + + def test_from_array(self) -> None: + """Test that overridden from_array works as expected.""" + + # Create a 5x5 DEM + data = np.ones((5, 5)) + transform = rio.transform.from_bounds(0, 0, 1, 1, 5, 5) + crs = CRS("EPSG:4326") + nodata = -9999 + vcrs = "EGM08" + dem = DEM.from_array(data=data, transform=transform, crs=crs, nodata=nodata, vcrs=vcrs) + + # Check output matches + assert isinstance(dem, DEM) + assert isinstance(dem.data, np.ma.masked_array) + assert np.array_equal(dem.data.data, np.ones((5, 5))) + assert dem.transform == transform + assert dem.crs == crs + assert dem.nodata == nodata + assert dem.vcrs == xdem.vcrs._vcrs_from_user_input(vcrs_input=vcrs) + + def test_from_array__vcrs(self) -> None: + """Test that overridden from_array rightly sets the vertical CRS.""" + + # Create a 5x5 DEM with a 2D CRS + transform = rio.transform.from_bounds(0, 0, 1, 1, 5, 5) + dem = DEM.from_array(data=np.ones((5, 5)), transform=transform, crs=CRS("EPSG:4326"), nodata=None, vcrs=None) + assert dem.vcrs is None + + # One with a 3D ellipsoid CRS + dem = DEM.from_array(data=np.ones((5, 5)), transform=transform, crs=CRS("EPSG:4979"), nodata=None, vcrs=None) + assert dem.vcrs == "Ellipsoid" + + # One with a 2D and the ellipsoid vertical CRS + dem = DEM.from_array( + data=np.ones((5, 5)), transform=transform, crs=CRS("EPSG:4326"), nodata=None, vcrs="Ellipsoid" + ) + assert dem.vcrs == "Ellipsoid" + + # One with a compound CRS + dem = DEM.from_array( + data=np.ones((5, 5)), transform=transform, crs=CRS("EPSG:4326+5773"), nodata=None, vcrs=None + ) + assert dem.vcrs == CRS("EPSG:5773") + + # One with a CRS and vertical CRS + dem = DEM.from_array( + data=np.ones((5, 5)), transform=transform, crs=CRS("EPSG:4326"), nodata=None, vcrs=CRS("EPSG:5773") + ) + assert dem.vcrs == CRS("EPSG:5773") + def test_copy(self) -> None: """ Test that the copy method works as expected for DEM. In particular @@ -97,7 +189,7 @@ def test_copy(self) -> None: # raster_attrs = ['bounds', 'count', 'crs', 'dtypes', 'height', 'indexes', 'nodata', # 'res', 'shape', 'transform', 'width'] # satimg_attrs = ['satellite', 'sensor', 'product', 'version', 'tile_name', 'datetime'] - # dem_attrs = ['vref', 'vref_grid', 'ccrs'] + # dem_attrs = ['vcrs', 'vcrs_grid', 'vcrs_name', 'ccrs'] # using list directly available in Class attrs = [at for at in _default_rio_attrs if at not in ["name", "dataset_mask", "driver"]] @@ -120,132 +212,129 @@ def test_copy(self) -> None: assert np.array_equal(r3.data, r2.data) - def test_set_vref(self) -> None: - """Tests to set the vertical reference""" + def test_set_vcrs(self) -> None: + """Tests to set the vertical CRS.""" - fn_img = xdem.examples.get_path("longyearbyen_ref_dem") - img = DEM(fn_img) + fn_dem = xdem.examples.get_path("longyearbyen_ref_dem") + dem = DEM(fn_dem) - # Check setting WGS84 - img.set_vref(vref_name="WGS84") - assert img.vref == "WGS84" - assert img.vref_grid is None + # -- Test 1: we check with names -- + + # Check setting ellipsoid + dem.set_vcrs(new_vcrs="Ellipsoid") + assert dem.vcrs_name is not None + assert "Ellipsoid (No vertical CRS)." in dem.vcrs_name + assert dem.vcrs_grid is None # Check setting EGM96 - img.set_vref(vref_name="EGM96") - assert img.vref == "EGM96" - assert img.vref_grid == "us_nga_egm96_15.tif" - # The grid argument should have priority over name and parse the right vref name - img.set_vref(vref_name="WGS84", vref_grid="us_nga_egm96_15.tif") - assert img.vref == "EGM96" + dem.set_vcrs(new_vcrs="EGM96") + assert dem.vcrs_name == "EGM96 height" + assert dem.vcrs_grid == "us_nga_egm96_15.tif" # Check setting EGM08 - img.set_vref(vref_name="EGM08") - assert img.vref == "EGM08" - assert img.vref_grid == "us_nga_egm08_25.tif" - # The grid argument should have priority over name and parse the right vref name - img.set_vref(vref_name="best ref in the entire world, or any string", vref_grid="us_nga_egm08_25.tif") - assert img.vref == "EGM08" + dem.set_vcrs(new_vcrs="EGM08") + assert dem.vcrs_name == "EGM2008 height" + assert dem.vcrs_grid == "us_nga_egm08_25.tif" + + # -- Test 2: we check with grids -- + dem.set_vcrs(new_vcrs="us_nga_egm96_15.tif") + assert dem.vcrs_name == "unknown using geoidgrids=us_nga_egm96_15.tif" + assert dem.vcrs_grid == "us_nga_egm96_15.tif" + + dem.set_vcrs(new_vcrs="us_nga_egm08_25.tif") + assert dem.vcrs_name == "unknown using geoidgrids=us_nga_egm08_25.tif" + assert dem.vcrs_grid == "us_nga_egm08_25.tif" # Check that other existing grids are well detected in the pyproj.datadir - # TODO: Figure out why CI cannot get the grids on Windows - if os.name != "nt": - img.set_vref(vref_grid="is_lmi_Icegeoid_ISN93.tif") - else: - with pytest.raises(ValueError): - img.set_vref(vref_grid="is_lmi_Icegeoid_ISN93.tif") + dem.set_vcrs(new_vcrs="is_lmi_Icegeoid_ISN93.tif") # Check that non-existing grids raise errors - with pytest.raises(ValueError): - img.set_vref(vref_grid="the best grid in the entire world, or any non-existing string") - - @pytest.mark.skip("This fails on Windows because the grids are not found") # type: ignore - def test_to_vref(self) -> None: - """Tests to convert vertical references""" - - # First, we use test points to test the vertical transform - # Let's start with Chile - lat = 43.70012234 - lng = -79.41629234 - z = 100 - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", module="pyproj") - # init is deprecated by - ellipsoid = pyproj.Proj(init="EPSG:4326") # WGS84 datum ellipsoid height - # EGM96 geoid in Chile, we expect ~30 m difference - geoid = pyproj.Proj(init="EPSG:4326", geoidgrids="us_nga_egm96_15.tif") - transformer = pyproj.Transformer.from_proj(ellipsoid, geoid) - z_out = transformer.transform(lng, lat, z)[2] - - # Check that the final elevation is finite, and higher than ellipsoid by less than 40 m (typical geoid in Chile) - assert np.logical_and.reduce((np.isfinite(z_out), np.greater(z_out, z), np.less(np.abs(z_out - z), 40))) - - # With the EGM2008 (catch warnings as this use of init is depecrated) - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", module="pyproj") - ellipsoid = pyproj.Proj(init="EPSG:4326") # WGS84 datum ellipsoid height - geoid = pyproj.Proj(init="EPSG:4326", geoidgrids="us_nga_egm08_25.tif") - transformer = pyproj.Transformer.from_proj(ellipsoid, geoid) - z_out = transformer.transform(lng, lat, z)[2] - - # Check final elevation is finite, higher than ellipsoid with less than 40 m difference (typical geoid in Chile) - assert np.logical_and.reduce((np.isfinite(z_out), np.greater(z_out, z), np.less(np.abs(z_out - z), 40))) - - # With GEOID2006 for Alaska - lat = 65 - lng = -140 - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", module="pyproj") - # init is deprecated by - ellipsoid = pyproj.Proj(init="EPSG:4326") # WGS84 datum ellipsoid height - geoid = pyproj.Proj(init="EPSG:4326", geoidgrids="us_noaa_geoid06_ak.tif") - transformer = pyproj.Transformer.from_proj(ellipsoid, geoid) - z_out = transformer.transform(lng, lat, z)[2] - - # Check that the final elevation is finite, lower than ellipsoid by less than 20 m (typical geoid in Alaska) - assert np.logical_and.reduce((np.isfinite(z_out), np.less(z_out, z), np.less(np.abs(z_out - z), 20))) - - # With ISN1993 for Iceland - lat = 65 - lng = -18 - # TODO: Figure out why CI cannot get the grids on Windows - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", module="pyproj") - # init is deprecated by - ellipsoid = pyproj.Proj(init="EPSG:4326") # WGS84 datum ellipsoid height - # Iceland, we expect a ~70m difference - geoid = pyproj.Proj(init="EPSG:4326", geoidgrids="is_lmi_Icegeoid_ISN93.tif") - transformer = pyproj.Transformer.from_proj(ellipsoid, geoid) - z_out = transformer.transform(lng, lat, z)[2] - - # Check that the final elevation is finite, lower than ellipsoid by less than 100 m (typical geoid in Iceland) - assert np.logical_and.reduce((np.isfinite(z_out), np.less(z_out, z), np.less(np.abs(z_out - z), 100))) - - # Check that the function does not run without a reference set - fn_img = xdem.examples.get_path("longyearbyen_ref_dem") - img = DEM(fn_img) - with pytest.raises(ValueError): - img.to_vref(vref_name="EGM96") - - # Check that the function properly runs with a reference set - img.set_vref(vref_name="WGS84") - mean_ellips = np.nanmean(img.data) - img.to_vref(vref_name="EGM96") - mean_geoid_96 = np.nanmean(img.data) - assert img.vref == "EGM96" - assert img.vref_grid == "us_nga_egm96_15.tif" - # Check that the geoid is lower than ellipsoid, less than 35 m difference (Svalbard) - assert np.greater(mean_ellips, mean_geoid_96) - assert np.less(np.abs(mean_ellips - mean_geoid_96), 35.0) - - # Check in the other direction - img = DEM(fn_img) - img.set_vref(vref_name="EGM96") - mean_geoid_96 = np.nanmean(img.data) - img.to_vref(vref_name="WGS84") - mean_ellips = np.nanmean(img.data) - assert img.vref == "WGS84" - assert img.vref_grid is None - # Check that the geoid is lower than ellipsoid, less than 35 m difference (Svalbard) - assert np.greater(mean_ellips, mean_geoid_96) - assert np.less(np.abs(mean_ellips - mean_geoid_96), 35.0) + with pytest.raises( + ValueError, + match="The provided grid 'the best grid' does not exist at https://cdn.proj.org/. " + "Provide an existing grid.", + ): + dem.set_vcrs(new_vcrs="the best grid") + + def test_to_vcrs(self) -> None: + """Tests the conversion of vertical CRS.""" + + fn_dem = xdem.examples.get_path("longyearbyen_ref_dem") + dem = DEM(fn_dem) + + # Reproject in WGS84 2D + dem = dem.reproject(dst_crs=4326) + dem_before_trans = dem.copy() + + # Set ellipsoid as vertical reference + dem.set_vcrs(new_vcrs="Ellipsoid") + ccrs_init = dem.ccrs + median_before = np.nanmean(dem) + # Transform to EGM96 geoid + dem.to_vcrs(dst_vcrs="EGM96") + median_after = np.nanmean(dem) + + # About 32 meters of difference in Svalbard between EGM96 geoid and ellipsoid + assert median_after - median_before == pytest.approx(-32, rel=0.1) + + # Check that the results are consistent with the operation done independently + ccrs_dest = xdem.vcrs._build_ccrs_from_crs_and_vcrs(dem.crs, xdem.vcrs._vcrs_from_user_input("EGM96")) + transformer = Transformer.from_crs(crs_from=ccrs_init, crs_to=ccrs_dest, always_xy=True) + + xx, yy = dem.coords() + x = xx[5, 5] + y = yy[5, 5] + z = dem_before_trans.data[5, 5] + z_out = transformer.transform(xx=x, yy=y, zz=z)[2] + + assert z_out == pytest.approx(dem.data[5, 5]) + + def test_to_vcrs__equal_warning(self) -> None: + """Test that DEM.to_vcrs() does not transform if both 3D CRS are equal.""" + + fn_dem = xdem.examples.get_path("longyearbyen_ref_dem") + dem = DEM(fn_dem) + + # With both inputs as names + dem.set_vcrs("EGM96") + with pytest.warns( + UserWarning, match="Source and destination vertical CRS are the same, " "skipping vertical transformation." + ): + dem.to_vcrs("EGM96") + + # With one input as name, the other as CRS + dem.set_vcrs("Ellipsoid") + with pytest.warns( + UserWarning, match="Source and destination vertical CRS are the same, " "skipping vertical transformation." + ): + dem.to_vcrs(CRS("EPSG:4979")) + + # Compare to manually-extracted shifts at specific coordinates for the geoid grids + egm96_chile = {"grid": "us_nga_egm96_15.tif", "lon": -68, "lat": -20, "shift": 42} + egm08_chile = {"grid": "us_nga_egm08_25.tif", "lon": -68, "lat": -20, "shift": 42} + geoid96_alaska = {"grid": "us_noaa_geoid06_ak.tif", "lon": -145, "lat": 62, "shift": 17} + isn93_iceland = {"grid": "is_lmi_Icegeoid_ISN93.tif", "lon": -18, "lat": 65, "shift": 68} + + @pytest.mark.parametrize("grid_shifts", [egm08_chile, egm08_chile, geoid96_alaska, isn93_iceland]) # type: ignore + def test_to_vcrs__grids(self, grid_shifts: dict[str, Any]) -> None: + """Tests grids to convert vertical CRS.""" + + # Using an arbitrary elevation of 100 m (no influence on the transformation) + dem = DEM.from_array( + data=np.array([[100]]), + transform=rio.transform.from_bounds( + grid_shifts["lon"], grid_shifts["lat"], grid_shifts["lon"] + 0.01, grid_shifts["lat"] + 0.01, 0.01, 0.01 + ), + crs=CRS.from_epsg(4326), + nodata=None, + ) + dem.set_vcrs("Ellipsoid") + + # Transform to the vertical CRS of the grid + dem.to_vcrs(grid_shifts["grid"]) + + # Compare the elevation difference + z_diff = 100 - dem.data[0, 0] + + # Check the shift is the one expect within 10% + assert z_diff == pytest.approx(grid_shifts["shift"], rel=0.1) diff --git a/tests/test_spatialstats.py b/tests/test_spatialstats.py index 9cffc05f..f9c4b73c 100644 --- a/tests/test_spatialstats.py +++ b/tests/test_spatialstats.py @@ -2,7 +2,6 @@ from __future__ import annotations import os -import time import warnings from typing import Any @@ -434,9 +433,9 @@ def test_sample_empirical_variogram_speed(self) -> None: subsample = 10 # First, run the xDEM wrapper function - t0 = time.time() + # t0 = time.time() df = xdem.spatialstats.sample_empirical_variogram(values=values, subsample=subsample, random_state=42) - t1 = time.time() + # t1 = time.time() # Second, do it manually with skgstat @@ -476,7 +475,7 @@ def test_sample_empirical_variogram_speed(self) -> None: # Index of valid values values_arr, mask_nodata = gu.raster.get_array_and_mask(values) - t3 = time.time() + # t3 = time.time() rems = skgstat.RasterEquidistantMetricSpace( coords=coords[~mask_nodata.ravel(), :], shape=shape, @@ -494,7 +493,7 @@ def test_sample_empirical_variogram_speed(self) -> None: bin_func=bin_func, maxlag=maxlag, ) - t4 = time.time() + # t4 = time.time() # Get bins, empirical variogram values, and bin count bins, exp = V.get_empirical(bin_center=False) @@ -508,20 +507,20 @@ def test_sample_empirical_variogram_speed(self) -> None: df2.drop(df2.tail(1).index, inplace=True) df2 = df2.astype({"exp": "float64", "err_exp": "float64", "lags": "float64", "count": "int64"}) - t2 = time.time() + # t2 = time.time() # Check if the two frames are equal pd.testing.assert_frame_equal(df, df2) # Check that the two ways are taking the same time with 50% margin - time_method_1 = t1 - t0 - time_method_2 = t2 - t1 - assert time_method_1 == pytest.approx(time_method_2, rel=0.5) + # time_method_1 = t1 - t0 + # time_method_2 = t2 - t1 + # assert time_method_1 == pytest.approx(time_method_2, rel=0.5) # Check that all this time is based on variogram sampling at about 70%, even with the smallest number of # samples of 10 - time_metricspace_variogram = t4 - t3 - assert time_metricspace_variogram == pytest.approx(time_method_2, rel=0.3) + # time_metricspace_variogram = t4 - t3 + # assert time_metricspace_variogram == pytest.approx(time_method_2, rel=0.3) @pytest.mark.parametrize( "subsample_method", ["pdist_point", "pdist_ring", "pdist_disk", "cdist_point"] diff --git a/tests/test_vcrs.py b/tests/test_vcrs.py new file mode 100644 index 00000000..d9f2c61f --- /dev/null +++ b/tests/test_vcrs.py @@ -0,0 +1,201 @@ +"""Tests for vertical CRS transformation tools.""" +from __future__ import annotations + +import pathlib +import re +from typing import Any + +import numpy as np +import pytest +from pyproj import CRS + +import xdem +import xdem.vcrs + + +class TestVCRS: + def test_parse_vcrs_name_from_product(self) -> None: + """Test parsing of vertical CRS name from DEM product name.""" + + # Check that the value for the key is returned by the function + for product in xdem.vcrs.vcrs_dem_products.keys(): + assert xdem.vcrs._parse_vcrs_name_from_product(product) == xdem.vcrs.vcrs_dem_products[product] + + # And that, otherwise, it's a None + assert xdem.vcrs._parse_vcrs_name_from_product("BESTDEM") is None + + # Expect outputs for the inputs + @pytest.mark.parametrize( + "input_output", + [ + (CRS("EPSG:4326"), None), + (CRS("EPSG:4979"), "Ellipsoid"), + (CRS("EPSG:4326+5773"), CRS("EPSG:5773")), + (CRS("EPSG:32610"), None), + (CRS("EPSG:32610").to_3d(), "Ellipsoid"), + ], + ) # type: ignore + def test_vcrs_from_crs(self, input_output: tuple[CRS, CRS]) -> None: + """Test the extraction of a vertical CRS from a CRS.""" + + input = input_output[0] + output = input_output[1] + + # Extract vertical CRS from CRS + vcrs = xdem.vcrs._vcrs_from_crs(crs=input) + + # Check that the result is as expected + if isinstance(output, CRS): + assert vcrs.equals(input_output[1]) + elif isinstance(output, str): + assert vcrs == "Ellipsoid" + else: + assert vcrs is None + + @pytest.mark.parametrize( + "vcrs_input", + [ + "EGM08", + "EGM96", + "us_noaa_geoid06_ak.tif", + pathlib.Path("is_lmi_Icegeoid_ISN93.tif"), + 3855, + CRS.from_epsg(5773), + ], + ) # type: ignore + def test_vcrs_from_user_input(self, vcrs_input: str | pathlib.Path | int | CRS) -> None: + """Tests the function _vcrs_from_user_input for varying user inputs, for which it will return a CRS.""" + + # Get user input + vcrs = xdem.dem._vcrs_from_user_input(vcrs_input) + + # Check output type + assert isinstance(vcrs, CRS) + assert vcrs.is_vertical + + @pytest.mark.parametrize( + "vcrs_input", ["Ellipsoid", "ellipsoid", "wgs84", 4326, 4979, CRS.from_epsg(4326), CRS.from_epsg(4979)] + ) # type: ignore + def test_vcrs_from_user_input__ellipsoid(self, vcrs_input: str | int) -> None: + """Tests the function _vcrs_from_user_input for inputs where it returns "Ellipsoid".""" + + # Get user input + vcrs = xdem.vcrs._vcrs_from_user_input(vcrs_input) + + # Check output type + assert vcrs == "Ellipsoid" + + def test_vcrs_from_user_input__errors(self) -> None: + """Tests errors of vcrs_from_user_input.""" + + # Check that an error is raised when the type is wrong + with pytest.raises(TypeError, match="New vertical CRS must be a string, path or VerticalCRS, received.*"): + xdem.vcrs._vcrs_from_user_input(np.zeros(1)) # type: ignore + + # Check that an error is raised if the CRS is not vertical + with pytest.raises( + ValueError, + match=re.escape( + "New vertical CRS must have a vertical axis, 'WGS 84 / UTM " + "zone 1N' does not (check with `CRS.is_vertical`)." + ), + ): + xdem.vcrs._vcrs_from_user_input(32601) + + # Check that a warning is raised if the CRS has other dimensions than vertical + with pytest.warns( + UserWarning, + match="New vertical CRS has a vertical dimension but also other components, " + "extracting the vertical reference only.", + ): + xdem.vcrs._vcrs_from_user_input(CRS("EPSG:4326+5773")) + + @pytest.mark.parametrize( + "grid", ["us_noaa_geoid06_ak.tif", "is_lmi_Icegeoid_ISN93.tif", "us_nga_egm08_25.tif", "us_nga_egm96_15.tif"] + ) # type: ignore + def test_build_vcrs_from_grid(self, grid: str) -> None: + """Test that vertical CRS are correctly built from grid""" + + # Build vertical CRS + vcrs = xdem.vcrs._build_vcrs_from_grid(grid=grid) + assert vcrs.is_vertical + + # Check that the explicit construction yields the same CRS as "the old init way" (see function description) + vcrs_oldway = xdem.vcrs._build_vcrs_from_grid(grid=grid, old_way=True) + assert vcrs.equals(vcrs_oldway) + + # Test for WGS84 in 2D and 3D, UTM, CompoundCRS, everything should work + @pytest.mark.parametrize( + "crs", [CRS("EPSG:4326"), CRS("EPSG:4979"), CRS("32610"), CRS("EPSG:4326+5773")] + ) # type: ignore + @pytest.mark.parametrize("vcrs_input", [CRS("EPSG:5773"), "is_lmi_Icegeoid_ISN93.tif", "EGM96"]) # type: ignore + def test_build_ccrs_from_crs_and_vcrs(self, crs: CRS, vcrs_input: CRS | str) -> None: + """Test the function build_ccrs_from_crs_and_vcrs.""" + + # Get the vertical CRS from user input + vcrs = xdem.vcrs._vcrs_from_user_input(vcrs_input=vcrs_input) + + # Build the compound CRS + + # For a 3D horizontal CRS, a condition based on pyproj version is needed + if len(crs.axis_info) > 2: + import pyproj + from packaging.version import Version + + # If the version is higher than 3.5.0, it should pass + if Version(pyproj.__version__) > Version("3.5.0"): + ccrs = xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=crs, vcrs=vcrs) + # Otherwise, it should raise an error + else: + with pytest.raises( + NotImplementedError, + match="pyproj >= 3.5.1 is required to demote a 3D CRS to 2D and be able to compound " + "with a new vertical CRS. Update your dependencies or pass the 2D source CRS " + "manually.", + ): + xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=crs, vcrs=vcrs) + return None + # If the CRS is 2D, it should pass + else: + ccrs = xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=crs, vcrs=vcrs) + + assert isinstance(ccrs, CRS) + assert ccrs.is_vertical + + def test_build_ccrs_from_crs_and_vcrs__errors(self) -> None: + """Test errors are correctly raised from the build_ccrs function.""" + + with pytest.raises( + ValueError, match="Invalid vcrs given. Must be a vertical " "CRS or the literal string 'Ellipsoid'." + ): + xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=CRS("EPSG:4326"), vcrs="NotAVerticalCRS") # type: ignore + + # Compare to manually-extracted shifts at specific coordinates for the geoid grids + egm96_chile = {"grid": "us_nga_egm96_15.tif", "lon": -68, "lat": -20, "shift": 42} + egm08_chile = {"grid": "us_nga_egm08_25.tif", "lon": -68, "lat": -20, "shift": 42} + geoid96_alaska = {"grid": "us_noaa_geoid06_ak.tif", "lon": -145, "lat": 62, "shift": 15} + isn93_iceland = {"grid": "is_lmi_Icegeoid_ISN93.tif", "lon": -18, "lat": 65, "shift": 68} + + @pytest.mark.parametrize("grid_shifts", [egm08_chile, egm08_chile, geoid96_alaska, isn93_iceland]) # type: ignore + def test_transform_zz(self, grid_shifts: dict[str, Any]) -> None: + """Tests grids to convert vertical CRS.""" + + # Using an arbitrary elevation of 100 m (no influence on the transformation) + zz = 100 + xx = grid_shifts["lon"] + yy = grid_shifts["lat"] + crs_from = CRS.from_epsg(4326) + ccrs_from = xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=crs_from, vcrs="Ellipsoid") + + # Build the compound CRS + vcrs_to = xdem.vcrs._vcrs_from_user_input(vcrs_input=grid_shifts["grid"]) + ccrs_to = xdem.vcrs._build_ccrs_from_crs_and_vcrs(crs=crs_from, vcrs=vcrs_to) + + # Apply the transformation + zz_trans = xdem.vcrs._transform_zz(crs_from=ccrs_from, crs_to=ccrs_to, xx=xx, yy=yy, zz=zz) + + # Compare the elevation difference + z_diff = 100 - zz_trans + + # Check the shift is the one expect within 10% + assert z_diff == pytest.approx(grid_shifts["shift"], rel=0.1) diff --git a/xdem/ddem.py b/xdem/ddem.py index d11888a1..a510f4b9 100644 --- a/xdem/ddem.py +++ b/xdem/ddem.py @@ -7,15 +7,15 @@ import geoutils as gu import numpy as np import shapely -from geoutils.raster import RasterType, get_array_and_mask +from geoutils.raster import Raster, RasterType, get_array_and_mask from rasterio.crs import CRS from rasterio.warp import Affine import xdem -from xdem._typing import NDArrayf +from xdem._typing import MArrayf, NDArrayf -class dDEM(xdem.dem.DEM): # pylint: disable=invalid-name +class dDEM(Raster): # type: ignore """A difference-DEM object.""" def __init__(self, raster: gu.Raster, start_time: np.datetime64, end_time: np.datetime64, error: Any | None = None): @@ -91,14 +91,14 @@ def time(self) -> np.timedelta64: @classmethod def from_array( cls: type[RasterType], - data: NDArrayf, + data: NDArrayf | MArrayf, transform: tuple[float, ...] | Affine, crs: CRS | int | None, start_time: np.datetime64, end_time: np.datetime64, - error: float = None, nodata: int | float | None = None, - ) -> dDEM: + error: float = None, + ) -> dDEM: # type: ignore """ Create a new dDEM object from an array. @@ -112,7 +112,7 @@ def from_array( :returns: A new dDEM instance. """ - return dDEM( + return cls( gu.Raster.from_array(data=data, transform=transform, crs=crs, nodata=nodata), start_time=start_time, end_time=end_time, diff --git a/xdem/dem.py b/xdem/dem.py index dc693cd3..d1dca259 100644 --- a/xdem/dem.py +++ b/xdem/dem.py @@ -1,76 +1,92 @@ """DEM class and functions.""" from __future__ import annotations -import json -import os -import subprocess +import pathlib import warnings -from typing import Any +from typing import Any, Literal -import pyproj import rasterio as rio +from affine import Affine from geoutils import SatelliteImage from geoutils.raster import RasterType -from pyproj import Transformer +from pyproj import CRS +from pyproj.crs import CompoundCRS, VerticalCRS -from xdem._typing import NDArrayf +from xdem._typing import MArrayf, NDArrayf +from xdem.vcrs import ( + _build_ccrs_from_crs_and_vcrs, + _grid_from_user_input, + _parse_vcrs_name_from_product, + _transform_zz, + _vcrs_from_crs, + _vcrs_from_user_input, +) +dem_attrs = ["_vcrs", "_vcrs_name", "_vcrs_grid"] -def parse_vref_from_product(product: str) -> str | None: - """ - - :param product: Product name (typically from satimg.parse_metadata_from_fn) - :return: vref_name: Vertical reference name +class DEM(SatelliteImage): # type: ignore + """ + The digital elevation model. + + The DEM has a single main attribute in addition to that inherited from :class:`geoutils.Raster`: + vcrs: :class:`pyproj.VerticalCRS` + Vertical coordinate reference system of the DEM. + + Other derivative attributes are: + vcrs_name: :class:`str` + Name of vertical CRS of the DEM. + vcrs_grid: :class:`str` + Grid path to the vertical CRS of the DEM. + ccrs: :class:`pyproj.CompoundCRS` + Compound vertical and horizontal CRS of the DEM. + + The attributes inherited from :class:`geoutils.Raster` are: + data: :class:`np.ndarray` + Data array of the DEM, with dimensions corresponding to (count, height, width). + transform: :class:`affine.Affine` + Geotransform of the DEM. + crs: :class:`pyproj.crs.CRS` + Coordinate reference system of the DEM. + nodata: :class:`int` or :class:`float` + Nodata value of the DEM. + + All other attributes are derivatives of those attributes, or read from the file on disk. + See the API for more details. """ - # Sources for defining vertical references: - # AW3D30: https://www.eorc.jaxa.jp/ALOS/en/aw3d30/aw3d30v11_format_e.pdf - # SRTMGL1: https://lpdaac.usgs.gov/documents/179/SRTM_User_Guide_V3.pdf - # SRTMv4.1: http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1 - # ASTGTM2/ASTGTM3: https://lpdaac.usgs.gov/documents/434/ASTGTM_User_Guide_V3.pdf - # NASADEM: https://lpdaac.usgs.gov/documents/592/NASADEM_User_Guide_V1.pdf, HGTS is ellipsoid, HGT is EGM96 geoid !! - # ArcticDEM (mosaic and strips): https://www.pgc.umn.edu/data/arcticdem/ - # REMA (mosaic and strips): https://www.pgc.umn.edu/data/rema/ - # TanDEM-X 90m global: https://geoservice.dlr.de/web/dataguide/tdm90/ - # COPERNICUS DEM: https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198 - - if product in ["ArcticDEM/REMA", "TDM1", "NASADEM-HGTS"]: - vref_name = "WGS84" - elif product in ["AW3D30", "SRTMv4.1", "SRTMGL1", "ASTGTM2", "NASADEM-HGT"]: - vref_name = "EGM96" - elif product in ["COPDEM"]: - vref_name = "EGM08" - else: - vref_name = None - - return vref_name - - -dem_attrs = ["vref", "vref_grid", "_ccrs"] - -class DEM(SatelliteImage): # type: ignore def __init__( self, filename_or_dataset: str | RasterType | rio.io.DatasetReader | rio.io.MemoryFile, - vref_name: str | None = None, - vref_grid: str | None = None, + vcrs: Literal["Ellipsoid"] + | Literal["EGM08"] + | Literal["EGM96"] + | VerticalCRS + | str + | pathlib.Path + | int + | None = None, silent: bool = True, **kwargs: Any, ) -> None: """ - Load digital elevation model data through the Raster class, parse additional attributes from filename or - metadata through the SatelliteImage class, and then parse vertical reference from DEM product name. - For manual input, only one of "vref", "vref_grid" or "ccrs" is necessary to set the vertical reference. + Instantiate a digital elevation model. + + The vertical reference of the DEM can be defined by passing the `vcrs` argument. + Otherwise, a vertical reference is tentatively parsed from the DEM product name. + + Inherits all attributes from the :class:`geoutils.Raster` and :class:`geoutils.SatelliteImage` classes. :param filename_or_dataset: The filename of the dataset. - :param vref_name: Vertical reference name - :param vref_grid: Vertical reference grid (any grid file in https://github.com/OSGeo/PROJ-data) - :param silent: Whether to display vertical reference setting - :param silent: boolean + :param vcrs: Vertical coordinate reference system either as a name ("WGS84", "EGM08", "EGM96"), + an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data). + :param silent: Whether to display vertical reference parsing. """ self.data: NDArrayf + self._vcrs: VerticalCRS | Literal["Ellipsoid"] | None = None + self._vcrs_name: str | None = None + self._vcrs_grid: str | None = None # If DEM is passed, simply point back to DEM if isinstance(filename_or_dataset, DEM): @@ -83,173 +99,189 @@ def __init__( warnings.filterwarnings("ignore", message="Parse metadata from file not implemented") super().__init__(filename_or_dataset, silent=silent, **kwargs) - # self.indexes can be None when data is not loaded through the Raster class + # Ensure DEM has only one band: self.indexes can be None when data is not loaded through the Raster class if self.indexes is not None and len(self.indexes) > 1: raise ValueError("DEM rasters should be composed of one band only") - # user input - self.vref = vref_name - self.vref_grid = vref_grid - self._ccrs = None + # If the CRS in the raster metadata has a 3rd dimension, could set it as a vertical reference + vcrs_from_crs = _vcrs_from_crs(CRS(self.crs)) + if vcrs_from_crs is not None: + # If something was also provided by the user, user takes precedence + # (we leave vcrs as it was for input) + if vcrs is not None: + # Raise a warning if the two are not the same + vcrs_user = _vcrs_from_user_input(vcrs) + if not vcrs_from_crs == vcrs_user: + warnings.warn( + "The CRS in the raster metadata already has a vertical component, " + "the user-input '{}' will override it.".format(vcrs) + ) + # Otherwise, use the one from the raster 3D CRS + else: + vcrs = vcrs_from_crs + + # If no vertical CRS was provided by the user or defined in the CRS + if vcrs is None: + vcrs = _parse_vcrs_name_from_product(self.product) - # trying to get vref from product name (priority to user input) - self.__parse_vref_from_fn(silent=silent) + # If a vertical reference was parsed or provided by user + if vcrs is not None: + self.set_vcrs(vcrs) def copy(self, new_array: NDArrayf | None = None) -> DEM: + """ + Copy the DEM, possibly updating the data array. + + :param new_array: New data array. + + :return: Copied DEM. + """ new_dem = super().copy(new_array=new_array) # type: ignore # The rest of attributes are immutable, including pyproj.CRS - # dem_attrs = ['vref','vref_grid','ccrs'] #taken outside of class for attrs in dem_attrs: setattr(new_dem, attrs, getattr(self, attrs)) return new_dem # type: ignore - def __parse_vref_from_fn(self, silent: bool = False) -> None: - """Attempts to pull vertical reference from product name identified by SatImg.""" - - if self.product is not None: - vref = parse_vref_from_product(self.product) - if vref is not None and self.vref is None: - if not silent: - print('From product name "' + str(self.product) + '": setting vertical reference as ' + str(vref)) - self.vref = vref - elif vref is not None and self.vref is not None: - if not silent: - print( - "Leaving user input of " - + str(self.vref) - + " for vertical reference despite reading " - + str(vref) - + " from product name" - ) - else: - if not silent: - print('Could not find a vertical reference based on product name: "' + str(self.product) + '"') + @classmethod + def from_array( + cls: type[DEM], + data: NDArrayf | MArrayf, + transform: tuple[float, ...] | Affine, + crs: CRS | int | None, + nodata: int | float | None = None, + vcrs: Literal["Ellipsoid"] + | Literal["EGM08"] + | Literal["EGM96"] + | str + | pathlib.Path + | VerticalCRS + | int + | None = None, + ) -> DEM: + """Create a DEM from a numpy array and the georeferencing information. + + :param data: Input array. + :param transform: Affine 2D transform. Either a tuple(x_res, 0.0, top_left_x, + 0.0, y_res, top_left_y) or an affine.Affine object. + :param crs: Coordinate reference system. Either a rasterio CRS, + or an EPSG integer. + :param nodata: Nodata value. + :param vcrs: Vertical coordinate reference system. + + :returns: DEM created from the provided array and georeferencing. + """ + # We first apply the from_array of the parent class + rast = SatelliteImage.from_array(data=data, transform=transform, crs=crs, nodata=nodata) + # Then add the vcrs to the class call (that builds on top of the parent class) + return cls(filename_or_dataset=rast, vcrs=vcrs) @property - def ccrs(self) -> pyproj.CRS: - """Set compound CRS, i.e. horizontal and vertical references""" + def vcrs(self) -> VerticalCRS | Literal["Ellipsoid"] | None: + """Vertical coordinate reference system of the DEM.""" - # Temporary fix for some CRS with proj < 7.2 - def get_crs(filepath: str) -> pyproj.CRS: - """Get the CRS of a raster with the given filepath.""" - info = subprocess.run( - ["gdalinfo", "-json", filepath], stdout=subprocess.PIPE, check=True, encoding="utf-8" - ).stdout + return self._vcrs - wkt_string = json.loads(info)["coordinateSystem"]["wkt"] + @property + def vcrs_grid(self) -> str | None: + """Grid path of vertical coordinate reference system of the DEM.""" - return pyproj.CRS.from_wkt(wkt_string) + return self._vcrs_grid - # Temporary fix to get all types of CRS - if pyproj.proj_version_str >= "7.2.0": - crs = self.crs - else: - crs = get_crs(self.filename) - - if self.vref == "WGS84": - # The WGS84 ellipsoid corresponds to no vertical reference in pyproj - self._ccrs = pyproj.CRS(crs) - elif self.vref_grid is not None: - # For other vrefs, keep same horizontal projection and add geoid grid (the "dirty" way: because init is so - # practical and still going to be used for a while) - # see https://gis.stackexchange.com/questions/352277/ - with warnings.catch_warnings(): - warnings.filterwarnings("ignore", module="pyproj") - self._ccrs = pyproj.Proj(init="EPSG:" + str(int(crs.to_epsg())), geoidgrids=self.vref_grid).crs + @property + def vcrs_name(self) -> str | None: + """Name of vertical coordinate reference system of the DEM.""" + + if self.vcrs is not None: + # If it is the ellipsoid + if isinstance(self.vcrs, str): + # Need to call CRS() here to make it work with rasterio.CRS... + vcrs_name = f"Ellipsoid (No vertical CRS). Datum: {CRS(self.crs).ellipsoid.name}." + # Otherwise, return the vertical reference name + else: + vcrs_name = self.vcrs.name else: - self._ccrs = None + vcrs_name = None - return self._ccrs + return vcrs_name - def set_vref(self, vref_name: str | None = None, vref_grid: str | None = None) -> None: + def set_vcrs( + self, + new_vcrs: Literal["Ellipsoid"] | Literal["EGM08"] | Literal["EGM96"] | str | pathlib.Path | VerticalCRS | int, + ) -> None: """ - Set vertical reference with a name or with a grid. + Set the vertical coordinate reference system of the DEM. - :param vref_name: Vertical reference name - :param vref_grid: Vertical reference grid (any grid file in https://github.com/OSGeo/PROJ-data) - - :return: + :param new_vcrs: Vertical coordinate reference system either as a name ("Ellipsoid", "EGM08", "EGM96"), + an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data). """ - # Using vref_name only for WGS84 ellipsoid or the EGM96/EGM08 geoids (used 99% of the time) - if isinstance(vref_grid, str): - - # Default behaviour: use grid if both name and grid are provided - if isinstance(vref_name, str): - print("Both a vertical reference name and vertical grid are provided: defaulting to using grid only.") - - if vref_grid == "us_nga_egm08_25.tif": - self.vref = "EGM08" - self.vref_grid = vref_grid - elif vref_grid == "us_nga_egm96_15.tif": - self.vref = "EGM96" - self.vref_grid = vref_grid - else: - if os.path.exists(os.path.join(pyproj.datadir.get_data_dir(), vref_grid)): - self.vref = "Unknown vertical reference name from: " + vref_grid - self.vref_grid = vref_grid - else: - raise ValueError( - "Grid not found in " + str(pyproj.datadir.get_data_dir()) + ": check if proj-data is " - "installed via conda-forge, the pyproj.datadir, and that you are using a grid available at " - "https://github.com/OSGeo/PROJ-data" - ) + # Get vertical CRS and set it and the grid + self._vcrs = _vcrs_from_user_input(vcrs_input=new_vcrs) + self._vcrs_grid = _grid_from_user_input(vcrs_input=new_vcrs) - # Otherwise, use name provided - elif isinstance(vref_name, str): - if vref_name == "WGS84": - self.vref_grid = None - self.vref = "WGS84" # WGS84 ellipsoid - elif vref_name == "EGM08": - self.vref_grid = "us_nga_egm08_25.tif" # EGM2008 at 2.5 minute resolution - self.vref = "EGM08" - elif vref_name == "EGM96": - self.vref_grid = "us_nga_egm96_15.tif" # EGM1996 at 15 minute resolution - self.vref = "EGM96" - else: - raise ValueError( - 'Vertical reference name must be either "WGS84", "EGM96" or "EGM08". Otherwise, provide' - " a geoid grid from PROJ DATA: https://github.com/OSGeo/PROJ-data" - ) + @property + def ccrs(self) -> CompoundCRS | CRS | None: + """Compound horizontal and vertical coordinate reference system of the DEM.""" - # Else, return an error + if self.vcrs is not None: + ccrs = _build_ccrs_from_crs_and_vcrs(crs=self.crs, vcrs=self.vcrs) + return ccrs else: - raise ValueError("Vertical reference name or vertical grid must be a string") - - def to_vref(self, vref_name: str = "EGM96", vref_grid: str | None = None) -> None: + return None + def to_vcrs( + self, + dst_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int, + src_vcrs: Literal["Ellipsoid", "EGM08", "EGM96"] | str | pathlib.Path | VerticalCRS | int | None = None, + ) -> None: """ - Convert between vertical references: ellipsoidal heights or geoid grids. + Convert the DEM to another vertical coordinate reference system. - :param vref_name: Vertical reference name - :param vref_grid: Vertical reference grid (any grid file in https://github.com/OSGeo/PROJ-data) + :param dst_vcrs: Destination vertical CRS. Either as a name ("WGS84", "EGM08", "EGM96"), + an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data) + :param src_vcrs: Force a source vertical CRS (uses metadata by default). Same formats as for `dst_vcrs`. :return: """ - # All transformations grids file are described here: https://github.com/OSGeo/PROJ-data - if self.vref is None and self.vref_grid is None: + if self.vcrs is None and src_vcrs is None: raise ValueError( - "The current DEM has not vertical reference: need to set one before attempting a conversion " - "towards another vertical reference." + "The current DEM has no vertical reference, define one with .set_vref() " + "or by passing `src_vcrs` to perform a conversion." ) - # Initial ccrs - ccrs_init = self.ccrs + # Initial Compound CRS (only exists if vertical CRS is not None, as checked above) + if src_vcrs is not None: + # Warn if a vertical CRS already existed for that DEM + if self.vcrs is not None: + warnings.warn( + category=UserWarning, + message="Overriding the vertical CRS of the DEM with the one provided in `src_vcrs`.", + ) + src_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=src_vcrs) + else: + src_ccrs = self.ccrs - # Destination crs: first, set the new reference (before calculation doesn't change anything, - # we need to update the data manually anyway) - self.set_vref(vref_name=vref_name, vref_grid=vref_grid) - ccrs_dest = self.ccrs + # New destination Compound CRS + dst_ccrs = _build_ccrs_from_crs_and_vcrs(self.crs, vcrs=_vcrs_from_user_input(vcrs_input=dst_vcrs)) - # Transform the grid - transformer = Transformer.from_crs(ccrs_init, ccrs_dest) + # If both compound CCRS are equal, do not run any transform + if src_ccrs.equals(dst_ccrs): + warnings.warn( + message="Source and destination vertical CRS are the same, skipping vertical transformation.", + category=UserWarning, + ) + return None + + # Transform elevation with new vertical CRS zz = self.data xx, yy = self.coords(offset="center") - zz_trans = transformer.transform(xx, yy, zz[0, :])[2] - zz[0, :] = zz_trans + zz_trans = _transform_zz(crs_from=src_ccrs, crs_to=dst_ccrs, xx=xx, yy=yy, zz=zz) + + # Update DEM + self._data = zz_trans.astype(self.dtypes[0]) # type: ignore - # Update raster - self.data = zz + # Update vcrs (which will update ccrs if called) + self.set_vcrs(new_vcrs=dst_vcrs) diff --git a/xdem/vcrs.py b/xdem/vcrs.py new file mode 100644 index 00000000..ab83c6b3 --- /dev/null +++ b/xdem/vcrs.py @@ -0,0 +1,333 @@ +"""Routines for vertical CRS transformation (fully based on pyproj).""" +from __future__ import annotations + +import http.client +import os +import pathlib +import warnings +from typing import Literal, TypedDict + +import pyproj +from pyproj import CRS +from pyproj.crs import BoundCRS, CompoundCRS, GeographicCRS, VerticalCRS +from pyproj.crs.coordinate_system import Ellipsoidal3DCS +from pyproj.crs.enums import Ellipsoidal3DCSAxis +from pyproj.transformer import TransformerGroup + +from xdem._typing import MArrayf, NDArrayf + +# Sources for defining vertical references: +# AW3D30: https://www.eorc.jaxa.jp/ALOS/en/aw3d30/aw3d30v11_format_e.pdf +# SRTMGL1: https://lpdaac.usgs.gov/documents/179/SRTM_User_Guide_V3.pdf +# SRTMv4.1: http://www.cgiar-csi.org/data/srtm-90m-digital-elevation-database-v4-1 +# ASTGTM2/ASTGTM3: https://lpdaac.usgs.gov/documents/434/ASTGTM_User_Guide_V3.pdf +# NASADEM: https://lpdaac.usgs.gov/documents/592/NASADEM_User_Guide_V1.pdf, HGTS is ellipsoid, HGT is EGM96 geoid !! +# ArcticDEM (mosaic and strips): https://www.pgc.umn.edu/data/arcticdem/ +# REMA (mosaic and strips): https://www.pgc.umn.edu/data/rema/ +# TanDEM-X 90m global: https://geoservice.dlr.de/web/dataguide/tdm90/ +# COPERNICUS DEM: https://spacedata.copernicus.eu/web/cscda/dataset-details?articleId=394198 +vcrs_dem_products = { + "ArcticDEM/REMA/EarthDEM": "Ellipsoid", + "TDM1": "Ellipsoid", + "NASADEM-HGTS": "Ellipsoid", + "AW3D30": "EGM96", + "SRTMv4.1": "EGM96", + "ASTGTM2": "EGM96", + "ASTGTM3": "EGM96", + "NASADEM-HGT": "EGM96", + "COPDEM": "EGM08", +} + + +def _parse_vcrs_name_from_product(product: str) -> str | None: + """ + Parse vertical CRS name from DEM product name. + + :param product: Product name (typically from satimg.parse_metadata_from_fn). + + :return: vcrs_name: Vertical CRS name. + """ + + if product in vcrs_dem_products.keys(): + vcrs_name = vcrs_dem_products[product] + else: + vcrs_name = None + + return vcrs_name + + +def _build_ccrs_from_crs_and_vcrs(crs: CRS, vcrs: CRS | Literal["Ellipsoid"]) -> CompoundCRS | CRS: + """ + Build a compound CRS from a horizontal CRS and a vertical CRS. + + :param crs: Horizontal CRS. + :param vcrs: Vertical CRS. + + :return: Compound CRS (horizontal + vertical). + """ + + # If a vertical CRS was passed, build a compound CRS with horizontal + vertical + # This requires transforming the horizontal CRS to 2D in case it was 3D + # Using CRS() because rasterio.CRS does not allow to call .name otherwise... + if isinstance(vcrs, CRS): + # If pyproj >= 3.5.1, we can use CRS.to_2d() + from packaging.version import Version + + if Version(pyproj.__version__) > Version("3.5.0"): + crs_from = CRS(crs).to_2d() + ccrs = CompoundCRS( + name="Horizontal: " + CRS(crs).name + "; Vertical: " + vcrs.name, + components=[crs_from, vcrs], + ) + # Otherwise, we have to raise an error if the horizontal CRS is already 3D + else: + crs_from = CRS(crs) + # If 3D + if len(crs_from.axis_info) > 2: + raise NotImplementedError( + "pyproj >= 3.5.1 is required to demote a 3D CRS to 2D and be able to compound " + "with a new vertical CRS. Update your dependencies or pass the 2D source CRS " + "manually." + ) + # If 2D + else: + ccrs = CompoundCRS( + name="Horizontal: " + CRS(crs).name + "; Vertical: " + vcrs.name, + components=[crs_from, vcrs], + ) + + # Else if "Ellipsoid" was passed, there is no vertical reference + # We still have to return the CRS in 3D + elif isinstance(vcrs, str) and vcrs.lower() == "ellipsoid": + ccrs = CRS(crs).to_3d() + else: + raise ValueError("Invalid vcrs given. Must be a vertical CRS or the literal string 'Ellipsoid'.") + + return ccrs + + +def _build_vcrs_from_grid(grid: str, old_way: bool = False) -> CompoundCRS: + """ + Build a compound CRS from a vertical CRS grid path. + + :param grid: Path to grid for vertical reference. + :param old_way: Whether to use the new or old way of building the compound CRS with pyproj (for testing purposes). + + :return: Compound CRS (horizontal + vertical). + """ + + if not os.path.exists(os.path.join(pyproj.datadir.get_data_dir(), grid)): + warnings.warn( + "Grid not found in " + + str(pyproj.datadir.get_data_dir()) + + ". Attempting to download from https://cdn.proj.org/..." + ) + from pyproj.sync import _download_resource_file + + try: + _download_resource_file( + file_url=os.path.join("https://cdn.proj.org/", grid), + short_name=grid, + directory=pyproj.datadir.get_data_dir(), + verbose=False, + ) + except http.client.InvalidURL: + raise ValueError( + "The provided grid '{}' does not exist at https://cdn.proj.org/. " + "Provide an existing grid.".format(grid) + ) + + # The old way: see https://gis.stackexchange.com/questions/352277/. + if old_way: + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", module="pyproj") + ccrs = pyproj.Proj(init="EPSG:4326", geoidgrids=grid).crs + bound_crs = ccrs.sub_crs_list[1] + + # The clean way + else: + # First, we build a bounds CRS (the vertical CRS relative to geographic) + vertical_crs = VerticalCRS( + name="unknown using geoidgrids=" + grid, datum='VDATUM["unknown using geoidgrids=' + grid + '"]' + ) + geographic3d_crs = GeographicCRS( + name="WGS 84", + ellipsoidal_cs=Ellipsoidal3DCS(axis=Ellipsoidal3DCSAxis.LATITUDE_LONGITUDE_HEIGHT), + ) + bound_crs = BoundCRS( + source_crs=vertical_crs, + target_crs=geographic3d_crs, + transformation={ + "$schema": "https://proj.org/schemas/v0.2/projjson.schema.json", + "type": "Transformation", + "name": "unknown to WGS84 ellipsoidal height", + "source_crs": vertical_crs.to_json_dict(), + "target_crs": geographic3d_crs.to_json_dict(), + "method": {"name": "GravityRelatedHeight to Geographic3D"}, + "parameters": [ + { + "name": "Geoid (height correction) model file", + "value": grid, + "id": {"authority": "EPSG", "code": 8666}, + } + ], + }, + ) + + return bound_crs + + +# Define types of common Vertical CRS dictionary +class VCRSMetaDict(TypedDict, total=False): + grid: str + epsg: int + + +_vcrs_meta: dict[str, VCRSMetaDict] = { + "EGM08": {"grid": "us_nga_egm08_25.tif", "epsg": 3855}, # EGM2008 at 2.5 minute resolution + "EGM96": {"grid": "us_nga_egm96_15.tif", "epsg": 5773}, # EGM1996 at 15 minute resolution +} + + +def _vcrs_from_crs(crs: CRS) -> CRS: + """Get the vertical CRS from a CRS.""" + + # Check if CRS is 3D + if len(crs.axis_info) > 2: + + # Check if CRS has a vertical compound + if any(subcrs.is_vertical for subcrs in crs.sub_crs_list): + # Then we get the first vertical CRS (should be only one anyway) + vcrs = [subcrs for subcrs in crs.sub_crs_list if subcrs.is_vertical][0] + # Otherwise, it's a 3D CRS based on an ellipsoid + else: + vcrs = "Ellipsoid" + # Otherwise, the CRS is 2D and there is no vertical CRS + else: + vcrs = None + + return vcrs + + +def _vcrs_from_user_input( + vcrs_input: Literal["Ellipsoid"] | Literal["EGM08"] | Literal["EGM96"] | str | pathlib.Path | CRS | int, +) -> VerticalCRS | BoundCRS | Literal["Ellipsoid"]: + """ + Parse vertical CRS from user input. + + :param vcrs_input: Vertical coordinate reference system either as a name ("Ellipsoid", "EGM08", "EGM96"), + an EPSG code or pyproj.crs.VerticalCRS, or a path to a PROJ grid file (https://github.com/OSGeo/PROJ-data). + + :return: Vertical CRS. + """ + + # Raise errors if input type is wrong (allow CRS instead of VerticalCRS for broader error messages below) + if not isinstance(vcrs_input, (str, pathlib.Path, CRS, int)): + raise TypeError(f"New vertical CRS must be a string, path or VerticalCRS, received {type(vcrs_input)}.") + + # If input is ellipsoid + if ( + (isinstance(vcrs_input, str) and (vcrs_input.lower() == "ellipsoid" or vcrs_input.upper() == "WGS84")) + or (isinstance(vcrs_input, int) and vcrs_input in [4326, 4979]) + or (isinstance(vcrs_input, CRS) and vcrs_input.to_epsg() in [4326, 4979]) + ): + return "Ellipsoid" + + # Define CRS in case EPSG or CRS was passed + if isinstance(vcrs_input, (int, CRS)): + if isinstance(vcrs_input, int): + vcrs = CRS.from_epsg(vcrs_input) + else: + vcrs = vcrs_input + + # Raise errors if the CRS constructed is not vertical or has other components + if isinstance(vcrs, CRS) and not vcrs.is_vertical: + raise ValueError( + "New vertical CRS must have a vertical axis, '{}' does not " + "(check with `CRS.is_vertical`).".format(vcrs.name) + ) + elif isinstance(vcrs, CRS) and vcrs.is_vertical and len(vcrs.axis_info) > 2: + warnings.warn( + "New vertical CRS has a vertical dimension but also other components, " + "extracting the vertical reference only." + ) + vcrs = _vcrs_from_crs(vcrs) + + # If a string was passed + else: + # If a name is passed, define CRS based on dict + if isinstance(vcrs_input, str) and vcrs_input.upper() in _vcrs_meta.keys(): + vcrs_meta = _vcrs_meta[vcrs_input] + vcrs = CRS.from_epsg(vcrs_meta["epsg"]) + # Otherwise, attempt to read a grid from the string + else: + if isinstance(vcrs_input, pathlib.Path): + grid = vcrs_input.name + else: + grid = vcrs_input + vcrs = _build_vcrs_from_grid(grid=grid) + + return vcrs + + +def _grid_from_user_input(vcrs_input: str | pathlib.Path | int | CRS) -> str | None: + + # If a grid or name was passed, get grid name + if isinstance(vcrs_input, (str, pathlib.Path)): + # If the string is within the supported names + if isinstance(vcrs_input, str) and vcrs_input in _vcrs_meta.keys(): + grid = _vcrs_meta[vcrs_input]["grid"] + # If it's a pathlib path + elif isinstance(vcrs_input, pathlib.Path): + grid = vcrs_input.name + # Or an ellipsoid + elif vcrs_input.lower() == "ellipsoid": + grid = None + # Or a string path + else: + grid = vcrs_input + # Otherwise, there is none + else: + grid = None + + return grid + + +def _transform_zz( + crs_from: CRS, crs_to: CRS, xx: NDArrayf, yy: NDArrayf, zz: MArrayf | NDArrayf | int | float +) -> MArrayf | NDArrayf | int | float: + """ + Transform elevation to a new 3D CRS. + + :param crs_from: Source CRS. + :param crs_to: Destination CRS. + :param xx: X coordinates. + :param yy: Y coordinates. + :param zz: Z coordinates. + + :return: Transformed Z coordinates. + """ + + # Find all possible transforms + with warnings.catch_warnings(): + warnings.filterwarnings("ignore", "Best transformation is not available") + trans_group = TransformerGroup(crs_from=crs_from, crs_to=crs_to, always_xy=True) + + # Download grid if best available is not on disk, download and re-initiate the object + if not trans_group.best_available: + trans_group.download_grids() + trans_group = TransformerGroup(crs_from=crs_from, crs_to=crs_to, always_xy=True) + + # If the best available grid is still not there, raise a warning + if not trans_group.best_available: + warnings.warn( + category=UserWarning, + message="Best available grid for transformation could not be downloaded, " + "applying the next best available (caution: might apply no transform at all).", + ) + transformer = trans_group.transformers[0] + + # Will preserve the mask of the masked-array since pyproj 3.4 + zz_trans = transformer.transform(xx, yy, zz)[2] + + return zz_trans