Vertical referencing

Vertical referencing#

xDEM supports the use of vertical coordinate reference systems (vertical CRSs) and vertical transformations for elevation data by conveniently wrapping PROJ pipelines through Pyproj in the DEM class.

Note

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.

For more reading on referencing for elevation data, see the Georeferencing intricacies guide page.

Quick use#

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.

To pass a vertical CRS argument, xDEM accepts string of the most commonly used ("EGM96", "EGM08" and "Ellipsoid"), any pyproj.crs.CRS objects and any PROJ grid name (available at https://cdn.proj.org/) which is automatically downloaded.

Hide the code for opening example data 
import xdem
import matplotlib.pyplot as plt

ref_dem = xdem.DEM(xdem.examples.get_path("longyearbyen_ref_dem"))

# Set current vertical CRS
ref_dem.set_vcrs("EGM96")
# Transform to a local reference system from https://cdn.proj.org/
trans_dem = ref_dem.to_vcrs("no_kv_arcgp-2006-sk.tif")

# Plot the elevation differences of the vertical transformation
(trans_dem - ref_dem).plot(cmap='RdYlBu', cbar_title="Elevation differences of\n vertical transform (m)")

/home/docs/checkouts/readthedocs.org/user_builds/xdem/checkouts/stable/xdem/vcrs.py:138: UserWarning: Grid not found in /home/docs/checkouts/readthedocs.org/user_builds/xdem/conda/stable/share/proj. Attempting to download from https://cdn.proj.org/...
  warnings.warn(
_images/3f5e8c6562bca4b7956a23c35a133f769a3261eb5e40db750eaa3a1df06db8f2.png

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

Automated vertical CRS detection#

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

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

  1. From the product name of the DEM, if parsed from the filename by the parse_sensor_metadata of geoutils.Raster.

See also

The Raster parent class can parse sensor metadata, see its documentation page.

# Open an ArcticDEM strip file, recognized as an ArcticDEM product
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

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'

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/stable/xdem/vcrs.py:138: UserWarning: Grid not found in /home/docs/checkouts/readthedocs.org/user_builds/xdem/conda/stable/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.

# Open a DEM and set its CRS
filename_dem = xdem.examples.get_path("longyearbyen_ref_dem")
dem = xdem.DEM(filename_dem, vcrs="Ellipsoid")
trans_dem = dem.to_vcrs("EGM96")
trans_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 returns a new DEM by default, but can also be done in-place. It vertically shifts each pixel value by the transformation at their coordinates:

import numpy as np

diff = trans_dem - dem
# Mean difference after transformation (about 30 m in Svalbard)
np.nanmean(diff)

-32.087208

# Plot the elevation differences due to the vertical transformation
diff.plot(cmap="RdYlBu", cbar_title="Elevation differences (m)")
_images/1d385861ca70b9d2486df1178219f329f1b19ec47f024360057f2281f7d41415.png