Global digital elevation model (GDEM) product generation by correcting ASTER GDEM elevation with ICESat-2 altimeter data

Advancements in scientific inquiry and practical applications have created a higher demand for the accuracy of global digital elevation models (GDEMs), especially for GDEMs whose main data source is optical imagery. To address this challenge, integrating GDEM and satellite laser altimeter data (global coverage and high-accuracy ranging) is an important research direction, in addition to the technological enhancement of the main data source. In this paper, we describe the datasets and algorithms used to generate a GDEM product (IC2-GDEM) by correcting ASTER GDEM (Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model) elevation data with ICESat-2 altimeter data. The algorithm scheme presents the details of the strategies used for the various challenges, such as the processing of DEM boundaries, the fusion of the different data, and the geographical layout of the satellite laser altimeter data. We used a high-accuracy global elevation control point dataset and multiple high-accuracy local DEMs as the validation data for a comprehensive assessment at the global scale. The results from the validation comparison show that the elevation accuracy of IC2-GDEM is evidently superior to that of the ASTER GDEM product: (1) the RMSE reduction ratio of the corrected GDEM elevation is between 16 % and 82 %, and the average reduction ratio is about 47 %; and (2) from the analysis of the different topographies and land covers, this error reduction is effective even in areas with high topographic relief (>15°) and high vegetation cover (>60 %). ASTER GDEM has been in use for more than a decade, and many historical datasets and models are based on its elevation data. IC2-GDEM facilitates seamless integration with these historical datasets, which is essential for longitudinal studies examining long-term environmental change, land use dynamics, and climate impacts. Meanwhile, IC2-GDEM can serve as a new complementary data source for existing DEMs (such as Copernicus DEM) mainly sourced from synthetic aperture radar (SAR) observation. By cross-validating qualities, filling data gaps, and conducting multi-scale analyses, it can lead to more reliable and comprehensive scientific discoveries, thereby improving the overall quality and reliability of Earth science research. The IC2-GDEM product is openly available at https://doi.org/10.11888/RemoteSen.tpdc.301229 (Xie et al., 2024).

1 Introduction

High-quality digital elevation model (DEM) data are of great significance for the field of Earth science research and scientific applications, including hydrological modelling, climate change research, natural hazard assessment, and ecosystem management. For example, high-quality DEMs can improve the accuracy of watershed delineation and flood risk assessment, promoting effective water resource management and disaster preparedness (Ao et al., 2024). In climate change research, high-quality DEMs are vital for the modelling of glacier dynamics and the sea level rise impact, thereby improving predictive accuracy and informing mitigation strategies (Fan et al., 2022; Cook et al., 2012; Tran et al., 2023; Tran et al., 2024a, b).

Airborne lidar, radar, and high-resolution imaging systems can produce DEMs with centimetre to sub-metre accuracy, but this kind of DEM is costly to derive. As a result, such DEMs are not extensively accessible at the global scale (Pham et al., 2018). Laser altimeters carried on satellite platforms (Schutz et al., 2005; Xie et al., 2021b; Dubayah et al., 2020; Martino et al., 2019; Tang et al., 2020), such as ICESat and ICESat-2, can provide accurate elevation information around the world, and their elevation information is freely available. The accuracy of the elevation information can attain a metre, sub-metre, or even centimetre level after data refinement (Li et al., 2023b; Zhu et al., 2022; Neuenschwander et al., 2020; Fernandez-Diaz et al., 2022; Li et al., 2021). Previous studies have refined ICESat and ICESat-2 altimeter data in the polar regions to produce and release polar DEM products (Shen et al., 2022; Fan et al., 2022). However, at present, it is still difficult to use satellite laser altimeter data to generate DEMs with a medium resolution at the global scale, due to the limitation of the spacing between ground tracks (Magruder et al., 2021). At the global scale, analysis of the spatial elevation at a medium resolution still relies on DEMs derived from other satellite technologies (optical imaging and synthetic aperture radar – SAR), but the accuracy of these kinds of DEM data is lower, at multiple metres or even more than 10 m (Meadows et al., 2024; Al-Areeq et al., 2023; Yue et al., 2017; Del Rosario González-Moradas et al., 2023; Hawker et al., 2019). This data characteristic is challenged by related scientific applications with a requirement for higher and higher elevation accuracy. Therefore, improving the accuracy of these kinds of DEM data has been widely focused on and studied by many scholars. In addition to technological improvement, integrating a DEM with satellite laser altimeter data is an important research direction. Scholars have also provided important guidance and reference to this field by using linear fitting, machine and deep learning, and other methods (Li et al., 2023c; Magruder et al., 2021; Hawker et al., 2022).

