Digital elevation model and orthophotographs of Greenland based on aerial photographs from 1978–1987

Abstract

Digital Elevation Models (DEMs) play a prominent role in glaciological studies for the mass balance of glaciers and ice sheets. By providing a time snapshot of glacier geometry, DEMs are crucial for most glacier evolution modelling studies, but are also important for cryospheric modelling in general. We present a historical medium-resolution DEM and orthophotographs that consistently cover the entire surroundings and margins of the Greenland Ice Sheet 1978–1987. About 3,500 aerial photographs of Greenland are combined with field surveyed geodetic ground control to produce a 25 m gridded DEM and a 2 m black-and-white digital orthophotograph. Supporting data consist of a reliability mask and a photo footprint coverage with recording dates. Through one internal and two external validation tests, this DEM shows an accuracy better than 10 m horizontally and 6 m vertically while the precision is better than 4 m. This dataset proved successful for topographical mapping and geodetic mass balance. Other uses include control and calibration of remotely sensed data such as imagery or InSAR velocity maps.

Design Type(s)	observation design
Measurement Type(s)	elevation
Technology Type(s)	objective lens
Factor Type(s)	Â
Sample Characteristic(s)	Greenland • glacier

Machine-accessible metadata file describing the reported data (ISA-Tab format)

Background & Summary

DEMs produced from aerial and satellite imagery are widely used in studies of geodetic glacier mass balance and landform mapping^1–10. In Greenland, photogrammetry has been used on terrestrial and aerial photographs for glaciological research and topographical mapping since the early 1930s and is still in use on aerial photographs today^10–14. Since the 1960s satellites have become an increasingly important platform for acquisition of stereoscopic imagery¹⁵. Gridded DEMs derived from optical imagery recorded by e.g., SPOT-5, Terra (ASTER) and WorldView satellites have been used to answer a range of outstanding research questions about the Greenlandic cryosphere^3–5. Nevertheless, high-latitude stereo-photogrammetric DEMs, are limited by the low-visible contrast of snow and ice surfaces that reduce the ability to resolve heights from an image stereopair^6,8. In Greenland, this limits the data acquisition to the coastal regions including the margin and outlets of the Greenland Ice Sheet (GrIS). Thus, to overcome this limitation, elevation data from other sources, such as radar altimetry, Synthetic Aperture Radar (SAR) and photoclinometry, are often combined with stereo-photogrammetry to produce a complete, but sometimes temporally inconsistent, elevation coverage of Greenland including the GrIS^7–9.

Here, we present a 25×25 m gridded DEM and 2 m resolution orthophotograph derived from ~3,500 vertical photographs (scale 1:150,000) acquired from 1978 to 1987 that cover all ice-free areas of Greenland and a significant part of the GrIS margin. In Fig. 1 the resulting footprints of the triangulated photos are shown by year of recording. The aerial campaigns were carried out by the Agency for Data Supply and Efficiency, SDFE (previously Geodetic Institute, then National Cadastre and Survey of Denmark, and Danish Geodata Agency), and ground control was surveyed and used for aero-triangulation of the photographs by the SDFE and the National Space Institute, Technical University of Denmark (DTU Space)^{14,16–18}. A challenge for stereo-photogrammetry in remote, Arctic regions is the scarcity and quality of the ground control⁶. The photographs used for our data set have previously been used for DEMs and mapping, but geodetic control has been of inconsistent quality^7,8,12,14. Here, we use the latest 2006 aero-triangulation with ground control provided by a subset of the GPS-based REFGR Greenlandic reference network supplemented with Doppler stations. Our technical validation comprises an assessment of the aero-triangulation, complemented by two validations using external data. We co-register DEM tiles and 50×50 km DEM blocks to ICESat laser altimetry in order to determine consistency and precision. Vertical accuracy is tested using Airborne Topographic Mapper (ATM) laser altimetry as reference data, which is also used as input to assign error to the orthophotographs.

Figure 1

Footprints of aero-triangulated photographs by year of recording.

