GEDTM30: global ensemble digital terrain model at 30 m and derived multiscale terrain variables

Training point data

As training points we consistently use GEDI and ICESat-2, two global exhaustive state-of-the-art sources of spaceborne LiDAR that also record terrain height (Dubayah et al., 2020; Markus et al., 2017). Global Ecosystem Dynamics Investigation (GEDI) is a full-waveform satellite LiDAR on the International Space Station, collecting 3D measurements over land (through recorded backscattered laser energy, i.e., waveforms) near-globally, covering from 51.6°N to 51.6°S and with a laser footprint diameter of about 25 m. To derive terrain height information, we use level 2 product, i.e., GEDI Level 2 Geolocated Elevation and Height Metrics products (GEDI02). The product provides waveform interpretation, including ground elevation, canopy top height, and relative height (RH) metrics. We first downloaded the GEDI data set covering the period between 2019.04.18 and 2023.03.16 from the Land Process Distributed Active Archive Center (LP DAAC). The original data set is stored in HDF5 format along GEDI’s orbital flight path. To facilitate processing, we restructured the data into 5° $\times$ 5° tile and applied a pre-filtering step. This filtering excluded low-quality GEDI points, i.e., GEDI shots with quality flag equal to 0, degradation flag larger than 0, and sensitivity smaller than 0.95. The final dataset comprised approximately five billion GEDI samples covering the latitudinal range of 51.6°N and 51.6°S.

The Ice, Cloud, and land Elevation Satellite-2 (ICESat-2) is the 2nd-generation of the laser altimeter ICESat mission, covering the Earth between 88°N and 88°S. ICESat-2 carries a photon-counting laser altimeter that initially aims to measure the elevation of ice sheets, glaciers, and sea ice, but can also be used to measure tree canopy height and land topography. For terrain modeling, we utilized the ATLAS/ICESat-2 L3A Land and Vegetation Height (ATL08) version 6 (Neuenschwander & Pitts, 2019). This data set provides height measurements along the track of both ground and canopy surfaces relative to the WGS84 ellipsoid, derived from the Advanced Topographic Laser Altimeter System (ATLAS). Initially, ATL08 processed canopy and ground elevations in fixed 100-m data segments; however, the release of version 6 in May 2023 introduced 20-m segments, offering improved resolution for terrain height (h_te_best_fit_20m) and canopy height (h_canopy_20m).

After data import, we combined ICESat-2 ATL08 and GEDI02 for terrain modeling to produce a complete consistent data set with the best estimates of terrain height. We specifically used the relative height at 95% (rh95) as the canopy height of GEDI02 in the counterparts to ATL08 canopy height (column h_canopy_20m) and terrain height (column elev_lowestmode) in the counterparts to ATL08 terrain height (h_te_best_fit_20m), despite the slight difference in shot diameter (20 m vs. 25 m). Both ICESat-2 ATL08 and GEDI02 represent elevations in the WGS84 vertical datum. In order to be compatible with existing digital elevation models, we applied a conversion from the WGS84 ellipsoid to the EGM2008 geoid as vertical reference for ICESat-2 ATL08 and GEDI02, using the 5 arc-degree resolution correction surface (Pavlis et al., 2012). Despite all these efforts to harmonize these data, we still observed discrepancies in terrain height between spaceborne LiDAR datasets and DSMs and many artifacts: for example, the presence of terrain samples over bodies of water, including lakes and oceans. Additionally, we noticed that the data set we prepared overrepresents bare-earth and underrepresents forest areas, which is suboptimal for object removal applications.

Cleaning and stratification of training data

To ensure a high-quality and representative training dataset, ICESat-2 and GEDI points must still be filtered and checked for anomalies. To remove any remaining anomalies, we overlaid the 5° × 5° tiled ICESat-2 and GEDI samples with GLO30, AW3D30, MERIT DEM, and JRC’s Global Surface Water (Pekel et al., 2016). The Global Surface Water dataset provides long-term water fraction measurements on a per-pixel basis from 1984 to 2021. We also excluded any samples in which the surface water fraction exceeded $(>)$20%. In addition, samples that showed significant discrepancies from DEMs were removed. Subsequently, for each tile, we computed the average standard deviation across all samples. The standard deviation for each sample was derived from four height values: the terrain height of the sample itself and the corresponding elevations from MERIT DEM, GLO30, and AW3D30. The mean of these standard deviations was used to establish a discrepancy threshold. Samples were removed if their absolute error, calculated as the difference between the MERIT DEM and the individual sample, exceeded two standard deviations. This method ensured that each tile had a reasonable discrepancy bound, reflective of the landform variability within the tile. To avoid excessively large discrepancy thresholds, we capped the bound at a maximum of 20 m in cases where the computed standard deviation was unusually high. After filtering, we randomly selected exactly 500,000 samples per tile to create a global sample pool. If the total number of filtered samples in a tile was below this threshold, all available samples were included. This approach ensured a balanced distribution of samples across global tiles before proceeding with stratification.

Following the initial filtering, the remaining samples from individual tiles were aggregated into a global dataset. To ensure a comprehensive representation of global terrain and aboveground objects, we implemented two stratification strategies: one for artificial objects and another for vegetation objects. Artificial objects include buildings, roads, and other man-made structures; Vegetation objects are forest, grassland, cropland, and also the green cover in urban parks, etc. The samples were initially separated using the binary mask that outlines the extent of human settlement (World Settlement Footprint, WSF2019) (Marconcini et al., 2021). The vegetation object samples were selected by locating pixels where WSF2019 is equal to 0, and the artificial object samples were selected by locating pixels where WSF2019 is equal to 1. Aftermath, we implemented the stratification on two sets of samples, respectively:

