Co-registration and Bias Corrections of Satellite Elevation Data Sets for Quantifying Glacier Thickness Change

There are an increasing number of digital elevation models (DEMs) available worldwide for deriving elevation differences over time, including vertical changes on glaciers. Most of these DEMs are heavily post-processed or merged, so that physical error modelling becomes difficult and statistical error modelling is required instead. We propose a three-step methodological framework for assessing and correcting DEMs to quantify glacier elevation changes: (i) remove DEM shifts, (ii) check for elevation-dependent biases , and (iii) check for higher-order, sensor-specific biases. A simple, analytic and robust method to co-register elevation data is presented in regions where stable terrain is either plentiful (case study New Zealand) or limited (case study Sval-bard). The method is demonstrated using the three global elevation data sets available to date, SRTM, ICESat and the ASTER GDEM, and with automatically generated DEMs from satellite stereo instruments of ASTER and SPOT5-HRS. After 3-D co-registration, significant biases related to elevation were found in some of the stereoscopic DEMs. Biases related to the satellite acquisition geometry (along/cross track) were detected at two frequencies in the automatically generated ASTER DEMs. The higher frequency bias seems to be related to satellite jitter, most apparent in the back-looking pass of the satellite. The origins of the more significant lower frequency bias is uncertain. ICESat-derived elevations are found to be the most consistent globally available elevation data set available so far. Before performing regional-scale glacier elevation change studies or mosaicking DEMs from multiple individual tiles (e.g. ASTER GDEM), we recommend to co-register all elevation data to ICESat as a global vertical reference system.