The present DEM and orthophotograph were produced to facilitate investigations of ice sheet and glacier behavior. This is reflected in the studies where the data set has already been used, for example, the data is used to determine the geodetic mass change in northwest and northeastern Greenland, and resolved the ice sheet dynamics back to the late-1970s and mid-1980s^19–21. Moreover, to better understand climate-related variability of glaciers in southeast Greenland, Bjørk et al.¹¹ used orthophotographs as reference for co-registration of a wide range of imagery spanning 80 years. Khan et al.²² determined elevation change and frontal positions of Helheim and Kangerlussuaq Glaciers using heights from the DEM and orthophotographs. Lea et al.²³ reconstructed Little Ice Age glacier geometry of Kangiata Nunaata Sermia in southwest Greenland. Finally, Kjeldsen et al.²⁴ used this data set to reconstruct the geodetic mass loss of the GrIS since the Little Ice Age.

The DEM and reliability mask are made available in GeoTIFF file format, while orthophotograph and footprint coverage are made available in JPEG 2000 (jp2) and ESRI shapefile file formats, respectively.

Methods

The materials used for our DEM and orthophotographs are the 1978–1987 aerial photos and the coordinate lists resulting from the aero-triangulation providing the geodetic control of the photographs. The super-wide-angle photographs were recorded at flying heights of 13,000 m using a WILD RC10 camera with a nominal focal length of 88 mm¹⁶. About 45% of the photographs were scanned on a photogrammetric film roll scanner for this work, the remainder having already been scanned on a photogrammetric film roll scanner for earlier mapping projects. Scan resolutions of 14 and 15 μm provide a nominal ground resolution of 2.1 and 2.25 m, respectively. The coordinate lists, image observations and camera calibration reports are available from the Agency for Data Supply and Efficiency¹⁶. The coordinate list contains the ~21,500 Greenland Reference Frame 1996 (GR96) coordinates and heights from the aero-triangulation of the ~3,500 photographs. We use the latest aero-triangulation from 2006 for control of the aerial photographs from which we derive our DEM and orthophotograph products¹⁸.

Aero-triangulation

The ground control points for the aero-triangulation are a subset of ~6,300 stations in the terrestrial triangulation network which had GR96 reference system coordinates assigned by GPS-based adjustment²⁵. GPS stations from the Greenland Reference Network (REFGR) were used in the network adjustment allowing for a better accuracy of coordinates than previous methods permitted^14,19,25. Coordinates in the triangulation network were recalculated, which improved accuracy to 25 cm in the denser part of the network, and 75 cm in the less dense parts^19,25. The NAVSAT/TRANSIT Doppler stations, an elder form of satellite navigation and positioning system, were recalculated to GR96 coordinates. In general, GPS stations have been deployed in the areas covered by 1985 photos in west Greenland including Inglefield Land in the northwest, while Doppler stations have been deployed in the rest (Fig. 2). The area covered by the 1978 aerial campaign of northern Greenland has little to no terrestrial triangulation network, thus in support of the 1978 aerial photo campaign, a network of Doppler stations was deployed.

Figure 2: Ground control used in the aero-triangulation of the photographs.

Triangulation points indicate ground control with coordinates and heights assigned by GPS or Doppler-based adjustment of the terrestrial triangulation network. Ground control outside of these networks is not connected to triangulation networks.

The adjustments resulted in a coordinate list, where stations had been assigned with new GR96 coordinates and heights, a subset (Fig. 2, Table 1) is used to provide geodetic ground control points for the aero-triangulation of the aerial photographs. This input list has been supplemented with special height control in the form of same height at lakes (the horizontal lake surface is used to obtain dz=0 observations) and zero height (h=0) at the coastline (Table 1).

Table 1 Number of stations, observations and redundancies

The aero-triangulation procedure was divided into in six regional parts (1–6), in which the geodetic control was not fixed in the adjustment. Here, we rename these six parts by geography for clarity, so that part 1 is west (W), part 2 is south (S), part 3 is southeast (SE), part 4 is east (E), part 5 is northeast (NE), and part 6 is the north (N). In order to minimize tensions in the overlapping border zones, the stations in the zones were re-adjusted, while maintaining stations outside the overlap remained fixed (Table 1). The resulting mean error on coordinates and the impact of the borders is discussed in the section Internal Validation.