Free global DEMs (GDEMs) with a medium resolution mainly include the Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model (ASTER GDEM) (Abrams et al., 2015), the Shuttle Radar Topography Mission (SRTM) DEMs (Farr et al., 2007), and the Copernicus GLO-30 DEM (Fahrland et al., 2022). The data sources of these DEMs can be divided into two main kinds: (1) optical imaging (for ASTER GDEM) and (2) SAR (for the SRTM DEMs and Copernicus GLO-30 DEM). Owing to the noise and anomalies resulting from the limitations inherent in optical imaging, the elevation quality of ASTER GDEM is typically deemed to be lower than that of the other GDEMs, for which the source is radar data (Meadows et al., 2024; Del Rosario González-Moradas et al., 2023; Purinton and Bookhagen, 2021), such as the SRTM DEMs and Copernicus GLO-30 DEM. In addition, a survey of some studies showed that the role of GDEMs for which the main source is optical imaging has not been completely replaced. For example, the ASTER GDEM elevation shows a higher quality than the SRTM elevation in some mountainous areas (Li et al., 2013; Yue et al., 2017) and also shows a better quality than the Copernicus GLO-30 DEM in some areas with steep topography (Del Rosario González-Moradas et al., 2023) or high vegetation cover (Okolie et al., 2024; Huang and Yu, 2024). Furthermore, some studies (Pham et al., 2018; Yamazaki et al., 2017; Crippen et al., 2016; Franks et al., 2020; Okolie and Smit, 2022) have reported that diversified choices and complementary use strategies for GDEMs are still favoured by many scholars and users. For instance, the data sources for the void-filling of the Copernicus GLO-30 DEM (currently the best available global DEM) consist of other DEMs, including ASTER GDEM, the SRTM DEMs, and several national DEMs (Del Rosario González-Moradas et al., 2023). Thus, obtaining more accurate ASTER GDEM elevation products has great significance for advancing the research of the GDEM-diversified choices and complementary use strategies. Researchers can use the enhanced ASTER GDEM in conjunction with other DEM products to cross-validate qualities, fill data gaps, and conduct multi-scale analyses (Del Rosario González-Moradas et al., 2023). This complementary use of multiple DEMs can lead to more reliable and comprehensive scientific discoveries, thereby improving the overall quality and reliability of geoscience research (Del Rosario González-Moradas et al., 2023; Ao et al., 2024). However, few previous studies have focused on integrating ASTER GDEM and satellite laser altimeter data to generate and release a new enhanced ASTER GDEM product at the global scale.

In this study, given this issue, we collected multiple sets of data with a global coverage, including ICESat-2 altimeter data, a global land cover product, a global vegetation index product, and ASTER GDEM Version 3 data (the most up-to-date version), to correct the ASTER GDEM elevation and build a corresponding corrected product to further explore and enhance the applicability of the GDEMs whose main data source is optical imagery. Taking into account the particularity of polar areas, i.e. the high variability of ice sheets and ice flow rates, the corrected product covers the global land areas but not the polar areas. Moreover, in our previous study, we presented a high-accuracy ICESat-2 elevation correction method (Li et al., 2021) and a DEM elevation correction model (Li et al., 2023c). In this study, based on these models, we further optimized and designed an automatic processing scheme for correcting ASTER GDEM elevation. The scheme does not require human intervention during the processing and is suitable for large-scale data processing. It introduces detailed strategies for the various challenges, such as the processing of DEM boundaries, the fusion of different data, and the geographical layout of the ICESat-2 altimeter data. The details of the scheme will provide a meaningful reference for related fields. Meanwhile, a high-accuracy global control point dataset (elevation RMSE: 0.5–3 m within different topographies) (Li et al., 2022; Xie et al., 2021b) and multiple local DEMs (LDEMs) with a high resolution and accuracy were used to validate the accuracy of the original and refined ASTER GDEMs, including comparison of the elevation accuracy in areas with different geolocations, topographic reliefs, and vegetation covers. The related analysis results will provide a beneficial supplement to the accuracy qualification of the latest ASTER GDEM with different geolocations, altitudes, topographies, and vegetation covers. Moreover, ASTER GDEM has been used for more than a decade, and many historical datasets and models are based on its elevation data. The release of the refined ASTER GDEM facilitates seamless integration with historical datasets, which is essential for longitudinal studies examining long-term environmental changes, land use dynamics, and climate impacts.

Throughout the remainder of this paper, we introduce the materials (Sect. 2), present the methodology (Sect. 3), and compare the accuracy of ASTER GDEM before and after elevation correction (Sect. 4). Finally, we draw our conclusions in Sect. 5.

2 Materials

2.1 Satellite laser altimeter data

The ICESat-2 satellite embarked on its mission in 2018 (Martino et al., 2019). This satellite operates in a near-polar orbit, skimming Earth at a low altitude of 496 km, featuring an orbital inclination of 92°, and completing its cycle every 91 d (Martino et al., 2019; Markus et al., 2017). The satellite (Markus et al., 2017) carries a new type of laser altimeter, i.e. a photon-counting laser altimeter, which observes global areas with three pairs of laser beams. These characteristics result in the observations in higher-latitude areas being denser than those in lower-latitude areas. Moreover, compared with previous altimeters, this altimeter only requires low laser transmission energy to observe a profile with a smaller laser footprint along the ground track.

The ATL08 product (Neuenschwander and Pitts, 2019) was used in this study, spanning a survey period from 2018 to 2022.

2.2 GDEM with a medium resolution

ASTER GDEM Version 3 (Abrams et al., 2020) was produced through the utilization of stereo-pair images captured by the ASTER instrument on board the Terra satellite. The DEM covers land areas between 83° N and 83° S, offering a spatial resolution of approximately 30 m. The whole DEM is split into over 20 000 files. The geo-boundary of each file is 1°×1°. Compared with other areas, ICESat-2 laser altimeter observations are denser in the polar areas, making these data particularly suitable for generating polar DEMs. Previous studies have refined ICESat and ICESat-2 altimeter data to produce and release polar DEM products. Moreover, ASTER GDEM does not fully cover polar areas, and correcting elevations in these regions is challenging due to the high variability in ice sheets and flow rates. Therefore, only the parts of the product within the land areas between 83° N and 60° S (except for the polar areas) were corrected.