Introduction
Applications of regional and global scale elevation products have increased substantially in geoscience.Surface elevation data are collected by many sensors using various techniques, and differencing between the multi-temporal elevation products is becoming a common method for monitoring surface changes, particularly of glaciers.The data are typically available as semi-continuous profiles or swaths of points, a network of points or a regular grid, the latter we will refer to as a Digital Elevation Model (DEM).There are three (nearly) global elevation products available to the public today.The Shuttle Radar Topography Mission (SRTM) flown in February 2000 provided continuous elevation data using interferometric SAR (InSAR) techniques (Farr et al., 2007).However, arctic coverage is limited since the mission only acquired data between 60 • N and 56 • S. The ICESat mission from 2003 to 2009 created the second nearly global dataset using space-borne Light Detection and Ranging (Lidar) (Zwally et al., 2002).The data are spatially limited to profiles of points rather than a continuous DEM, and arctic coverage is denser than at mid and low latitudes due to the polar orbiting strategy of the satellite.The third nearly global elevation dataset is the newly released ASTER GDEM based upon a composition of automatically generated DEMs from Advanced Spaceborne Emission and Reflection radiometer (ASTER) stereo scenes acquired from 2000-present (METI/NASA/USGS, 2009).In all of these datasets, errors and biases may persist from sensor instabilities, limitations of the techniques, bad surveying conditions on the ground and post-processing artifacts.The errors occur at a range of scales that directly affect measurement accuracy and precision, and increases the significance level a glacier thickness change requires for adequate detection through elevation differencing.
Today, the high temporal and spatial availability of spaceborne elevation data in remote areas where glaciers are present increases the capability to quantify glacier changes.Some studies use the data sets as they are, without correcting for biases between them (e.g.Rignot et al., 2003;Sund et al., 2009;Muskett et al., 2009) which may lead to flawed estimates of glacier volume changes or false-detection of surgetype behaviour (Berthier, 2010).The consequences of uncorrected biases in the previously named and other studies is not known to us.
Many of the data sets available to researchers today and those tested in this study are the result of second-level processing.This means that the conversion procedures between the original data acquisition (i.e.laser return waveforms, radar interferograms or stereo-imagery) to final elevation data is difficult to access and complicated and thus errors can not be anymore easily physically determined or modeled based upon the original transformation equations and acquisition parameters.Therefore, we use statistical approaches to analyze errors and to determine potential bias corrections.Even if physical modelling of errors might be preferable, an advantage of the statistical error modelling approach is that universal methods can be developed that may be widely applicable to different types of elevation data and irrespective of the sensor systems and processing procedures applied.
The most important correction is to co-register the two elevation data products such that the pixels of each DEM represent the same location and area on the Earth surface.This may be accomplished by minimizing stable terrain elevation residuals using a 2-D linear regression, e.g.iteratively shifting one DEM to the other within a ±5 pixel window (Rodriguez et al., 2006;Berthier et al., 2007;Howat et al., 2008).Other studies have corrected DEMs using single or multiple linear regression corrections between elevation and the location and terrain parameters (Gorokhovich and Voustianiouk, 2006;Bolch et al., 2008;Peduzzi et al., 2010).In terrestrial and airborne laser scanning, 3-D Least Squares Matching (LSM) is used to minimize the Euclidean distances between the points of point clouds, often allowing not only for shifts but also for rotations and scales between the two or more datasets (Gruen and Akca, 2005;Miller et al., 2009).Another commonly applied correction to DEMs is an elevation dependent bias adjustment (Berthier et al., 2004(Berthier et al., , 2007(Berthier et al., , 2010;;Paul and Haeberli, 2008;Kääb, 2008) which may have significant implications for glacier elevation changes because glaciers spread a range of altitudes which define their ablation and accumulation areas.Biases have also been found associated with the satellite acquisition geometry, specifically to satellite attitude parameters which may be significant enough to warrant a correction (Berthier et al., 2007).This type of correction will only apply to those data products where it is significant; e.g.satellite stereoscopic DEMs.

Objectives and case study locations
The motivation behind this study is to address the accuracy of comparisons between the globally available elevation data sets with particular attention towards detecting glacier thickness changes.This involves classifying and understanding the errors, and especially biases, associated with each of the data products and to suggest corrections that may improve the accuracy and precision of the differences.The first objective is to present a simple and effective universal method to co-register elevation products without the need for specialized software and with a high degree of automation.We argue that this method should be used as a first step in elevation comparison due to the varying location accuracies of the different sensors.In a second step, after centering the two data products to each other, analysis of remaining anomalies is compiled to detect both linear and non-linear biases, and to determine which errors require correction and how they affect glacier thickness changes.In contrast to the first step, the universal 3-D co-registration, the procedures applied in the latter steps are highly dependent on the sensor type and post-processing used for the elevation data.We will therefore only show examples for these secondary adjustments using ASTER satellite stereo as a scenario.
Two sites are chosen for this study.The first is the midlatitude high alpine region of the southern Alps in New Zealand.The region is chosen because of its alpine glacier characteristics, high elevation range, and availability of stable terrain from which to exemplify the biases and derive corrections.In this case-study, SRTM, ICESat, ASTER GDEM, and automatically generated ASTER DEMs from the US Geological Survey Land Processes Distributed Active Archive Center (LPDAAC) are compared.The second site is the high Arctic alpine region of Svalbard where ground control is limited to nunataks and along the strandflats.Automatically generated DEMs from ASTER and SPOT5-HRS are used in combination with ICESat and an aero-photogrammetrically derived DEM.

Stereoscopic DEMs
Stereoscopic DEMs are generated using photogrammetric principles.In along-track stereo, a parallax measurement gives the difference between the projected stereo rays of the same object onto the Earth's ellipsoid and can be converted to height provided the observer positions, the sensor pointing angles and camera parameters are known (Lillesand et al., 2004).Examples of satellite stereoscopic geometries are nadir and backward looking sensors (e.g.ASTER), forward and backward looking sensors (e.g.SPOT-5 HRS), forward, nadir and backward looking sensors (e.g.ALOS PRISM), or sensors that can be freely rotated to any stereo geometry The Cryosphere, 5, 271-290, 2011 www.the-cryosphere.net/5/271/2011/(e.g.Ikonos, WorldView, Pleiades), all of which are examples of pushbroom systems (line scanners).
Satellite stereoscopy is slightly more complicated than traditional photogrammetry from aerial frame imagery due to the typical pushbroom acquisition strategies and to the greater effect of Earth's rotation and curvature from the higher flying height of satellites (Toutin, 2004;Kääb, 2005).Image orientation may be solved from Ground Control Points (GCP) and a satellite orbital model (Toutin, 2004) that is implemented in commonly available software like PCI Geomatica ® .Automated approaches are becoming more common for deriving the relative and/or absolute orientation of stereo images using direct measurements of the satellite's attitude and position (i.e.pointing information, auxiliary and ancillary data) (for more details, see Schenk, 1999).The latter is the approach for both satellite stereo DEMs used in the this study: the ASTER DEMs produced by LPDAAC using the SilcAst software (product AST14) (Fujisada et al., 2005) and the SPOT5-HRS DEMs (Bouillon et al., 2006;Korona et al., 2009), as for instance available through the IPY SPIRIT (SPOT 5 stereoscopic survey of Polar Ice: reference Images and Topography) program.
The stereoscopic ASTER instrument, in orbit since 1999 aboard the Terra platform, contains a nadir and backward VNIR sensor (0.76-0.86 µm) separated by ≈30 • corresponding to a B/H ratio of 0.6 (ERSDAC, 2005;Toutin, 2008).The ground swath is 60 km while the image and reported DEM ground resolution is 15 and 30 m, respectively.The HRS instrument, aboard the SPOT5 satellite since 2002, contains forward and backward panchromatic sensors (0.48-0.7 µm), both 20 • from nadir providing a B/H ratio of 0.8 (Berthier and Toutin, 2008).The 120 km ground swath is twice as large as ASTER, with a ground pixel resolution of 10 m across track and 5 m along track, and a final DEM resolution of 40 m (Korona et al., 2009).
Errors associated with stereoscopic DEMs are related to the errors in the orientations of the stereo-scenes, either from GCP-based solutions or direct on-board determination, and to the ability of the matching algorithms to locate the corresponding points on two or more images.Errors in the parallax determination are both due to imperfect matching procedures and to imperfect image quality such as from lack of visible contrast, cloud cover, shadows and topographic distortions.Errors related to the parallax matching often result in blunders rather than bias, whereas errors related to the image orientation will typically induce bias.ASTER DEM uncertainty is reported to be typically within 15-60 m RMSE in the vertical depending upon terrain type (Toutin, 2002(Toutin, , 2008;;Kääb et al., 2002;Hirano et al., 2003;Kääb, 2005;Fujisada et al., 2005) and between 15 and 50 m horizontally (Fujisada et al., 2005;Iwasaki and Fujisada, 2005).SPOT5 uncertainty is reported to be between 10-25 m vertically (Berthier and Toutin, 2008;Korona et al., 2009) and greater than 15 m in the horizontal (Bouillon et al., 2006;Berthier and Toutin, 2008).In relationship to pushbroom sensors (e.g.ASTER and SPOT5 HRS), it has been shown that variation in the satellite's attitude induces biases within the raw images acquired as well as final DEMs produced (Leprince et al., 2007;Berthier et al., 2007).
The ASTER GDEM (the most recently released nearly global elevation product available) is a compilation of automatically generated ASTER SilcAst DEMs where the procedures for scene merging are not clearly defined.The validation summary reports, "that the ASTER GDEM contains residual anomalies and artifacts that most certainly degrade its overall accuracy, represent barriers to effective utilization of the GDEM for certain applications, and give the product a distinctly blemished appearance in certain renditions" (METI/NASA/USGS, 2009).The sources for the artifacts are residual cloud blemishes and the algorithm used to compile and generate the final DEM, the latter which is of most significance.Nonetheless, METI/NASA released this product publicly as an experimental/research grade product in hope of deriving a better Global DEM in the future.

Interferometric SRTM DEM
The Shuttle Radar Topography Mission (SRTM), launched in February 2000, mapped the Earth from 60 • N to 56 • S using single-pass synthetic aperture radar (SAR) interferometry (Farr et al., 2007).SAR interferometry uses the phase differences between two radar images acquired with a small base-to-height ratio.These phase differences are the photogrammetric equivalent to a "parallax" measurement allowing retrieval of topography (Rosen et al., 2000).We use the SRTM3, V2 without void filling (NASA et al., 2002).Many glacier elevation change studies have used this as a base dataset to compare to both newer and older data products (Rignot et al., 2003;Berthier et al., 2004;Larsen et al., 2007;Schiefer et al., 2007;Paul and Haeberli, 2008).Typically reported vertical uncertainties of the dataset are ≈ ±10 m which is lower than the mission standards of ±16 m (Rodriguez et al., 2006).However, vertical biases are present due to instability of the sensor and/or platform (Rabus et al., 2003), and elevation-dependent biases have also been shown due to penetration of the C-band radar waves (center frequency at 5.3 GHz) into snow and ice (Rignot et al., 2001;Berthier et al., 2006).Rignot et al. (2001) determined that the phase center of the C-band signal return was 1 to 10 m into the surface depending upon the snow conditions (i.e.dry vs. wet) in Greenland and Alaska.In Svalbard, the volumetric phase center of the C-band varied from ≈1 to 5 m along a profile from ablation to firn zones, respectively (Müller, 2011).Corrections for depth penetration are hardly used for the SRTM data, and is extremely difficult to correct for as knowledge of the snow conditions at the time of acquisition is required yet hardly available.We do not consider radar wave penetration in this study.

ICES at Lidar profiles
In 2003, the NASA Ice, Cloud, and land Elevation Satellite (ICESat) was launched with the Geoscience Laser Altimeter System (GLAS) acquiring elevation measurements in a 40-70 m elliptical footprint every 170 m (Zwally et al., 2002).ICESat obtained global coverage of elevations along profiles with a denser track sampling in the arctic due to the polar orbit.The rapid failure of the first laser invoked a curtailed orbital acquisition program.Nonetheless, the GLAS lasers operated for the following five years, collecting nearly two billion elevation point measurements before the last laser failed in November 2009.The altimeter has proven to be accurate to within ±15 cm over flat deserts (Fricker et al., 2005), and crossover track differences over low sloped glaciers on the order of ±1 m (Brenner et al., 2007;Moholdt et al., 2010a).ICESat products are freely available from NSIDC (www.nsidc.org),and are the third global elevation product publicly available and tested in this study.ICESat has been extremely successful for glacier applications in terms of elevation changes (Howat et al., 2008;Pritchard et al., 2009;Moholdt et al., 2010b) but also for determining the accuracy of newer satellite products (Korona et al., 2009) and older topographic maps (Nuth et al., 2010).Two second level products are available, GLA06 for smooth ice sheets and GLA14 for rougher terrain surfaces.The products vary by the number of Gaussian peak fits used to determine the maximum return-echo amplitude, maximum 2 and 6 respectively.Mean differences between the two products is ≈0.15 m though variations up to ±3 m occur (Kääb, 2008;Nuth, 2011).ICE-Sat release 531 is used for this study; the GLA14 products (Zwally et al., 2010b) are used for analysis of stable terrain whereas analysis of ice is using the GLA06 product (Zwally et al., 2010a).

Post-processing
The global elevation data sets used here are the result of the combination and post-processing of individual original data tiles, in particular SRTM (Rabus et al., 2003) and ASTER GDEM (METI/NASA/USGS, 2009).Among others, these procedures include vertical merging of overlapping elevations and horizontal mosaicking.Furthermore, the independent elevation data strips (in the SRTM and ASTER GDEM) may contain horizontal and vertical shifts of varying directions in which the merging process will entrain errors related to the varying mis-registrations.These steps make the original sensor geometry inaccessible and thus prevent the physical modelling of errors and error propagation.Similar problems arise also for other elevation data sets such as from airborne laser scanning or aerophotogrammetry, but usually at much lower levels if proper strip overlaps/adjustments and aerotriangulation procedures are applied.

Methods
To minimize the significance level a glacier elevation change requires for detection, we seek to analyze elevation differences on terrain assumed to be stable (e.g. off glacier) for 3 potential biases 1. the geo-location of the data (x, y, and z matrices), 2. an elevation dependent bias, and 3. biases related to the acquisition geometry of the data.
We analyze each bias individually and present solutions for each of these iteratively, rather than combining all into one full regression or co-registration adjustment (e.g.Gorokhovich and Voustianiouk, 2006;Racoviteanu et al., 2007;Peduzzi et al., 2010).The reason for that is to be able to follow and understand individual error terms, and to decide individually on their correction.Furthermore, it will become clear why a multiple regression based upon some combination of these terrain parameters will be significant, though such a correction may not be geometrically appropriate (e.g., see Peduzzi et al., 2010).
Elevation differences are calculated by re-sampling the spatial resolution of one of the DEMs to the other, or in cases involving ICESat, interpolating the DEM at the estimated centroid of the ICESat footprint.Bi-linear interpolation is used in both cases.Determining which DEM should be re-sampled to the other is subjective and will vary for each study.However, this decision should be well considered as differences in the corrections may occur depending upon whether one samples to the larger pixel size or vice versa (Paul, 2008).It could be worthwhile to check corrections by re-sampling in each direction to determine such influences.In the case studies presented here, the oldest DEM is re-sampled to the newest DEM.We use the population of assumed stable terrain elevation differences to analyze the quality of the comparison.Glacier and water pixels/points are removed using land and glacier polygon masks.In New Zealand, the glacier masks are derived from ASTER imagery (Gjermundsen, 2007) while the ocean and lake boundaries were downloaded freely from GADM database of Global Administrative Areas (http://www.gadm.org) (GADM, 2010).In Svalbard, the glacier masks are a part of the new digital Svalbard glacier atlas which is soon to be released by the Norwegian Polar Institute to the GLIMS dataset.The ocean is masked using data from the Norwegian Polar Institute mapping department.

A universal co-registration correction
Two DEMs of the same terrain surface that are not perfectly aligned experience a characteristic relationship between elevation differences and the direction of the terrain (aspect) that is precisely related to the x-y-shift vector between them.

The relationship between elevation error and aspect has been
The Cryosphere, 5, 271-290, 2011 described previously (Schiefer et al., 2007;Van Niel et al., 2008;Gorokhovich and Voustianiouk, 2006;Peduzzi et al., 2010), although corrections applied in the latter two studies were not analytical but based upon multiple regression adjustments to elevation.Gorokhovich and Voustianiouk (2006) showed the significance of the correlation between elevation differences and aspects on large slopes but overlooked the underlying cause as described in Kääb (2005).
The simplicity of this relationship and detection of unaligned DEMs lies in the similarity of elevation differences with the hillshade of the terrain (Fig. 1), a function that is based upon terrain slope and aspect.The correction of the mis-alignment requires a more detailed derivation.Figure 2 shows a schematic drawing and a real example where one DEM is shifted to the second.Resulting elevation differences (dh) are larger on steeper slopes due to the relationship of the magnitude (a) of the shift vector and the elevation errors to the tangent of the slope of the terrain (α): Additionally, dh are positive on eastern slopes and negative on western slopes, exemplifying the relationship to terrain aspect ( ).Because is usually defined circular from the north (azimuth), the direction of the shift can be modeled using a cosine of the difference between and the horizontal directional component of the shift vector.Combining this with the relation described by Eq. ( 1) derives the full analytical solution by relating the elevation differences to the elevation derivatives slope and aspect (Kääb, 2005): where dh is the individual elevation difference, a is the magnitude of the horizontal shift, b is the direction of the shift vector, α is the terrain slope, is the terrain aspect and dh is the overall elevation bias between the two elevation data sets.Slope and aspect can be calculated by any standard GIS  or mathematical software, and different approaches exist depending upon application.In this case, the finite difference method is more appropriate then the D8 method (Wilson and Gallant, 2000).To remove the error dependency on slope due to an x-y shift, we normalize the vertical deviations by dividing by the tangent of slope at that pixel.This produces a clean sinusoidal relationship between elevation difference and aspect (Fig. 2).The transformation of Eq. ( 2 where Three cosine parameters (a, b and c) are solved using least squares minimization where the amplitude of the cosine (a) is directly the magnitude of the shift vector, b is the direction of the shift vector and c is the mean bias between the DEMs divided by the mean slope tangent of the terrain (see Fig. 2).Because the solution to this actually analytical relationship is solved using the terrain which is not an analytical surface, the first solution may not be the final solution and iteration of the process is required to arrive at an ultimate solution.We choose to stop the iteration after the improvement of the standard deviation is less than 2% or if the magnitude of the solved shift vector is less than 0.5 m.The final correction is applied to the corner coordinates of the un-registered DEM by solving the x-and y-components of the shift vector from the magnitude (a) and direction (b).The mean bias determined by inverting Eq. ( 4) is added to the DEM using an estimate of the mean slope of the terrain used to solve Eq. ( 3).The internal consistency of the universal co-registration correction can be assessed providing three elevation data sets (e.g.Z 1 , Z 2 and Z 3 ) are available.Correction vectors can be solved between each of the datasets and the residual between the triangulation of these vectors is an estimate of the accuracy of co-registration.As an example: where − −− → Z 1 Z 2 is the correction vector from DEMs Z 1 to Z 2 etc., but may also be the elevation difference matrices themselves.Equation ( 5) states that the correction vector from Z 1 to Z 3 should equal the sum of correction vectors from Z 1 to Z 2 and from Z 2 to Z 3 .The residual of the triangulated shift vectors is the remaining un-removed shift between the three datasets and thus represents the internal consistency of the datasets.

Elevation dependent correction
An elevation dependent bias can for instance result from an uneven spatial distribution of the GCPs in the x-y-z-planes which leads to a poorly resolved stereo orientation further causing a distortion of the z-scale in the measurement of parallaxes.In these cases, either a linear or polynomial relationship between the elevation differences and elevation have been used to adjust the DEMs; e.g.: where Z is elevation, κ and τ are the regression parameters and n is the order of the polynomial (e.g. 1 for linear).The range of previously applied linear corrections varies from 1 to 40 m per 1000 m (Berthier et al., 2007(Berthier et al., , 2010;;Kääb, 2008).Figure 3 shows an example between a 2003 and 2002 ASTER DEM (as described in Sect.5) where a significant elevation dependent bias is apparent, which leads to a correction of ≈10 m per 1000 m (Table 2).An elevation dependent bias is also suggested to exist within the SRTM over non-glaciated terrain (Berthier et 2006,2007) and has been corrected for in some studies (Surazakov and Aizen, 2006;Schiefer et al., 2007), though this bias may also be the result of varying resolutions (see Paul, 2008).The SRTM is also expected to contain some bias due to penetration of the C and X Band radar waves into snow and ice, which has been suggested to be up to 10 m (Rignot et al., 2001).Either way, an elevation dependent bias is extremely significant for estimating glacier volume changes because a glacier and its mass balance varies predominantly with elevation, and thus a bias with elevation either from failure of the z-scale or from radar wave penetration into snow/ice will directly affect the measurement and interpretation of either glacier thinning or thickening.Linear bias with elevation causes either over-or under-estimated elevation changes of a shrinking glacier depending upon whether the bias stems from the newest or oldest topography, respectively (Berthier et al., 2006).

Along/cross track corrections
While the above co-registration and elevation-dependent bias are in principle universal to all types of elevation data, additional individual error characteristics apply according to the sensor type and method used for DEM generation.Along/Cross track biases are the errors associated with the satellite geometry, and may only be relevant for satellite stereoscopic DEMs.Few studies demonstrate that such along/cross track error exists.Leprince et al. (2007) showed that an along track pattern with a frequency of 11-12 cycles per scene existed within the geo-location of pixels of an ASTER scene, corresponding to the 11-12 tie points where the Terra satellite acquires specific attitude information (ERSDAC, 2007).They relate this bias specifically to the under-sampling of the pitch information.Berthier et al. (2007) find elevation biases of a SPOT5 cross-track stereo DEM of up to 12 m which they can reproduce using the highly sampled attitude measurements, specifically the roll in this example.To analyze these errors, we rotate the coordinate system from X-and Y -to cross (X track ) and along (A track ) track directions, respectively, using a preliminary along track angle (θ) estimated from the two corners of available data in the scene: Bias adjustments, if required, are fitted to these parameters using higher order polynomials, as described in the following sections.Section 5.2 provides an example of this bias and a correction using polynomials.
Errors related to the acquisition geometry are not restricted to stereo elevation data, but may also be present in interferometric DEMs.Height errors in InSAR generated DEMs generally derive from phase noise, atmospheric distortions and the imaging geometry (Knöpfle et al., 1998).In terms of geometry, the baseline length, along track position and platform height will all induce elevation errors within InSAR generated DEMs (Farr et al., 2007).

Error propagation
Errors within elevation data, whether a DEM or individual points, are commonly estimated by comparing to independently acquired GCPs, generally of a much higher accuracy than that of the elevation source being tested.The quantification of this error, assuming the GCPs are absolutely correct, typically uses 2 measures of statistical spread of the residuals, the Root Mean Square Error (RMSE) or the standard deviation (σ ), assuming Gaussian distribution of the residuals (randomness).However, if the mean difference of the residuals does not equal zero, then the RMSE is not a proper estimate of the statistical error distribution, and the mean and standard deviation should be reported (Li, 1988;Fisher, 1998).In this study, we do not use GCPs for accuracy determination, but rather create a residual population of the difference between two independent data sources over stable terrain.These residuals represents the relative errors between elevation data sets, rather than absolute.
Standard principles of error propagation are used for estimating errors between two DEMs (Burrough et al., 1998).For example, if one DEM has a random error, σ 1 , and the second DEM, σ 2 , then the resulting error of a statistically independent elevation difference point or pixel is defined as: However, elevation data, especially DEMs contain a degree of spatial autocorrelation that should be accounted for.The adapted error equation is then: where r is the correlation between σ 1 and σ 2 (Burrough et al., 1998;Etzelmüller, 2000).Determination of r requires semi-variogram analysis and advanced statistical procedures (Bretherton et al., 1999;Rolstad et al., 2009).When analyzing and quantifying glacier elevation changes, not just the spread of elevation changes is desired but rather the mean of the elevation changes over a particular area, e.g. a glacier or glacier zone.The standard error equation about the mean is defined (Davis, 2002), where n is the number of measurements.Two approaches to apply this equation to autocorrelated datasets are to use ε as defined in Eq. ( 10) or to use ε as defined in Eq. ( 9) and define n as the amount of un-correlated measurements.In the latter approach, some studies have assumed an autocorrelation distance of 0.1 km (Koblet et al., 2010), 0.5 km (Berthier et al., 2010) or 1 km (Nuth et al., 2007;Kääb, 2008).
www.the-cryosphere.net/5/271/2011/The Cryosphere, 5, 271-290, 2011 5 Results I: New Zealand case study The terrain of New Zealand varies from the high alpine characteristics of the Southern Alps to the low gentle slopes towards the southeast of the mountain chain.In this case study, comparisons of the global data sets cover both alpine and low-slope topographies while individual ASTER DEM comparisons are localized on four glaciers around Mt. Cook: Franz Josef, Fox, Tasman and Murchison glaciers.Franz Josef and Fox glacier are located on the west side of the mountain divide and consequently receive large amounts of accumulation due to the large eastwest precipitation gradient (Fitzharris et al., 1999) and experience high magnitudes of ablation (Anderson et al., 2006).These glaciers are generally quite steep with rather short response times (Oerlemans et al., 2005).Tasman and Murchison glaciers are located on the east of the divide and experience much less inter-annual variability of accumulation and ablation and have debris covered tongues (Kirkbride, 1995).

Global data sets
The SRTM DEM, ICESat and the ASTER GDEM seem to be well oriented to each other at the regional scale.The universal correction of ASTER GDEM to ICESat and the SRTM both resulted in a 10 m shift to the Northwest direction.The shift between SRTM and ICEsat was less than a meter.A triangulation of these three shift vectors resulted in error residuals less than 0.3 m in the x and y directions, and ≈1.5 m in the z-direction.However in this case, the solutions for the shifts are completely dependent upon the size of the area chosen for analysis due to the post-processed merging of individual tiles.The σ of vertical differences between the SRTM and ICESat is ≈5 m whereas σ of the vertical differences between ASTER GDEM with SRTM and ICESat is twice that (≈10 m) after co-registration.

Individual ASTER DEMs
This section refers to individual ASTER DEMs as computed by LPDAAC using the SilcAst software and onboardderived orientation parameters only (no GCPs; ASTER product AST14).We compare each possible combination of the data in Table 1 producing 10 differential DEMs.The population is first filtered using a 3 σ filter to remove the largest outliers.For each DEM pair, three potential adjustments are applied iteratively using the population of stable terrain difference pixels: 1. Co-register the DEMs using Eq.(3).Practically, we solve for the parameters (a, b and c) iteratively until the improvement of σ is less than 2%.We convert each iteration of the magnitude (a) and angle (b) of the shift vector into its x-and y-components using trigonometry and sum up the iterations to determine the final adjustment that is applied either on the corner coordinates or the x and y vectors of one of the DEMs.
2. Search and adjust for any elevation dependent bias.We use a robust regression of the elevation differences with elevation to solve Eq. ( 6) which is then used to correct one of the DEMs.
3. Search for any bias related to the acquisition geometry of the ASTER satellite.Here we search for biases that occur in the along and cross track directions of the satellite overpass.Higher order polynomials (6th to 8th order) are then fitted to the elevation differences with along/cross track directions which is used to adjust one of the DEMs.
Table 2 shows the results for each DEM pair before any adjustments are applied and after each correction is applied iteratively.In total, the three corrections improve the σ of stable terrain from 8-69%.The most significant improvement is obtained through co-registration of the DEM pair.Each individual ASTER DEM has a unique linear x-, y-and z-bias to the SRTM, independent of any other scenes.The direction of the shift is not uniform for all scenes which has important repercussions on the quality of the algorithms used to create the ASTER GDEM.
Figure 4 shows the processing sequence for differencing the 2006 ASTER DEM to the SRTM which we take to be the reference surface.About 20% of the ASTER scene is covered by semi-transparent clouds that result in erroneously high elevations in the DEM, which are removed by 3 σ filters on the elevation differences.Figure 4a and b shows the original elevation differences and their relationship with aspect which results in a shift vector of 30 m to the northeast.The elevation-dependent bias is not significant enough to warrant a correction (see Table 2).After co-registration (Fig. 4c), a visible pattern remains related to the along/cross track directions (Fig. 4d and e).We fit 8th and 6th order polynomials to the differences in the along/cross track direction, respectively, and adjust first along track before recalculating the cross track correction.The two adjustments applied to the ASTER DEM (1st -co-registration, 2nd -Along/Cross track) resulted in a 35% and 6% improvement Table 2. New Zealand Southern Alps.ASTER DEM and SRTM difference statistics on stable terrain (σ ) for the original population of elevation differences after adjusting for the mean and for each correction applied, in sequence.The parameter solutions for the corrections are given for both the co-registration and the elevation bias correction.The improvement of the standard deviation is the total improvement of all three corrections.The units for all parameters are meters except for κ which has the units meters per 1000 m.The most significant elevation dependent bias corrections occur in the 2001 and 2003 scenes where the corrections are as much as 10 m per 1000 m.In these scenes, the ocean covers ≈30% making the potential distribution of automatically generated tie points not uniformly distributed in space.Whether this refinement is performed within the DEM generation is not completely known to us, though it may provide an explanation to why these scenes contain large elevation scale distortions.Alternatively, the 2002 and 2006 scenes do not contain any ocean or significant distortions.However, an elevation dependent bias may be confounded with biases related to the sampling resolution (Paul, 2008).We test this on the 2001 and 2003 scenes by taking the inverse differences, i.e. resampling the ASTER to the 90 m SRTM DEM, and find that the slope of the elevation dependent bias does not change significantly though it is slightly smaller.This may be an indication that the ASTER 30 m products are actually representative for a pixel size more similar to the SRTM resolution of 90 m than the nominal one of 30 m.

Difference Original
In 9 out of 10 cases of Table 2, the along/cross track corrections improved σ of elevation difference residuals.However, in scenes that are missing a significant proportion of data (e.g., the 2001 and 2003 scenes contain ≈30% ocean), the spatial distribution of terrain differences are not uniformly distributed in the along and cross track directions.Therefore, the cross track direction is under sampled when the along track correction is completely sampled, and vice versa.This leads to overcompensation in either end of the along/cross track corrections, and therefore both along and cross track corrections are not applied in all cases shown in Table 2.
The noteworthy example for the along/cross track corrections is the difference between the 2006 and 2002 ASTER SilcAst DEMs.A shift of ≈45 m in a NE-SW direction is observed and corrected (Fig. 5a).An elevation dependent bias showed not more than ≈1 m per 1000 m which is not significant enough for correcting.Slight along/cross track biases are present up to ±5 m that are corrected using a 6th order polynomial (Fig. 5b and c).The post-correction pattern of elevation difference (Fig. 5d) reveal linear cross track features that run along track of the flight path of the ASTER scene.These features are similar to those discovered by Leprince et al. (2007), which they relate to jitter of the instrument and under-sampling of the sensor attitude information in the along-track direction (specifically, the pitch).The geometric correction of the ASTER pixels relies on a lattice of 12 by 11 points along/cross track, respectively, where precise satellite attitude measurements are acquired.A linear interpolation is used for geolocating all pixels in between these lattice points (ERSDAC, 2007).The number of cycles apparent in the mean vertical differences along track (Fig. 4d) is ≈10-12 cycles.The vertical amplitude of these variations is ±2 m.We choose not to correct for this bias as it is below the significance level of our dataset.However, if a very precise DEM is available to difference to (e.g.laser scanning), these higher frequency bias corrections will probably be above the significance level.The previous example shows that along/cross track biases exist within the ASTER SilcAst DEMs and that corrections can be applied with relatively good confidence.We find that along/cross track bias occur at 2 frequencies.A lower frequency pattern is the most significant with 2-3 cycles within an ASTER scene and an amplitude of up to 5 m.The cause for this lower frequency bias is unclear to us at the present time.Leprince et al. (2007) who also found jitter did not observe the lower frequency bias.In contrast to our data (i.e.LPDAAC), they use their own sensor model involving GCPs.The higher frequency bias has ≈10-12 cycles per scene and an amplitude of 2 m.The visibility of the higher frequency pattern confirms the lower frequency bias correction.The unrecorded pitch variations which are the hypothesized cause of the higher frequency bias occur independently for each scene acquisition.They are integrated into the DEM creation, most likely during the back-looking pass of the satellite because small variations (jitter) in the back-looking pitch cause slight variations in the looking angle directly affecting the vertical component of the parallax estimates.In this case, the unrecorded pitch variations of both stereo-pairs seem to have been in opposite directions and overlay each other constructively (i.e.added to each other) as otherwise the vertical variations would vanish (i.e.destructive superposition).
For detecting glacier elevation changes, co-registering the two DEMs to each other is the most important correction.Without co-registration, the elevation dependent and along/cross track bias corrections would not be as clearly visible and definable.The accuracy of the co-registration is shown in Table 3 as the 10 triangulated residuals (Eq.5) from the 5 datasets (SRTM + 4 ASTER DEMs).The coregistration solution has an internal horizontal accuracy of at least 1/3 of an ASTER DEM pixel (30 m), though often 1/10 of a pixel.The nominal vertical accuracy lies around 1-2 m, though 4-5 m in worst case scenario.
Table 3.The universal co-registration vector residuals and the RSS (Root Sum of Squares) of the total vector mean and standard deviations of the elevation change residuals as solved through triangulation of three DEMs.dz and σ dz are the mean and standard deviation of the triple vertical difference in the DEMs.These estimates represent an internal coherency between the three datasets that reflect the residual shift resulting from uncertainties in the solution of the universal co-registration correction.2000-2002-2003 8.5 −3.5 −4.6 10.3 −4.9 4.9 2000-2001-2002 −8.3 6.4 4.9 11.5 5.4 5.5

The ASTER GDEM
The statistics presented in Sect.5.1 about the ASTER GDEM are similar to those from the validation summary (METI/NASA/USGS, 2009) with biases of up to 10 m and RMSE of 5-50 m.Analysis of the spatial distribution of the elevation differences between SRTM and the GDEM reveal large-scale linear features (Fig. 6) which are highly related to the number of images used in the GDEM for a specific location (METI/NASA/USGS, 2009).Simple inspection of the false hillshade (as shown in Fig. 1) of the elevation differences between the GDEM and SRTM reveal the multiple directional shifts within the product.In addition, consistent bias persists over distances of 10-20 km with multiple slightly sinusoidal patterns of amplitudes of up to 10-20 m.Furthermore, constructive superposition of the high frequency bias is easily visible within the low slope regions as denoted by the white ellipses in Fig. 6.Therefore, the bias of individual ASTER scenes is incorporated into the GDEM and future compilations of a GDEM will benefit greatly from first co-registration of the individual tiles, and second removal of both the high and low frequency biases apparent in the individual ASTER SilcAst DEMs.In terms of glaciological research, the GDEM may be an appropriate data source for deriving area-altitude distributions of glaciers, which can be useful for volume change estimation using ICESat (Moholdt et al., 2010b) or for providing elevation input data required for spatial mass balance modelling.

Glacier elevation changes
The detection of glacier elevation changes is dependent upon both glacier characteristics and data quality.Assuming no bias and a precision of ±15 m for each ASTER SilcAst DEM (i.e.σ of off glacier terrain differences presented in Sect.5.2), the uncertainty associated with an individual difference pixel is ±21 m (using Eq. 9) for a single year DEM difference and ±3.5 m yr −1 for a 6 year difference.For the glaciers on the east side of the divide, elevation change rates of the tongue have been ≈ −1 m yr −1 between 1890-1964/1971 (Skinner, 1964;Hochstein et al., 1995) and more recently up to −4.5 m yr −1 from 1986-1990 (Quincey and Glasser, 2009).This is at the limit of what can be detected given the data in Table 1.In the 2000-2006 time epoch, the here estimated frontal thinning of the glaciers on the east side of the divide is between 1 and 4 m yr −1 (Fig. 7).The elevation changes also suggest that Murchison glacier experiences more rapid frontal thinning than the Tasman glacier.
On the west side of the divide, a different story persists where Fox and Franz Josef glacier experience slight thinning (1-2 m yr −1 ) at the highest elevations and thickening (5-10 m yr −1 ) at the glacier fronts.This is consistent with recent measurements of velocity covering the same time epoch (Herman et al., 2011).

Results II: Svalbard case study
The archipelago of Svalbard contains ≈34 000 km 2 of glaciers, about 60% of the land area.The availability of stable terrain is limited to nunatak areas between the glaciers and the strandflats around the coastline (Hisdal, 1985).A 2003 ASTER SilcAst DEM is tested against a 2008 SPOT5-HRS DEM from the IPY-SPIRIT Project (Korona et al., 2009), a 1990 aerophotogrammetric DEM from the Norwegian Polar Institute (description and accuracy of the dataset can be found in Nuth et al., 2007Nuth et al., , 2010) ) and 2003-2008 ICE-Sat data (Table 4).The 1990 dataset is partially incomplete with a missing strip over the center of the ASTER scene.This has few repercussions besides the along/cross track adjustments described in Sect.6.3.The landform characteristics within the ASTER scene is ≈65% glacier, 10% stable terrain and 25% ocean.

Universal co-registration correction
The four datasets in Table 4 allow the derivation of 6 shift vectors (Table 5).The aerophotogrammetric DEM and ICE-Sat (DI ) resulted in the smallest shift vector (≈3 m) and an RMSE (3.6 m) of stable terrain after two iterations.We expect the aerophotogrammetric DEM to be of the highest quality and accuracy, thus the impressive coherence with ICE-Sat further confirms previously published ICESat horizontal and vertical accuracies (Fricker et al., 2005;Luthcke et al., 2005;Magruder et al., 2005;Shuman et al., 2006;Brenner et al., 2007).For the other 5 comparisons, the SPOT5-HRS DEM compared better than the ASTER, with a shift vector The solution, SD and SI , of ≈20 m ( 1 2 pixel) and an RMSE of 8 and 5 m, respectively.All 3 shifts for the ASTER SilcAst (AD, AS, AI ) resulted in vector magnitudes of 80-100 m, or ≈2-3 times the pixel size.Figure 8 shows the vertical differences before and after co-registering the ASTER SilcAst DEM to the SPOT5-HRS DEM.The final fit solution was obtained after 3 iterations as opposed to 2 iterations common for all the ASTER DEMs tested in New Zealand.We additionally tested the two DEMs generated from the ASTER scenes acquired directly before and after the acquisition of the scene in Fig. 8 (i.e.same flight path).The shift vectors for these were all in the same direction and magnitude (not shown here).
The shift vector magnitudes for ASTER (2-3 pixels) is much larger than that of SPOT (0.5 pixels) in reference to ICESat and the aerophotogrammetric DEM.This reflects the more accurate satellite positioning and sensor pointing information of the SPOT5-HRS sensor as compared to ASTER.The elevation difference RMSE of the ASTER SilcAst products are double (≈20 m) those from the SPOT comparison to the aerophotogrammetric DEM or ICESat.This mainly reflects the different spatial image resolution, but presumably also the different stereo configuration (forward-backward) of the SPOT5-HRS sensor with a base-to-height ratio of 0.8 that provides a more precise parallax measurement than the nadir-backward configuration of ASTER (base-to-height ratio of 0.6).The results in Table 5 suggest that the cross-track positioning is less accurate than the along-track positioning.
Despite the spatial limitation of ICESat to ascending and descending tracks, it may still be used as a reference for any relative DEM, given a large enough distribution of points over stable nunataks.Schenk et al. (2005) showed the feasibility of using ICESat as ground control for historical vertical imagery and complimentary aerophotogrammetric DEMs by selecting visible nunatak areas and minimizing the vertical deviations of these areas through a 2-D regression.Figure 9 shows the first iteration for the comparison of ASTER to the aerophotogrammetric DEM (AD) and to ICESat (AI ).The sinusoidal relationship in both graphs are similar, though the variation in the relationship between AI is much larger due to the smaller sample size (less than 600 pts) of available stable terrain ICESat footprints (Table 5).
The internal consistency of the universal co-registration correction and the coherence between data as tested by triangulating the shift vectors is presented in Table 5.From the 4 elevation products and 6 shift vectors available between them, four error vectors can be calculated (Table 6).The lowest errors occur between the combinations SDI and ASD with horizontal positioning errors of less than 5 m whereas larger errors of ≈10 m occurs in the combinations ADI and AIS (Table 6).This difference is mainly caused by one less www.the-cryosphere.net/5/271/2011/The Cryosphere, 5, 271-290, 2011 significantly defined shift vector, AI .Figure 9 shows that the ASTER to ICESat comparison is noisier due to a smaller sample size (≈600) and spatial distribution of stable terrain elevation points.The solution to Eq. ( 3) is weaker than other solutions involving ICESat; for example, SI contain more than 6000 stable terrain elevation differences and the distribution of these differences with aspect are a lot more uniform than that of AI (Fig. 9b).Nonetheless, despite the limited number of points used for AI , the correction to ICESat was still as precise as 1/3 an ASTER pixel (Table 6).

Along track bias correction
Visual inspection of Fig. 8 reveal suspicious elevation difference trends, particularly the positive changes along the southwest part of the scene and negative changes in the northeast.Conveniently, ICESat acquisitions from the same year as each of the DEMs with a repeat track from 10 October 2003 and 3 March 2008 is available.This repeat track has a cross track separation of ≈15 m and is similar to the along track direction of the ASTER satellite overpass.The comparisons between the ASTER and SPOT5-HRS DEMs with the ICE-Sat profiles are shown in Fig. 10.Despite the extreme limitation of stable terrain throughout the ASTER scene, we detect an along track bias of up to ±10 m between the ASTER and SPOT5 DEMs shown in Fig. 10a.The differences between the ASTER and the 2003 ICESat track is similar to the along track bias (Fig. 10d).Along track biases are not present between the ICESat track and the SPOT5 DEM.The slightly negative mean difference is probably a summer melt signal, especially significant in the first 5 km of the profile which ascend the front of Storbreen glacier with an apparent 5 month loss of ≈2-2.5 m (Fig. 10d).After correction, the 2008-2003 differences between the DEMs are similar to the ICESat repeat track differences (Fig. 10e).This example proves the significance and feasibility of such corrections to the ASTER DEMs, even where stable terrain is limited (less than 10% of the scene) and distributed unevenly over the scene.Corrections of the along/cross track biases seem to remove most of the spatially visible trends between the ASTER and SPOT (Fig. 11).We use the same along/cross track corrections to difference the 2003 ASTER DEM with the 1990 aerophotogrammetric DEM because the missing strip in the 1990 data may weaken the significance of along/cross track biases.The mean bias between the adjusted 2003 ASTER DEM and the 1990 DEM (−0.7 m) is therefore corrected.

Glacier elevation changes
Svalbard glaciers, as opposed to New Zealand glaciers have much lower rates of ablation and accumulation.The elevation change rates of the previous decades are typically between −3 and +1 m yr −1 (Nuth et al., 2010;Moholdt et al., 2010b).We denote a number of significant anomalies and glacier trends within

Conclusions and perspectives
The aims of this study were to detect, analyze and statistically correct the various errors and biases that exist within publicly available terrain elevation products.We present a simple and robust co-registration method for DEM pairs using the elevation difference residuals and the elevation derivatives of slope and aspect.The method represents the complete analytical solution of a 3-D shift vector between two DEMs.The solution to Eq. ( 3) returns statistically significant results for situations when full continuous surface residuals are available but also when stable terrain is limited to less than 10% of the scene and in comparisons between a DEM and the spatially limited ICESat elevation tracks.By triangulating the co-registration residuals between three elevation data sets, we estimate an internal precision of at least 1/3 but up to 1/10 of an ASTER or SPOT5 pixel in the horizontal and between 1-4 m vertically.The co-registration accuracy increases with availability of stable terrain.In this study, ≈600 difference points between ICESat and ASTER effectively co-registered the data products to at least 1/3 of a pixel.The improvement in σ of elevation residuals through co-registration amounted to 5-70% depending upon the magnitude of the shift vector.We suggest that co-registration be tested and, if necessary, performed whenever elevation differencing is used for estimation of glacier changes, and other terrain changes.The magnitude of the bias induced by not co-registering is directly related to the direction and magnitude of the shift with the direction and slope of the glacier surface.That implies that for very flat glaciers a correction effect might be small, but also that the correction effect for steeper glaciers might far exceed the signal intended to detect.Unless there is a perfectly random distribution of (glacier) slopes and aspects within a study area, omitting to correct a significant shift will not only result in an increased RMSE of the elevation differences, but induce a systematic vertical bias, which can easily be estimated from Eq. (3) for a given shift magnitude and direction.Our co-registration method is advantageous to the method of minimizing the elevation differences by iteratively shifting one DEM to the other because it only requires 2-3 iterations as opposed to more than 20 iterations for the latter method.The two methods theoretically result in the same co-registration parameters, given use of the same measure of minimization (i.e.RMSE, σ , r 2 etc.).However, the latter method requires interpolation and pixel down-sampling to return sub-pixel adjustements, whereas our method is independent of pixel size.In addition, our method directly provides error estimates on the co-registration parameters of magnitude and direction (a and b in Eq. 3), which can help distinguish flawed solutions and help detect potential multiple shifts between the data within the automated co-registration technique.
In this study, large elevation dependent biases occurred within the ASTER DEMs that covered less than 70% of the land surface.This may imply that the spatial and elevational distribution of automatically generated tie points affects the tuning of the stereo model within the automated process.It is difficult to determine whether the SRTM has a significant elevation dependent bias; all tests were not as convincing as Fig. 3.An elevation dependent bias caused by penetration of the SRTM C band radar is however much more significant.Determination of this type of bias is out of the scope of this paper.More research should certainly be focused on for example, comparing glacier DEMs created at roughly the same time as the SRTM to analyze the magnitude of this bias.
Significant along/cross track biases are specifically found within the ASTER SilcAst DEMs.Longer frequency along track biases contain amplitudes as large as ±10 m which we adjust using 6-8th-order polynomials.A higher frequency bias is detected with ≈10-12 cycles per scene which may be related to the under-sampled pitch of the backward looking sensor, similar to those found with the nadir looking camera (Leprince et al., 2007).The amplitude of this bias is 1-2 m, which we regard as under the significant limitation of our statistical adjustments.It is important to note, that, since every ASTER SilcAst DEM individually is affected by these high-frequency variations, a differential DEM might contain in the best case a destructive superposition of these variations (i.e.error elimination), or in the worst case a constructive superposition (i.e.error maximization).A prime example is the ASTER GDEM where constructive superposition of jitter is clearly visible despite the compilation from numerous ASTER DEMs (Fig. 6).As a main conclusion from our study, we suggest a methodological framework (Fig. 12) for whenever DEM (or elevation) comparison is to be performed for glacier research.The first and most important step is to co-register the elevation data products.Our method can easily be implemented in free or standard geoinformation systems, spreadsheet softwares (i.e.excel), or standard programming environments such as MATLAB or IDL.The only functionalities necessary are: computation of DEM differences, slope and aspect; simple DEM attribute algebra (here dh / tan(α) ); curve-fitting including fitting of sines or cosines; DEM shifting.If no curve fitting functionality is available, the necessary shift parameters can easily be estimated from a scatter plot as shown in Fig. 2. The method can be fully automated.The correction of any further, secondary, biases is dependent on the individual sensor systems and DEM post-processing procedures.However, it should be noted that these biases can easily mimic real glaciological processes such as surges or mass-balance variations with altitude.
We found the ICESat-derived elevations to be the most consistent globally available elevation data set available so far.It could be used as reference to register DEMs to in any regional-scale study.This would lead to a consistent global reference frame for glacier elevation change studies.As a consequence, we recommend for instance, to consider within a new compilation of the ASTER GDEM to reference any individual ASTER DEM to ICESat elevations before merging these individual DEMs to the global data set.A similar procedure, at least for testing, might be appropriate for other ongoing or future global DEM projects such as TanDEM-X or SPOT5-HRS.

Fig. 1 .
Fig. 1.The elevation differences before co-registration (left) between ASTER DEMs in 2006 and 2002 from New Zealand (described in Sect.5.2 and shown in Fig. 5) are remarkably similar to the hillshade of the DEMs (right).The location of the subsetted region is depicted in the 2006 ASTER image (center).

Fig. 2 .
Fig. 2. Top: 2-D scheme of elevation differences induced by a DEM shift.Bottom: The scatter of elevation differences between 2 DEMs showing the relationship between the vertical deviations normalized by the slope tangent (y-axis) and terrain aspect (x-axis).The example is the DEM differences between the 2002 and 2003 DEM shown in Sect. 5.The equation for the solved sinusoidal curve fit is shown along with the three unknown solution parameters, a, b and c.

Fig. 3 .
Fig. 3. Example of elevation differences between 2003 and 2002 ASTER DEMs from Sect.5.2 before and after applying an elevation dependent bias correction using a 3rd order polynomial.The two DEMs were first co-registered before checking for an elevation dependent bias.Glacier masks are indicated by black outlines.

Fig. 4 .
Fig. 4. The processing sequence applied to the 2006 ASTER DEM as compared to the SRTM DEM.(a) Shows the original elevation differences while (b) shows the first iteration relationship to aspect and the cosine-fit solution.(c) Shows the elevation differences after co-registration and (d) and (e) shows the residual relationship to the along track and cross track directions, respectively, with the polynomial correction.(f) Shows the elevation differences after the along and cross track corrections.(g) Shows the final histograms of the three elevation difference maps of (a), (c) and (f).All glaciers, lakes, ocean and outliers due to clouds are masked out in black.
Fig. 5. (a), (b) and (c) show the first three corrections applied between two ASTER SilcAst DEMs from 2006 and 2002.(d) Shows the resulting elevation differences with a plot of the mean elevation differences along track.The linear cross-track features that run along track seem to have an amplitude of 1-2 m in the vertical direction.These fluctuations are thought to be induced by unrecorded pitch variations of the satellite, jitter.

Fig. 6 .
Fig. 6.Vertical differences between SRTM and the ASTER GDEM (a).Elevation changes exemplified for along (b) and cross track (c) profiles.Significant bias of 10-15 m is visible in (a), (b) and (c) which is highly related to the number of ASTER DEMs used to compile the GDEM.Also, the bias is consistent for lengths of up to 20-40 km as seen in (b) and (c).Apparent jitter is denoted in the difference image.

Fig. 7 .
Fig. 7. Elevation changes from 2000 to 2006 on the Fox, Franz Josef, Tasman and Murchison glaciers.Only six year changes are shown due to their level of significance.Fox and Franz Josef glaciers show slight thinning at upper elevations and thickening at the glacier fronts, while the Tasman and Murchison glaciers experience slight thinning that is at the limits for significant detection using these datasets.

Fig. 8 .
Fig. 8. Elevation differences before and after co-registration of the ASTER DEM to the SPOT5-HRS DEM.The final shift is ≈2.5 ASTER pixels to the west-north-west.The green lines are the ICE-Sat tracks.The glacier area is shaded with a transparent gray to emphasize the stable terrain differences.

Fig. 9 .
Fig. 9.The first iteration of the co-registration between the ASTER DEM and the aerophotogrammetric DEM (a) and the ASTER DEM and ICESat (b).

Fig. 10 .Fig. 11 .
Fig. 10.Elevation differences between the SPOT5-HRS and ASTER DEMs over the non-glaciated regions (a) and the glaciated regions (b).The graph inset in (a) shows the along track bias measured from the stable terrain and the 6th order polynomial correction.The green line is an ICESat repeat track from 10 October 2003 and 2 March 2008.The elevation profile from the 2003 ICESat track is shown in (c), The differences between DEMs and the ICESat track closest in time to the DEM acquisition is compared (d).The difference between the ASTER DEM and ICESat is similar to the bias correction as determined between the ASTER and SPOT DEMs.The elevation changes between 2008 and 2003 are shown in (e) before and after correcting for the along track ASTER bias.ICESat to ICESat differences are made by a simple along track interpolation as the cross track separation was not greater than 15 m, which is well within the footprint size of ICESat.

Fig. 12 .
Fig. 12.A suggested methodology for comparing DEMs or elevation products for glacier change detection.

Table 1 .
New Zealand elevation data type, date acquisition, resolution, and scene ID.

Table 4 .
Data sources used in Svalbard in Sect.6), their acquisition date(s) and resolution.

Table 5 .
Shift vectors between the 4 data types in Svalbard as tested in Sect.6.1.X, Y and Z are the 3 components of the full co-registration adjustment vector between the datasets in meters and the RMSE is calculated after correction.

Table 6 .
Error vectors revealed through triangulation.All units are meters.