In the aero-triangulation, there is a good overall redundancy per station, although a few stations have low redundancy, thus the ability to detect observation errors is very good which makes the triangulation reliable.

Coordinates and heights used for the aero-triangulation were assigned a priori errors in order to weight the triangulation¹⁸. Image observations were assigned an error of 10–20 μm, while coordinates were assigned an error of 8 m. The much larger assigned error on coordinates relative to those of the survey network reflects the accuracy with which the stations (ground control) can be found in the photographs, and not the accuracy of the surveyed terrestrial network¹⁴. GPS heights have been assigned errors of 0.5 m, while Doppler heights range between 1 and 7 m. This is due to multiple sources of height information for these stations. The remaining triangulated heights are assigned an error of 3.5 m. Special height control, i.e., zero heights and same heights, have been assigned errors of 0.9 and 5 m, respectively.

Variance estimation is shown in Table 2, where we can assess the relationship between assigned a priori and a posteriori error estimates. Values close to one indicate agreement, however this is also dependent on the degrees of freedom in the adjustment (not shown here), which is reflected in the result. There is no weight normalization as assigned errors are based on empirical knowledge and experience, and therefore maintained. The mono errors in the regions E to N (parts 4 to 6) are smaller due to the higher a priori mean error of the mono observations. Weights are reduced for about 100 stations using thresholds of 3–4σ, to decrease the influence of large errors on the triangulation.

Table 2 Variance estimation

The ground control forming the basis of this aero-triangulation is the product of various instruments and methods. Using fixed GPS stations, DTU/SDFE recalculated the terrestrial triangulation network. Stations determined by the Doppler technique were also recalculated. Although the aero-triangulation is reliable (internally consistent) and well executed, it is also very dependent on the quality and density of the ground control.

DEM and orthophotographs

The aero-triangulation was set up with BAE System’s digital photogrammetric application SOCET SET 5.6 for DEM and orthophotograph generation. Projects were set up in SOCET SET for each UTM zone covering Greenland (zones 19–27). Coordinates were transformed for each UTM zone and heights were transformed WGS 84 ellipsoid heights using KmsTrans2012 (ref. 26). The heights in the 25×25 m DEM grid were measured by the NGATE (Next Generation Automated Terrain Extraction) module of SOCET SET 5.6 (ref. 27). Two strategies were used for automated collection of the elevations: the standard adaptive (ngate.strategy) and a low contrast (ngate_low_sp.strategy). The latter strategy is better at sampling heights where contrast is low in the photographs. Settings were set for highest precision. The two grids were merged to one grid with priority for elevations with the highest Figure-Of-Merit (FOM), a value assigned to each height indicating how reliably it is measured. Heights and FOM values were processed to GeoTIFF files cropped along latitude and longitudes. Heights are not average heights for each 25×25 m DEM post, but rather they should be considered spot heights, i.e., the height of the 2×2 m pixel at the post²⁷.

The orthophotographs are sampled using the most nadir strategy and bilinear interpolation. There is a 100 m seamline feathering along the seamlines and the photographs have been radiometrically balanced using a custom dynamic range adjustment. Geospatial Data Abstraction Library (GDAL) was used to post process the orthophotographs into jp2 format²⁸.

Data Records

DEM, orthophotograph, and supporting data are archived at the National Oceanic and Atmospheric Administration (NOAA) National Centers for Environmental Information (NCEI), with access details provided in the Data Citation 1. The file formats used are GeoTIFF, jp2 (JPEG 2000), and ESRI shapefile. In addition, both raster and vector formats have an ESRI projection file. Examples of the data products are shown in Fig. 3.

Figure 3: Examples of the data products from the head of Nuup Kangerlua (Godthåbsfjorden).

(a) Orthophotograph, (b) hillshade DEM, and (c) reliability mask where white=measured heights, and black=interpolated or outside of boundary.

Digital elevation model

