Common Pitfall in Elevation differencing
Contents
Common Pitfall in Elevation differencing¶
Data from different sources¶
The lidar data we have used so far are acquired from the same source (QSI). This made the DEM differencing to be quiet straightforward because the crs, resolution and extent are the same. At times, the elevation data might be from diffeerence sources. It will thus require extra effort of resampling and reprojection to accurately derive snow depth. Let’s do DEM differencing again but using snow-on and snow-off DEM from different sources. The snow-on data is an 0.5m resolution raster acquired by QSI over GrandMesa in 2020. A 3m DEM data acquired by ASO in 2016 will be used as the snow-off DEM.
As a first step, we have to downsample the QSI data to the resolution of the ASO.
import os
import numpy as np
#plotting packages
import matplotlib.pyplot as plt
import seaborn as sns
import hvplot.xarray
#geospatial packages
import geopandas as gpd #for vector data
import xarray as xr
import rioxarray #for raster data
import rasterstats as rs
%%bash
# Retrieve a copy of data files used in this tutorial from Zenodo.org:
# Re-running this cell will not re-download things if they already exist
cd /tmp
wget -q -nc -O input.zip "https://github.com/snowex-hackweek/tutorial-data/releases/download/SnowEx_2022_lidar/input.zip"
unzip -q -n input.zip
!gdal_translate -tr 3, 3 /tmp/input/QSI_0.5m_GrandMesa13Feb_Winter_DEM.tif ./output/QSI_3m_GrandMesa13Feb_Winter_DEM.tif
#check the raster metadata
!gdalinfo ./output/QSI_3m_GrandMesa13Feb_Winter_DEM.tif
Driver: GTiff/GeoTIFF
Files: ./output/QSI_3m_GrandMesa13Feb_Winter_DEM.tif
Size is 7621, 3574
Coordinate System is:
PROJCRS["unnamed",
BASEGEOGCRS["Unknown datum based upon the GRS 1980 ellipsoid",
DATUM["Not specified (based on GRS 1980 ellipsoid)",
ELLIPSOID["GRS 1980",6378137,298.257222101004,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4019]],
CONVERSION["Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-111,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",0.9996,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",500000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (737387.500000000000000,4330284.500000000000000)
Pixel Size = (3.000000000000000,-3.000000000000000)
Metadata:
AREA_OR_POINT=Area
Image Structure Metadata:
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 737387.500, 4330284.500) (108d15'19.09"W, 39d 5'21.85"N)
Lower Left ( 737387.500, 4319562.500) (108d15'32.54"W, 38d59'34.42"N)
Upper Right ( 760250.500, 4330284.500) (107d59'28.62"W, 39d 4'58.38"N)
Lower Right ( 760250.500, 4319562.500) (107d59'43.36"W, 38d59'11.03"N)
Center ( 748819.000, 4324923.500) (108d 7'30.87"W, 39d 2'16.69"N)
Band 1 Block=7621x1 Type=Float32, ColorInterp=Gray
NoData Value=-3.4028235e+38
#read the data as arrays
snow_on_dem = rioxarray.open_rasterio('./output/QSI_3m_GrandMesa13Feb_Winter_DEM.tif', masked =True)
snow_on_dem.hvplot.image(x= 'x', y = 'y', rasterize = True, cmap = 'viridis')
#read the snow-off DEMs
snow_free_dem = rioxarray.open_rasterio('/tmp/input/ASO_3M_GrandMesa_Summer_DEM.tif',
masked=True)
snow_free_dem.hvplot.image(x= 'x', y = 'y', rasterize = True, cmap = 'viridis')
The extent and crs of the snow-off DEM is different from the snow-on. Let’s match the snow-off DEM with the snow-on DEM
snow_free_dem_match = snow_free_dem.rio.reproject_match(snow_on_dem)
snow_free_dem_match.hvplot.image(x= 'x', y = 'y', rasterize = True, cmap = 'viridis')
snow_depth = snow_on_dem - snow_free_dem_match
snow_depth.hvplot.image(x= 'x', y = 'y', rasterize = True, cmap = 'viridis')
About 17 m snow depth…. What’s the issue?
!gdalsrsinfo -o proj4 /tmp/input/ASO_3M_GrandMesa_Summer_DEM.tif
+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs
!gdalsrsinfo -o proj4 ./output/QSI_3m_GrandMesa13Feb_Winter_DEM.tif
+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs
The snow-on and snow-free DEM are referenced to different datum, so we need to do vertical datum transformation. Matching the vertical datum of data from different sources is important for elevation data.
The vertical datum of ASO is in ellipsoid (WGS 84). QSI is in UTM Zone 12 North. Horizontal Datum used is NAD83 (2011). NAD 83 is a 3-D reference system, it also serves as the horizontal control datum for the United States, Canada, Greenland and Central America using the Geodetic Reference System 1980 (GRS 80) ellipsoid, an earth-centered reference ellipsoid. See this link for more reading. The vertical datum is NAVD88 (Geoid12b).
Let’s transform ASO to same datum as QSI.
ls /tmp/input/egm96_15.gtx
/tmp/input/egm96_15.gtx
#coordinate transformation
!gdalwarp -s_srs "+proj=utm +zone=13 +datum=WGS84 +units=m +no_defs" -t_srs "+proj=utm +zone=12 +ellps=GRS80 +units=m +no_defs +geoidgrids=/tmp/input/egm96_15.gtx" /tmp/input/ASO_3M_GrandMesa_Summer_DEM.tif ./output/ASO_3M_GrandMesa_Summer_DEM_geoid_96.tif
#read the snow-off DEMs
snow_free_dem = rioxarray.open_rasterio('./output/ASO_3M_GrandMesa_Summer_DEM_geoid_96.tif',
masked=True)
#reproject the ASO DEM to match the QSI DEM
snow_free_dem_match = snow_free_dem.rio.reproject_match(snow_on_dem)
snow_depth = snow_on_dem - snow_free_dem_match
snow_depth = snow_depth.where((snow_depth > 0) & (snow_depth < 3))
snow_depth.hvplot.image(x= 'x', y = 'y', rasterize = True, cmap = 'viridis')
Important notes on vertical datum transformation¶
GDAL uses the PROJ library to carry out CRS transformations. PROJ uses pre-defined grids for datum transformations. Vertical datum transformations are defined using a GTX file. In the example above, when we did the transformation to EGM96 vertical CRS, it used the parameters supplied in the egm96_15.gtx file. Many such grid files are included when you download and install GDAL. But there are other grids which are very large and are not included by default. If you want to convert to a vertical datum whose grid files are not included by default, you need to download them separately and copy them to appropriate directory on your system. Learn more about PROJ Resource files.
One such vertical CRS is EGM2008 (EPSG:3855). Say you wanted to transform WGS84 heights to this new and more precise geoid model. You can run a command such as below
gdalwarp -s_srs "+proj=longlat +datum=WGS84 +no_def" -t_srs "+proj=longlat +datum=WGS84 +no_defs +geoidgrids=egm08_25.gtx" /tmp/input/ASO_3M_GrandMesa_Summer_DEM.tif ./output/ASO_3M_GrandMesa_Summer_DEM_geoid_08.tif
But if you run it, you may get an error such as below
ERROR 1: Cannot open egm08_25.gtx.
This is because the EGM2008 grid is very large and is not included in many GDAL installations. To fix this, you can download the grid file from http://download.osgeo.org/proj/vdatum/
Copy the egm08_25.gtx file to the PROJ resource directory on your computer. The location of this directory will depend on your platform and the installation method. Some common paths are as below
Windows: C:\OsGeo4W64\share\proj\
Mac: /Library/Frameworks/PROJ.framework/Resources/proj/
Linux: /usr/share/proj/
Coregistion¶
Another potential source of error for lidar differencing is coregristration. If the differenced value still doesn’t make sense after matching the crs and extent, it might be a misregistration issue.
Available Data for Lidar Remote Sensing¶
Various lidar data has been acquired to monitory snow over the western US.
2013 - 2019: ASO data (Available on NSIDC)
2020 ASO data: Here are some .zip files to access ASO 2020 data. The folder contains snow depth and SWE data products
2021 Quantumn data over selected areas. The lidar data for 2020 and 2021 covers the location shown in the map below and include:
DEMs (SnowOn and SnowFree)
DSMs (SnowOn and SnowFree)
Vegetation Height (SnowOn and SnowFree)
Snow Depth
Intensity Images (SnowOn and SnowFree)
2022: Riegl flight over MCS
References:
Earthdatascience
Spatial thought. https://spatialthoughts.com/2019/10/26/convert-between-orthometric-and-ellipsoidal-elevations-using-gdal/