The Reference Elevation Model of Antarctica (REMA) is the
first continental-scale digital elevation model (DEM) at a resolution of
less than 10 m. REMA is created from stereophotogrammetry with submeter
resolution optical, commercial satellite imagery. The higher spatial and
radiometric resolutions of this imagery enable high-quality surface
extraction over the low-contrast ice sheet surface. The DEMs are registered
to satellite radar and laser altimetry and are mosaicked to provide a
continuous surface covering nearly 95 % the entire continent. The mosaic
includes an error estimate and a time stamp, enabling change measurement.
Typical elevation errors are less than 1 m, as validated by the
comparison to airborne laser altimetry. REMA provides a powerful new
resource for Antarctic science and provides a proof of concept for
generating accurate high-resolution repeat topography at continental
scales.

Ice sheet surface elevation is among the most fundamental datasets in
glaciology. Often, investigations aimed at, for example, quantifying mass
balance or ice flow modeling are limited by the spatial and temporal
resolution and accuracy of surface elevation measurements. The polar regions
have particularly poor topographic data due to their remoteness and the
latitudinal limits of global datasets, such as the Shuttle Radar Topography
Mission (SRTM, limited to of 60∘ N and 56∘ S). For most
of Antarctica, continuous grids of surface elevation, generally termed
digital elevation models (DEMs), have been limited to spatial resolutions of
> 500 m and/or vertical errors reaching tens of meters or more
(e.g., DiMarzio et al., 2007; Griggs and Bamber, 2009; Cook et al., 2012;
Fretwell et al., 2013; Helm et al., 2014; Slater et al., 2018). This limits
their utility for geodetic applications, such as rectifying satellite
imagery. Openly available global DEMs, including the 30 m ASTER GDEM
(https://asterweb.jpl.nasa.gov/gdem.asp, last access: 13 February 2019) and
recently released 90 m TanDEM-X DEM
(https://geoservice.dlr.de/web/dataguide/tdm90, last access: 13 February 2019)
have large errors over ice sheet interiors and, in the case of the latter,
include a several-meter bias due to penetration of the X band into firn
(Wessel et al., 2018). Further, these DEMs do not have definitive time
stamping, limiting their use for elevation change measurements. Since
existing DEMs were mostly constructed from satellite ranging data (radar and
lidar) for which errors increase with surface slope, errors tend to be the
largest in areas of more complex terrain, such as the coasts, mountain
ranges, and outlet glacier interiors (Bamber and Gomez-Dans, 2005). These
areas, which include the Antarctic Peninsula and the Amundsen Sea outlet
glaciers, are also the areas where some of the largest changes in ice sheet
dynamics and mass balance are occurring. Measurements of ice shelf thinning
also require high precision because 1 m of thinning equates to only
∼0.1 m of surface elevation change. Complex patterns of
change, such as from ice shelf basal melting and subglacial hydrology, are
not often observable from current DEMs.

Figure 1Maps of coverage of individual digital elevation models (DEMs)
produced from stereoscopic submeter imagery for the REMA project, with color
indicating the number of repeats, for (a) all data, (b) DEMs that passed
visual quality inspection (note regional decrease in repeat coverage due to
change in procedure), and quality-controlled DEMs with registrations within
acceptable criteria from (c) CryoSat-2 and (d) the ICESat GLAS 2-D campaign.

High-precision elevations and elevation changes are obtainable from airborne
laser altimetry but are only available along narrow (hundreds of meters) swaths over
limited areas of the ice sheet, and at infrequent intervals in time. While
ICESat-2, launched on 15 September 2018, will greatly increase the density and
coverage of altimeter observations, it will not provide the continuous
surface required for modeling and data processing applications. Further, a
precise and time-stamped reference DEM will greatly benefit satellite
altimeter missions by providing a validation surface, a basis for data
filtering, and slope corrections for radar ranges and offset ground tracks.

The objective of the Reference Elevation Model of Antarctica (REMA) project
is to provide a continuous, time-stamped reference surface that has a resolution 1 to
2 orders of magnitude higher than currently available (Fig. 1)
and with absolute uncertainties of 1 m or less, depending on the
availability of ground control, and relative uncertainties (i.e., precision)
of decimeters. With REMA, therefore, any past or future point observation of
elevation provides a measurement of elevation change. Further, REMA may
provide corrections for a wide range of remote-sensing processing
activities, such as image orthorectification and interferometry, provide
constraints for geodynamic and ice flow modeling, mapping of grounding
lines, studies in surface processes, and field logistics planning.

REMA was constructed from stereophotogrammetric DEMs extracted from pairs of submeter-resolution commercial satellite
imagery and vertically registered to radar and laser altimetry data. Here we
describe the source datasets and the algorithms used to build REMA, as well
as a validation of the final product using airborne altimetry from multiple
sources.

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.1 DEM processing

DEMs were generated from scene pairs using the open-source and fully
automated SETSM software package (Noh and Howat, 2017) on the Blue Waters
supercomputer at the National Center for Supercomputing Applications (NCSA).
Prior to DEM processing, WorldView 1 and WorldView 2 images were destriped using the
wv_correct function within the Ames Stereo Pipeline software
package (Shean et al., 2016). SETSM DEMs commonly have artifacts at scene
edges due to lack of constraint on triangulated irregular network (TIN)
generation and neighbor-based filtering at the scene boundaries. These
artifacts appear as unrealistically high relief extending tens to hundreds of
pixels from the scene edge. To detect and remove boundary artifacts, we
computed the average slopes over square 21-pixel kernels. All pixels within
a square 13-pixel radius of those with a mean slope greater than 1 were then
removed. Enclosed gaps were then filled, so that only gaps touching the
scene edge remain. Isolated clusters of less than 1000 pixels were then
removed. A convex hull algorithm that includes concavity across data gaps
was then applied to the remaining data to define the cropped scene boundary.

We additionally filtered each scene DEM for erroneous surfaces resulting
from clouds or opaque shadows using the density of successful matches in the
DEM extraction processes as given in the match tag file. We derive a match
point density field by calculating the fraction of successful matches within
square 21-pixel kernels. Pixels are then filtered if the match point density
is below 0.9.

The filtered scene DEMs were then merged with adjoining scenes to form DEM
strips comprising the overlapping area of the original stereopair image
strips, performing three-dimensional coregistration using the iterative
least-squares method of Nuth and Kääb (2011) and Levinson et al. (2013) and
with distance-weighted averaging over the overlapping areas. Extensive
erroneous surfaces due to clouds or water bodies, for example, will cause errors in
coregistration. Therefore, the scene was not merged if the root mean square
of the residual differences in elevation between the overlapping area of the
coregistered scenes was greater than 1 m. In this case, the strip was broken
into separate segments and was treated as separate DEMs during the global
mosaicking step described in Sect. 3. We note that the coregistration
procedure may not provide correct horizontal offsets in extremely flat, or
uniformly sloping, terrain because the small range in slopes and aspects may
not yield a confident regression. We could not identify such cases, however,
suggesting that there is enough surface variation at these high resolutions
(2–8 m) for the method to be successful.

From the archive of all imagery collected over the Antarctic continent as of
July 2017, and with a cloud cover classification of 20 % or less, we
processed 79 362 individual strip pairs to create 187 585 DEM strip
segments, with 66 401 of these covering West Antarctica and mountainous
areas of East Antarctica processed to a resolution of 2 m and the remaining
areas to
a resolution of 8 m (Fig. 1a). The lower resolution over regions of,
generally, flat ice sheet surface was chosen to save computational costs.
This equates to 122 567 288 km2 of total coverage, including repeat
coverage, and coverage of 13 987 485 km2 (or 98 %) of the continent,
including islands located above 60∘ S. The largest gap
occurs over the “pole hole”, south of approximately 88∘ S,
with smaller gaps, mostly on occluded sides of mountains and in areas of
persistent clouds such as the Antarctic Peninsula. These gaps will receive
priority tasking in the future.

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 km2 of quality-controlled DEMs cover an area of 13 567 969 km2,
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).

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 km2 is covered by CryoSat-2 registered DEMs, or
29 901 958 km2 including repeat coverage (Fig. 1c), with registered
ICESat DEMs covering 4 897 600 km2, including 8 739 128 km2 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 -0.39±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.

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.

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.

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.

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.

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.

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.

Stereophotogrammetry from high-resolution commercial satellite imagery has
enabled the first elevation mapping of nearly an entire continent at a
horizontal resolution of less than 10 m and with a vertical error of less than
1 m. The construction of REMA demonstrates the highly complementary
characteristics of satellite altimetry, either from laser or radar, and
stereo DEMs; altimetry provides highly accurate, but relatively sparse,
control points to which the stereo DEMs provide a continuous surface of
similar precision but lower accuracy. The combination of the two provides an
effective method for maximizing resolution, coverage, and accuracy.

Its geographic location, the flatness of the ice sheet and lack of
vegetation all make Antarctica the easiest case for application of these
methods. Polar orbiting satellites, with little competing demand for imagery,
provide the most abundant data at the poles. The flatness and lack of
vegetation simplify registration to satellite altimetry and ambiguities in
the canopy versus ground height. These complications will need to be
considered when expanding these methods to lower latitudes.

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.

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.

This work was supported by US National Science Foundation Office of Polar
Programs grants 1543501 to Ian M. Howat and 1559691 to Paul Morin. High-performance computing
resources were provided by the National Center for Supercomputer
Applications at the University of Illinois.

The Reference Elevation Model of Antarctica (REMA) is the first continental-scale terrain map at less than 10 m resolution, and the first with a time stamp, enabling measurements of elevation change. REMA is constructed from over 300 000 individual stereoscopic elevation models (DEMs) extracted from submeter-resolution satellite imagery. REMA is vertically registered to satellite altimetry, resulting in errors of less than 1 m over most of its area and relative uncertainties of decimeters.

The Reference Elevation Model of Antarctica (REMA) is the first continental-scale terrain map at...