The 25×25 m gridded DEM heights are stored in 32 bit floating point GeoTIFF file format. UTM zones, predominant year of photography, and file size considerations have been used to subset the DEM into multiple files. GeoTIFF file sizes have been kept under 4 Gigabytes for user friendliness. The structure of the DEM filenames is aerodem_*year*_*utm zone*_*subset*. tif, where aerodem indicates DEM data, *year* is the predominant year of photography in the file, *utm zone* indicates native projection of the file, and *subset* indicates if DEM has been cropped due to file size considerations (Table 3).

Table 3 Filename structure of the DEM GeoTIFF files

Reliability mask

The 25×25 m gridded reliability information is stored as 8 bit unsigned integer GeoTIFFs. The naming convention is identical to the DEM names in Table 3, but prefixed with rm. The Figure Of Merit (FOM) is a numerical value with a range of 0–100 assigned by the terrain extraction process. It may indicate one of three things for a given post measurement: It may be an error flag value, indicating that the automatic measurement was questionable; it may indicate a successful or good measurement; it may be an edit flag value, indicating the type of editing that was used (lake-filled, interpolated, etc.). As shown in Table 4 (available online only), FOM numbers greater than or equal to 40 indicate successful automatic correlation. These large FOMs are proportional to the correlation coefficient, so the larger the number, the better the measurement²⁹. In summary, values in the range of 2–21 have been interpolated; values in the range 22–38 are manually edited or LIDAR points, which do not occur in the reliability mask; the value 39 is the largest value assigned to posts that did not automatically correlate; and, finally, values of 40–99 are assigned to posts which did automatically correlate with increasing quality of correlation.

Table 4 Definition of FOM values

Digital orthophotograph

The 2 m resolution digital orthophotographs are stored in 8 bit JP2/JPEG 2000 (JP2ECW) file format. The naming convention of files containing the orthophotographs is identical to that of the DEM file names, with the exception that aerodem is replaced with g150 in the file name.

Photograph ground footprint coverage

The ~3,500 photograph footprints are stored as a single ESRI polygon shapefile with reference UTM 24 / WGS 84. The dates are stored for each footprint polygon in format yyyymmdd.

Technical Validation

In the absence of other measures of accuracy, the RMS result of a bundle block adjustment of aerial photographs is sometimes reported as an estimate of the accuracy of a photogrammetric model and its DEM product^1,21,30. Motyka et al.² and Khan et al.²⁰ reported the RMS result, but subsequently used a vertical accuracy estimate based on testing of the DEM product using ATM laser altimetry as reference data. Kjær et al.¹⁹ used the same method, while Howat et al.⁸ applied the same method but used ICESat laser altimetry as reference to validate the GIMP DEM.

We describe and assess the precision, or model fit, of the bundle block adjustment by analyzing the spatial prediction mean error. Instead of a relatively few number of checkpoints to determine the accuracy, we use external reference data in the form of laser altimetry to validate our aerophotogrammetric DEM. We co-register entire DEM tiles and then 50×50 km DEM blocks to ICESat satellite laser altimetry to determine horizontal and vertical co-registration vectors, and RMSE. Last, we test the DEM against ATM laser altimetry to test vertical accuracy—a result we also use as input into our estimate of the horizontal error of our orthophotographs.

Internal validation

The RMS result is a measure of model fit or precision, since the covariance matrix is an expression of the geometry of the survey network and the measurement accuracy. This assumes that the functional and stochastic models are correct: the latter assumes error to be normally distributed, observational gross errors to be down weighted or eliminated using blunder detection, and that variance estimation of groupings of observations is acceptable. There should be adequate redundancy for testing the reliability of the adjustment (Methods section).

We assess the quality and properties of the photogrammetric adjustment by examining maps of the mean horizontal and height errors (see Fig. 4). Error on coordinates and heights are found from the diagonal elements of the a posteriori covariance matrices resulting from the bundle block adjustment.

Figure 4: A posteriori mean errors from the aero-triangulation.

Mean errors for both ground control points (GCP) and tie points are shown. (a) horizontal mean errors, (b) mean errors on height. Plot of the result files. Modified after Engsager et al.¹⁸.

