The Greenland Ice Mapping Project ( GIMP ) land classification and surface elevation datasets

Introduction Conclusions References


Introduction
The objective of the Greenland Ice sheet Mapping Project (GIMP) is to establish benchmark data sets for observing ice sheet change.Such data sets include ice-sheet-wide ice velocity (Joughin et al., 2010;Moon et al., 2012) and surface elevation maps, as well as time series of ice velocity and elevation for selected areas of rapid, ongoing change.Production of these data requires spatial classification of ice-covered and ice-free surfaces for horizontal and vertical co-registration of data through subtraction of offsets over ice-free (i.e., stable) terrain.Processing of these data also requires a digital elevation model (DEM) at a resolution commensurate with the resolution of the imagery being processed.The resolution of synthetic aperture radar interferometer (InSAR) ice flow measurements, for example, are from 10 to 100 m.This resolution is 1 to 2 orders of magnitude finer than available DEMs for the ice sheet Additionally, a comprehensive mapping of the ice edge from data collected over a narrow time window provides a benchmark measurement for change detection.For this reason, GIMP produced a classification mask at the highest possible spatial resolution that could be achieved with widely available imagery.
Here we present, in the following order, the construction of three geospatial products for the GIMP project: (1) a complete, 15 m resolution image mosaic, (2) ice-covered and icefree terrain classification masks, also posted to 15 m resolution, and (3) a complete, altimeter-registered, digital elevation model posted at 30 m.

Projection and grid
We use a polar stereographic projection centered on Greenland, with an origin at 90 • N, 45 • W, a standard parallel of 70 • N and a reference to the WGS84 ellipsoid.This projection has the European Petroleum Survey Group (EPSG) code 3413 and is the National Snow and Ice Data Center's (NSIDC) standard planar projection for the Arctic.Software for converting between this projection, UTM and geographic coordinates are widely available.Most of the remote sensing data used here were obtained in UTM projection and were converted to our polar stereographic projection using the Geospatial Data Abstraction Library (GDAL) open source package (http://gdal.org).All mosaics are re-gridded, through bilinear interpolation, to a master grid with an origin at 59.1996 • N, 55.7983 • W (Fig. 1).The grid is divided into 36 tiles of 6 rows by 6 columns, each with dimensions of 249.3 by 450 km.These dimensions were selected because, first, they are divisible by 15 m, which is the resolution of Landsat-7 Enhanced Thematic Mapper Plus (ETM+) band-8 (panchromatic) imagery and, thus, is the base-level resolution we adopt and, second, they yield a 15 m resolution, 8-bit (integer) image of slightly less than 500 MB uncompressed.

