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 gravityrelated, 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), ora
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 thevcrs
argument),The setting method
set_vcrs()
to explicitly set the vertical CRS of aDEM
,The transformation method
to_vcrs()
to explicitly transform the vertical CRS of aDEM
.
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 3dimensional,
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 [ellipsoidalvertical]:
 Lat[north]: Geodetic latitude (degree)
 Lon[east]: Geodetic longitude (degree)
 [up]: Gravityrelated 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]: Gravityrelated 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 bygeoutils.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 

Ellipsoid 

Ellipsoid 

Ellipsoid 

Ellipsoid 

Ellipsoid 

EGM96 

EGM96 

EGM96 

EGM96 

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 reset 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]: Gravityrelated 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/stable/xdem/vcrs.py:120: 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]: Gravityrelated 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]: Gravityrelated 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]: Gravityrelated 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 areaofinterest 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")
trans_dem = dem.to_vcrs("EGM96")
trans_dem.vcrs
<Vertical CRS: EPSG:5773>
Name: EGM96 height
Axis Info [vertical]:
 H[up]: Gravityrelated 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 inplace. 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)")