CrevDEM

Tools for extracting crevasse location and volume from high-resolution digital elevation models.

Note Our paper using CrevDEM to assess crevasse change across the Greenland Ice Sheet is now in review. Find the preprint on EarthArxiv.

CrevDEM is a Python package for extracting crevasse location and volume from high-resolution Digital Elevation Models (DEMs) of glaciers and ice sheets provided by the ArcticDEM and REMA projects.

Warning These tools are designed for scientific purposes only. ArcticDEM and REMA strips are not suitable for detecting metre-scale and snow-covered crevasses, and should not be used for field hazard assessment.

The functions within provide a complete workflow for:

1. Automatically downloading and clipping 2 m strips of the version 4.1 ArcticDEM (Porter et al. 2022) and Reference Elevation Model of Antarctica (REMA; Howat et al. 2022).
2. DEM geoid correction and filtering of bare rock and proglacial sea-ice/mélange.
3. Extracting crevasse presence following the method outlined by Chudley et al. (in prep).

The principle of crevasse extraction is based around Black Top Hat filtering of a detrended surface (Kodde et al. 2007). The optimal kernel size can be determined quantiatively through the use of variogram analysis: a notebook is provided to aid with this.

Cite

Please cite the source paper when using CrevDEM:

Chudley, T. R., Howat, I. M., King, M. D., and MacKie, E. (in review) An increase in crevasses across accelerating Greenland Ice Sheet margins. Preprint doi:10.31223/X58099

As always when using ArcticDEM and REMA products, please cite the datasets appropriately and acknowledge the PGC.

Installation

Install CrevDEM

After downloading, crevdem can be installed from the top-level directory via pip install .:

git clone https://github.com/trchudley/crevdem
cd crevdem
pip install .

CrevDEM has the following dependencies. When installing into a conda environment, it may be beneficial to install these before the pip install.

• rioxarray
• Rasterio
• Shapely
• NumPy
• OpenCV

The variogram analysis notebook requires additional packages, including scikit-gstat. If you would rather not install these yourself using conda or similar, you can use pip install .[variogram] during the initial install.

Install supplementary datasets

Supplementary datasets are required to be available locally to complete geoid correction and filtering of non-glacial regions: specifically, the BedMachine Greenland v5 or BedMachine Antarctica v3 respectively (Morlighem et al. 2022a, 2022b) and, for Greenland, the GrIMP ice mask (Howat, 2017) for bedrock filtering. File or directory paths will be requested in the relevant functions.

These can be downloaded from the NSIDC manually (Greenland BedMachine, Antarctic BedMachine, Greenland surface mask) but for conveninence, download scripts are provided in the supp_data directory. The BedMachine download scripts are provided by the NSIDC, and require an Earthdata user account and password to be provided. The files will be downloaded into the directory the scripts are run.

cd supp_data
python download_bedmachine_greenland_v5.py
python download_grimp_2015_15m.py
python download_bedmachine_antarctica_v3.py

Usage

The sections below briefly outline the purpose of user-exposed functions available through the package. In order to see them in action, Jupyter Notebooks are provided in the ./notebooks directory in order to provide an introduction into the use of CrevDEM.

Information on the required and optional input variable for individual functions can be accessed through Python's help() function, e.g. help(crevdem.load_aws).

Loading DEM strips from the cloud

Note An updated and expanded version of these tools has been released as an independent Python package for searching, download, and preprocessing of both ArticDEM and REMA strips. See pdemtools for more information..

load_aws() - Returns the selected ArcticDEM/REMA strip, downloaded from the relevant AWS bucket, as an xarray DataArray suitable for further processing by crevdem. Option to filter to bounds and bitmask. 2 m DEM strips are large in size and loading remotely from AWS may take some time.

load_local() - Loads the desired ArcticDEM/REMA DEM strip, from local filepaths, as an xarray DataArray suitable for further processing by crevdem. Option to filter to bounds and bitmask.

Preprocessing the DEM strips

mask_bedrock() - Returns a bedrock-masked DEM. Can either provide your own mask (as a DataArray) using the mask variable (where bedrock = 0/False and ice/ocean = 1/True), or provide the path to a directory containing the GrIMP 15 m classification mask using the grimp_mask_dir variable.

geoid_correct() - Returns a geoid-corrected DEM. Can provide either your own geoid (as a DataArray) using the geoid variable, or the filepath to an appropriate BedMachine dataset using the bedmachine_fpath variable.

mask_melange() - Returns a DEM with mélange/ocean region, as identified by get_melange_mask() function, filtered out. If no likely sea level is identified, returns the original DEM. DEM must be geoid-corrected.