2.3 Auxiliary data

A global land cover product and a global vegetation index product were used to generate the evaluation attribute set for the GDEM elevation. In this study, the global land cover product was FROM-GLC10, for which the resolution is about 10 m (Gong et al., 2019), and the global vegetation index product was GFCC30TC, which has a resolution of about 30 m (Sexton et al., 2013; García-Álvarez and Lara Hinojosa, 2022). These two data sources were matched with ASTER GDEM by resampling.

2.4 Validation data

Two different kinds of validation data were used in this study, according to the type of survey platform (i.e. satellite and airborne platforms). We adopted this strategy for validating ASTER GDEM before and after more comprehensive elevation correction in order to reduce the impact of biased validation. The first was a high-accuracy global elevation control point dataset (HAGECPD) (Li et al., 2022). The source of the HAGECPD data is ICESat altimeter data with a survey time between 2003 and 2009. In HAGECPD, the elevation control points cover the global land area between 83° N and 60° S. The resolution of these elevation control points is more than a hectometre along the ground track, and the distance separating the ground tracks typically measures around ∼7 km. Figure 1 illustrates the elevation control points distributed across a grid of 1°×1°. The RMSE of these elevation control points is about 0.5 to 3 m. The second source of validation data was LDEM data from around the world, for which the spatial distribution of the data is depicted in Fig. 2. These LDEM data were generated from lidar data collected via an airborne platform or images captured by unoccupied aerial vehicles in conjunction with lidar data. Compared with the first type of validation data, these LDEM data have high resolution and accuracy. The detailed characteristics of the LDEM data are listed in Table 1. These LDEM data were collected according to the following criteria: (1) the data were freely available, (2) the topography and land cover of these data were diverse, and (3) the geolocations of these data were diverse around the world.

Figure 1 Distribution of the elevation control points at 1°×1° (redrawn from Li et al., 2022).

Figure 2 Distribution of the LDEMs with high resolution and accuracy.

Table 1 Main data characteristics of the LDEMs with high resolution and accuracy.

Compared with the second type of validation data, HAGECPD has a lower data density and accuracy but a global coverage. The HAGECPD data were used to validate the accuracy of the corrected ASTER GDEM with the different geographical locations and altitudes at the global scale. The LDEM data were used to validate the accuracy of the corrected ASTER GDEM with the different topographies and vegetation cover as they have a high data density and accuracy and their areas have diverse topographies and vegetation cover.

3 Methodology

The processing flow for correcting GDEM elevation, as shown in Fig. 3, includes four main parts: (1) DEM pre-possessing, (2) construction of the elevation evaluation attribute set, (3) generation of the elevation deviation, and (4) DEM elevation correction.

Figure 3 Main flowchart for GDEM correction.

3.1 GDEM pre-possessing

The elevations of ASTER GDEM were corrected one by one for each 1°×1° DEM file. In the correction, the following challenges needed to be considered: (1) there can be insufficient ICESat-2 altimeter data in some DEM files, especially in the land margin areas with a low proportion of effective land area, which can reduce the performance of the correction model within these areas; (2) elevation continuity between adjacent DEM elevations should be ensured after correction; and (3) DEM boundary pixels lack enough adjacent information, which can affect the accuracy of the evaluation attribute set of these boundary pixels. To address these challenges, it was necessary to expand the area around the DEM file to be corrected. We adopted a strategy to expand the boundaries of the central DEM (the processing DEM) by using data from neighbouring DEMs. As shown in Fig. 3a, the expanded area includes two types: the modelling area and the data buffer area. The modelling area was used to address challenges (1) and (2), and the data buffer data were used to address challenge (3). The calculations for the modelling and data buffer areas were as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{c}{\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}\\ {\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}}\end{array}\right]=R_{\text{ref}}\cdot \left(\left[\begin{array}{c}{r}_{\text{DEM}}^{\text{ul}}\\ c_{\text{DEM}}^{\text{ul}}\end{array}\right]-k\cdot L\right),\text{(2)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{c}{\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}}\\ {\mathrm{\Lambda }}_{\mathrm{m}}^{\text{max}}\end{array}\right]=R_{\text{ref}}\cdot \left(\left[\begin{array}{c}{r}_{\text{DEM}}^{\text{dr}}\\ c_{\text{DEM}}^{\text{dr}}\end{array}\right]+k\cdot L\right),\text{(3)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{c}{\mathrm{\Phi }}_{\mathrm{b}}^{\text{min}}\\ {\mathrm{\Lambda }}_{\mathrm{b}}^{\text{max}}\end{array}\right]=R_{\text{ref}}\cdot \left(\left[\begin{array}{c}{r}_{\text{DEM}}^{\text{dr}}\\ c_{\text{DEM}}^{\text{dr}}\end{array}\right]-k\cdot L-\mathrm{\Delta }L\right),\text{(4)}& {\displaystyle }& {\displaystyle }\left[\begin{array}{c}{\mathrm{\Phi }}_{\mathrm{b}}^{\text{min}}\\ {\mathrm{\Lambda }}_{\mathrm{b}}^{\text{max}}\end{array}\right]=R_{\text{ref}}\cdot \left(\left[\begin{array}{c}{r}_{\text{DEM}}^{\text{dr}}\\ c_{\text{DEM}}^{\text{dr}}\end{array}\right]+k\cdot L+\mathrm{\Delta }L\right),\text{(5)}& {\displaystyle }& {\displaystyle }L=max\left(\left(r_{\text{DEM}}^{\text{dr}}-r_{\text{DEM}}^{\text{ul}}\right),\left(c_{\text{DEM}}^{\text{dr}}-c_{\text{DEM}}^{\text{ul}}\right)\right),\end{array}$$