• Stratification of samples for vegetation objects removal: We stratified the vegetation object samples by canopy height in a 5 m interval until 30 m. There are seven strata, 0–5 m, 6–10 m, 11–15 m, 16–20 m, 20–25 m, 25–30 m, and above 30 m. The canopy height obtained from the LiDAR samples from h_canopy_20m or rh95 from ICESat-2 or GEDI respectively. In addition, to detect vegetation in various landforms, we further stratified the samples by Iwahashi & Pike’s (2007) landform classification. This layer is a 1 km spatial resolution raster data, presenting relief classes which are classified using an unsupervised nested-means algorithms and a three part geometric signature. The nested partitioning applies gradient, connectivity, and textures for three rounds from eight partitions to total 16 partitions, where half of the steeper land pixels feeds into a round of classification. In total, this two-tier stratification resulted in 102 nested strata (7 canopy height classes $\times$16 landform classes). We randomly selected approximately 45,000 samples per stratum, yielding a total of approximately five million samples for vegetation object removal.

• Stratification of samples for artificial object removal: To remove artificial objects, samples were stratified according to building height using a 5 m interval up to 30 m, forming seven strata, 0–5 m, 6–10 m, 11–15 m, 16–20 m, 20–25 m, 25–30 m, and above 30 m. Building height data was extracted from the Global Human Settlement-Building Height dataset (GHS-BUILT-H R2022A) (Schiavina et al., 2022). This data set extracts the building height from the filtering of a composite of DEM and the filtering of satellite imagery using linear regression, providing global coverage at a 100-m resolution. For each building height stratum, we randomly selected 285,000 samples, ensuring a total of approximately two million samples for artificial-object removal.

During data processing, we identified a critical gap in the sample pool: the underrepresentation of mountain summits. These regions were either not measured by global LiDAR datasets or filtered out during the data cleaning process. Given the limitations of decision tree-based models, which extrapolate poorly beyond the observed data range, missing high-altitude samples could lead to inaccurate predictions in extreme elevation areas. To address this issue, we incorporated additional summit samples from OpenStreetMap (OSM) (OpenStreetMap Contributors, 2017) and the ICESat-2 and GEDI’s samples on the glaciers. The samples on glacier and permanent ice were identified using the University of Marylands Global Land Analysis & Discovery (GLAD)’s Global Land Cover and Land Use Change (GLCLUC) data set (Potapov et al., 2022), specifically selecting the points that are classified as permanent ice. Since altitude values in OpenStreetMap often differ from DEM-derived elevations, we replaced OSM summit heights by the average of the corresponding values from AW3D30 and GLO30 to ensure consistency with our model. Most of the OSM-derived summits were located in the Himalayas, the South Island of New Zealand, the Alps, and the Andes. Ultimately, we incorporated approximately 3,700 summit samples from OSM and randomly selected 50,000 additional LiDAR samples from glacier and ice-covered regions, enhancing the overall representativeness of high-altitude terrain and permanent ice of high latitude region in the dataset.

Covariates

To model global terrain height at a 1 arc-s grid spacing, we prepared an extensive set of covariates that capture essential information about terrain and surface objects. Two primary DSMs, GLO30 and AW3D30, were chosen as the state-of-the-art global open elevation data sets. However, both required preprocessing to ensure consistency in global modeling. For example, GLO30 provides no data for Azerbaijan, so we filled this gap using AW3D30 data. In contrast, AW3D30 contains sporadic but significant anomalies, including unrealistic craters and spikes, likely introduced during data processing. To correct for these inconsistencies, we derived the 95th percentile of absolute error between AW3D30 and GLO30 using the globally stratified sample set (as described in “Cleaning and stratification of training data”). This threshold was then used to identify and mask noisy pixels in AW3D30, replacing them with the corresponding values of GLO30. In addition, we incorporated global auxiliary data sets that significantly impact the removal of vegetation and artificial objects. Table 1 summarizes the covariates, including variable type, grid spacing, temporal range, and reference sources. The auxiliary layers can be roughly categorized into three thematic groups:

• Vegetation object modeling: To estimate vegetation object height, we included models of tree cover and canopy height. Specifically, we selected the GLAD Global Tree Cover 2010 dataset, derived from annual Landsat 7 ETM+ composites spanning 2000–2012 (Hansen et al., 2013). This dataset provides a continuous canopy density measure ranging from 0 to 100 and has been widely used in machine learning-based digital terrain model (DTM) applications (Hawker et al., 2022; Hengl et al., 2020). For canopy height estimation, we incorporated two of the most detailed global models: the GLAD UMD Canopy Height Model (Potapov et al., 2021) and the ETH Canopy Height Model (Lang et al., 2023). The GLAD UMD model integrates multiple regional models through local fitting, while the ETH model employs an ensemble of convolutional neural networks (CNNs) trained under a sparse supervision framework. Despite both models being trained on GEDI LiDAR data, they exhibit systematic biases, where GLAD UMD tends to underestimate canopy height, whereas ETH tends to overestimate it. In addition to the canopy height model, we also applied sobel_filter with a 3 $\times$ 3 window from Whitebox Workflow (Lindsay, 2016) to detect the edge of canopy and added them as covariates.