Landsat 7 ETM+/RADARSAT-1 image mosaic
Our first objective is to assemble an ice-sheet-wide imagery mosaic to be used for mapping and land surface classification at the highest possible spatial resolution and within as narrow a time window as possible to enable change detection.South of ∼ 81.2 • N, we use Landsat 7 ETM+ imagery orthorectified and distributed by the U.S. Geological Survey (http://glovis.usgs.gov/).Using 1 August 2000 as a target date, we selected imagery from July and August, as close in time as possible to 1 August for the years, in preferential order, of 2000, 1999, 2001 and 2002.All imagery was automatically filtered for clouds using the algorithm presented in Luo et al. (2008), adapted to Landsat by Oreopoulos et al. (2011) and visually inspected for quality.In some cases additional manual cloud masking was required.In order to increase the consistency of the grayscale between images, each digital number image was converted to reflectance, including corrections for sun angle and distance using the parameters provided in the metadata.Multispectral bands 1 through 4 were pan-sharpened to 15 m posting using band-8 and a simple and fast additive method in which the band-8 image was down-sampled to 30 m and differenced from each multispectral band.The difference image was then up-sampled to 15 m using bilinear interpolation and added to the band-8 image.An example of pan-sharpening is given in Fig. 2.
The pan-sharpened reflectance images were then regridded via cubic convolution and mosaicked to the reference grid.Where images overlapped, the pixel that was closest in time to the target date of 1 August 2000 was selected.No edge feathering was applied.The mosaicked images were then converted back to a byte precision digital number by linearly scaling the reflectance values to the global minimum and maximum for each band (Fig 1 .).
The USGS employs two levels of geo-registration processing for their imagery (see http://landsat.usgs.gov/Landsat_Processing_Details.php).First, Standard Terrain Correction (Level 1T) incorporates both ground control points and a DEM for terrain corrections.Geodetic accuracy depends on the accuracy of the ground control and the quality of the DEM and is better than 90 m.Imagery covering the periphery and margin of the ice sheet, where features are visible on the surface, are processed to L1T.For L1T imagery, the root mean square of the residual between the geo-location model and the ground control are provided in the imagery metadata and are typically on the order of several meters.Second, Systematic Correction (Level 1G) uses only the satellite ephemeris for geo-location, providing a 1σ geometric accuracy within 250 m.Scenes over the featureless interior of the ice sheet are typically processed to L1G.
North of the maximum extent of Landsat we include synthetic aperture radar amplitude imagery mosaics acquired between October and December 2000 by the RADARSAT-1 satellite.These data were produced by the Applied Physics Laboratory at the University of Washington as part of GIMP (Joughin et al., 2010).The data are distributed at 20 m resolution and were up-sampled through bilinear interpolation to 15 m to match the resolution of Landsat band-8.We merged the RADARSAT and Landsat band-8 imagery by applying a stretch to the RADARSAT image so that the histograms of both data sets match where they overlap.As with Landsat, the primary source of geolocation error in the RADARSAT imagery is error in the DEM used for terrain correction and these errors are similar in magnitude to the Landsat mosaic (Moon and Joughin , 2008).
The final image mosaic (Fig. 1) is distributed in tiles, with one image for each band, plus an index image in which each pixel gives the index number of its corresponding source image in an accompanying metadata file.The metadata file lists each Landsat scene identification number (scene ID) used in the mosaic for that tile, the acquisition time and the rootmean-square control point registration error where available.The original scene ID, acquisition date and geo-location error for any pixel in an image can thus be obtained using the index image.

Land classification masks
Land classification masks are needed for co-registration of repeat imagery and elevation data, as ice surfaces can change with time while areas of exposed bedrock provide control.Further, the accurate delineation of ice boundaries provides a benchmark for measuring future ice margin changes.Landsat-7 ETM+ data are commonly used for mapping snow and ice, either manually, by tracing the margin with a computer mouse directly on the imagery, or automatically, from multispectral classification techniques (e.g., Rastner et al., 2012).Automatic methods are far more efficient and are effective for ice and snow that is free of surface debris.However, the drawbacks of automated, multispectral classification methods are that (1) they cannot differentiate between seasonal/ephemeral snow cover and glacial ice, (2) they fail at marine margins when dense packs of icebergs and sea ice are present, (3) much of the marginal ice of the Greenland Ice Sheet and surrounding glaciers is debris-covered and (4) Landsat does not cover the most northern regions of the ice sheet.For these reasons, we abandoned multispectral mapping methods in favor of manual digitization of the panchromatic and pan-sharpened multispectral image mo-saic presented in Sect.3.Even with manual methods, the ice margin can be difficult to locate visually in areas of abundant debris and snow cover.Margins of debris-covered ice were identified by breaks in surface slope, emerging melt water streams, color differences and the presence of small melt water ponds typical of debris-covered glaciers.Similarly, glaciers were differentiated from perennial snowfields by visible crevassing, surface moraines and the existence of a visible toe.Snowfields without these features were not classified as glaciers.Using the same method, we also digitized the coastline to produce an ocean mask, with the null of the ice and ocean masks being ice-free terrain (including freshwater lakes).
Uncertainty in these classification masks arise from three sources of error: (1) image pixel resolution, (2) image georegistration and (3) erroneous selection or non-selection of pixels (i.e., mapping error).All error sources are expected to vary randomly in space, although there is likely a systematic component of error source (2) over distances equivalent to the size of a single image (e.g., 185 km for Landsat 7) due to errors in the registration model used to orthorectify the image, which typically is on the order of ±5 m, or one third of a pixel for L1T-processed imagery.
Error source (1) contributes a random error of one pixel for each ice boundary pixel.The position of any point of the ice margin has an uncertainty of 21 m while the total error for a given area of ice is then (8N ) 1/2 x 2 , where N is the number of boundary pixels and x is the pixel posting in meters.
Erroneous selection or non-selection of pixels can be due to debris cover, shadows and misidentification by the operator, as well as the ambiguity of delineating an ice boundary at glacier fronts ending in packs of icebergs.Without ground control, delineation of the ice edge in areas of debris cover, terminal moraines and persistent snow cover is subjective.These errors are difficult to quantify.We estimated uncertainties due to ambiguity in the ice edge and operator error by comparing mappings done by three different operators over the same area.On average, each operator identified 24.21 km (1614 pixels) of ice margin over the common area, with a 660 m (44 pixels) difference between the maximum and minimum mappings, giving an estimated error of ±3 %, which is similar to other comparisons (Paul et al., 2013).This error, however, is expected to vary widely by particular location and size of area considered.
Initial versions of the GIMP classification mask have been used and analyzed in two studies.Rastner et al. (2012) compared the version 1.1 GIMP classification to their own, semiautomated delineation of peripheral glaciers and ice caps, which also utilized Landsat 7 data.They found an overall difference in classified area of 6 %.This difference was mostly due to misclassification of debris-covered margin in GIMP.That study incorporated the GIMP classification into their data set for far northern regions, and their combined map has been included in the global Randolph Glacier Inventory (Pfeffer et al., 2014).Citterio and Ahlstrøm (2013) compared the version 1.2 GIMP classification to glacier outlines mapped from aerial photography in the 1980s and were able to measure local changes in margin positions between the data sets.They also detected some classification errors.Errors detected in both of these studies have been corrected in the current version 2.0 of the mask, along with additional quality control by our team.Both the ice and ocean classification masks were used in the production of the digital elevation model, described next.

Digital elevation model
The quality of data over most of the Greenland Ice Sheet in global elevation data sets, is too poor to be of use for glaciological applications.The standard DEM used in glaciological studies was created from a combination of satellite radar altimeter and aerial photogrammetry (Bamber et al., 2001) with a posting of up to 1 km.This DEM was enhanced to 625 m posting through photoclinometry by Scambos and Haran (2002).While these DEMs are accurate to a few meters over the relatively flat interior of the ice sheet, they have poor resolution over the steeper margins and higher-relief periphery.
Our objective is to enhance DEM resolution and accuracy, particularly over the ice sheet margin and periphery, by integrating high-quality photogrammetric topography data into the existing low-resolution DEM and registering the DEM to elevations acquired by the Geoscience Laser Altimeter System (GLAS) aboard the Ice, Cloud, and land Elevation Satellite (ICESat).Our approach follows the schematic shown in Fig. 3.We focus on generating a continuous surface and we ignore temporal changes in ice elevation, which are over 100 m near the fronts of some rapidly retreating glaciers, and produce a DEM that approximates the mean elevation over the ICESat era (2003)(2004)(2005)(2006)(2007)(2008)(2009).We first present each input data set and then describe the procedure for merging them, followed by a description of errors and artifacts in the resulting DEM.

ICESat GLAS
All data are referenced to elevations obtained from ICESat GLAS between 2003 and 2009.We use the 633 products of the GLA12 release corrected for time-varying elevation biases, as estimated based on apparent variation of the mean sea surface height (Shepherd, 2012; see their supporting online material).Poor-quality returns were removed using techniques developed for elevation-change estimation that identify the best-quality returns based on parameters that describe the shape and amplitude of the returned laser pulse (Shepherd et al., 2012).Elevations were corrected for detector saturation, and the time-varying bias correction should remove offsets associated with campaign-to-campaign variations in the shape of the transmitted pulse (Borsa et al., 2014).Elevations calculated in this way should be accurate to better than 0.1 m, or two orders of magnitude smaller than the expected DEM uncertainty.

Photo-enhanced Bamber (PEB) DEM
The most widely used DEM for the entire ice sheet is that presented in Bamber et al. (2001), created from a combination of radar altimeter and stereo-photogrammetric data from the mid 1990s.These data were validated against airborne altimeter data, also from the mid-1990s, with a reported, icesheet-wide 1σ error of ±7 m and errors of several hundred meters at the coasts.This DEM was subsequently enhanced through photoclinometry with Advanced Very High Resolution Radiometer (AVHRR) imagery (Scambos and Haran, The PEB DEM was provided by the NSIDC in a spherical Lambert azimuthal projection at a posting of approximately 627 m.These data were re-gridded to EPSG 3413 and upsampled to 30 m posting using bilinear interpolation.The regridded data were then co-registered to the ICESat GLAS point cloud using an iterative, 3-D conformal transformation (Noh and Howat, 2014).This procedure results in residuals between the DEM surface and ICESat point cloud with a normal distribution and a mean of 0. Co-registration was performed on 25 by 25 km tiles with 5 km of overlap.The co-registered tiles were then mosaicked with linear distanceweighted edge feathering.The root mean square (RMS) of the residuals between the PEB DEM and the ICESat point cloud following co-registration are given in Table 1.The total RMS error of ±21.8 m is nearly three times higher than reported by Bamber et al. (2001) and Scambos and Haran (2002), likely due to the more extensive sampling by ICESat relative to the airborne altimetry used in those studies, especially over ice-free terrain, where errors are much higher.The RMS errors over the interior ice sheet are more consistent with reported errors.

GDEM V2
The Global Digital Elevation Model (GDEM) is a global, 30 m posted DEM produced by the Ministry of Economy, Trade and Industry (METI) of Japan and the United States National Aeronautics and Space Administration (NASA) (Slater et al., 2011).The GDEM is created by average-stacking individual stereo-photogrammetric DEMs acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) between 2000 and 2010.Following an initial release in 2009, Version 2 was released in October 2011.The GDEM is distributed in 1 • × 1 • tiles in geographic projection.The distribution includes metadata giving the number of individual AST14DEM granules that were stacked to obtain each posted elevation.No information, however, is given regarding which scenes were used, so the time period of elevation measurements cannot be determined directly.
GDEM data quality is poor over much of the ice sheet owing to low-contrast surfaces on snow and ice (Fig. 4b).Additionally, artifacts due to shadows, clouds and blunders in the automated matching algorithm are abundant over all terrains.Following reprojection and gridding of the GDEM Version 2 to the GIMP grid, we applied a pyramiding standard deviation filter in which the DEM is smoothed to progressively finer resolutions and differenced from the nativeresolution DEM.Pixels with differences exceeding 2.5σ of the mean are discarded.Since ice-covered terrain is substantially smoother than ice-free terrain, we apply this filter separately to the two land classifications, using the masks derived in Sect. 4. Following automated filtering, we manually masked blunders visible on a hillshade image of the DEM.These procedures removed nearly all data from above 1600 m elevation, which is approximately the average mass balance equilibrium line altitude (Fig. 4c).Following filtering and masking, GDEM covers 30 % of Greenland's total area and respectively 92 and 19 % of its total ice-free and ice-covered terrain.For this project, we obtained all available SPIRT DEM products over Greenland.Each DEM was reprojected to EPSG 3413 and the WGS-84 ellipsoid and up-sampled to 30 m.As advised in Korona et al. (2009), we use version 2 of each DEM and mask out all interpolated pixels.We then applied the same filtering and masking procedure as used for the GDEM.
Each individual SPIRIT DEM was then co-registered to overlapping regions of the filtered GDEM using the 3-D conformal transformation (Noh and Howat, 2014).This provided a consistent registration between the SPIRIT and GDEM data sets to facilitate merging.Each individual SPIRIT DEM was then stacked into a single mosaic by taking the median elevation at each pixel, keeping track of the number of individual measurements.The resulting filtered SPIRIT mosaic (Fig. 4d) covers 10 % of Greenland's total area and respectively 24 and 8 % of its total ice-free and ice-covered terrain.The most continuous coverage is along the southwestern and southern coasts, with approximately 50 % of the land and ice area covered in each tile (Table 1), or most of the land and ice area below 1500 m elevation.

CNES mean sea surface height
Stereo-photogrammetric methods typically cannot resolve open-water surfaces due to the lack of features, so that these surfaces are usually interpolated from the shoreline.This and the presence of icebergs result in spurious sea surface heights in stereo-photogrammetric DEMs.To ensure correct sea surface heights, we apply the ocean mask derived in Sect. 4 to the final DEM and replace those ocean surfaces with the CLS11 mean sea surface height product from the Centre National d'Etudes Spatiales (CNES).The CNES CLS11 is the 16-year mean of surface height measurements made by several satellite radar altimeters and gridded to 1/3 of a degree (Schaeffer et al., 2012).We reproject these data to EPSG 3413 and up-sample them to the 30 m GIMP grid using bilinear interpolation.

Data merging
Following co-registration and stacking, the SPIRIT DEM mosaic was differenced from the GDEM and the differences were extrapolated across the grid using an inversedistance interpolation.The extrapolated difference map was then added to the SPIRIT stack.The GDEM and SPIRIT DEMs were then merged under the following conditions at each pixel: 1.If there was a GDEM value but no SPIRIT value, the pixel is assigned the GDEM value.
2. If there was a SPIRIT value but no GDEM value, the pixel is assigned the corrected SPIRIT value.
3. If there were both GDEM and SPIRIT values and the pixel is over ice-free terrain, the pixel is assigned GDEM value.This is due to the GDEM's higher spatial resolution.
4. If there were both GDEM and SPIRIT values and the pixel is over ice, the pixel is assigned the average of the GDEM and SPOT values, weighted by the Nnumber of observations, where Nequals 1 for GDEM plus the number of individual SPIRIT DEMs used in the stack described in Sect.5.4.
The merged GDEM and SPIRIT DEM (merged G&S) was then co-registered to the ICESat GLAS point cloud using the 3-D conformal transformation (Noh and Howat, 2014).The RMS validation errors of the merged G&S DEM are given in Table 1.To assess the improvement in validation score provided by the higher-resolution data, Table 1 also gives the RMS errors for the PEB DEM exclusive to areas of overlap with the merged G&S DEM.On average, the merged G&S DEM improves the validation score by a factor of 8 over the ICESat-registered PEB DEM.
To combine the merged G&S DEM and PEB DEM, the PEB DEM was first adjusted by differencing it from the merged G&S DEM and interpolating the differences across areas of no data in the merged G&S DEM.The difference was then added to the PEB DEM and the two DEMs were combined using the following rules at each pixel: 1.If there was a merged G&S DEM value, the pixel is assigned the merged G&S value.
2. If there was no merged G&S DEM value, the pixel is assigned the adjusted PEB DEM value.
An ocean mask (see Sect. 4.) is then applied and those pixels are replaced with the CLS11 sea surface heights, as described in Sect.5.5.The final GIMP DEM thus provides an altimeter-registered, relief-enhanced version of the PEB DEM.The enhancement is most pronounced over regions of high relief on the margin and periphery of the ice sheet (Fig. 5).Notably, whereas outlet glaciers are not clearly defined in the PEB DEM, the GIMP DEM resolves outlet glacier termini and fjord walls in detail (see example in Fig. 5).

Errors and artifacts
The ICESat validation errors of the completed GIMP DEM are mapped in Fig. 6 and listed by tile in Table 2.The overall RMS of the differences between the GIMP DEM and ICESat elevation is ±9.1 m, which is less than half that of the ICESat-registered PEB.The error on ice-free terrain (±18.3 m) is over twice that of ice-covered terrain (±8.5 m), which is to be expected considering the higher relief at the ice-free margin.We note that an unknown amount of this error can be attributed to differences in the geometries of the ICESat footprint, which has a typical diameter of 70 m, and  the DEM pixels.The effect of this difference will increase with slope.Additionally, over ice, some amount of the validation error can be attributed to temporal variations in surface elevation, ranging from decimeters over the interior to tens of meters over rapidly thinning outlet glaciers.Besides ice thinning, the advection of crevasses and other surficial expressions with ice flow contributes an unknown error.These validation errors should, therefore, be viewed as an upper bound for the true standard data error.The largest validation errors exist for the most northern regions, for which few high-resolution data exist and coverage is mostly from the PEB DEM.Higher errors, exceeding ±20 m, are also found in areas of extreme relief, such as the Geikie Peninsula (tile 4,2 and 5,2), where gaps in highresolution data coverage exist over steep mountain glaciers and ice caps.
To examine the influence of slope on validation error, we plot RMS and mean and median GIMP DEM and ICESat differences by slope in Fig. 7.As expected, RMS errors are the smallest for the flattest surfaces (e.g., the interior ice sheet), increasing with slope to a peak of ±24 m at 2 • .RMS error then decreases to ±13 m for 5 • slopes before increasing again.The peak in RMS error at 2 • slopes corresponds roughly with the equilibrium line of the ice sheet and, thus, the boundary between the merged G&S and PEB DEM.Errors in both the PEB and merged G&S DEMs result in spurious, step-like transitions between the two (Fig. 8a).This effect results in the continuous zone of large errors running along the southeast margin, which is especially steep (Fig. 6).A positive peak in mean and median errors of 2.1 and 0.9 m, respectively, at 1.1 • shows a positive bias in the GIMP DEM relative to ICESat over the area just inland of the snow line.Over steeper terrain, this bias becomes increasingly negative (i.e., the GIMP DEM reports increasingly lower elevations than ICESat), from −0.5 m at 5 • to −1.5 at 25 • .Since most of the coverage of surfaces with 1 • slopes are from the PEB DEM, the positive bias could be explained by either slopedependent errors in the PEB or thinning between the PEB and ICESat epochs, with neither effect completely compensated  during co-registration.The cause of the negative bias over steeper terrains is unknown.Since these biases are spatially variable and are small (< 10 %) relative to the random error, we do not correct for them.
Where merged G&S coverage exists above the snow line, the apparent surface is much rougher, with pitting resulting from blunders in the surface matching procedure used to gen- ∼ 10 m tall rises in the surface visible here as large bumps in the center of the glacier and clustered around its margins.Also visible is a region of dense icebergs, also known as mélange, at the glacier front in the bottom right corner.This mélange region was not masked due to retreat of the glacier between the time the ice front was mapped and the elevation data were acquired.erate the DEMs (Fig. 8b).These roughness features typically have amplitudes of several meters.
Rapid ice thinning and front retreat also cause DEM artifacts.Many fast-moving outlet glaciers thinned by tens of meters, reaching over 100 m in some cases, during the data collection period.This thinning causes offsets between DEM surfaces acquired at different times and, when stacked, can result in spurious offsets and discontinuities in the surface.An example of this effect is shown in Fig. 9. Additionally, ice-front retreat between date of the imagery used in construction of the ice cover mask and DEM data acquisition causes incomplete masking of the ocean boundary.For outlet glaciers, this often means that areas of dense icebergs remain in the DEM.This is also shown in the example in Fig. 9.

Conclusions
We have presented three novel and complete high-resolution geographic data sets for the Greenland Ice Sheet designed to serve as benchmarks for change detection, control for image orthorectification and other topography or surface classification-dependent processing and constraints for ice sheet models.A second phase of the GIMP project began in 2013.This second phase will include improvements to the GIMP DEM through the inclusion of sub-meter resolution www.the-cryosphere.net/8/1509/2014/The Cryosphere, 8, 1509-1518, 2014 DEMs acquired by the Worldview series of satellites and from airborne laser altimetry collected through NASA's Operation IceBridge.A major goal of the next phase of DEM improvements will be to provide date stamping for each pixel in the DEM so that it may be used for change detection.Additionally, complete, updated image mosaics and masks will be constructed from imagery acquired by the newly launched Landsat 8 satellite, providing an opportunity for ice-sheetwide change detection.
As with all data sets produced as part of the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) program, all GIMP data sets will be available online and at no cost through the NASA Distributed Active Archive Center at the NSIDC (www.nsidc.org/data/measures).Prior to distribution at the NSIDC, betaversion data sets are available from the Glacier Dynamics Research Group at the Byrd Polar Research Center (bprc.osu.edu/GDG/data.php).In both cases, data access requires registration and acceptance of a data use agreement.Announcements of updates will be sent to all registered users.

Figure 2 .
Figure 2. Example subset of a natural color (bands 1, 2 and 3) Landsat-7 ETM+ image at (a) original 30 m resolution and (b) after pansharpening to 15 m using band-8.This example highlights the improvement in resolution of common ice sheet features including a glacier calving front, medial moraine and supraglacial melt water stream channels.

Figure 3 .
Figure 3. Schematic of the approach used to the produce the GIMP DEM from the three source data sets.

Figure 4 .
Figure 4. Grayscale representations of the input digital elevation models used to create the GIMP DEM, including the (a) photo-enhanced Bamber (PEB) DEM, (b) GDEM2, (c) filtered and masked GDEM2 and (d) SPIRT mosaic.

5. 4
Photogrammetrically derived DEMs over Greenland were produced from images acquired in 2007 and 2008 as part of the Satellite Pour l'Observation de la Terre (SPOT-5) stereoscopic survey of Polar Ice: Reference Images and Topographies (SPIRIT) program.A description of data set production and validation is given inKorona et al. (2009).The SPIRIT DEM is distributed in UTM projection and referenced to the EGM96 Geoid and posted at 40 m.Two versions of each DEM, processed with different correlation parameters, are provided, along with data quality and interpolation masks.Korona et al. (2009) reports a slightly better precision and accuracy of SPIRIT DEM (< ±5 m) over ASTER DEMs based on validation experiments with ICESat.For this project, we obtained all available SPIRT DEM products over Greenland.Each DEM was reprojected to EPSG 3413 and the WGS-84 ellipsoid and up-sampled to 30 m.As advised inKorona et al. (2009), we use version 2

Figure 5 .
Figure 5. Hillshade representations of tile 1,2 (see Fig. 1) for the (a) PEB and (b) GIMP DEM showing the improvement in resolution of the margin and ice-free periphery.Enlargements over the front and fjord of Jakobshavn Isbrae for (c) PEB and (d) the GIMP DEM highlight the improvement in resolving outlet glacier termini and fjords.

Figure 6 .
Figure 6.Map of the root mean square (RMS) of differences between the GIMP DEM and ICESat elevations.Note the logarithmic color scale.Individual point differences along ICESat tracks are linearly interpolated to a 1 km by 1 km map grid in polar stereographic projection.

Figure 7 .
Figure 7. Plots of surface slope, in degrees from horizontal, versus the (a) GIMP DEM/ ICESat RMS validation error and (b) the (blue curve) mean and (red curve) median difference between the GIMP DEM and ICESat elevations.

Figure 8 .
Figure 8. Hillshade representations of examples of common artifacts in the GIMP DEM.Boxes denote the location of enlargements.(a) Example of an offset in the boundary between the merged GDEM/SPIRIT mosaic and the PEB DEM in a particularly steep section of the ice sheet margin from tile 0,2.The apparent "cliff" resulting from this offset is up to 5 m high.(b) Example of rough ice sheet surfaces resulting from blunders in stereo-photogrammetric DEM extraction over relatively featureless terrain.The smoother areas are PEB DEM data coverage.The roughness amplitude is less than 3 m.

Figure 9 .
Figure 9. Hillshade representation of the GIMP DEM for the terminus of Kangerdlugssuaq glacier, East Greenland (tile 4,2) showing artifacts created by rapid elevation change.In this case, ice thinning of ∼ 100 m over the data collection period results in spurious,

Table 1 .
Coverage and ICESat validation statistics for the data sets used in the GIMP DEM.Tile boundaries are delineated in Fig. 1.
* Statistics only for areas of overlap with the merged GDEM and SPIRIT DEM.

Table 2 .
Land classification and ICESat validation statistics for each GIMP DEM tile.Tile boundaries are delineated in Fig. 1.