where $[r_{\text{DEM}}^{\text{ul}},c_{\text{DEM}}^{\text{ul}}{]}^{\text{T}}$ and $[r_{\text{DEM}}^{\text{ul}},c_{\text{DEM}}^{\text{ul}}{]}^{\text{T}}$, respectively, represent the pixel coordinates of the left vertex and bottom-right pixel coordinates of the DEM to be corrected; R_ref represents the transformation matrix of pixel coordinates to geographic coordinates for the DEM to be corrected, which can be obtained directly from the DEM file; and k and ΔL represent the pixel expansion coefficient and pixel expansion increment, respectively. The initial setting of k was 0.1. ΔL was set to 10, which could ensure that the boundary pixels of the modelling area had sufficient adjacent information.

3.2 Generation of the GDEM elevation evaluation attribute set

The sources for the attribute set included ASTER GDEM (elevation and quality data), the land cover data, and the vegetation cover data. However, there were inconsistencies in the resolution of these data sources. Before constructing the attribute set, it was necessary to match each ASTER GDEM pixel with the other data based on the criterion of the nearest geographical location.

After matching, each pixel $p_i({\mathit{\lambda }}_i,{\mathit{\phi }}_i,h_i,q_i,{\mathit{\eta }}_i,{\mathit{\zeta }}_i)$ (called an attribute pixel) contained the geographical location (longitude λ_i and latitude φ_i), elevation h_i, elevation quality q_i, vegetation cover index ζ_i, and land cover category η_i. Then, for each DEM of ASTER GDEM, the corresponding attribute pixel sets P_b and P_m were extracted according to the boundaries of the modelling and data buffer areas, respectively, as shown in the following formulas:

$$\begin{array}{}\text{(6)}& {\displaystyle }{P}_{\mathrm{b}}& {\displaystyle }=\bigcup _{i=\mathrm{1}}^N{p}_i\left({\mathit{\lambda }}_i\in \left[{\mathrm{\Lambda }}_{\mathrm{b}}^{\text{min}},{\mathrm{\Lambda }}_{\mathrm{b}}^{\text{max}}\right],{\mathit{\phi }}_i\in \left[{\mathrm{\Phi }}_{\mathrm{b}}^{\text{min}},{\mathrm{\Phi }}_{\mathrm{b}}^{\text{max}}\right]\right),\text{(7)}& {\displaystyle }{P}_{\mathrm{m}}& {\displaystyle }=\bigcup _{i=\mathrm{1}}^N{p}_i\left({\mathit{\lambda }}_i\in \left[{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}},{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{max}}\right],{\mathit{\phi }}_i\in \left[{\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}},{\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}\right]\right),\end{array}$$

where N represents the total pixel number of the DEM to be processed and its neighbouring DEM. According to Earth's surface information, $P_{\mathrm{b}}^{\text{nw}}$ (P_b within the water-free area) and $P_{\mathrm{m}}^{\text{nw}}$ (P_m within the water-free area) were extracted, as shown in the following formulas:

$$\begin{array}{}\text{(8)}& {\displaystyle }{P}_{\mathrm{b}}^{\text{nw}}& {\displaystyle }=P_{\mathrm{b}}\cap P_{\text{nw}},\text{(9)}& {\displaystyle }{P}_{\mathrm{m}}^{\text{nw}}& {\displaystyle }=P_{\mathrm{m}}\cap P_{\mathrm{b}}^{\text{nw}},\text{(10)}& {\displaystyle }{P}_{\text{nw}}& {\displaystyle }=\bigcup _{i=\mathrm{1}}^N{p}_i(h_i\ne \mathrm{0},{\mathit{\eta }}_i\ne \mathrm{60},{\mathit{\zeta }}_i\ne \mathrm{200}),\end{array}$$

where P_nw represents the set of attribute pixels within the water-free area of the DEM to be processed and its neighbouring DEM. Then, according to our previous study (Li et al., 2023c), the attribute pixels in $P_{\mathrm{b}}^{\text{nw}}$ were used to generate the evaluation features for the topography, and those in $P_{\mathrm{m}}^{\text{nw}}$ were used to generate the other evaluation features. The evaluation attribute set in the model includes topography, surface coverage, spatial distribution, and data source quality (Li et al., 2023c).

Table 2 Characteristics of the input and output data.

3.3 Generation of the GDEM elevation deviation

