REMA was constructed from stereoscopic imagery collected by four commercially operated satellites: WorldView-1, WorldView-2, WorldView-3, and GeoEye-1, launched in 2007, 2009, 2014, and 2008, respectively (Table 1). These satellites are operated by DigitalGlobe Inc. and their images are distributed via the Polar Geospatial Center (PGC) under a scientific use licensing agreement with the US National Geospatial-Intelligence Agency (NGA). While the images are only available to US federally funded investigators, derived products, including DEMs, may be openly distributed. The push-broom sensors aboard these polar orbiting satellites provide optical imagery with pixel ground resolutions of less than 0.5 m in the panchromatic band. Their camera pointing capabilities allow them to obtain overlapping images from different look angles, yielding convergence angles between image pairs that are appropriate for stereophotogrammetric DEM extraction. Using only the rational polynomial coefficients (RPCs) derived from satellite positioning, these DEMs may have translational errors (biases) of several meters. These can be reduced through ground control registration to a point-to-point (relative) error of 20 cm (Noh and Howat, 2015; Shean et al., 2016), which is comparable to the uncertainty of airborne lidar. Importantly, unlike other common stereo-capable imagers, such as from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), the high spatial and radiometric resolutions of these imagers enable high-quality elevation extraction over low-contrast surfaces such as snow cover and ice sheet interiors/accumulation zones (Noh and Howat, 2015; Shean et al., 2016).

The satellite imagers targeted swaths of the ground surface in each orbital pass, creating image strips that are up to 200 km long and between 11 and 17 km wide, depending on the sensor (Table 1). To ease data handling, the data provider divided each strip into approximately square subsets, or scenes, with ∼10 % overlap, prior to delivery. Pairs of strips covering the same ground area were selected for use as DEM stereopairs if they had a convergence angle greater than 10^∘ (Hasewaga et al., 2000), a difference in sun elevation angles of less than 10^∘ and a time separation of less than 10 days to reduce the likelihood of errors due to surface change, such as snowfall or ice motion. Pairs of images from the same sensor and different sensors were used. Through off-nadir camera pointing, we were able to successfully obtain DEMs over flat surfaces up to approximately 88^∘ S. All REMA products are in polar stereographic projection, with a central meridian of 0^∘ and standard latitude of −71^∘ S and are referenced to a WGS84 ellipsoid. A flow chart of the REMA workflow is provided in Fig. S1 of the Supplement.

2.2 DEM strip quality control and registration

Hillshade representation images were generated for each DEM strip segment and these were visually inspected and classified based on visual quality (i.e., lack of erroneous surfaces due to clouds, shadows). Such erroneous surfaces appear as chaotic textures in the hillshade image that contrast with the actual topography. DEMs were accepted if no erroneous surfaces were identified in the hillshade image, manually edited to mask erroneous surfaces where errors were small and isolated, or rejected if errors were too extensive to be edited. Of the 187 585 strip DEM segments, 130 386 (69 %) were visually inspected and classified. The remaining 31 % of strips were not visually inspected because we switched from inspecting every strip to only inspecting strips needed for a single mosaic coverage part way through the quality control process. This resulted in fewer inspected strips for regions inspected after this change in procedure. Of the visually inspected strips, 43 550 (33 %) passed quality control, with 19 971 (15 %) requiring manual masking. In total, the 55 491 482 km² of quality-controlled DEMs cover an area of 13 567 969 km², or 95 % of the continent (Fig. 1b).