In Fig. 4 several patterns are apparent. First, the mean error increase with distance to ground control, i.e., the errors increase towards the interior of the ice sheet and ice caps away from the geodetic control found in the ice-free areas. The borders between the fixed and free DEM adjustment blocks are clearly seen due to the readjustment in the overlapping areas. The horizontal error map clearly shows that a section in southeast Greenland has no horizontal control, but only zero height control (Figs 2,4a).This can be recognized by the large horizontal mean errors. The color scale has been saturated above 10 m for visualization purposes. Out of the total of ~21,500 stations, 183 (0.9% of all stations) have height errors greater than 10 m. Errors larger than 10 m have a maximum at 35 m with a mean error of 14 m. In addition, 138 stations (0.6% of all stations) have horizontal errors greater than 10 m, with a maximum of 19 m and a mean of 12 m.

External validation: DEM co-registration to ICESat laser altimetry

This section describes the external validation of the aerophotogrammetric DEM tiles and 50×50 km blocks of the DEM around Greenland using ICESat laser altimetry, and is an expanded validation of that described in Kjeldsen et al.²⁴. We use ICESat³¹ laser altimetry for our validation as it is one of the most consistent global elevation products available³². Despite the lack of full spatial coverage, each ICESat footprint of ~60–90 m returns an elevation estimate related to the histogram peak of elevation in each footprint and are often at the decimeter precision on smooth flat surfaces^33,34, and often better than 1–4 meters on steeper topography depending upon surface roughness^3,32,35,36. For this particular study and validation, we use ICESat altimetry from the GLA12 Release 31 (ref. 31) product using the WGS 84 datum with ellipsoid heights in the native UTM projections that span from zones 19–27.

Our validation methodology is based upon co-registration methods described by Kääb³⁷ and Nuth & Kääb³² that relate the 3-D co-registration vector between two elevation surfaces to terrain slope (α) and aspect (ψ). The co-registration parameters are determined by robust least squares minimizations of stable terrain elevation changes between the DEMs and ICESat (dh) using the following equation:

$$dh=a⋠\cos (b−\mathrm{ψ})⋠\tan (\mathrm{Î\pm })+c$$

Where a and b is the magnitude and direction of the horizontal co-registration vector, respectively, and c is the mean vertical bias between the two elevation data sources. This method uses the relationship that mis-registrations between elevation data has with terrain characteristics. It however requires a sufficient sample of elevation differences on stable terrain (i.e., non-glacier) that contains some degree of terrain slope as the flat surfaces are unsolvable and ideally a uniform distribution of terrain aspects is available.

We perform the co-registration at two different scales, the first on each aero-photogrammetric DEM tile (23 tiles) and the second on a 50×50 km grid over all the DEMs. The rectangular 50 km block size is chosen in order to have enough ICESat footprint elevations to perform the co-registration, even for those DEMs that are more than 70–80% covered by glacier ice. A 50 km block covers about the same area as the stereo coverage of a 2×4 block of photos. All slopes less than 5 degrees are removed and a curvature filter is applied to remove regions where resolution variation between the datasets may cause spurious elevation differences. The significance of the co-registration solution is strongly dependent upon the number and distribution of stable terrain elevation change points, and solutions with less than 200 points are statistically susceptible to weak solutions.

The precision and accuracy of the co-registration depends strongly on the sample of input elevation differences, the terrain characteristics, and the resolution of the input data. Previous analysis has revealed that the method is capable to co-register to ~10% of the pixel size (resolution)^32,38. In our study, the limiting factor is the ICESat footprint size which varies from about 60 to 90 m, and thus we predict a horizontal co-registration precision of 6–9 meters. We test this by running a minimum Monte Carlo simulation by selecting randomly a sample of elevation differences to determine the co-registration solutions for each of the grid points. These tests revealed a Root Mean Square Error (RMSE) of 3 and 1 m for the horizontal and vertical co-registration parameters, respectively.