The accuracy of ICESat-2 altimeter data significantly exceeds that of ASTER GDEM data. ICESat-2 altimeter data can be used as the source of the accurate elevation references (i.e. seed points) for ASTER GDEM to generate the training samples. For the training samples, the seed points can be expected to be of better quality, followed by a more even distribution, and of greater quantity. However, there is a contradiction between these indicators. At the global scale, it was necessary to take comprehensive account of these indicators in each DEM area of ASTER GDEM. For this comprehensive consideration, in our prior study, we presented a method for deriving accurate elevations from ICESat-2 altimeter data with the comprehensive evaluation labels. We then analysed the performance of each evaluation label and discussed a balanced strategy between the quality and distribution or quantity of the extracted elevations (Li et al., 2021). Here, based on our previous study, we extracted two types of seed points from the ICESat-2 altimeter data through different combinations of evaluation labels. The first type of seed point (FYSP) was extracted by a comprehensive evaluation label to ensure as high an elevation quality as possible. The second type of seed point (SYSP) was extracted by partial assessment labels (outliers and atmospheric assessments) to ensure good elevation quality and good quantity and distribution. According to the two types of seed points, there are two strategies for adjustment of the adoption: (1) the adjustment strategy of the seed point number (ASSPN) and (2) the adjustment strategy of the seed point distribution (ASSPD).

For ASSPN, it was necessary to consider two situations for each 1°×1° DEM: (1) the quantity of ICESat-2 altimeter data increases with the increase in latitude, and (2) the effective land area is different at the same latitude. To this end, the type of seed point to adopt was related to the maximum latitude ${\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}$ and the effective land proportion of the modelling area r_land, as shown in the following equations. If the number of FYSP was higher than N_S, FYSP was adopted; otherwise, SYSP was adopted.

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{N}_{\mathrm{S}}=r_{\text{land}}\cdot \cos \left({\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}\right)\cdot N_{\mathrm{o}}\text{(12)}& {\displaystyle }& {\displaystyle }{r}_{\text{land}}={\displaystyle \frac{\left|P_{\mathrm{m}}^{\text{nw}}\right|}{\left|P_{\mathrm{m}}\right|}}\end{array}$$

N_o represents the original number threshold. ASSPD can be expected to obtain a better spatial distribution of seed points and is a supplement to ASSPN. ASSPD is used for a situation in which ASSPN adopts FYSP as SYSP has a good quantity and distribution. ASSPD quantifies the distribution of the first type of seed point by dividing each modelling area into sub-regions, and then it uses the second type of seed point to replace the regions with a low number of seed points. As shown in the following formulas, each modelling area was evenly divided into sub-regions (30×30 was adopted in this study) to generate a set P_I where a statistic results from the number of seed points in each sub-region.

$$\begin{array}{}\text{(13)}& {\displaystyle }& {\displaystyle }{P}_{\mathrm{I}}=\bigcup _{i=\mathrm{1}}^{n_{\mathrm{I}}}\left({\mathrm{\Gamma }}_i\right)\text{(14)}& {\displaystyle }& {\displaystyle }{\mathrm{\Gamma }}_i=\sum _{j=\mathrm{1}}^n∥{\displaystyle \frac{{\mathit{\gamma }}_j^i}{\left|{\mathit{\gamma }}_j^i\right|}}∥\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_j^i={\mathit{\tau }}_j\left({\mathit{\lambda }}_j\in \left[{\mathrm{\Lambda }}_i^{\text{min}},{\mathrm{\Lambda }}_i^{\text{max}}\right),{\mathit{\phi }}_j\in \left[{\mathrm{\Phi }}_i^{\text{min}},{\mathrm{\Phi }}_i^{\text{max}}\right)\right)\text{(16)}& {\displaystyle }& {\displaystyle }{\mathrm{\Lambda }}_i^{\text{min}}={\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}}+(i-\mathrm{1})\cdot {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{max}}-{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}}}{\sqrt{n_{\mathrm{I}}}}}\text{(17)}& {\displaystyle }& {\displaystyle }{\mathrm{\Lambda }}_i^{\text{max}}={\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}}+i\cdot {\displaystyle \frac{{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{max}}-{\mathrm{\Lambda }}_{\mathrm{m}}^{\text{min}}}{\sqrt{n_{\mathrm{I}}}}}\text{(18)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_i^{\text{min}}={\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}}+(i-\mathrm{1})\cdot {\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}-{\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}}}{\sqrt{n_{\mathrm{I}}}}}\text{(19)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_i^{\text{max}}={\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}}+i\cdot {\displaystyle \frac{{\mathrm{\Phi }}_{\mathrm{m}}^{\text{max}}-{\mathrm{\Phi }}_{\mathrm{m}}^{\text{min}}}{\sqrt{n_{\mathrm{I}}}}}\end{array}$$

After obtaining P_I, the median N_median (except for zeros, thus reducing the influence of the blank regions of seed points) was calculated. The parts P_add from which the number of seed points was less than half of N_median were then identified, as shown below. Finally, in the sub-regions belonging to P_add, SYSP was replaced with FYSP.

$$\begin{array}{}\text{(20)}& {\displaystyle }& {\displaystyle }{N}_{\text{median}}=M\left(P_{\mathrm{I}}\right({\mathrm{\Gamma }}_i\ne \mathrm{0},i=\mathrm{1},\mathrm{2},\mathrm{\dots },N_{\mathrm{I}}\left)\right)\text{(21)}& {\displaystyle }& {\displaystyle }{P}_{\text{add}}\in \left\{P_{\mathrm{I}}<{\displaystyle \frac{N_{\text{median}}}{\mathrm{2}}}\right\}\end{array}$$

M(⋅) represents the median calculation.

Finally, all of the collected seed points matched the ASTER GDEM elevations within the modelling area, based on the criterion of the nearest geographical location, and the matched elevations were used to obtain the elevation deviations of ASTER GDEM.

3.4 GDEM elevation correction