All strips were vertically registered to altimetry point clouds obtained from CryoSat-2 radar and ICESat-1 Geoscience Laser Altimeter System (GLAS) 2-D campaign (25 November to 17 December 2008). The ICESat-1 GLAS covers all areas of Antarctica north of 86^∘, with that limit, known as the pole hole, due to orbital constraints. We use version 34 of the Level 2 GLAH12 altimetry data distributed by the National Snow and Ice Data Center (Zwally et al., 2014). The ground footprint of each laser shot has a diameter of approximately 70 m and an accuracy over flat surfaces of ±0.1 m with small variations due to snow surface properties (Shuman et al., 2006). Launched in April, 2010, CryoSat-2 carries the KU-band SAR/Interferometric Radar Altimeter (SIRAL) instrument with along- and across-track resolutions of 450 m and 1 km, respectively, in its higher-resolution interferometric (SARIn) mode. CryoSat-2 registration points were obtained from the point of closest arrival (POCA) locations in SARIn mode derived using a slope-threshold retracker (Gray et al., 2017). We use only the SARIn mode data and not the low-resolution-mode (LRM) measurements because we did not feel confident that, over the scale of a DEM strip, the slope-driven error in LRM elevations would reliably average to zero. Each DEM strip was smoothed and down-sampled to a 32 m grid spacing, filtered to remove rough terrain, and then interpolated to the CryoSat-2 SARIn-mode point cloud locations. For CryoSat-2 registrations, we estimated the linear temporal trend in the surface height from the time series of all points within each DEM, so that each altimetric point measurement would provide an estimate of the surface height at the time of DEM acquisition. We did not apply a similar time-dependent correction to the ICESat-1 data because the time span between the altimetry measurements and the DEM was much larger. Further, we only use ICESat-1 data in the absence of quality CryoSat-2 SARIn mode data, which is predominantly in the slow-flowing interior of the ice sheet where changes in surface elevation are expected to be less than the resolution of repeat surface height observations on sub-decadal timescales (e.g., Helm et al., 2014).

Figure 2Plot of vertical bias (i.e., the mean of residuals) between REMA strip DEMs and overlapping point clouds from ICESat-1 laser altimetry and CryoSat-2 radar altimetry. Only strips with standard deviations in residuals less than 0.25 m are plotted. The solid line is unity, shifted by the mean difference among biases (0.39 cm).

Download

The median difference between the DEMs and the altimeter point clouds provides an estimate of the DEM's vertical bias. For CryoSat-2 data, only vertical bias corrections with a 1σ uncertainty of less than 0.1 m and residuals with a standard deviation of less than 1 m were used in mosaicking. For ICESat-1, we impose a lower maximum threshold in the standard deviation of the residuals of 0.35 m because such strips were mostly used over the flatter interior terrain of CryoSat-2's LRM coverage. A total of 6 679 897 km² is covered by CryoSat-2 registered DEMs, or 29 901 958 km² including repeat coverage (Fig. 1c), with registered ICESat DEMs covering 4 897 600 km², including 8 739 128 km² of repeat coverage (Fig. 1d).

Strips with both CryoSat-2 and ICESat-1 registration within the bias correction uncertainty thresholds allow for an estimate of the biases in CryoSat-2 height estimates due to the penetration of microwaves into the snow and firn layer (i.e., the penetration depth), or biases due to the retracking algorithm (i.e., where the retracker identifies a point on the leading edge of the waveform that does not correspond perfectly to the surface). Such biases are assumed negligible for the 1064 nm wavelength pulse of ICESat-1's laser altimeter and, therefore, the difference between the ICESat-1 and CryoSat-2 bias corrections should give an estimate of the CryoSat-2 bias. Figure 2 plots the vertical bias corrections from ICESat-1 and CryoSat-2 for 227 strips for which standard deviations of residuals were less than 0.25 cm. These strips were distributed across the entire area of CryoSat-2 SARIn coverage and, therefore, the mean difference between CryoSat-2 and ICESat-1 bias corrections should not be sensitive to local variability in surface elevation change over the period between the two missions. The mean difference between the two corrections is $-\mathrm{0.39}\pm \mathrm{0.35}$ m. We expect the bias in the CryoSat-2 data to depend on surface density and surface slope (Wang et al., 2015), but we do not have a straightforward way of predicting the bias, and we did not find a clear spatial or elevational dependence of the CryoSat-2–ICESat differences. Therefore, we added a uniform value of 0.39 m to the CryoSat-2-registered heights, regardless of the location of the strips and the surface type.

Quality-controlled strip DEMs were mosaicked into 100 km by 100 km tiles with a 1 km wide buffer on each side to enable coregistration and feathering between tiles. For each tile, strips with altimetry registration were added first, in order of ascending vertical error, with a linear distance-weighted edge feather applied to the strip boundaries. The error value at each pixel is the standard error from the residuals of the registration to altimetry, and the date stamp is the day of DEM acquisition. The ±0.35 m errors in bias for CryoSat-2 registered tiles were not included in this error estimate. In areas where edges of strips have been feathered, the error and date stamp are averaged with the same weighting as the elevation. Once all registered strips were added, unregistered strips were added to fill gaps and are coregistered to the existing registered data in the mosaic. Each quality-controlled, unregistered strip that overlaps a data gap was tested for the precision of three-dimensional coregistration, using the Nuth and Kääb (2011) algorithm, with the strip with the smallest coregistration error, defined as the root mean square of the elevation difference between the mosaic and the coregistered DEM, selected to fill the gap with the coregistration offset applied. Again, a distance-weighted feathering was applied to smooth strip edges.