For co-registration of each entire DEM tile, all horizontal and vertical adjustments are less than 5 meters, which means that the DEMs are accurate to the precision of our method applied to ICESat data. For the denser grid, co-registration parameters are generally less than 15 meters horizontally and less than 10 m vertically. At the 1σ confidence level, the aerophotogrammetric DEM contains an accuracy of 10 m horizontally and 6 m vertically while the precision is better than 4 m (Fig. 5). The largest displacements between the DEMs and ICESat occur for those solutions that use rather small sample sizes (<200 elevation difference points), often those sections with the largest percent glacier cover or along the coast (Fig. 6). Although the large horizontal (and vertical) predicted mean error seen in the southeast (Fig. 4a) seems to be detected, it is questionable how confidently we can say this is due to the reliability of our method. On the southeast coast the test and reference data are a poor match, which is reflected in the small number of samples per DEM blocks. The entire southeast coast is a narrow strip of land (i.e., stable terrain), between the ice margin and the coast, and this is compounded by increasing spacing between ICESat tracks with decreasing latitude³⁹.

Figure 5: Histograms of the horizontal (a) and vertical (b) co-registration displacements for each 50 km×50 km grid cell show that the aero-photogrammetric DEM compilation is generally accurate to within 10 m horizontally and 6 m vertically with a precision greater than 4 m (1σ confidence level) (c).

The red bars show the fraction of displacements determined from 200 elevation difference samples or less.

Figure 6: Map of the horizontal (a) and vertical (b) components of the co-registration vectors between 50 km by 50 km sections of the aerophotogrammetric DEM compilation and ICESat laser altimetry.

(c) The RMSE of stable terrain differences after adjusting for the 3D mis-registration.

We plot the direction and magnitude of the co-registration vectors in Fig. 7. Vectors greater than 25 m plot along the coast or the ice margin, and are characterized by small sample sizes. Regional systematic error is apparent everywhere, but smaller and far less pronounced in the west and the south, where the ground control has the highest density. As density decrease, the ground control becomes more sensitive to the ability to identify them in the photographs and error on the coordinates, resulting in the larger magnitude and more regionalized pattern.

Figure 7: Magnitude and direction of the co-registration.

There is some spatial consistency of the vertical adjustments between the aerophotogrammetric DEM and ICESat, which is likely to be related to the density of the original input ground control that is used to constrain the aerotriangulation during the adjustment of the photogrammetric model.

As an additional verification of our validation methods we use two SPOT5-HRS DEMs generated from satellite stereo images acquired in 2008 and 2014 during the IPY-SPIRIT campaign³. The DEMs are centered on Kangerlussuaq (1981–2008) and Daugaard Jensen Glacier (1987–2014) on the central east coast of Greenland (IDs: GES_08-013_Kangerdlugssuaq_Glacier & GES_14-012_DaugaardJensen), each scene covers more than 10000 km² with about 3000 km² ice-free terrain. For validation, we co-register The IPY-SPIRIT DEMs to ICESat and to the G150 DEM. While the G150 products are generally well aligned to the ICESat framework, the IPY-SPOT DEM of Daugaard Jensen Glacier has a misalignment of about a DEM pixel (40 m), which is both captured by co-registration to the G150 products as well as to ICESat. The co-registration vectors are then triangulated (vector sum) to produce a residual between the three datasets that is not larger than 4.5 meters horizontally and 0.25 meters vertically. This test strongly verifies the precision of our co-registration methods combined with the precision of the datasets used here, both aerial and satellite based.

In summary, the aerophotogrammetric DEM is consistent with the entire ICESat acquisition, at least to an accuracy of ±6 m (1σ), which is possible from the native resolution of the input data, constrained by the ICESat footprint size. In terms of precision of the aerophotogrammetric DEMs, our comparisons on stable terrain (that contain steep slopes) show elevation difference standard deviations for each individual DEM tile and for our 50 km co-registration grid consistently around 4 meters, or less than one fifth of the native DEM resolution. This is the maximum conservative estimate since our comparisons are based upon ICESat altimetry, and thus these estimates are the combined precision of the ICESat and aerophotogrammetric DEMs. Nonetheless, there is ample evidence that these data products are not worse than the estimates provided by this external validation, and thus these quality indicators are therefore conservative estimates which may be significantly smaller when compared with more precise validation data.

External validation: Vertical accuracy using ATM