For each DEM file of ASTER GDEM, a random forest regression algorithm was used to train its elevation correction model. The random forest regression function (fitrensemble) of the MATLAB platform was directly used for its processing efficiency and compatibility. The function can be viewed and downloaded from https://ww2.mathworks.cn/help/stats/fitrensemble.html (last access: 1 June 2024). The function method was selected for bootstrap aggregation (bagging and random forest) (Breiman, 2001). For this method, we adopted the recommended or default setting values for the parameter selection, including the tree number (100) and learners (tree). Moreover, when the number of training samples is too low, the area (i.e. the modelling area) needs to be expanded for selecting training samples in order to ensure that there are sufficient training samples. If this number is less than 100, it is necessary to expand the modelling area by adjusting k, as shown in the following formula, and then repeating as in Fig. 3a, i.e. “DEM prepossessing”.

$$\begin{array}{}\text{(22)}& k=min\left(\mathit{\right\{}\mathrm{0.1},\mathrm{0.2},\mathrm{\dots }\mathrm{1}\mathit{\}},L+\mathrm{2}\cdot (L-\mathrm{\Delta }L\left)\right)\end{array}$$

After obtaining the model, all the DEM elevations (i.e. z components) were corrected by the model, and then a new DEM file was generated with the same format.

4 Results and discussion

4.1 Release of the corrected ASTER DEM product

Based on the presented scheme, we corrected the ASTER GDEM elevation with ICESat-2 altimeter data and then stored the corrected ASTER GDEM product (IC2-GDEM) in GeoTIFF format (.tif) with the same projection and datum as ASTER GDEM (Xie et al., 2024). The resolution of IC2-GDEM is about 30 m (grid). The survey date of its main source data is between 2000 and 2013, and its elevation correction source is ICESat-2 laser altimeter data with a survey date from 2018 to 2022. More details about the data source of IC2-GDEM are listed in Table 2. The IC2-GDEM product has been openly released via the National Tibetan Plateau Data Center (DOI: 10.11888/RemoteSen.tpdc.301229).

Figure 4 shows the distribution of the IC2-GDEM product. Figure 5 shows the average of elevation corrections of IC2-GDEM relative to ASTER GDEM within the global 1°×1° grid. IC2-GDEM covers the global land area between 83° N and 60° S, except for the polar regions. Under this coverage, approximately 99.98 % of the ASTER GDEM elevation has been corrected, and the remaining ASTER GDEM elevation has not been corrected, mainly because of the lack of sufficient ICESat-2 seed points. That is, the ICESat-2 seed points did not satisfy the constraints of the presented methodology. The uncorrected ASTER GDEM elevation areas are mainly parts of islands and reefs where the geolocation is near low-latitude regions.

Figure 4 The ASTER GDEM product after ICESat-2 altimeter data correction.

Figure 5 The average of elevation corrections of IC2-GDEM relative to ASTER GDEM within the global 1°×1° grid.

4.2 Validation by the high-quality elevation control point dataset at the global scale

By using the HAGECPD data (survey date: 2003–2009), we validated and then contrasted the accuracy of the original and corrected ASTER GDEM elevations within various continents, including Asia, Africa, South America, Oceania, North America, and Europe. Figure 6 provides the statistical results.

Figure 6 Comparison of the accuracies of the original and corrected ASTER GDEM elevations within different continents.

The results from the comparison indicate that the accuracy of the corrected ASTER GDEM elevation is higher than that of the original GDEM elevation. In the different continents, the RMSEs of the corrected ASTER GDEM elevation are about 2 to 4.5 m, and those of the original ASTER GDEM elevation are about 8 to 10 m. After the elevation correction, the errors of the ASTER GDEM elevation are reduced by more than 45 %, and the maximum reduction ratio exceeds 70 %. There are two main reasons for this inconsistency in elevation accuracy improvement among the different continents. The first is that there are significant differences in the quality of ASTER GDEM across the continents. The second is that the topographic relief and the land covers across the continents are different. These influencing factors are also important evaluation attributes in the correction model of DEM elevation. Moreover, there is an obvious difference in the corrected ASTER GDEM elevation errors for the different continents. The elevation errors for the continents in order from lowest to highest are for Oceania, South America, Africa, North America, Asia, and Europe. The order of the original ASTER GDEM elevation errors is similar to the above order. To further evaluate the correction performance for the various geolocations, we scrutinized the accuracy of the original and corrected ASTER GDEM elevations in different longitude and latitude regions. The results from the scrutiny are displayed in Figs. 7 and 8.

Figure 7 Comparison of the errors of the original and corrected ASTER GDEM elevations within different longitude regions.

Figure 8 Comparison of the errors of the original and corrected ASTER GDEM elevations within different latitude regions.

From Fig. 7, in the different longitude regions, it is apparent that, compared with the original ASTER GDEM elevation, the systematic error and standard deviation after ASTER GDEM elevation correction are significantly reduced. In particular, for the original ASTER GDEM elevation, there is a large difference in the systematic error. The systematic error of the original ASTER GDEM elevation in the low-longitude regions is smaller than that in the high-longitude regions. The systematic error after ASTER GDEM elevation correction is close to zero in the different longitude regions.