• Artificial object modeling: To model artificial objects, we included datasets that capture both the extent and height of built structures. The World Settlement Footprint 2019 (WSF2019) dataset was used to delineate human settlements. WSF2019, which also plays a role in sample stratification (“Cleaning and stratification of training data”), provides high-precision building outlines derived from Sentinel-1 and Sentinel-2 imagery in 2019. To estimate building height, we incorporated two global datasets: GHS-BUILT-H R2022A and 3D-GloBFP (Che et al., 2024). The 3D-GloBFP integrates remote sensing data and building morphology using an extreme gradient boosting (XGBoost) model. It is a vector dataset that provides building footprints with height attributes. We then rasterized 3D-GloBFP to 100 m resolution for alignment with GHS-BUILT-H R2022A. Although both datasets offer global coverage, we observed gaps in artificial object representation within 3D-GloBFP. To mitigate this issue, we gap-filled missing buillding in rasterized 3D-GloBFP layers using corresponding GHS-BUILT-H R2022A values.

• Ancillary remote sensing and topographic data: In addition to DSMs and structural data, we incorporated optical remote sensing indices and topographic derivatives to refine terrain modeling. Long-term Normalized Difference Vegetation Index (NDVI), Normalized Difference Wetness Index (NDWI), and Near-Infrared (NIR) reflectance were included for the periods 2006–2010 and 2011–2015. These temporal windows align with the acquisition periods of AW3D30 and GLO30, ensuring consistency in object removal. The indices were derived from Landsat Analysis Ready Data (GLAD ARD) (Potapov et al., 2020), a harmonized dataset integrating imagery from Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI/TIRS. For each time period, we removed cloud from the scenes, aggregated the available images and computed the 5th, 50th, and 95th percentiles for each index to capture the distribution of surface conditions over time. Additionally, we incorporated topographic information from ETOPO2022 (MacFerrin et al., 2024), an updated global topographic-bathymetric dataset at a 15-arc-s resolution. ETOPO2022 integrates various data sources, including bare-earth DEMs, sea and lake bathymetry, and ice-sheet bed elevation. Its land topography component is derived from aggregated FABDEM and GLO30, making it suitable for global DTM applications. From this data set, we extracted slope information to provide a gentle steepness.

Global-to-local transfer learning modeling for DTM prediction

A DEM is a general term that encompasses both DSMs and DTMs. A DSM represents the Earth’s surface, including all above-ground objects such as buildings and trees. In contrast, a DTM represents “bare earth”, without natural and human-made objects (Guth et al., 2021). The primary distinction between DSMs and DTMs lies in their application: DSMs are often used for telecommunication or other line-of-sight planning, while DTMs are used for simulating natural processes, such as hydrological modeling, landform classification, and geological structure analysis. In this study, we focus on modeling the Earth’s terrain, employing the best available methods to remove tree canopy and other surface objects from the DSM. In other words, we implemented a direct modeling of the terrain from DSMs to a DTM instead of e.g., fitting a model to estimate canopy height, then derive canopy height from a DSM. This approach is different from previous global DTM mapping research, which focuses on object detection and removal to minimize the alteration of the bare earth surface (Dusseau, Zobel & Schwalm, 2023; Hawker et al., 2022; Pronk et al., 2024; Uhe et al., 2025; Yamazaki et al., 2017). In contrast, our direct modeling approach adopts Random Forest (RF) to merge multiple DSMs as base input maps, and purposively add layers that represent or reflect the aboveground objects and terrain complexity. In addition, we adopt a framework to train a global-to-local model in order to address the global diversity of land cover and landform. Global-to-local modeling is a two-step transfer learning training framework: the global model trains with representative, globally stratified sampling, and the local model inherits the global model, training with additional locally stratified sampling. This framework not only preserves the generalization of the global model, but also reinforces the local features to achieve better performance (Fig. 2).

We first trained a global RF model using 10% random global samples for hyperparameter tuning, reducing the processing time and keeping the model representative of the whole data set. The remaining 90% samples were used in cross-validation. After model fine-tuning, all samples were then used to train a global RF model with 100 decision trees. This model was used as a starting point for training locally enhanced models with additional global LiDAR terrain height samples (transfer learning), selected according to 5° × 5° tiles. For each tile, additional samples, derived according to our global modeling stratification strategy (see “Cleaning and stratification of training data”), were used for hyperparameter tuning and training of 100 additional decision trees. Locally enhanced models, with 200 decision trees in total, were trained for all tiles that met the minimal training sample requirement ( $\sim$average 10 points in each vegetation stratum); otherwise, the global model was used to run DTM predictions. Our predictions estimated the mean and standard deviation of DTM elevation considering all individual decision trees of the final RF models, implemented with the library scikit-learn (an extend version of the class RandomForestRegressor).

Lastly, the prediction outputs were filtered by a adaptive filter (Gallant, 2011) and bilateral filter (Tomasi & Manduchi, 1998) to remove speckle noise from the data fusion. The filters are also used in FABDEM and MERIT DEM (Hawker et al., 2022; Yamazaki et al., 2017). We used the functions adaptive_filter, with 5 $\times$ 5 neighborhood size, and the function bilateral_filter, setting the standard deviation distance to 2 and the standard deviation intensity to 3. The functions are implemented in Whitebox Workflows and the filters were only applied on the tree and building cover pixels (Global 2010 Tree Cover is not 0 and WSF2019 is not 0), in order to avoid potential artifacts in non-aboveground-object areas.

Model evaluation

We validated our models and predictions using independent reference datasets and using visual inspection and evaluation metrics. The quality of the DTM is not only the accuracy of pixel values to the reference, but also the consistency of the terrain surface. Therefore, the validation focuses on not only the accuracy metrics, but also the cross section, spatial visualization on DTM and hillshading. We incorporated independent source of validation datasets including:

• High quality LiDAR datasets: A collection of reference LiDAR DTMs (Guth et al., 2023) used by Bielski et al. (2024), Guth et al. (2024) for the evaluation of global DEMs. The reference DTM locations are mainly in Europe, North and South America, and completely under-represent Asia and Africa.

• Global distributed Global Navigation Satellite System (GNSS) stations: Altitude data from GNSS stations collected from Nevada Geodetic Laboratory (geodesy.unr.edu/velocities/midas.IGS14.txt). The data set is global coverage and is also used to validate the Global Digital Elevation Merged Model 2024 (GDEMM2024) (Ince, Abrykosov & Förste, 2024).

In addition to the validation dataset, three metrics are employed to measure the prediction performance: the root mean square error (RMSE), the mean absolute error (MAE), and their mean error (ME):

(1) $$RMSE=\sqrt{\frac{1}{N}{\mathrm{\Sigma }}_{i=1}^N{(\hat{y}\text{-}y)}^2},$$

(2) $$MAE=\frac{1}{N}{\mathrm{\Sigma }}_{i=1}^N|\hat{y}\text{-}y|,$$

(3) $$ME=\frac{1}{N}{\mathrm{\Sigma }}_{i=1}^N(\hat{y}\text{-}y).$$

RMSE is more sensitive to outliers, whereas MAE measures individual differences equally weighted; ME measures the prediction bias, indicating whether the model overestimates or underestimates the height of the terrain. In addition to point-to-point metrics, for comparing with raster reference DTMs about preserving DTM structure, we employed structural similarity index measurement (SSIM) (Wang & Bovik, 2009; Wang et al., 2004).

Pyramid representation of DTM