get_melange_mask() - Returns a mask of mélange/ocean regions of a DEM, using sea level as returned by the get_sea_level() function. Input DEM must be geoid-corrected. In returned mask, land/ice is True and ocean is False.

get_sea_level() - Returns estimated sea level following method of Shiggins et al. (2023). If no candidate sea level is identified, None is returned. Input DEM must be geoid-corrected.

Extracting crevasse presence

find() - Returns crevasse depths, batch processed from input DEM strip. Parameters default to Chudley et al. generic workflow for Greenland marine margins, but can be modified. This function is a wrapper for the detrend, bth_filter, threshold_depth, interpolate_surface, and calc_depth functions.

• detrend() - Returns a detrended DEM DataArray using a large gaussian filter. Standard deviation size should be >> the features of interest (in the default crevdem settings, the gauss_std to be 3* the range).

• bth() - Returns a black-top-hat-filtered DEM DataArray from the (detrended) DEM DataArray. Kernel diameter is set following the range distance.

• threshold_depth() - Returns crevasse mask (crevasse = 1; not crevasse = 0) DataArray from BTH-filtered DataArray. Mask is filtered to the threshold BTH value, which is set to 1 metre in the default workflow.

• interpolate_surface() - Returns a 'crevasse-filled' DEM from the original DEM and crevasse mask, using the GDAL FillNodata algorithm (inverse distance weighting) to fill crevasse-masked regions. Smoothing iterations are applied to smooth out artefacts.

• calc_depth() - Returns final crevasse depth, calculated from the raw DEM and the filled DEM.

Improvements

The tool is presented as-is, but requests/contributions to functionality are welcome (thomas.r.chudley@durham.ac.uk). Avenues for future work include the following:

• Currently, the only automated input mask is the GrIMP mask. However, it would be relatively trivial to include an automated BedMachine mask option for use in Antarctica (and Greenland). The resolution is large (150 m), so some automated interpolation would probably be necessary.
• The default PGC bitmask has a habit of leaving in some cloud blunders. Ian Howat has implemented an additional filter (in Matlab) to catch these remnant cloud blunders in the latest ArcticDEM and REMA mosaic tools. This could be rewritten in Python to filter out these blunders.
• The currently chosen Guassian and BTH filters (as implemented in OpenCV) return Nan if NaNs are present within the kernel, leading to loss of data at glacier margins. It is worth searching for or developing alternative implementations to account for NaN values. astropy already has such a function for Gaussian filters but is very slow over large datasets. A test numba.stencil implementation was also slow (>6 m vs <1 s for OpenCV).

References

Chudley, T. R., et al. (in review). An increase in crevasses across accelerating Greenland Ice Sheet margins. Preprint: https://doi.org/10.31223/X58099

Howat, I. (2017). MEaSUREs Greenland Ice Mapping Project (GIMP) Land Ice and Ocean Classification Mask, Version 1 [Data Set]. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/B8X58MQBFUPA

Howat, I., et al. (2022). The Reference Elevation Model of Antarctica – Strips, Version 4.1. Harvard Dataverse https://doi.org/10.7910/DVN/X7NDNY

Kodde, M. P., et al. (2007). Automatic glacier surface analysis from airborne laser scanning. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 36(3), 221–226.

Morlighem, M. et al. (2022). MEaSUREs BedMachine Antarctica, Version 3 [Data Set]. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/FPSU0V1MWUB6

Morlighem, M. et al. (2022). IceBridge BedMachine Greenland, Version 5 [Data Set]. NASA National Snow and Ice Data Center Distributed Active Archive Center. https://doi.org/10.5067/GMEVBWFLWA7X

Porter, C., et al. (2022). ArcticDEM - Strips, Version 4.1. Harvard Dataverse. https://doi.org/10.7910/DVN/OHHUKH

Shiggins, et al. (2023). Automated ArcticDEM iceberg detection tool: insights into area and volume distributions, and their potential application to satellite imagery and modelling of glacier–iceberg–ocean systems, The Cryosphere, 17, 15–32, https://doi.org/10.5194/tc-17-15-2023

Acknowledgements

ArcticDEM: DEMs are provided by the Polar Geospatial Center under NSF-OPP awards 1043681, 1559691, and 1542736.

REMA: DEMs are provided by the Byrd Polar and Climate Research Center and the Polar Geospatial Center under NSF-OPP awards 1543501, 1810976, 1542736, 1559691, 1043681, 1541332, 0753663, 1548562, 1238993 and NASA award NNX10AN61G. Computer time provided through a Blue Waters Innovation Initiative. DEMs produced using data from Maxar.