Similar to the analysis results in Fig. 7, in the different latitude regions, it is also apparent that the systematic error and standard deviation after ASTER GDEM elevation correction are significantly reduced compared with the original ASTER GDEM elevation, as displayed in Fig. 8. Notably, the data density of the ICESat-2 laser altimeter data varies greatly in different latitude regions, and the data density of the ICESat-2 laser altimeter data in the low-latitude regions is generally less than that in the high-latitude regions (Neuenschwander and Pitts, 2019; Markus et al., 2017). However, the errors of the ASTER GDEM elevation corrections in low-latitude regions (especially near the Equator) are no larger than those in other latitude regions, which indicates that the correction results are less affected by the difference in the data density of the ICESat-2 laser altimeter data. There are two main reasons for the observed phenomenon in DEM elevation correction. First, the DEM elevation correction model generalizes well across various spatial distributions (Li et al., 2023c). Second, multiple strategies have been designed for extracting seed points from ICESat-2 laser altimeter data during the generation of GDEM elevation deviation. These strategies ensure an even distribution of seed points, thereby mitigating the impact of varying densities on the ICESat-2 laser altimeter data. Moreover, the errors of the ASTER GDEM elevation corrections are larger in high-latitude regions, especially in the Northern Hemisphere. The primary cause of this phenomenon is that the original ASTER GDEM elevation errors in the high-latitude regions are significantly larger than those in the other regions.

4.3 Validation by LDEMs

By using the LDEM data from around the world (as shown in Fig. 2), we validated the errors of the ASTER GDEM elevation before and after correction and then qualified the error reduction ratio of the corrected ASTER GDEM elevation. Table 3 displays the results.

Table 3 Comparison of the errors of the original and corrected ASTER GDEM elevations in the different LDEM areas.

From Table 3, it can be seen that the errors of the ASTER GDEM elevation after correction are significantly reduced. For all of the validation areas, the average RMSE of the original ASTER GDEM is 9.38 m and that of the corrected ASTER GDEM is 3.37 m, and the reduction ratio of the elevation error of the corrected ASTER GDEM is 64.05 % when compared with the original ASTER GDEM. For the individual validation areas, the RMSEs of the original ASTER GDEM are between 3.82 and 19.37 m, and those of the corrected ASTER GDEM are between 1.63 and 12.51 m. Compared with the original ASTER GDEM, the reduction ratio of the elevation error of the corrected ASTER GDEM ranges from 16.14 % to 82.49 %, corresponding to an average of 47.66 %.

Figure 9 provides a more detailed evaluation of the error distribution of the original and corrected ASTER GDEM elevations. In each plane, the red and orange histograms represent the error distribution of the original and corrected ASTER GDEMs. From Fig. 9, it can be seen that the error distributions of the corrected ASTER GDEM elevation show better symmetry at zero error than those of the original ASTER GDEM elevation. This phenomenon indicates that the systematic error of the original ASTER GDEM elevation has been corrected well. Meanwhile, the statistical dispersion of the error distribution of the corrected ASTER GDEM elevation is significantly smaller than that of the original ASTER GDEM elevation.

Figure 9 Comparison of the original and corrected ASTER GDEM elevation errors.

Moreover, there are differences between these statistical descriptions. This phenomenon is mainly caused by differences in the topography and Earth's surface. To confirm this, we combined other data, including LDEM data, vegetation cover index data (GFCC30TC), and land cover data (FROM-GLC10), to evaluate the influence of these two factors. Results from the analysis are provided in Figs. 10 and 11. For terrain, we used the common topographic indicator – slope – to quantify the topographic relief. The slope was obtained through ArcGIS processing. After the slope was obtained, we first divided it into different slope intervals and then calculated the errors after ASTER GDEM elevation correction in each interval. Meanwhile, the errors of the original ASTER GDEM elevation were added to the statistical analysis to further supplement the investigation of the correction results. To examine the different surfaces of Earth, we first divided it into urban, bare land, and vegetation areas according to the land cover data and then carried out statistical analysis of the errors of the original and corrected ASTER GDEM elevation in each area. We then carried out sub-division and statistical analysis of the vegetation areas through the vegetation cover data.

From Fig. 10, it is apparent that, in the different slope regions, the errors of the corrected GDEM elevation are significantly reduced when compared with those of the original GDEM elevation, which indicates that the corrected results have a better elevation quality in areas with different topographic reliefs. Meanwhile, this kind of error decrease shows a weakening trend with the increase in the slope. There are two main reasons for this trend. The first reason is that the quality of the data from the ICESat-2 laser altimeter, which were used as the seed points, generally decreases as the slope increases. The second reason is that the accuracy of the original ASTER GDEM elevation also decreases as the slope increases, as shown in Fig. 10.

Figure 10 Comparison of the original and corrected ASTER GDEM elevation errors with different slopes.

Figure 11 Comparison of the original and corrected ASTER GDEM elevation errors with different land covers.

From Fig. 11, in the three different Earth surface categories, it can be seen that the errors of the corrected GDEM elevation are all significantly smaller than those of the original GDEM elevation. For the different land covers, the corrected GDEM elevation in the bare land areas shows the best elevation correction. This means that the DEM elevation correction model can reduce the impact caused by the inherent noise of optical data. Compared with the other two types of ground objects, the bare land area is not shielded by ground objects, which is conducive to optical observation of ground texture (Li et al., 2023c; Pham et al., 2018) and laser detection of ground elevation (Li et al., 2023a; Neuenschwander and Pitts, 2019). Therefore, there is little need to consider the influence of ground objects on terrain elevation in bare land areas. For the other two types of land cover, the corrected GDEM elevation quality in the areas with low vegetation cover (less than 20 % of the vegetation cover index) is comparable to that in urban areas, but with the increase in vegetation cover, the difference in the corrected GDEM elevation quality between the two types of land cover gradually increases. Moreover, relative to the original GDEM elevation, the quality improvement of the corrected GDEM elevation shows a reduced trend with the increase in vegetation cover. Similar to the analysis of topographic relief, there are two main reasons for this phenomenon: the quality of the ICESat-2 laser altimeter data and the accuracy of the original ASTER GDEM elevation both decrease as the vegetation cover increases.