If CryoSat-2 registered data were available within a tile, those data were used, and any ICESat-1-registered data were ignored. If neither CryoSat-2 or ICESat registered data were available, the quality-controlled strip with the most coverage of the tile was added first and served as a relative reference. Unregistered strips were then coregistered to the mosaic and added as described above. Figure 3 shows the distribution of tiles registered to CryoSat-2, ICESat-1, or alignment to neighbors.

Figure 3Map of the registration data source for each 100 km by 100 km REMA tile.

Download

Tiles around the edge of the ice sheet and within the zone of CryoSat-2 SARIn mode collection were mostly registered to contemporaneous CryoSat-2 altimetry, with the exception of coastal tiles with too little land surface or with extensive crevassing that prevented successful altimetry registration. Most of the interior tiles were registered to ICESat-1 and therefore have a nominal date stamp of late December 2008, although little or no secular surface elevation change is expected in these regions on sub-decadal timescales (e.g., Helm et al., 2014). Some tiles that were missing registration, and thus registered through alignment to neighboring tiles, are found around the pole hole and along a narrow zone to its northeast. In most cases, the lack of registration was caused by a registration error larger than the thresholds defined in Sect. 2.3, likely because of the extreme off-nadir angles required for the satellites to acquire stereo imagery in the far south. Tile edges were feathered to smooth any offsets.

Figure 4Maps of REMA (a) elevation error, obtained from the root mean square of the differences in elevation between the DEM and altimetry data following registration, or the differences among co-registered DEMs in the case of alignment (note the logarithmic color scale), and (b) date stamp obtained from the date of image acquisition.

Download

