Why pDEMtools?#
This page reproduces the ‘Statement of Need’ of the Journal of Open Source Software manuscript for
pdemtools.
ArcticDEM and REMA are high-resolution, time-stamped 2-metre-resolution DEMs of the polar regions provided by the Polar Geospatial Center (PGC). They are extracted by applying stereo auto-correlation techniques (Noh & Howat, 2017) to pairs of submetre Maxar satellite imagery. The data includes Worldview-1, Worldview-2, Worldview-3, and GeoEye-1, beginning in 2007 (ArcticDEM) or 2009 (REMA) and ongoing to the present day. Products are available as tens of thousands of time-stamped ‘strips’ (Porter et al. 2022; Howat et al. 2022a) constructed from individual scene pairs, or as a single mosaic (Porter et al. 2023; Howat et al. 2022a) compiled from the combined stack of strips. Strips allow users to perform change detection by comparing data from different seasons or years, whilst mosaics provide a consistent and comprehensive product over the entire polar regions.
As Earth Science has moved into the ‘big data’ era, increasing amounts of Arctic- and Antarctic-focused resources are available as public, cloud-optimised datasets. New approaches are providing Python tools to act as combined API and processing tools, such as icepyx (Scheick et al. 2023) or pypromice (How et al. 2023). From 2022 (ArcticDEM v4.1 and REMA v2), the PGC DEM products are hosted as Cloud Optimised GeoTIFFs (CoGs) in a SpatioTemporal Asset Catalog (STAC), a standardised structure for cataloguing spatiotemporal data. However, the PGC STAC is currently a ‘static’ rather than ‘dynamic’ STAC, which means there is no convenient Application Programming Interface (API) for searching the datasets in response to user queries. This limits the ability of users to programmatically interact with ArcticDEM and REMA data in a quick and efficient manner. The pdemtools package has two aims: the first is to provide a Python-focussed alternative for searching and downloading ArcticDEM and REMA data, emulating dynamic STAC query tools such as pystac; whilst the second is to provide commonly used processing functions specific to the needs of ArcticDEM and REMA users (a focus on ice sheet and cryosphere work), as well as the particular strengths of the ArcticDEM and REMA datasets (high-resolution and multitemporal).
The pdemtools search() tool and load module allow for convenient access to the ArcticDEM and REMA datasets. Mosaics can be downloaded from a one-line load.mosaic() function, whilst the search() function allows for convenient filtering of a locally downloading ArcticDEM/REMA strip index according to variables such as date, region of interest, spatial coverage, temporal baseline, source sensors, accuracy, and cross-track data. The results of searches are returned as a geopandas dataframe, and can be downloaded using the load.from_search() function. Elevation models are returned as xarray DataArrays Hoyer _et al. 2017 with geospatial metadata via the rioxarray extension - a standard format for storing and processing n-dimensional geospatial data within the geospatial Python community. By utilising standardised formats, the aim is to allow the user to quickly move beyond pdemtools into their own analysis in whatever format they desire, be that xarray, numpy or dask datasets, DEM analysis Python packages such as xdem for advanced coregistration or richdem for flow analysis, or exporting to geospatial file formats for analysis beyond Python.
After download, there exist a number of (pre-)processing steps that are near universally common in topographic analyses. These include geoid-correction, co-registration of time-series data, and/or the construction of terrain parameters such as hillshade, slope, aspect, and curvature. pdemtools contains pre-built functions to perform these processing steps, as well as further functionality specific to ArcticDEM and REMA use cases. For instance, we include functions to quickly extract the EIGEN-6C4 geoid Foerste et al. 2014 and Greenland/Antarctic bedrock masks directly from local versions of the Greenland and Antarctic BedMachine datasets (Morlighem et al, 2022a, 2022b), reprojecting and resampling the data to match the target DEM. Options for ingesting user-provided mask and geoid data are also provided. Additionally, partial derivatives of the surface used to calculate terrain parameters (∂z/∂x, ∂z/∂y, ∂2z/∂x2, ∂2z/∂y2, ∂2z/∂x∂y) are calculated following Florinsky (2009), as opposed to more common methods such as Zevenbergen and Thorne (1987). The newer approach computes partial derivatives of elevation based on fitting a third-order polynomial, by the least-squares approach, to a 5 $\times$ 5 window as opposed to the more common 3 $\times$ 3 window. This is more appropriate for high-resolution DEMs: curvature over a 10 m window for the 2 m resolution ArcticDEM/REMA strips will lead to a local denoising effect that limits the impact of noise common in high-resolution photogrammetric products. These methods are also adapted into a co-registration routine, which otherwise follows the commonly used approach of Nuth & Kääb (2011).
We aim to grow pdemtools by implementing new methods developed by the ArcticDEM and REMA research community. For instance, we currently include sea-level-filtering and iceberg detection routines outlined by @shiggins_automated_2023, and invite community contributions or requests of other routines that will be of use to users of pdemtools. Ongoing research projects making use of pdemtools are applying ArcticDEM and REMA data to the mapping of crevasses, ice cliff heights, and subglacial lakes, as well as the initiation of ice sheet models. It has also been used within training exercises at the 2024 Polar Geospatial Center Data Workshop, contributing to a growing international network of pdemtools users.