We test the DEM vertical accuracy with Airborne Topographic Mapper (ATM) laser altimetry from the period 1994–2014 (ref. 40). The accuracy of the reference data should be at least three times better than the data being evaluated⁴¹, and ATM data are therefore well suited for this purpose, as the accuracy and horizontal resolution of the ATM altimetry data are orders of magnitude better than our DEM⁴². ATM flights are flown to capture elevations on ice sheets and glaciers, thus for glaciological applications the test is often spatially close to the object of interest. The method captures the mean, random error from heights and the error induced from the horizontal displacement⁴³. Kjær et al.¹⁹ found a clear relationship between slope and error in the northwestern part of the present data set, which made it possible to model vertical accuracy on ice from ATM testing on ice-free terrain. We crop the ATM data to ice-free terrain using the land coverage of GEOGREEN2 (ref. 12) with a 50 m buffer and extract elevation differences using bilinear interpolation. The resulting raw elevation differences have a Greenland-wide 1σ=28.7 m. We remove blunders and outliers using an 86.1 m threshold (~3σ) on the elevation differences. Elevation differences are then filtered using the Figure-Of-Merit greater than or equal to 40, hence we assign error to measured pixels. Finally, we also remove slopes greater than 20 degrees to assess the impact of slope-induced error. In order to assess the spatial variability of error and not bias our analysis towards areas with a higher density of flightlines, we also assign the elevation differences to a 10×10 km grid (>20 samples) shown in Table 5.

Table 5 Elevation accuracy from spot height checks, resampling to 10×10 km grid, and mean horizontal displacement magnitude from the ICESat co-registration

The coverage of the 1985 aerial campaign (including the 1987 photos from the west) coincide with the best ground control, while in the southeast, which was flown in 1981, the terrain is steep and rugged and more slope induced error can be expected. Therefore, we assign a vertical accuracy by predominant year of photography.

Indeed, the impact from slope-filtering on SD_ATM (ATM SD) is greatest on the 1981 and 1987 regions. Thus, our data suggest that this variability is best explained by a combination of horizontal displacements and the more steep and rugged terrain. In Table 5, it is also shown, that even though the ground control has similar quality in the 1978 region as for 1981 and 1987, the error can be small. The difference is due to a much less sloping terrain in the north, thereby reducing the effect of slope-induced error. Arguably, this could also be explained by the high density of samples around the Thule Base¹⁹; however, the statistics for the 10×10 km grid are roughly the same as for the spot heights, thus underlining that the difference is largely caused by terrain differences. We use the SD_ATM as input to our horizontal accuracy of the orthophotographs and consider our result as conservative, due to the ~3σ outlier threshold and no slope corrections. SD_ATM has a spread of filtered elevation differences of 1σ=5–10 m.

Accuracy of orthophotographs

We calculate horizontal mean and maximum mean error using the vertical accuracy found from our ATM validation (SD_ATM)^11,44. The gridded DEM is TIN interpolated, but the merging of the two DEMs produced using different strategies, occasionally produces artefacts in the interpolation. We do not account for interpolation error and assign the larger maximum mean errors (6–13 m) to our orthophotographs as a conservative estimate, shown in Table 6.

Table 6 Horizontal mean and maximum mean error assigned to the orthophotographs.

Completeness

The ability to resolve height from a stereo image-pair is dependent on high-visible contrast, continuous and unambiguous surface textures⁶. This makes it challenging to produce photogrammetric DEMs of glaciers and ice sheets in high-latitude regions due to the low-contrast surface of ice and snow on glaciers, and the shadows created by mountains and valleys^6,8.

We find the completeness of two land coverages, glaciers and ice sheets (Ice/snow) and ice-free terrain/bedrock (Ice-free), by calculating the percentage of successfully resolved heights to the number of posts in each land coverage. Thus, completeness is the percentage of measurements with FOM values in the 40–99 range relative to FOM values in the 2–99 range (Table 4 (available online only)). We use GEOGREEN2 (ref. 12) map data to mask our two land coverage classes so that the coastline makes up the outer boundary and the edge of the innermost strips of aerial photographs makes up the boundary in the interior. The ICE coverage of GEOGREEN2 (ref. 12) is then used to differentiate between Ice/snow and Ice-free terrain/bedrock. The results are shown in Table 7, where it is evident that the ability to resolve heights is impacted by snow and topography. Visual inspection of the reliability masks reveal that low contrast on snow makes it difficult to resolve heights in the interior of ice sheet and glaciers, and in particular the G150 DEM products in the east (1987) and southeast (1981) are affected by this (e.g., Fig. 8). The topography is also more rugged in the east and southeast, with numerous deeply incised valleys and nunataks, which is also reflected in the results for the ice-free terrain.

