Skip to content

Latest commit

 

History

History
240 lines (186 loc) · 8.77 KB

vertical_ref.md

File metadata and controls

240 lines (186 loc) · 8.77 KB
file_format jupytext kernelspec
mystnb
formats text_representation
md:myst
extension format_name
.md
myst
display_name language name
xdem-env
python
xdem

(vertical-ref)=

Vertical referencing

xDEM supports the use of vertical coordinate reference systems (vertical CRSs) and vertical transformations for DEMs by conveniently wrapping PROJ pipelines through Pyproj in the {class}~xdem.DEM class.

**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 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<xdem.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<pyproj.crs.CRS> with only a vertical axis for a CRS based on geoid heights (e.g., the EGM96 geoid).

In practice, a {class}pyproj.CRS<pyproj.crs.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.
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<xdem.DEM.reproject>` and {func}`DEM.to_vcrs<xdem.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,
: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")

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)

# The vertical CRS is extracted automatically
dem.vcrs

:tags: [remove-cell]

import os
os.remove("mydem_with3dcrs.tif")
  1. From the {attr}~xdem.DEM.product name of the DEM, if parsed from the filename by {class}geoutils.SatelliteImage.
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).

: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")

# 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

:tags: [remove-cell]

os.remove("SETSM_WV03_20151101_104001001327F500_104001001312DE00_seg2_2m_v3.0_dem.tif")

Currently recognized DEM products:

   :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",
# Set a geoid vertical CRS based on a string
dem.set_vcrs("EGM96")
dem.vcrs

# Set a vertical CRS extended from the ellipsoid of the DEM's CRS
dem.set_vcrs("Ellipsoid")
dem.vcrs

**No need to download the grid!** This is done automatically during the setting operation, if the grid does not already exist locally.

# 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,
# 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.
# 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.

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`.

# 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:

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)