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 DEM class.

Important

A DEM already possesses a crs attribute that defines its 2- or 3D CRS, inherited from 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 Automated vertical CRS detection 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 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 pyproj.CRS with only a vertical axis for a CRS based on geoid heights (e.g., the EGM96 geoid).

In practice, a pyproj.CRS with only a vertical axis is either a BoundCRS (when created from a grid) or a 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 DEM that implicitly tries to set the vertical CRS from the metadata (or explicitly through the vcrs argument),
The setting method set_vcrs() to explicitly set the vertical CRS of a DEM,
The transformation method to_vcrs() to explicitly transform the vertical CRS of a DEM.

Caution

As of now, Rasterio 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 DEM.reproject and DEM.to_vcrs, respectively.

Automated vertical CRS detection#

During instantiation of a DEM, the vertical CRS vcrs is tentatively set with the following priority order:

From the crs of the DEM, if 3-dimensional,

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)

<Compound CRS: EPSG:9707>
Name: WGS 84 + EGM96 height
Axis Info [ellipsoidal|vertical]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
- [up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich
Sub CRS:
- WGS 84
- EGM96 height

# The vertical CRS is extracted automatically
dem.vcrs

<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
- [up]: Gravity-related height (metre)
Area of Use:
- undefined
Datum: EGM96 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

From the product name of the DEM, if parsed from the filename by geoutils.SatelliteImage.

# 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

'Ellipsoid'

Currently recognized DEM products:

DEM	Vertical CRS
ArcticDEM	Ellipsoid
REMA	Ellipsoid
EarthDEM	Ellipsoid
TanDEM-X global DEM	Ellipsoid
NASADEM-HGTS	Ellipsoid
NASADEM-HGT	EGM96
ALOS World 3D	EGM96
SRTM v4.1	EGM96
ASTER GDEM v2-3	EGM96
Copernicus DEM	EGM08

If your DEM does not have a .vcrs after instantiation, none of those steps worked. You can define the vertical CRS explicitly during DEM instantiation with the vcrs argument or with set_vcrs(), with user inputs described below.

Setting a vertical CRS with convenient user inputs#

The vertical CRS of a DEM can be set or re-set manually at any point using set_vcrs().

The vcrs argument, consistent across the three functions DEM, set_vcrs() and 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

<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: EGM96 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

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

'Ellipsoid'

Any PROJ grid name available at 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.

# Set a geoid vertical CRS based on a grid
dem.set_vcrs("us_noaa_geoid06_ak.tif")
dem.vcrs

/home/docs/checkouts/readthedocs.org/user_builds/xdem/checkouts/latest/xdem/vcrs.py:120: UserWarning: Grid not found in /home/docs/checkouts/readthedocs.org/user_builds/xdem/conda/latest/share/proj. Attempting to download from https://cdn.proj.org/...
  warnings.warn(

<Bound CRS: {"$schema": "https://proj.org/schemas/v0.2/projjso ...>
Name: unknown using geoidgrids=us_noaa_geoid06_ak.tif
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- undefined
Coordinate Operation:
- name: unknown to WGS84 ellipsoidal height
- method: GravityRelatedHeight to Geographic3D
Datum: None
- Ellipsoid: undefined
- Prime Meridian: undefined
Source CRS: unknown using geoidgrids=us_noaa_geoid06_ak.tif

Any EPSG code as int,

# Set a geoid vertical CRS based on an EPSG code
dem.set_vcrs(5773)
dem.vcrs

<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: EGM96 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

Any 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

<Vertical CRS: EPSG:3855>
Name: EGM2008 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: EGM2008 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

Transforming to another vertical CRS#

To transform a DEM to a different vertical CRS, 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 pyproj.transformer.TransformerGroup, considering the area-of-interest as the DEM extent 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

<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
- H[up]: Gravity-related height (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: EGM96 geoid
- Ellipsoid: undefined
- Prime Meridian: undefined

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)

32.087208