Table 7 Completeness defined as the percentage of successful height measurements of the total possible in the photograph coverage

Figure 8: Decadal elevation change calculated using the G150 DEM and IPY-SPIRIT SPOT5-HRS³ DEM products.

(a) Kangerlussuaq Glacier 1981–2008 and (b) Dauggaard-Jensen Glacier 1987–2014. Note the elevation difference legend has been saturated at −30 m, and elevation difference transects are plotted in the insert show actual values. Analysis of elevation change on Kangerlussuaq Glacier can be found in Khan et al.²² and Kjeldsen et al.²⁴.

Summary of technical validation

For the internal validation we noted that the adjustment is reliable and consistent, and the spatial distribution of the mean error as expected reflects the dependence on accuracy and density of ground control and the adjustment strategy.

Our co-registration with ICESat laser altimetry shows, that at a 1σ confidence level the aerophotogrammetric DEM has an accuracy of 10 m horizontally and 6 m vertically, while the precision is better than 4 m. The spot height check corroborates this pattern with errors SD_ATM in the range of 5–10 m.

The quality of the ground control can be roughly be divided into two regions: The northwest, west and south, where coordinates and heights for a dense triangulation network have been determined using GPS based triangulation; in the southeast, east and northeast, the Doppler technique provides the coordinates and is supplemented with heights from various sources, and supported by local GPS and Doppler based triangulation networks. This pattern of deployed ground control is largely reflected in both of our external validations. Co-registration to ICESat altimetry reveals that the pattern of direction and magnitude of the co-registration parameters coincides with the spatial distribution of the ground control.

Variability in slope and ground control coincides with the reported error of the DEM obtained from comparison with ATM data. The range of the vertical accuracy is small (5–10 m), and comparable to the DEM co-registraton result. Therefore, we assign an overall accuracy of 10 m horizontally and 6 m vertically to the DEM with a precision better than 4 m, found from the DEM registration to ICESat data.

Usage Notes

Due to the merging of two elevation grids collected with two different strategies, interpolated heights are not always consistent. For elevation change analysis, we recommend using the reliability mask as a filter and treat elevations with Figure-Of-Merit (FOM) values less than 40 in the grid as outliers.

The ESRI shapefile containing photograph footprints with recording dates is recommend for assigning dates to the observed elevations.

Orthophotographs may occasionally look ’liquid’ due to localized poor interpolation. In case of uncertainty, the reliability mask can be used to ascertain measured point density in the area of interest combined with the DEM to assess orthophotograph reliability.

Perspectives

This circum-Greenland DEM product presents an important baseline geometric dataset for Greenland and the Greenland Ice Sheet. A number of studies have already shown the historic value of the dataset^{11,19–24}, especially for estimating the recent past changes of the ice sheet and will also prove important for future studies of geometry changes and geodetic mass balance. Figure 8 shows a comparison with two DEMs acquired from the IPY-SPIRIT campaign³ on the central east coast of Greenland. The decadal changes are clearly visible both on the larger outlet glaciers, but also on many of the smaller glaciers and ice fields surrounding. Despite limitations of photogrammetric products on surfaces with low visual contrast that fail to achieve accurate elevations at the highest positions of the glacier, the consistency between the three data sets is remarkable, stressing the precision and accuracy of both the G150 DEM and the SPOT5 satellite DEM products but also the ability to use ICESat as a reference framework. In summary, the G150 DEM presents the first high resolution, systematic and consistent, terrain products covering the entire Greenland coastline and therefore is an important baseline dataset for future scientific discoveries. It further provides an important time stamp from the 1980s for the frontal geometry of the entire GrIS, a useful and essential product for past and future glacier volume change estimations and accurate modelling constrain.