Finally, we applied a coastline mask using the British Antarctic Survey (BAS) land–ice classification polygons (https://add.data.bas.ac.uk/, last access: 13 February 2019). Since this coastline is of a lower resolution (tens to hundreds of meters) and does not precisely match REMA in several areas, we masked as water all surfaces within 800 m of the coastline that were less than 2 m from the local mean sea level. Improving the delineation of the coastline is an objective of future versions of REMA.

The REMA mosaic includes a vertical error estimate, based on altimetry registration; coregistration for tiles aligned to neighbors, as described above; and grids of the day of data acquisition (Fig. 4). The 68th and 90th percentiles of errors are 0.63 and 1.00 m, respectively. Errors are highest in rougher terrain, such as the Antarctic Peninsula and Transantarctic Mountains. Higher errors also exist in zones of extensive crevassing along the coast and for tiles without control that are registered through alignment, and for which errors are thus propagated from the neighbors. The smallest errors are in the interior of the CryoSat-2 SARIn mask.

The mean date for REMA is 9 May 2015 with a standard deviation of 432 days. The mosaicking procedure resulted in no systematic distribution of date by acquisition time, but younger data tend to cover the higher latitudes, while older data cover areas of high science interest, as a result of long-term targeting. Our method of DEM registration to CryoSat-2 altimetry, described in Sect. 2.3, accounts for differences in time between the altimetry and DEM acquisitions, yielding temporal constraints on elevation for rapidly changing coasts and areas of fast flow. Even though many of the interior DEMs were registered to ICESat-1 data from late 2008, we retain the strip acquisition time in the date stamp as time-dependent changes in these regions are expected to be small relative to the data error (e.g., Helm et al., 2014). Areas of local change, such as over subglacial lakes, should be small enough so as not to substantially affect tile registration. Caution, however, should be used when assessing changes in tiles registered to ICESat-1. Tiles that are registered through neighbor alignment are given the weighted mean day of the data in the neighboring buffers.

We additionally filled gaps in the 100 m and coarser versions of the REMA mosaic using existing lower-resolution DEMs. For the Antarctic Peninsula area, we use the 100 m positing edited ASTER GDEM mosaic by Cook et al. (2012). For the rest of the continent, we use Helm et al. (2014). We filled the mosaic by reprojecting and linearly regridding the lower-resolution DEM to REMA, differencing the regridded fill DEM with REMA in areas of data overlap and then adjusting the fill DEM by the difference. The gaps are then filled with the adjusted fill DEM data. This process did cause artifacts in high-relief areas and where errors along the strip edges are propagated into the gaps by interpolation. These will be corrected or removed in future REMA versions as additional imagers are collected to fill gaps. Example hillshade representations of the REMA mosaic DEM that demonstrate the resolution are provided in Figs. S2, S3, S4, and S5.

Two common artifacts in the REMA DEM are noisy surfaces due to opaque shadows, typically on the south sides of mountains, and repeated horizontally offset surface features resulting from ice motion among stacked DEMs in the mosaic. Shadow artifacts can occur in both the strip and mosaic DEMs and appear as rough surfaces over the area of the shadow. Shadows reduce the confidence of the stereopair matching algorithm within the DEM-generation software, resulting in noisy surfaces. Examples of shadow artifacts are provided in Fig. S6. Ghosting only occurs in the mosaic, as it results from stacking multiple overlapping DEMs, and is most commonly found on fast-moving glaciers and ice shelves where crevasses and rifts advect with the ice. Ghosting artifacts appear as repeated offset features that may fade in and out due to the feathering applied in the mosaicking procedure. Examples of ghosting for ice shelf rifts are provided in Figs. S3 and S7.

Figure 5Validation results for DEM strips. Histograms of the 68th and 90th percentiles of the absolute value of vertical residuals, or linear error, among each of three Operation IceBridge lidar systems and REMA DEM strips after registration.

Download

We provide an independent validation of the REMA strip DEMs and mosaic through comparison to airborne lidar altimetry acquired by the U.S. National Aeronautics and Space Administration's Operation IceBridge (OIB) between 2009 and 2017. Data from three different lidar systems are available: the Airborne Topographic Mapper (ATM), the Land, Vegetation, and Ice Sensor (LVIS), and the ICECAP laser altimeter system. The ATM is a conically scanning, 531 nm, 5 kHz lidar with a nominal footprint size of 1 m and a single shot accuracy of ±0.1 m (Martin et al., 2012). The LVIS system is a high-altitude, 1064 nm, 500 Hz scanning lidar with a 20–25 m footprint and a vertical accuracy similar to that of ATM (Hofton et al., 2008). The ICECAP laser altimeter operates at 905 nm with a footprint resolution of 25 m along track by 1 m across track and an accuracy of 0.12 m (Young et al., 2015). All data were obtained from the National Snow and Ice Data Center (https://nsidc.org/, last access: 5 November 2018). For ATM we used the Level 1B elevation data product, while we used Level 2 geolocated elevation products for LVIS and ICECAP.

Figure 6Validation results for mosaic tiles. Histograms of the medians and linear errors at the 68th and 90th percentiles (LE68 and LE90) obtained from the differences between each REMA tile and the three NASA OIB airborne laser altimeters. The altimeter elevations are subtracted from the REMA elevations, so that a positive median of residuals (top plots) indicates the REMA surface is higher than the altimeter surface. Vertical red dashes are the median values.

Download

All lidar data were provided in geographic coordinates referenced to the WGS84 ellipsoid and were converted to the REMA polar stereographic projection. We selected lidar data collected within 18 months of the REMA strip acquisition date or mosaic date stamp. Strips were then three-dimensionally registered to each lidar point cloud, with the vertical residuals providing an estimate of precision. Histograms of the 68th and 90th percentiles of the absolute values of vertical residuals, or the linear error, LE, between each lidar system and the coregistered strips are shown in Fig. 5. The medians of the 68th percentiles are 0.44, 0.48, and 0.52 m for ATM, LVIS, and ICECAP airborne lidars, respectively, and 0.84, 0.98, and 1.01 for the 90th percentiles. These values are similar to those found by Shean et al. (2016) using DEMs created from the same imagery as REMA using ASP software, from Summit Camp, Greenland. They are, however, larger than the ∼0.3 m found in comparisons between 2 m data and ATM over coastal Greenland, which is likely due to a combination of resolution and the larger, less rigorously quality-controlled lidar and DEM datasets used here. Examination of outliers reveals that errors are often due to clouds and other errors in the various lidar datasets, as well as the DEMs. Thus, our data support the finding of Noh and Howat (2015) that the DEMs constructed from cloud-free imagery with adequate illumination and appropriate base–height ratio are of internal accuracy comparable (i.e., among locations for a single DEM) to that of available airborne lidar data.

Figure 7Median of elevation differences by mosaic tile between REMA and each of the three NASA Operation IceBridge lidar systems. Only measurements collected less than 1 year apart are used.

Download

For tiles, which are registered to satellite altimetry through the mosaicking process described above, we linearly interpolate the REMA grid to the center coordinates of each overlapping lidar data point collected within 1 year of the REMA data and differenced the interpolated REMA elevation from the lidar elevation. We then obtained the medians of the differences of all points within each tile, as well as the 68th and 90th percentiles of their absolute values (termed the linear error, or LE68 and LE90 for the respective percentiles). Histograms of these values are shown in Fig. 6, and the medians and root mean square of the residuals are mapped in Figs. 7 and 8. REMA elevations are, at the median, 0.06 and 0.47 m higher than ATM and LVIS, and 0.16 m lower than ICECAP elevations, while the LE68 values are 1.04, 1.19, and 0.77 m, and the LE90 values are 1.78, 1.74, and 1.25 m for ATM, LVIS, and ICECAP, respectively. The lower error values for the ICECAP data would be expected due to the typically lower sloped terrain of East Antarctica, where these data are collected. We find no significant relationship, however, between slope and error. The median difference and root-mean-square error values mapped in Figs. 7 and 8 are largely consistent with those given by the CryoSat-2 and ICESat-1 registration errors in Fig. 4a, with the largest errors found over areas of crevassing and rifting on the coasts, in the high mountains of the Antarctic Peninsula, and over tiles registered through alignment, such as around the pole hole. As with these strip comparisons, the comparison with tiles also reveals errors in the lidar datasets, likely caused by clouds and aircraft positioning errors.

Figure 8Root mean square of point elevation differences by mosaic tile between REMA and each of the three NASA Operation IceBridge lidar systems. Only measurements collected less than 1 year apart are used.

Download

The REMA datasets include all individual DEM strips (described in Sect. 2.3) and the mosaic in 100 km by 100 km tiles (Sect. 3), all as 32-bit floating point raster files in GeoTIFF format. The strip DEMs are either 2 or 8 m resolution, depending on the region (Fig. 9), and include a metadata text file giving the version, projection, and processing information. No ground control or altimetry registration is applied to the strip DEMs in the current (version 1) release. Version 1 includes 66 401, 2 m and 121 184, 8 m strip DEMs, totaling 45 TB uncompressed.

The mosaic is 8 m resolution everywhere and is registered to satellite altimetry data as described in Sect. 3. Each mosaic tile includes error estimate and date files, also in GeoTIFF format, as described above. The error file has 32-bit floating point precision, whereas the date file has 16-bit integer precision in units of days since 1 January 2000. No void filling is applied to the 8 m tiles. Version 1 includes 1524 mosaic tiles, totaling 1.0 TB uncompressed. In addition to the 8 m tiles, reduced-resolution resampled versions are provided at 100 m, 200 m, and 1 km resolutions. The reduced-resolution datasets have an alternate filled version. Previews of the non-annotated and annotated complete mosaic DEM hillshade representation maps provided by the PGC are shown in Figs. S8 and S9. Full-resolution versions of the maps are available at http://maps.apps.pgc.umn.edu/id/2364 and http://maps.apps.pgc.umn.edu/id/2365 (last access: 13 February 2019).

Figure 9Map showing coverage of 2 and 8 m resolution of DEM strips in the REMA version 1 release.

Download

All of the REMA products described above are openly available from the U.S. Polar Geospatial Center at https://www.pgc.umn.edu/data/rema/ (Howat et al., 2019). Imagery used to produce the REMA DEMs is available to US federally funded researchers through the Polar Geospatial Center by request.

The supplement related to this article is available online at: https://doi.org/10.5194/tc-13-665-2019-supplement.

IMH conceived this study, authored the filtering, registration and mosaicking software, and performed the validation analysis. CP managed the data production workflow. BES performed the CryoSat-2 registration. MJN authored the DEM extraction software and assisted with data production. PM assisted in data access and production. All authors assisted in paper preparation.

The authors declare that they have no conflict of interest.