The absence of LDEM data in Asia and Africa (difficult to obtain for free) is a noteworthy consideration. The lack of such data in these continents could introduce uncertainties, as the validation relies on the available LDEM data. To reduce the impact of these uncertainties, we also used the other satellite laser altimeter data with global coverage, i.e. ICESat altimeter data, as another source of the validation data in this study. The resolution and accuracy of the ICESat altimeter data may not match those of the LDEM data. This means that the validation of finer topographic details in Asia and Africa may not be fully accounted for. Despite this gap, the similar validation accuracy between the two data is advantageous, particularly given the average reduction ratios close to 47 %. This similarity suggests that ICESat altimeter data can serve as a reliable alternative and a beneficial supplement for areas lacking LDEM data.

4.4 Limitations and future work

From the above analysis, while IC2-GDEM demonstrates significant improvements in elevation accuracy, its quality of elevation correction is notably limited in areas with steep slopes or a high vegetation cover index. For research or applications with high-quality elevation requirements, this limitation can be quantified and identified through topographic relief calculations or in combination with vegetation cover data. Moreover, the temporal differences between the ICESat-2 data survey and ASTER GDEM collection may lead to elevation inconsistencies, especially in extremely dynamic landscapes (e.g. coastal erosion areas). To this end, for the validation areas, we selected multiple areas located near the coast to further evaluate the performance of DEM elevation correction in this kind of landscape, including the coast of the Netherlands, the coast of Australia, the islands of New Zealand, the western coast of the United States, and the island of Hawaii, as shown in Fig. 2. Validation results (Fig. 9) show that the corrected DEM elevations still exhibit significant accuracy improvements. This phenomenon is due to the temporal discrepancy between the datasets being generally smaller than the ASTER GDEM elevation error (about 10 m). To further mitigate these inconsistencies, two main strategies are proposed. First, using repeated ICESat-2 observations to analyse areas with potential elevation changes can assess the application potential of IC2-GDEM in dynamic landscapes. Second, integrating IC2-GDEM with other elevation data collected closer to the ASTER GDEM collection time can reduce analysis deviations in dynamic landscapes. In future work, we will identify areas with elevation changes using the above strategies and integrate other GDEMs to further enhance IC2-GDEM elevation quality in dynamic landscapes.

5 Data availability

All the data used in this study are open-source. IC2-GDEM is available at https://doi.org/10.11888/RemoteSen.tpdc.301229 (Xie et al., 2024). The ICESat-2 ATL08 product is available at https://doi.org/10.5067/ATLAS/ATL08.005 (Neuenschwander et al., 2021). The ASTER GDEM product is available at https://doi.org/10.5067/ASTER/ASTGTM.003 (NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team, 2019). The global land cover data (FROM-GLC10) are available at https://data-starcloud.pcl.ac.cn/zh/resource/1 (Gong et al., 2019). The vegetation cover index data (GFCC30TC v003) are available at https://doi.org/10.5067/MEaSUREs/GFCC/GFCC30TC.003 (Townshend, 2016). The HAGECPD data are available at https://doi.org/10.11888/Geogra.tpdc.271727 (Xie et al., 2021a). The high-resolution and accurate LDEM data are available at https://data.linz.govt.nz/data (Land Information New Zealand, 2023), https://coast.noaa.gov/dataviewer/#/lidar/search (NOAA Office for Coastal Management, 2023), https://www.ahn.nl/ahn-viewer (Ahn Viewer, 2023), https://www.ign.gob.ar/NuestrasActividades/Geodesia/ModeloDigitalElevaciones/Mapa (Instituto Geográfico Nacional, 2023) and https://doi.org/10.26186/89644 (Geoscience Australia, 2015).

6 Conclusion

In this study, we introduce a new open-source dataset named IC2-GDEM. IC2-GDEM is generated by directly refining the ASTER GDEM elevation with ICESat-2 altimeter data, FROM-GLC10 data, and GFCC30TC data. This strategy provides a cost-effective way of improving elevation accuracy, which is especially beneficial for developing regions with limited resources for high-accuracy DEMs. After validation at the global scale, IC2-GDEM presents superior elevation quality and application potentials: (1) it is less affected by the difference in the data density of the ICESat-2 laser altimeter data; (2) there is an obvious improvement in the elevation quality for the different topographies, even in mountainous areas with higher slopes; (3) a clear improvement in elevation quality was observed across various land covers, including areas with dense vegetation cover; (4) it is expected to promote seamless integration with the historical datasets of ASTER GDEM for longitudinal studies of long-term environmental changes, land use dynamics, and climate impacts; and (5) it will be a new complementary data source for other GDEMs, such as SRTM and Copernicus DEM, to cross-validate qualities, fill data gaps, and conduct multi-scale analyses.