Multiscale terrain variables, also called terrain attributes, were derived from the GEDTM30, which is crucial to modeling soil, water, and biological processes (Behrens et al., 2018; Kulp & Strauss, 2018; Lozbenev et al., 2021; Moore, Grayson & Ladson, 1991). The multiscale process comprises two steps: (1) the DTM is down-scaled to coarse resolution, and (2) the original layer and the down-scaled DTM layers run the terrain variables derivation independently (Fig. 3). In our research, the objective is to produce the terrain variables in grid spacing of 30, 60, 120, 240, 480 and 960 m, increasing scale by a factor of 2. We ran the derivation using Whitebox Workflows version 1.3.4 (WbW) in Python (https://www.whiteboxgeo.com/whitebox-workflows-for-python/). We set identical arguments for each function on different resolutions to make all derivation consistent. For resampling from coarser to finer resolution, we recommend using cubic splines interpolation. Note that we have initially tested derivation in SAGA GIS, GRASS GIS and GDAL (all provide extensive libraries for terrain variables derivation), however, Whitebox Workflows have shown to be fastest and easiest to implement on large grids.

Terrain variables

There are two types of terrain variables: regional and local (Olaya, 2009b; Shary, 1995; Wilson, 2012; Wilson & Burrough, 1999). The total number of information that can be derived from a relief model is extensive and probably counts in hundreds of variables. We selected 15 fundamental terrain variables that represent measures of local topographic position, curvature, hydrology, incoming solar light and shadow. Local terrain variables include hillshade, slope, openness, curvature and topographic position such as deviation from mean elevation and spherical standard deviation of normals. They are introduced below:

• Hillshade: Hillshade, also called shaded relief, is a terrain visualization technique that accounts for the positioning of a theoretical light source and the local slope and aspect of the elevation surface (Horn, 1981). In this research, we derived mutlidirectional hillshade, which integrates hillshade values calculated from multiple sources of illumination. We used the function MultidirectionalHillshade and set altitude to 30 degrees, and Z-factor to 1 as default.

• Slope in degree: The slope gradient, measured in degrees, measures the steepness of the terrain surface in each grid cell in a DEM. The algorithm used to calculate the slope used the third-order bivariate Taylor polynomial method described by Florinsky (2016). We used the function slope, which applies a polynomial fit of the elevations within the 5 × 5 neighborhood surrounding each cell, to mitigate against the outlier elevations.

• Openness: The dominance or enclosure of a location on an irregular surface. It has two viewer perspectives, positive and negative openness. The mean value of all zenith angles gives positive openness, while the mean nadir value gives negative openness (Yokoyama, Shirasawa & Pike, 2002). In our research, we used the function openness applied the formula from Yokoyama, Shirasawa & Pike (2002) to derive both positive and negative openness. We set search distance to three grid cells.

• Spherical standard deviation of normals: The measure of the angular dispersion of the surface within a local neighborhood of a specified size. It indicates the complexity, texture, and roughness of the surface shape. Grohmann, Smith & Riccomini (2010) discovers that the vector dispersion, a related measure of the angular dispersion, increases monotonically with scale. The scale relation can therefore be estimated by isolating the amount of surface complexity associated with specific scale ranges: large scale reflects the texture; small scale reflects the complexity. We used the function spherical_std_dev_of_normals and the local neighborhood size was set to 3 in our research.

• Difference from mean elevation: The calculation of the difference between the elevation of each DEM pixel and the mean elevation within a moving window of a size. It is also a measure of local topographic position. In our research, we used the function difference_from_mean_elevation and set both filterx and filtery to 3.

• Maximal curvature: The curvature of a principal section with the highest value of curvature at each DEM pixel. The positive values correspond to ridge positions while negative values are indicative of closed depression (Florinsky, 2016). It is measured in unit m⁻¹. We used the function maximal_curvature, applied the formula from Shary, Sharaya & Mitusov (2002), with a polynomial fit of elevations within the 5 $\times$ 5 neighborhood surrounding each cell. Log-transform is adopted.

• Minimal curvature: The curvature of a principal section with the lowest value of curvature at each DEM pixel. The positive values correspond to hills while negative values are indicative of valley positions. It is measured in unit m⁻¹. We used the function minimal_curvature, applied the formula from Shary (1995), Shary, Sharaya & Mitusov (2002), with a polynomial fit of the elevations within the 5 $\times$ 5 neighborhood surrounding each cell. Log-transform is adopted.

• Profile curvature: Profile, or vertical curvature (unit: m⁻¹) is the curvature of the rate of change in slope along a flow line. Positive values correspond to flow acceleration, while negative values are indicative of flow deceleration. It is the curvature of a normal section having a common tangent line with a slope line on each DEM surface. We used the function profile_curvature, applied the formula from Shary (1995), Shary, Sharaya & Mitusov (2002), with a polynomial fit of the elevations within the 5 $\times$ 5 neighborhood surrounding each cell. Log-transform is adopted.

• Tangential curvature: Tangential, or horizontal curvature (unit: m⁻¹) is the curvature of an inclined plane perpendicular to both the direction of flow and the surface. The positive values correspond to flow divergence while negative values are indicative of flow convergence (Wilson & Gallant, 2000). We used the function tangential_curvature, applied the formula from Shary (1995), Shary, Sharaya & Mitusov (2002), with a polynomial fit of the elevations within the 5 $\times$ 5 neighborhood surrounding each cell. Log-transform is adopted.

• Ring curvature: The product of horizontal excess and vertical excess. Ring curvature (unit: m⁻²) has values equal to or greater than zero, and high values indicates the greater potential of the flow motion. It is used to measure flow line twisting (Shary, 1995). We used the function ring_curvature, applied the formula from Shary (1995), Shary, Sharaya & Mitusov (2002), with a polynomial fit of the elevations within the 5 $\times$ 5 neighborhood surrounding each cell. Log-transform is adopted.

• Shape index: The shape indicator. It is independent of the size i.e., the amount of curvature, as distinct from the type of curvature. The shape index represents the intuitive notion of local shape, with positive values indicative of convex landforms, negative values corresponding to concave landforms (Koenderink & Van Doorn, 1992). We used the function shape_index applied the formula from Koenderink & Van Doorn (1992), with a polynomial fit of the elevations within the 5 $\times$ 5 neighborhood surrounding each cell.

The regional terrain variables consider other parts of the DEM or the entire scene in addition to the exact point where they should be calculated (Olaya, 2009a). We selected mainly hydrological variables that require flow enforcement to ensure that interior topographic depressions and flat areas have been removed from the input DEM. The depression-breaching algorithm used for flow enforcement in this work was implemented based on Lindsay (2020). The approach uses a least-cost path analysis to identify the breach channel that connects pit cells (i.e., grid cells for which there is no lower neighbor) to some distant lower cell. We used the function BreachDepressionsLeastCost, which combines breaching at first and filling at the second. We set the maximal search distance to 100 pixels, and used an increment value of 0.001 to connect input DTM pixels within depressions to lower cells. Any depressions that could not be resolved using a 100-cell-long breach channel were subsequently resolved using depression filling. After preprocessing the DTM, we produced specific catchment area, slope length and steepness factor, and topographic wetness index:

• Specific catchment area: The upstream catchment area of a unit contour length. In our research, we computed SCA using Qin et al.’s (2007) method. It relates the degree of flow dispersion from a grid cell to the local maximum downslope gradient. We used the function QinFlowAccumulation and set upper-bound exponent and slope to 10 and 45 respectively.

• Slope length and steepness factor: One of the components used to estimate the physical potential for sheet and rill erosion in upland catchments (Moore, Grayson & Ladson, 1991), so-called length-slope factor (LS factor) or sediment transport index. The inputs of LS factor are slope in degree and SCA which are computed in advance and described above. We used the function sediment_transport_index and set SCA exponent to 0.4 and slope exponent to 1.3. Log-transform is adopted for the result.

• Topographic wetness index: The index describes the propensity for a site to be saturated to the surface given its contributing area and local slope characteristics. It is commonly used in the TOPMODEL rainfall-runoff framework (Beven, 1997). The inputs of TWI are slope in degree and SCA which are computed in advance. We use the function WetnessIndex. Log-transform is adopted.

The methods that were adopted in the derivation functions and their details are documented in, and all right are reserved ©WhiteboxTools user manual (Lindsay, 2018).

Preparation of land mask

We created an ensemble land mask based on GLO30, AW3D30 and ESA World Cover. The three products so far have the most complete land coverage, and they are complimenting each other in creating the most inclusive land mask. The World Cover of ESA is a global land cover product for 2020 and 2021 at 10 m resolution (Zanaga et al., 2022) and contains the finest land detail, but inland water is defined as the same class as ocean. Masking the water would end up losing the inland water pixels. GLO30 and AW3D30 are DSMs (described in “Covariates”). Both cover the whole period of 2006–2015, being possible to capture the emergent islands or submerged coastline due to sea level rise in this period. However, these DSMs define the ocean as sea level (altitude $=0$) which could also happen at inland area and lead to losing inland area pixels where its elevation is equal to 0. Finding the union of “not water” for the ESA World Cover and “not at sea level” for GLO30 and AW3D30 is the solution to create the only ocean mask land mask.

To produce a consistent land mask we aggregated ESA world cover to 1 arc s using mode as resampling method in GDAL to match the grid spacing of GLO30 and AW3D30. If the pixel reaches the condition where ESA World Cover is not water, or GLO30 is not 0, or AW3D30 is not 0, the pixel would be identified as land, and set to value 100. The other pixels are then set to 0 and saved as no data. In addition, the definition of the coastline could be controversial to countries, and a pixel can also not be recovered once it has been masked. Therefore, we decided to apply a buffer to the only-ocean-mask land mask. The buffer was created by aggregating the 1-arc-s mask to 30-arc-s grid spacing by maximum, to extend the land $\sim$1 km at the coastline. Subsequently, the 30-arc-s mask was then down-scaled back to 1-arc-s grid spacing by cubic spline. Ultimately, we filtered all the pixels with value (not 0) to 100 and defined 100 as the land pixel, and defined 0 as the ocean pixel. The land mask was used in the production of GEDTM30 and also in the production of 15 derived terrain variables.

Tiling system and projection

We implemented data parallelization followed by the defined Equi7 tiling system (Bauer-Marschallinger, Sabel & Wagner, 2014). Equi7 is a continent-based projection system that minimizes data oversampling over global land surface. However, if we would implement the original Equi7 tiling system to derive variables without any overlap, final mosaics would eventually show tile-boundary effects. To resolve this issue, a spatial overlap is required to create seamless global maps which are artifact-free and hydrologically connected. The previous article published by Amatulli et al. (2022) proposed an irregular tiling system, which retains the entire drainage basin in each tile to preserve the lateral and longitudinal connectivity. Such an irregular tilling system contains instances like partially included drainage basins across multiple tiles, or the tile could not be enlarged beyond the maximum number of grid cells. The irregular tiling system used by Amatulli et al. (2022) is in our opinion not generalizable to multiscale terrain variables derivation. Instead, we implemented an extended tiling system with consistent tiles optimized for hydrological modeling. We considered three aspects of parallelization:

• Size of tiles;

• Overlap of tiles;

• Software implementation (aiming at optimized computing and use of RAM) by using S3 (Amazon’s Simple Storage Service) and InfiniBand to access large number of tiles in parallel;

After initial testing, we decided that the 600 km by 600 km Equi7 tiles could be ideal as the tiling basis for derivation. To derive the local terrain variables we padded 1.92 km on each Equi7 tile, which preserves at least two pixels from all sides in the coarsest resolution (960 m) input DTM; Regional terrain variables instead need to consider the whole scene of the input data, so larger overlap is required. We discovered that padding with 16%, which is about 102 km in 1 arc s input DTM, is sufficient to mitigate the boundary effect. Although the extended tile caused around 1.78 times of data operation toward size and time, it ensures the hydrological connectivity persist until the boundary of the Equi7 tile.

In practice, the input GEDTM30 was first reprojected onto the Equi7 projected coordinate system (geographic to projected coordinate system) by cubic spline resampling approach, and cropped in extended tiles for the regional terrain variables. Subsequently, each tile was aggregated to multiscale resolution in 60, 120, 240, 480 and 960 m by average using GDAL. The original 30 m tile applied a Gaussian filter to reduce noise. Aftermath, the DTM tiles ran through the regional derivation function sequentially, with the breaching and depression-filling procedure applied before the hydrological analysis. Aftermath, the processed local and regional variables were truncated back to 600 km by 600 km with 960 m padding. Table 2 summarizes the overall information of the 15 terrain variables mentioned in “Terrain Variables” with tile overlapping distance.

The conversion between the Equi7 and WGS84 coordinate systems requires special consideration near the 180th meridian. The Equi7 grid is continent-based, adhering to the boundaries of continental landmasses, whereas the WGS84 system uses a global reference ellipsoid with a defined longitudinal boundary at the 180th meridian. During reprojection, if an Equi7 tile spans the 180th meridian, the reprojected tile in WGS84, similar to a polygon, would contain the few vertices approximate to 180^∘ and also some approximates to −180^∘ on the other side. Creating a bounding box from these reprojected tiles for cropping input GEDTM30, the longitudinal range would extend across the entire world from −180^∘ to 180^∘. Consequently, the attempt to read the huge bounding box of the 1-arc-s input DTM results in dramatically slowing down the process or even excessive memory usage. To address this issue, a dedicated set of irregular tiles in W84 was generated and reprojected into the Equi7 system. These irregular tiles are strategically located to mitigate boundary effects in variables and ensure comprehensive data coverage. The issue is particularly pronounced in regions such as the Chukotka Autonomous Okrug in Russia, the North Island of New Zealand, and several Pacific islands. Additionally, the Equi7 Oceania tile system fails to adequately cover specific Pacific territories, including Nouvelle-Calédonie, Yap (Micronesia), Palmerston Island, etc. To resolve this, additional irregular tiles are created in the EPSG:4326 coordinate system, ensuring seamless projection transformations without exceeding the 180th meridian boundary and providing full coverage of these regions. Figure 4 illustrates the tiling system, including regular and irregular tiles that resolve the issues of the 180th meridian and missing islands.

Global DTM production result

Figure 5 shows the distribution of the tiles whether local models were produced. As a result, the global-to-local model is used to predict 1,157 land tiles in 5^∘ ×5^∘ grid. Ultimately, 1,019 local models have been fitted, where 138 tiles do not have sufficient samples for the local models. The red tiles were only predicted by global models due to insufficient amount of training data. The background is the land mask that was specially created to buffer the coastline. It also shows the promising result that covers small Pacific islands such as Nauru and French Polynesia and covers until the northernmost point of land.

Assessment of global-to-local modeling

The analysis conducted for all reference DTMs demonstrates a reduction in the average RMSE from 4.14 to 3.86 m, representing an approximate 9% improvement when comparing the DTM generated using only the global model to that produced by the global-to-local model. Figure 6 illustrates an visual example of the differences between the predictions and the reference DTM. Assessment of the global-to-local model demonstrates an improvement in local accuracy while maintaining a consistent terrain representation. To evaluate the impact of incorporating additional local samples, we generated predictions using both the global model and the global-to-local (locally enhanced) model for an area in Trentino, Italy. Two tiles were produced and compared against the reference DTM. In summary, this shows that the locally enhanced model exhibits improved accuracy. The upper histogram in Fig. 6 shows that the locally enhanced model reduces both error and bias compared to the global model alone. The lower map highlights how additional local samples contribute to improving accuracy while maintaining the consistency of the terrain.

Accuracy comparison with source DSMs

We further evaluated GEDTM30 against source DSMs (GLO30 and AW3D30) to assess its ability to remove vegetation and artificial objects while improving accuracy. First, we used a multidirectional hillshade (WhiteboxTools QGIS Plugin, 45-degree altitude) to visually inspect the efficiency of object removal. Figure 7A illustrates the removal of forest and vegetation objects in Departamento de Artigas, Uruguay. The red circle and arrows highlight a tree patch visible in the Google Satellite image, which is effectively removed in GEDTM30. Furthermore, GEDTM30 integrates the roughness characteristics of GLO30 and AW3D30. Figure 7B demonstrates the removal of artificial objects in São Paulo, Brazil, where buildings are effectively removed while topographic characteristics are preserved.

Next, Fig. 8 presents a comparative evaluation of object removal against reference DTMs. We evaluated GEDTM30, GLO30 and AW3D30 using cross-sectional analysis, visual inspection, and error distribution analysis. The left panel (cross section) of Fig. 8 displays terrain roughness consistency and noise levels between DEMs. The middle panel provides a full-scene visualization of object removal and modeled terrain consistency, in order to inspect occasional model-induced hallucinations arose due to optimization in machine learning methods (Guth et al., 2024). The right panel presents error distributions via histograms, highlighting the effectiveness of object removal. For removal of vegetation objects, we extracted the pixel that falls on the canopy cover 10–50% and the canopy cover $>50$% for comparison. For the removal of artificial objects, we mask the pixels where WSF2019 classified areas as non-artificial objects. Figure 8A focuses on forest/vegetation object removal in Departamento de Artigas, Uruguay and Fig. 8B focuses on building/artificial object removal in São Paulo, Brazil.

Third, GEDTM30 is validated with reference DTMs using evaluation metrics, RMSE, MAE, and ME. We compared with the source DSMs, GLO30 and AW3D30 at the visual inspected sites (Fig. 8), São Paulo for artificial removal and Departamento de Artiga, Urugauy for vegetation removal. Table 3 summarizes the accuracy metrics grouped by inspection sites and land use and cover (LULC). The results show that GEDTM30 has the best RMSE, MAE and ME of all groups, and has approximately a 30% improvement in accuracy of GLO30 in the build-up and tree cover $>50$%.

Accuracy comparison with other global DTMs

GEDTM30 was also compared with other global DTMs, including MERIT DEM, FABDEM, and FathomDEM (Hawker et al., 2022; Uhe et al., 2025; Yamazaki et al., 2017), across diverse landforms. We selected three representative study areas: Trentino, Italy (Fig. 9), Chincoteague, the United States, (Fig. 10), and Vestfold, Norway, (Fig. 11), representing mountainous region, flood plain, and boreal forest. The first column of these figures is the visualization of the reference DTM and the global DTMs. The second column is the difference (reference DTM-global DTM), and the third column is the mutlidirectional hillshade visualization.

In Trentino (Fig. 9), all DTMs have performed better in the valley surrounded by mountains, which is relatively flat. As the slope increases, MERIT DEM has a significantly higher terrain estimate than other two models, and FABDEM also overestimate the terrain height/underestimate the tree height significantly in the steep area in the scene. GEDTM30 and FathomDEM have a rather balanced result. However, GEDTM30 slightly overestimates the steep area and underestimates the high attitude area in the bottom right corner and top left corner of the scene. On the right side of Fig. 9, mutlidirectional hillshade shows GEDTM30 is slightly noisier/softer than FABDEM and FathomDEM as the same 1-arc-s grid spacing DTMs. Finally, FABDEM appears to show pit holes on the slope, locating at the left and upper side of the valley, which could be due to the issues in the height removal algorithm used.

In Chincoteague (Fig. 10), a flat coastal plain, we incorporated land-use/land cover map change from GLAD for the period 2000–2021 (Potapov et al., 2022) to collate the difference map. The primary land use and land cover here is the build-up area (cyan), the stable crop field (orange), and the stable tree cover (green). FABDEM and GEDTM30 have a more significant error in stable tree cover than in built-up area and stable crop field. GEDTM30 does not remove some tree patches, leading to an overestimate of the terrain height, whereas FABDEM has a more extreme error than GEDTM30. It appears that FathomDEM performs the best in object removal and has the smallest error. Ultimately, visual inspection of all DTMs using mutlidirectional hillshade, shows that FathomDEM appears to be the most accurate containing the least artifacts, followed by MERIT DEM, and then GEDTM30 and FABDEM.

Vestfold, Norway (Fig. 9) is a boreal forest region with rugged terrain. The difference against the reference DTM is generally smaller than in Trentino (Fig. 9), but both GEDTM30 and FABDEM tend to underestimate terrain height, whereas MERIT overestimates elevations. FathomDEM hence has a rather balanced result. Overall, the high error region, same as Trentino, locates at steeper area. Fathom remains the smallest error among the target DTMs. The terrain roughness in GEDTM remains the greatest however.

Table 4 summarizes the accuracy metrics of the above three sites. The results show regardless of the locations, show the consensus of better accuracy in FathomDEM, followed by GEDTM30 and then FABDEM and MERIT DEM. In addition, to better estimate the model performance on different land cover, we aggregated the pixels from all the reference DTM sites and grouped them by land use and land cover. Table 5 summarizes the result of the evaluation metrics on various stable LULC classes aggregated from GLAD Global LULC Change (Potapov et al., 2020). As a result, FathomDEM performs overall the best in cropland, built-up, forest areas, whereas FABDEM seems to have a lower error in bare earth and sparse grassland. GEDTM30, FathomDEM, and FABDEM all maintain mean errors (ME) under 1 m for most of the LULCs, except for dense short vegetation in GEDTM30. On average, Fathom has $\sim 25$% smaller RMSE than GEDTM30.

Beyond remote sensing-derived DTMs, we further validated models using independent, globally distributed GNSS station data. Figure 12 illustrates the distribution of the GNSS stations and the comparison of heights of GNSS stations and target DTMs. Figure 12A is the distribution of the GNSS stations and the elevation of each record in orthometric height. Figure 12B presents a histogram of GNSS-DTM elevation differences (GNSS height-DTM). The legend indicates accuracy assessments among the global DTMs. GEDTM30 achieves the best mean at 7.34 m, standard deviation at 7.77 m, and RMSE at 10.69 m when compared to GNSS stations. However, systematic discrepancies persist between remote sensing DTMs and GNSS measurements, as previously reported by Ince, Abrykosov & Förste (2024) in their validation.

Software testing and optimization

This section summarizes the results of optimization test and the estimated computation time for global terrain variables derivation. The test was carried out by running the 60 m grid spacing derivation for all tiles (operation details in “Optimizing of terrain variables derivation”). The processing time for each tile and process is recorded. Figure 14 illustrates the average individual and cumulative processing time of terrain variables derivation in 60 m grid spacing. According to the Fig. 14, breaching depression takes the longest average processing time at 194 s per tile and the topographic wetness index takes the shortest to complete at the average processing time of 5 s. The complete computational workflow averages 742 s per 600 km by 600 km Equi7 tile at 60 m grid spacing. Extrapolating from this, a 30 m resolution grid spacing is estimated to complete within two days, using 1,344 CPUs on 14 servers equipped with 1 TB of RAM each and connected to S3 storage servers with InfiniBand. Using this infrastructure, the possible update for the multi-resolution variables can be completed in less than 60 h.

Gallery of derived terrain variables

A total of 15 GEDTM30-derived terrain variables were generated, covering aspects of topographic position, light and shadow, landform characteristics, and hydrology. Figure 15 shows an example of all variables at a grid spacing of 1 arc s in Virunga National Park, located at the intersection of Rwanda, Uganda, and the Democratic Republic of the Congo (DRC). The variables include positive openness, negative openness, slope in degree, profile curvautre, tangential curvature, mutlidirectional hillshade, minimal curvature, maximal curvature, shape index, specific catchment aream topograhpic wetness index, ring curvauture, spherical standard deviation of normals (SphericalStdDevOfNormals), difference from mean elev (DiffFromMeanElev), and LS factor. These variables were generated using functions within Whitebox Workflows, with additional derivation details provided in the Materials and methods section.

Figure 16 additionally illustrates the pyramid representation of the topographic wetness index (TWI) at different grid spacing (30, 60, 120, 240, 480, and 960 m) in the Jhuoshuei River Basin, central Taiwan. The representation of the TWI pyramid in the alluvial plain highlights how the channel size and the distinctiveness of the cells evolve with increasing resolution. In the flat plain region shown in Fig. 16, the 30 m grid spacing TWI clearly differentiates between wet and dry pixels. As the size of the grid increases, the primary channel becomes more prominent and by the coarsest resolution (960 m), the flat plain exhibits a nearly uniform wetness distribution. This multiscale TWI representation effectively captures a gradient of wetness, where some cells remain consistently dry, others transition from dry to wet at coarser resolutions, and some are persistently wet across all scales.

Consistency check

Figure 17 presents an assessment of spatial consistency across both processing tile and watershed boundaries, as well as a comparison with a leading global hydrological product, Hydrography90m (Amatulli et al., 2022). Two large drainage basins were selected for visual inspection: Lake Victoria in Africa (A) and Paraná Basin in South America (B). Watershed delineations are based on HydroBASINS Level 3 (Lehner & Grill, 2013), overlaid with the Equi7 processing tile grid in the upper panels of Fig 17. The lower panels display the TWI from GEDTM30 at the pixel level, alongside the corresponding TWI from Hydrography90m. The GEDTM30 TWI exhibits finer spatial detail due to its higher resolution, while maintaining a consistent spatial pattern with Hydrography90m. In addition, GEDTM30 preserves continuity not only across HydroBASINS watershed boundaries but also across processing tile boundaries.