Available online at www.sciencedirect.

com

Global and Planetary Change 59 (2007) 110 125

www.elsevier.com/locate/gloplacha

Evaluating digital elevation models for glaciologic applications:

An example from Nevado Coropuna, Peruvian Andes
Adina E. Racoviteanu a,b,, William F. Manley b , Yves Arnaud c , Mark W. Williams a,b
a
Department of Geography, University of Colorado, CB 260, Boulder CO 80309, USA
b
Institute of Arctic and Alpine Research, University of Colorado, CB 450, Boulder CO 80309, USA
c
Institut de Recherche pour le Dveloppement, UR GREAT ICE, Laboratoire de Glaciologie et Gophysique de l'Environnement,
38402 Saint Martin d'Hres Cedex, France
Available online 10 January 2007

Abstract

This paper evaluates the suitability of readily available elevation data derived from recent sensors the Advanced Spaceborne
Thermal Emission and Reflection Radiometer (ASTER) and the Shuttle Radar Topography Mission (SRTM) for glaciological
applications. The study area is Nevado Coropuna (6426 m), situated in Cordillera Ampato of Southern Peru. The glaciated area was
82.6 km2 in 1962, based on aerial photography. We estimate the glacier area to be ca. 60.8 km2 in 2000, based on analysis of the
ASTER L1B scene.
We used two 1:50,000 topographic maps constructed from 1955 aerial photography to create a digital elevation model with
30 m resolution, which we used as a reference dataset. Of the various interpolation techniques examined, the TOPOGRID
algorithm was found to be superior to other techniques, and yielded a DEM with a vertical accuracy of 14.7 m. The 1955 DEM
was compared to the SRTM DEM (2000) and ASTER DEM (2001) on a cell-by-cell basis. Steps included: validating the DEM's
against field GPS survey points on rock areas; visualization techniques such as shaded relief and contour maps; quantifying errors
(bias) in each DEM; correlating vertical differences between various DEM's with topographic characteristics (elevation, slope and
aspect) and subtracting DEM elevations on a cell-by-cell basis.
The RMS error of the SRTM DEM with respect to GPS points on non-glaciated areas was 23 m. The ASTER DEM had a RMS
error of 61 m with respect to GPS points and displayed 200300 m horizontal offsets and elevation spikes on the glaciated area
when compared to the DEM from topographic data.
Cell-by-cell comparison of SRTM and ASTER-derived elevations with topographic data showed ablation at the toes of the
glaciers (25 m to 75 m surface lowering) and an apparent thickening at the summits. The mean altitude difference on glaciated
area (SRTM minus topographic DEM) was 5 m, pointing towards a lowering of the glacier surface during the period 19552000.
Spurious values on the glacier surface in the ASTER DEM affected the analysis and thus prevented us from quantifying the glacier
changes based on the ASTER data.
2006 Elsevier B.V. All rights reserved.

Keywords: Digital Elevation Model (DEM); Geographic Information Systems (GIS); Andes; Glaciology; DEM generation; ASTER; SRTM

1. Introduction
Corresponding author. Department of Geography, University of
Colorado, CB 260, Boulder CO 80309, USA. Tel.: +1 303 492 5546; Digital elevation models (DEM's) are beginning to
fax: +1 303 492 6388. see wide use in glaciological applications. Some studies
E-mail address: [email protected] (A.E. Racoviteanu). have used DEM's to extract components of glacier
0921-8181/$ - see front matter 2006 Elsevier B.V. All rights reserved.
doi:10.1016/j.gloplacha.2006.11.036
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 111

topography (slope and aspect), which were then com- historical research in glaciology. In such areas, digitized
bined with satellite images to map glacier areas (Klein elevation contours from old topographic maps still con-
and Isacks, 1996; Duncan et al., 1998; Sidjak, 1999; stitute a ready source of historical data on glacier eleva-
Kb et al., 2002a; Paul et al., 2002). In addition, tion and area. However, there is no established
DEM's have been used as tools to derive hypsometry interpolation method especially suitable for creating
maps at different time steps and to quantify vertical continuous elevation data from these topographic maps
surface changes on glaciers in remote areas, as indirect for accurate representation of glacier terrain. The ac-
measurements of mass balance (Khalsa et al., 2004; curacy of various techniques to construct DEM's from
Berthier et al., 2004). digitized contour data has been addressed in GIS liter-
Several studies have explored ways to assess glacier ature, such as Burrough and McDonnell (1998) and
mass balance and volumetric change by using time Wood and Fisher (1993), but the glaciological commu-
series of digital elevation data. For example, Etzelmller nity has yet to agree on a suitable interpolation method.
(2000), Etzelmller and Bjrnsson (2000) and Etzel- For instance, Etzelmller and Bjrnsson (2000) used an
mller et al. (1993) discussed GIS techniques to Inverse Distance Weighted (IDW) interpolator to create
quantify changes in elevation, terrain roughness, glacier a continuous surface from radar profile lines on a gla-
hypsometry and flow patterns using grid-based DEM's. cier. Other authors (e.g. Mennis and Fountain, 2001)
Rentsch et al. (1990), Vignon et al. (2003) and Rivera chose a spline interpolation for representation of glacier
and Casassa (1999) estimated changes in glacier volume and sub-glacier topography. Alternatively, Gratton et al.
and mass balance based on reference elevations from (1990) chose a Triangular Irregular Network (TIN) de-
topographic data. rived from digitized contours to represent rugged glacier
The availability of new remote sensing platforms topography at the Columbia Icefield. The choice of
with high resolution, global coverage and low costs interpolation method depends on terrain topography and
provide the potential to calculate glacier mass balances the type of data analysis needed. So far, only a few
in remote areas with little existing glacial information, glaciological studies (e.g. Cogley and Jung-Rothen-
such as the Andes of South America. Two of the new hausler, 2004) provided a careful quantitative evaluation
sensors are the Advanced Spaceborne Thermal Emis- of interpolation accuracy over glaciated area using
sion and Reflection Radiometer (ASTER) sensor and ground data. At present we do not know how sensitive
the Shuttle Radar Topography Mission (SRTM). The glacial mass balance calculations are to the type of
ASTER sensor acquires simultaneous stereo images interpolation method used to create glacier elevation
from different directions, suitable for generation of surfaces.
DEM's. Glaciologists would like to evaluate changes in Here we assess the suitability of readily available
glacial mass balance over time by comparing changes in SRTM and ASTER datasets for mass balance studies at
DEM properties acquired at different times (Etzelmller, a remote mountain area in the Peruvian Andes. The
2000). Efforts are being undertaken to provide accuracy SRTM DEM's released by USGS and the ASTER
assessments of the new DEM's from SRTM and ASTER DEM's generated at the EROS Data Center constitute
imagery. For instance, ASTER-derived DEM's (30 m information in the public domain available to all
resolution) have been validated at several sites (e.g. glaciologists. Our objectives are: 1) to evaluate the
Welch et al., 1998; Lang and Welch, 1999; Kb, 2002; suitability of various interpolation techniques to con-
Hirano et al., 2003). The recently released SRTM-3 struct DEM's for glaciological studies; 2) to assess
(90 m resolution) datasets have been evaluated only at a elevation differences between DEM's from satellite data
few study sites, on non-glaciated terrain (e.g. Falorni and the DEM from topographic data, 3) to identify the
et al., 2003; Rabus et al., 2003). Therefore, it is unclear spatial distribution of these errors with respect to topo-
if the data available to glaciologists from these new graphic characteristics (elevation, slope and aspect) and
remote sensing instruments provide sufficient spatial 4) to ultimately distinguish a glacier signal from multi-
and temporal resolution to detect a glacial signal without temporal DEM's.
extensive calibration.
Combining these new satellite-derived DEM's with 2. Study area
DEM's constructed from topographic maps has the
potential to extend the time series of glacial change over Our study area is Nevado Coropuna (Fig. 1), situated
many decades. This ability to construct glacial mass in the volcanic Cordillera Ampato in Southern Peru
balances over many decades is particularly well-suited (15241551S latitude and 71517300W longi-
to remote areas such as the Andes where there is little tude) (Fig. 2). The Ampato range consists of 93 glaciers,
112 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

Fig. 1. Nevado Coropuna (6426 msl), Peruvian Andes.

with an estimated average glacier thickness of 35 m mode on both rock and glaciated areas (Fig. 2). The light
and a total glaciated area of 146.73 km2 based on 1962 and portable Garmin was preferred over the bigger
aerial photography (Ames et al., 1989). Nevado Trimble Pathfinder unit. However, a disadvantage of the
Coropuna is the highest peak in Cordillera Ampato, navigation-mode GPS is that it cannot be differentially
ranging from 4600 m at the base (Lake Pallacocha) to corrected. The accuracy of the Garmin unit was tested at
over 6400 m at the main summit, with gentle sloping Green Lakes Valley Long-Term Ecological Research
lava flows and glaciated terrain. There have been no (LTER) study site in Colorado (Ackerman et al., 2001).
comprehensive field measurements of glacial properties Horizontal errors of measurements taken with the
on Nevado Coropuna until recently. The only glacier Garmin unit were within 3.9 m of the differentially
mapping was carried out by Ames et al. (1989) who corrected data obtained with the Trimble Pathfinder
reported a glaciated area of 82.6 km2 for Nevado (Ackerman et al., 2001). As a rule of thumb, we consider
Coropuna based on planimetric analysis of 1962 aerial the vertical accuracy of the Garmin unit to be about 1.5
photography. This manuscript fills this gap and com- times bigger than horizontal errors ( 10 m). The GPS
plements two ice-core-drilling expeditions conducted on elevations referenced to the WGS84 ellipsoid were
Coropuna in 2003 by l'Institut de Recherche pour le converted to orthometric heights (heights above the
Dveloppement, France (GREAT ICE project) and the EGM96 geoid) by subtracting geoid heights calculated
Ice Core Paleoclimatology Research Group at the Byrd for each point based on latitude and longitude (Rapp,
Polar Research Center, Ohio State University. Results 1996).
from these paleoclimatic studies will provide an isotopic
record for the region, which can help understand the 3.2. Construction of the DEM from topographic data
climatic variability in the region and assess present-day
glacier fluctuations in the area. Two 1:50,000 topographic maps, constructed from
1955 aerial photography by Instituto Geogrfico Nacio-
3. Methods nal (IGN) of Peru were needed to cover the study site.
The maps used Provisional South American Datum of
3.1. Field data collection 1956 for Peru, and elevations were referenced to mean
sea level. The maps were scanned and georeferenced
In June and August 2003, GPS points were obtained based on UTM grids with a positional accuracy
using a handheld Garmin Etrex GPS unit in navigation (calculated as root mean square error in the X and Y
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 113

Fig. 2. Location map of the study area. The ASTER Level 1A image from July 2001 is draped over a shaded relief map of the topographic DEM. Also
shown are GPS transects surveyed in the field (red dots).

coordinates) of 4 m. Contour lines with 25 m spacing and McDonnell, 1998). TOPOGRID interpolates direct-
were digitized on screen, and attributed the corre- ly from the contour lines by determining areas of steep-
sponding elevation values read from the topographic est slope and generating terrain morphology. Ancillary
map. Additional GIS layers digitized from the topo- hydrologic data (streams and lakes) are used to define
graphic maps included lakes, streams, spot heights and drainage based on the ANUDEM algorithm for hydro-
the 1955 snowline. logic modeling described by Hutchinson (1988).
We examined common interpolation routines to Triangulated Irregular Network (TIN) data structures
create continuous data from the digitized contours. are terrain models represented by continuous triangular
These included: Inverse Distance Weighted (IDW), facets that store elevation at irregularly spaced nodes
Splines (TOPOGRID) and Triangulated Irregular Net- (Burrough and McDonnell, 1998).
work (TINs). The IDW method estimates the Z value of
an unknown point based on a distance-weighted average 3.3. SRTM and ASTER datasets
of elevation points within a neighborhood (Burrough
and McDonnell, 1998). Spline techniques use a piece- The Shuttle Radar Topography Mission (SRTM)
wise function to fit a curve through all the data points. acquired data in February 2000, from which digital
The TOPOGRID algorithm is a more sophisticated elevation models are created (Rabus et al., 2003). Pre-
spline technique (thin plate spline) that fits a smoothing liminary elevation datasets with 90 m resolution
surface through the data points to minimize artifacts (SRTM-3) were recently released for South America.
(excessively high or low spurious values) (Burrough An elevation dataset (1-degree latitude by 1-degree
114 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

longitude) was obtained for the study area, referenced to to examine the representation of topography in each
UTM projection and resampled to 30 m resolution. DEM.
Elevations are in meters, referenced to the EGS84 Difference maps were constructed by subtracting the
EGM96 geoid (USGS, 2003). DEM from topographic data from both the ASTER and
Two ASTER scenes acquired by along-track stereo SRTM-derived DEM's on a cell-by-cell basis. We
channel (3), with nadir (3n) and aft-viewing (3b) orien- examined correlations between vertical differences and
tations (Kb, 2002) were obtained from the Land topographic characteristics (elevation, slope and aspect).
Processes DAAC at EROS Data Center: one Level 1 B Errors on the non-glaciated areas (bias) were quanti-
scene from October 2000 and one Level 1A scene from fied by performing trend surface analyses on the dif-
July 2001. The cloud-free 2001 ASTER scene, shown in ference maps. After removing the bias, we examined the
Fig. 2, was used to extract a DEM using automated remaining elevation differences on glaciated areas to
stereo auto-correlation procedures at the USGS EROS distinguish a glacier signal using histograms, summary
Data Center. Ground control points (GCP's) were statistics and color maps of the height differences.
required to obtain an absolute ASTER DEM (where
locations are fitted to UTM coordinate system and 4. Results and discussion
elevations referenced to mean sea level) (Hirano et al.,
2003). Eight GCP's were digitized from the topographic 4.1. Topographic interpolation results
maps at river crossings, spot elevations and road inter-
sections and identified on the 3n and 3b bands of the An examination of the RMSEz values for DEM's
Aster image, following the protocol of Hirano et al. derived from topographic data (Table 1) shows that no
(2003) and Khalsa et al. (2004). The resulting ASTER- interpolation method performed perfectly. The vertical
derived DEM had 30 m post spacing. accuracy of the DEM created with the TOPOGRID
Various ASTER scenes were evaluated to find the algorithm was 14.7 m based on 61 spot elevations. The
one that provided the best glacial extent, based on image other interpolation methods yielded RMSEz values that
contrast and minimal snow coverage. We used the ranged from ~21 to 24 m, which is 3040% greater
ASTER L1B scene obtained in October 2000 (end of the than the TOPOGRID algorithm (Table 1). A one-way
dry season in the Andes) to delimitate the glacier out- analysis of variance (ANOVA) test showed that there
line. An unsupervised ISODATA clustering classifica- was a significant difference in the DEM's created with
tion (Aniya et al., 1996; Paul, 2001) was performed various interpolation methods at the 0.1 significance
using ASTER VNIR channels (1, 2 and 3) to delimitate level ( p-value = 0.07). The RMSEz of 14.7 m using the
the ice extent. The resulting raster image was converted TOPOGRID algorithm is only slightly bigger than half
to polygon coverage and was visually checked to ensure of the contour interval (25 m), which is considered an
correspondence with glaciated areas on the color com- acceptable vertical accuracy for DEM's derived from
posite image. topographic maps (Cogley and Jung-Rothenhausler,
2004).
3.4. DEM validation and comparison All DEM's constructed from contours lines display
terracing effects due to denser sampling along the
We focused on evaluating errors in the vertical co- contour lines, because points closer to the contour lines
ordinate (Z), estimated as root mean square errors are interpolated using the same elevation values
(RMSEz). The Z coordinate is the only unconstrained (Burrough and McDonnell, 1998). The terracing effect
value, since X and Y coordinates were used to locate is most visible on flat surfaces where contours are spread
corresponding grid cells in all DEM's. Moreover, ele- apart, and is most severe when using local interpolator
vation is the coordinate of interest in glaciological
applications because changes in surface elevation over Table 1
time can be an indicator of mass balance changes Evaluation of different interpolation methods used to construct DEM's
(Etzelmller, 2000). The RMSEz of the various inter- from the topographic maps
polated methods was calculated with respect to evenly Interpolation RMSEz spot elevations Terracing Other artifacts
distributed spot elevations digitized from topographic method (m)
maps. The RMSEz of the SRTM and ASTER DEM's TOPOGRID 14.7 Light Cones
was calculated with respect to GPS points from non- IDW 24.2 Severe No
glaciated areas. Visualization techniques (shaded relief TIN 22.0 No Triangulation
SPLINE 21.1 Moderate No
maps, elevation contours and slope maps) were used
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 115

Fig. 3. The effect of interpolation methods on representation of terrain topography at a subsection of the study area. a) original contour lines, with
25 m interval; b) shaded relief map of the DEM created with the IDW method; c) shaded relief map of the TIN data structure.

methods such as IDW (Fig. 3b). Etzelmller and large differences in glacier surface representation by
Bjrnsson (2000) used the IDW interpolation method DEM's as a function of interpolation algorithm used.
to create a continuous surface of glacier thickness. This Based on minimizing both RMSEz and artifacts
systematic terracing artifact was also reported in other (terracing and triangulation), we chose the DEM created
glacier studies that used tension splines for surface with the TOPOGRID algorithm (denoted as TOPO
representation (e.g. Mennis and Fountain, 2001). DEM) as the 1955 reference elevation dataset.
Terracing is known to affect subsequent calculations of
topographic characteristics (slope, aspect and profile 4.2. Accuracy assessment for the SRTM and ASTER
curvature) (Wilson and Gallant, 2000) which are of derived DEM's
interest glaciological applications. For our study area,
with gentle sloping terrain, the DEM created with the SRTM elevations and ASTER elevations were
TOPOGRID algorithm yielded the smoothest surface. checked against 64 GPS points from non-glaciated
Minimal terracing is detected in this DEM, and appears terrain. We focused on non-glaciated terrain to validate
as spikes on the histogram of elevation values (Fig. 4a). the DEM's because elevation changes might have oc-
Other glaciological studies (e.g. Gratton et al., 1990) curred on the glacier between 1955 and 2000/2001. We
preferred TIN's because of their advantage of capturing present the elevation differences of the DEM's as
complex terrain variations, accurately representing RMSEz relative to the GPS points, and not the absolute
ridges and streams, and reducing data redundancy on vertical accuracy. The RMSEz of the SRTM DEM
flat terrain (Burrough and McDonnell, 1998). For relative to the GPS points was 23.4 m. Since the vertical
Coropuna, the TIN structure introduced noticeable tri- accuracy of GPS points is 10 m (cf. Section 3.1), this
angular discretization on the gentle sloping lava flows gives an absolute vertical accuracy of 23.4 m 10 m.
and smooth glacier surface (Fig. 3c), which we con- The specified SRTM accuracy standard is 16 m for
sidered unacceptable. Our results show that there are global coverage (Rabus et al., 2003). The frequency
116 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

Fig. 4. Histograms of elevation values for the three DEM's analyzed. a) DEM created with TOPOGRID algorithm (TOPO DEM), b) SRTM DEM and
c) ASTER DEM. Spikes on the elevation histograms represent elevation values along the contour lines used in the interpolation and point to the
terracing effect.

histogram of SRTM-derived elevations (Fig. 4b) in- absolute vertical accuracy of 61.2 m 10 m, which is
dicates a normal distribution, with a few anomalously bigger than the specified accuracy of 750 m for abso-
high values (spikes). Water bodies are not well defined lute ASTER DEM's (Lang and Welch, 1999). ASTER-
and appear noisy or rough due to low radar back- derived elevations are consistently higher than the GPS
scatter (USGS, 2003). Height differences between points, and the magnitude of the vertical differences
SRTM elevations and GPS elevations tend to be between ASTER and GPS increases with elevation
randomly distributed (Fig. 5a), with SRTM elevations along the GPS transects (Fig. 5b). Such vertical errors
being both lower and higher than the GPS elevations. were reported in other studies. Kb (2002) found an
The comparison of elevations from the ASTER DEM overall accuracy of 60 m RMSEz when an absolute
with GPS points shows both a large RMSEz and a ASTER DEM was compared to a reference DEM on
vertical bias. The RMSEz of the ASTER DEM relative complex mountain terrain. However, better accuracy
to the GPS points is 61.2 m. This corresponds to an ( 18 m RMSEz) was found at a section of moderate
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 117

Fig. 5. Plots of height differences between the DEM's from satellite data and GPS elevation along GPS transects on non-glaciated (black diamonds)
and glaciated areas (grey triangles). a) SRTM elevations minus GPS elevations; b) ASTER elevations minus GPS elevations.

topography in the same study (Kb, 2002), suggesting based on 15 control points identified at lakes, stream
that errors in the ASTER DEM's tend to increase in crossings, and noticeable terrain features such as ridges.
rugged mountain terrain. Additional validations of the SRTM and ASTER-
Some errors come from some noise due to banding derived DEM's were performed by comparing their
in the ASTER L1A scene, visible in the DEM as spikes elevations with the reference 1955 DEM from topographic
on the elevation histogram (Fig. 4c). Comparison of data on a cell-by-cell basis. Subtracting the reference
contours derived from the ASTER DEM with contours DEM from the SRTM DEM yielded a mean difference
from topographic data revealed positional offsets as of 1.8 m and a standard deviation of 15.7 m. The range of
much as 300 m in X and 200 m in Y. These offsets were vertical differences was 113 m/+121 m, with the largest
not consistent throughout the study area, pointing to- differences occurring on non-glaciates areas, at valley
wards a distortion in the ASTER DEM in the X and Y bottoms and sharp moraine ridges, as well as on a few flat
coordinates. Such horizontal offsets have been observed areas where interpolation from contour lines produced
at other study areas with high relief (Dwyer, LP DAAC erroneous values (either spikes or sinks).
Project Scientist, USGS EROS Data Center, personal Subtracting the reference DEM from the ASTER
communication). To correct the offsets, the ASTER DEM yielded a mean difference of 80.5 m and
DEM was fitted to the georeferenced ASTER L1A standard deviation of 28.1 m. The range of differences
scene using a second order polynomial transformation was 86 m/ +500 m. The large positive differences
118 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

Fig. 6. Color maps of height differences (in meters) between the SRTM elevations and topographic elevations shown as 3D perspectives. a) 2000
SRTM DEM (before trend removal) minus 1955 TOPO DEM. The NE-SW spatial trend in elevation differences is visible. b) 2000 SRTM DEM (after
trend removal) minus 1955 TOPO DEM. Areas depicted in white represent NODATA in the SRTM DEM. Also shown is the glacier extent obtained
by classification of the October 2000 L1B ASTER scene (black line).

of +500 m come from spikes of erroneous elevation band 3b (Kb et al., 2002b) or low contrast in the
values that occur on glaciated summits and sharp ridges. ASTER scenes over snow and ice, causing failure in the
Such spikes of up to 500 m in ASTER DEM's were noted image-matching process (Toutin, 2002). While distor-
at other areas on sharp peaks (Kb et al., 2002b). tions in the ASTER DEM may be due to lack of adequate
Elevation waves with about 200300 m amplitude on ground control points, in other studies in mountainous
low contrast glacier areas were also reported by the terrain, introducing more GCP's did not significantly
USGS EROS Data Center (Wessels, U.S. Geological remove this effect (Kb et al., 2002b).
Survey, Alaska Science Center, personal communica- Height differences between the SRTM/ASTER
tion). These elevation errors are due to either steep DEM's and the reference DEM shown on color maps
northern slopes which are missed by the back-looking in Figs. 6a and 7a, indicate a systematic bias in both
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 119

Fig. 7. Color maps of elevation differences (in meters) between ASTER DEM and the TOPO DEM's, in meters. a) 2001 ASTER DEM (before trend
removal) minus 1955 TOPO DEM, with the NW-SE spatial trend visible; b) 2001 ASTER DEM (after trend removal) minus the 1955 TOPO DEM.
Black lines represent the 2000 glacier extent derived from the ASTER scene.

difference maps, with positive residuals on non-glaciated gression on X and Y coordinates and is an inclined
areas in the north and negative residuals in lower valleys surface of the form:
in the south. Similar biases in residuals were noted in
other comparisons of ASTER DEM's with IGN topo- f f x; yg a0 a1 X a2 Y ;
graphic data (e.g. Vignon et al., 2003) as well as
comparison of ASTER DEM's with photogrammetric- where a0 is the intercept and a1 and a2 are the slopes
derived DEM's (Kb et al., 2002b). We modeled the (Burrough and McDonnell, 1998).
variation of residuals over the non-glaciated area for the The bias of elevation differences between the SRTM
two datasets by fitting various polynomial surfaces DEM and the reference DEM is a tilted surface, oriented
through the residuals. The best fit in terms of R-square towards the NNE (5.42), which dips at a rate of 1.9 m
was obtained by a first order polynomial, suggesting that vertical per 1 km northing, and has a range of 21 m to
the magnitude of the residuals increases linearly with 20 m across the DEM. For the ASTER minus TOPO
location. The polynomial is derived by multiple re- DEM, the bias is a tilted surface oriented towards the
120 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

Fig. 8. Frequency histograms of elevation differences between the DEM's on non-glaciated (grey lines) vs. glaciated areas (black lines), after the trend
removal. a) SRTM DEM minus TOPO DEM; b) ASTER DEM minus TOPO DEM.

NNW (349.30), which dips at a rate of 2 m vertical per or low values) and they do not affect subsequent
1 km northing, and has a range of 47 m to 96 m across analysis of the glaciated areas.
the DEM. After the trend removal, we examined the effect of
Once the trend was removed, the mean statistics for slope and aspect on the vertical differences between the
the elevation differences on non-glaciated areas yielded: DEM's. Correlations with slope yielded a coefficient
(Pearson's r) of 0.54 for SRTM minus reference DEM
SRTM DEM minus TOPO DEM (Fig.8a): mean = 0, and 0.69 for ASTER minus the reference DEM. The
std. deviation = 9.5 m plots of vertical differences with respect to slope (Fig.9)
ASTER DEM minus TOPO DEM's (Fig. 8b): show that elevation errors in the SRTM and ASTER
mean = 0, std. deviation = 20.6 m. DEM's tend to increase with slope. On slopes less than
45, there is almost no difference in SRTM-derived
The histograms of elevation differences on non- elevations and the topographic DEM, but corresponding
glaciated areas are close to normally distributed ASTER elevations are consistently higher than the ref-
(Fig. 8 ab). Large standard deviations on non- erence DEM. For the SRTM DEM, elevation errors of
glaciated areas point to artifacts in the DEM's (high up to 25/+ 50 m occur on 6065 slopes. For the
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 121

Fig. 9. Vertical differences between DEMs from satellite data and

DEM from topographic data as a function of terrain slope: SRTM
minus TOPO DEM (black squares) and ASTER minus TOPO DEM
(grey dots). Largest vertical differences between SRTM/ASTER and
the topographic DEM occur on steepest slopes.

ASTER DEM, elevation errors greater than 100 m occur

on steep slopes (6077) and correspond to the spike
artifacts in the ASTER DEM. These results are con-
sistent with trends noted in ASTER-derived DEM's. For
instance, Kb (2002) and Kb et al. (2002b) found
large vertical differences in the ASTER DEM's com-
pared to reference DEM's from topographic data at
steep slopes.
Vertical differences between the SRTM DEM and the
reference DEM (Fig. 10a) do not depend on slope
aspect. The mean elevation differences between SRTM
and reference DEM range from 5 m on S-facing Fig. 10. Radar charts of vertical differences between the DEM's as a
slopes to 1 m on NNE-facing slopes. The ASTER function of aspect. a) SRTM DEM minus TOPO DEM; b) ASTER
DEM displays bigger mean vertical differences DEM minus TOPO DEM.
ranging from 3 m on W-facing slopes to + 19 m on
E-facing aspects (Fig. 10b). We expected bigger eleva- elevations acquired on the glaciated area yielded a
tion errors on N-facing slopes, which are normally RMSEz of 27 m, which is 4 m bigger than on non-
missed by the back-looking band 3b (Kb et al., glaciated areas. SRTM elevations are both higher and
2002b). For instance, decreased vertical accuracy of lower than the GPS points on the glacier (Fig. 5a). The
ASTER DEM's on northern slopes was reported by RMSEz of ASTER elevations with respect to GPS
Kb (2002) and Kb et al. (2002b). Our results show points on the glacier was 98 m, and residuals increase
that mean errors tend to occur on aspects between 0 with altitude (Fig 5b). This large vertical bias
and 180, not only on N-facing slopes. This suggests
that the large vertical errors in the ASTER DEM
cannot be entirely explained by the back-looking Table 2
Statistics summary of map differences for the glaciated areas vs. non-
ASTER channel. glaciated areas after trend removal
Statistics SRTM TOPO DEM ASTER TOPO DEM
4.3. Glacier signal from the SRTM DEM (m) (m)
Glaciated Non-glaciated Glaciated Non-glaciated
We checked elevations from the SRTM and ASTER
DEM's against GPS points on glaciated areas. The Mean 5 0.0 28.5 0.0
Std. deviation 15.8 9.5 26 20.5
difference between SRTM elevations and 56 GPS
122 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

implies that the ASTER DEM, created with a limited glaciated area (Fig. 6b) show ablation at the toes of the
number of GCP's, is not suitable for glaciological glaciers ( 25 m to 75 m surface lowering) along with
interpretation. an apparent thickening at the summits (2550 m).
To quantify the glacier signal from the SRTM DEM, Similar comparisons of ASTER data to topographic data
we examined the mean elevation differences (SRTM in Cordillera Blanca (Peru) revealed a loss of altitude of
minus topographic DEM) on glaciated areas after re- as much as 23 m at the glacier toes (Vignon et al.,
moving the NNE-SSW spatial trend. Once the trend was 2003). Ablation in the lower parts of the glaciers (via ice
removed, the elevation differences on the glaciated area melting and sublimation) was also observed from field
were negatively skewed (Fig. 8a), with a mean of 5 m measurements in other tropical glaciers (Kaser et al.,
and a standard deviation of 15.8 m (Table 2). We con- 1990; Kaser, 1999).
sider the remaining mean difference of 5 m 15.8 m as Thickening in the accumulation zone of the glaciers
a signal of glacier thinning (95% confidence interval). is a less common trend and was observed in some
Average height differences between the SRTM and mountain glaciers around the world during the 1961
topographic DEM on the glaciated area increase with 1997 time period (Dyurgerov and Meier, 2000). How-
altitude, with a correlation coefficient (Pearson's r) of ever, in the climatic context of Coropuna, an average
0.62 (Fig. 11a). Cell-by-cell comparison of elevations thickening of 2550 m firn in 50 yr, or 0.5 m/yr would
from SRTM data with topographic DEM within the represent a mean increase in precipitation of 250

Fig. 11. Vertical differences between the DEM's as a function of altitude on the glaciated area averaged in a 150 150 m neighborhood (grey
symbols). Also shown are mean trends calculated as average difference in 2 m elevation bands (black dots). a) SRTM DEM minus TOPO DEM;
b) ASTER DEM minus TOPO DEM.
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 123

500 mm water equivalent. At the col, our results agree few decades. Ames et al. (1989) reported a total
with field data from the ice core drilling of June 2003, glaciated area of 146.7 km2 based on 1962 aerial
which also point to an accumulation of 0.51 m firn / photography. The total glaciated area in the Ampato
yr in the col (Ginot, Laboratoire de Glaciologie et range was estimated to be 105 km2 based on Landsat
Geophysique de l'Environnement, Grenoble, personal TM imagery (Morales-Arnao, 1999). This corre-
communication). However, at the summits, the in- sponds to a retreat of 27% in Cordillera Ampato
crease of 0.51 m firn/yr from the DEM comparison is from 1962 to the end of the 20th century. Glacial
24 times bigger than ice core results (.26 m firn/yr) retreat in Peru has also been observed in other areas,
(Ginot, Laboratoire de Glaciologie et Geophysique especially in Cordillera Blanca (Kaser et al., 1990;
de l'Environnement, personal communication). In the Hasternath and Ames, 1995; Georges, 2004).
upper part of the glaciated area the noise may be too
high to be able to infer a positive change of 2550 m 5. Conclusions and further applications
in altitude.
The elevation differences between ASTER eleva- Using DEM's derived from topographic and satellite
tions and reference DEM on glaciated area are data at different steps in time holds potential for glacier
positively skewed (Fig. 8b), with a mean of 28.5 m analysis. DEM's constructed from old topographic data
and a standard deviation of 26 m. Comparison of GPS still constitute a valid elevation dataset for comparison
points with corresponding ASTER elevations on with more recent DEM's for glaciology purposes. Here
glaciated areas (Fig. 5b) shows that the ASTER we created a DEM from 1:50,000 topographic data for
DEM is too high on glaciated terrain, with a RMSEz Nevado Coropuna and tested different interpolation
error of 98.3 m with respect to GPS points. The techniques. Based on RMSEz and visual analysis, the
ASTER DEM is also systematically higher than the TOPOGRID algorithm was found to be superior to the
topographic DEM (Fig. 11b). However, an examina- other techniques examined, with the smallest RMSEz
tion of cell-by-cell differences between the ASTER error and least interpolation artifacts.
DEM and the reference DEM (Fig. 7b) shows Error analyses were performed on all DEM's to
negative residuals ( 50 to 25 m) in the ablation characterize the bias present in the various DEM's. We
areas of the southern glaciers. While surface lowering removed the spatial bias to distinguish a glacier signal.
at the glacier toes is consistent with results from the We found that the SRTM dataset with a RMSEz of
SRTM DEM, we could not quantify the glacier signal 23.4 m 10 m was suitable for glaciological applica-
from the ASTER DEM due to the altitudinal bias and tions after some calibration. However, in areas of
the large elevation spikes on the glacier surface, rugged terrain, the SRTM resolution (90 m) was not
which are affecting the mean statistics. We suspect that sufficient to accurately represent the topography.
the large elevation spikes were due to lack of contrast Comparison of the 2000 SRTM DEM with the DEM
over the glacier in the ASTER image. The VNIR and from 1955 topographic data points to an average
SWIR gain level settings, which are based on sun thinning of 5 m on the glacier surface, with a
angle, can be optimized for snow targets to provide significant lowering of the glacier surface at the glacier
maximum contrast over ice and snow (Raup et al., toes and an apparent accumulation on the summits. We
2000). However, the ASTER images currently avail- attribute the vertical differences of more than 25 m at the
able for Coropuna were not acquired with these settings summits to possible errors in either the SRTM data or
and therefore provided little contrasts over the the topographic data at higher elevations and steeper
glacierized surface. slopes. While lowering of glacier surface at the toes was
visible in the ASTER DEM, large elevation errors and
4.4. Changes in glacier extent and volume altitudinal bias did not allow quantifying a glacier signal
from the ASTER data.
For 1962, Ames et al. (1989) reported a glaciated In conclusion, the analysis of multi-temporal DEM's
area of 82.6 km 2 on Nevado Coropuna based on to quantify glacier changes is extremely sensitive to the
planimetric analysis of 1962 aerial photography. quality and spatial resolution of the DEM's. For our
Based on the ASTER L1B scene from October study of glacier change using DEM's on Nevado
2000, we obtained a glacier area of 60.8 km 2 , Coropuna, we found that several steps were necessary:
which represents a loss of 26% in glacier area from referencing all elevation data to the same vertical datum;
1962 to 2000. Our results are consistent with glacier evaluation of DEM differences in non-glaciated areas;
retreat observed in Cordillera Ampato during the last testing the DEM's against field GPS survey points;
124 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125

visualization techniques such as shaded relief, slope Dyurgerov, M.B., Meier, M.F., 2000. Twentieth century climate
angle and comparison of contours; removing the biases change: evidence from small glaciers. PNAS 97 (4), 14061411.
Etzelmller, B., 2000. On the quantification of surface changes using
in the elevation datasets. grid-based digital elevation models (DEMs). Transactions in GIS 4
Future steps to minimize large error differences (2), 129143.
occurring in DEM's derived from satellite data include Etzelmller, B., Bjrnsson, H., 2000. Map analysis techniques for
filtering and smoothing of the DEM's (Toutin, 2001, glaciological applications. International Journal of Geographical
Information Science 14 (6), 567581.
2002; Hirano et al., 2003). These techniques may
Etzelmller, B., Vatne, G., Odegard, R.S., Sollid, J.L., 1993. Mass-
help to better distinguish and quantify glacier surface balance and changes of surface slope, crevasse and flow pattern of
changes. Erikbreen, Northern Spitsbergen an application of a geographical
information-system (GIS). Polar Research 12 (2), 131146.
Acknowledgments Falorni, G., Istanbulluoglu, E., Bras, R.L., 2003. A method for
improving SRTM DEMs in high-relief terrain. EOS Transactions
AGU 84 (46).
The research was supported by the NSF/IGERT Georges, C., 2004. 20th century glacier fluctuations in the Cordillera
sponsored graduate training program entitled Carbon, Blanca (Per). Arctic, Antarctic, and Alpine Research 36 (1),
Climate and Society Initiative (CCSI) at University of 100107.
Colorado. The French program ACI Observation de la Gratton, D.J., Howarth, P.J., Marceau, D.J., 1990. Combining DEM
Terre supported the use of remote sensing for parameters with Landsat MSS and TM imagery in a GIS for
mountain glacier characterization. IEEE Transactions on Geosci-
glaciological studies. The Niwot Ridge LTER project ence and Remote Sensing 28 (4), 766769.
provided logistical support (NSF DEB 9810218). We Hasternath, S., Ames, A., 1995. Recession of Yanamaery glacier in
thank the GLIMS team for facilitating free access to the Cordillera Blanca, Per during the 20th century. Journal of
ASTER data through the NASA ESE Pathfinder project; Glaciology 41 (137), 191196.
USGS EROS data center for provinding the ASTER and Hirano, A., Welch, R., Lang, H., 2003. Mapping from ASTER stereo
image data: DEM validation and accuracy assessment. ISPRS
SRTM data; Bernard Pouyaud and Bernard Francou Journal of Photogrammetry and Remote Sensing 57 (56),
(UR GREAT ICE, IRD France) for supporting the 356370.
participation of the main author in the ice core drilling Hutchinson, M.F., 1988. Calculation of hydrologically sound digital
expedition on Coropuna; the Glaciology Unit at elevation models. Proceedings of Third International Symposium
INRENA, Huaraz, for providing helpful assistance; the on Spatial Data Handling, Sydney.
Kb, A., 2002. Monitoring high-mountain terrain deformation from
UNAVCO facility in Boulder for assisting with GPS repeated air- and spaceborne optical data: examples using digital
corrections. aerial imagery and ASTER data. ISPRS Journal of Photogram-
metry and Remote Sensing 57 (12), 3952.
References Kb, A., Paul, F., Maisch, M., Hoelzle, M., Haeberli, W., 2002a. The
new remote sensing-derived Swiss glacier inventory: II. First
Ackerman, T., Erickson, T., Williams, M.W., 2001. Combining GIS results. Annals of Glaciology 34, 362366.
and GPS to improve our understanding of the spatial distribution of Kb, A., Huggel, C., Paul, F., Wessels, R., Raup, B., Kieffer, H.,
snow water equivalence (SWE). Proceedings of ESRI USERS Kargel, J., 2002b. Glacier monitoring from ASTER imagery:
Conference, San Diego, CA. accuracy and applications. Proceedings of EARSeL LIS-SIG
Ames, A., Muoz, G., Verstegui, J., Zamora, M., Zapata, M., 1989. Workshop, Berne, March 1113, 2002.
Glacier Inventory of Peru. Part I. Hidrandina S.A. Unit of Kaser, G., 1999. A review of modern fluctuations of tropical glaciers.
Glaciology and Hydrology. Huaraz-Per. Global and Planetary Change 22, 93103.
Aniya, M., Sato, H., Naruse, R., Skvarca, P.C.G., 1996. The use of Kaser, G., Ames, A., Zapata, M., 1990. Glacier fluctuations and
satellite and airborne imagery to inventory outlet glaciers of the climate in the Cordillera Blanca, Per. Annals of Glaciology 14,
Southern Patagonia Icefield, South America. Photogrammetric 136140.
Engineering and Remote Sensing 62, 13611369. Khalsa, S.J.S., Dyurgerov, M., Khomova, T.B.R., Barry, R.G., 2004.
Berthier, E., Arnaud, Y., Baratoux, D., Vincent, C., Rmy, F., 2004. Space-based mapping of glacier changes using ASTER and GIS
Recent rapid thinning of the Mer de Glace glacier derived from tools. IEEE Transactions on Geoscience and Remote Sensing, 42
satellite optical images. Geophysical Research Letters 31, L17401. (190), 21772183.
doi:10.1029/2004GL020706. Klein, A.G., Isacks, B.L., 1996. Glaciers: tracking change in the
Burrough, P.A., McDonnell, R.A., 1998. Principles of Geographic central Andes Mountains. GIS World 9, 4852.
Information Systems. Oxford University Press, Oxford. Lang, H.R., Welch, R., 1999. Algorithm Theoretical Basis Document
Cogley, J.G., Jung-Rothenhausler, F., 2004. Uncertainty in digital for ASTER Digital Elevation Models, Version 3.0. Jet Propulsion
elevation models of Axel Heiber Island, Arctic Canada. Arctic, Laboratory, Pasadena, CA. 69 pp.
Antarctic and Alpine Research 36 (2), 249260. Mennis, J.L., Fountain, A.G., 2001. A spatio-temporal GIS database
Duncan, C.C., Klein, A.J., Masek, J.G., Isacks, B.L., 1998. for monitoring alpine glacier change. Photogrammetric Engineer-
Comparison of late Pleistocene and modern glacier extents in ing and Remote Sensing 67 (8), 967975.
central Nepal based on digital elevation data and satellite imagery. Morales-Arnao, B., 1999. Glaciers of Peru. In: Williams Jr., R.S.,
Quaternary Research 49, 241254. Ferrigno, J.G. (Eds.), Satellite image atlas of glaciers of the world:
 A.E. Racoviteanu et al. / Global and Planetary Change 59 (2007) 110125 125

South America. United States Geological Survey Professional Paper Sidjak, R., 1999. Glacier mapping of the Illecillewaet icefield, British
1386-I. United States Government Printing Office, Washington. Columbia, Canada, using Landsat TM and digital elevation data.
Paul, F., 2001. Evaluation of different methods for glacier mapping International Journal of Remote Sensing 20 (2), 273284.
using Landsat TM. Proceedings of EARSeL Workshop on Remote Toutin, T., 2001. Elevation modeling from satellite visible and infrared
Sensing of Land Ice and Snow, Dresden, June 16 -17, 2000, vol. 1, (VIR) data. International Journal of Remote Sensing 22 (6),
pp. 239245. 10971125.
Paul, F., Kb, A., Maisch, M., Kellenberger, T., Haeberli, W., 2002. Toutin, T., 2002. Three-dimensional mapping with ASTER stereo data
The new remote- sensing-derived Swiss glacier inventory: I. in rugged topography. IEEE Transactions on Geoscience and
methods. Annals of Glaciology 34, 355361. Remote Sensing 40 (10), 22412247.
Rabus, B., Eineder, M., Roth, A., Bamler, R., 2003. The shuttle radar USGS, 2003. SRTM documentation. Available at: ftp://edcsgs9.cr.usgs.
topography missiona new class of digital elevation models gov/pub/data/srtm/Documentation/ (accessed May 16, 2003).
acquired by spaceborne radar. ISPRS Journal of Photogrammetry Vignon, F., Arnaud, Y., Kaser, G., 2003. Quantification of glacier
and Remote Sensing 57 (4), 241262. volume change using topographic and ASTER DEMs: a case study
Rapp, R.H., 1996. Use of potential coefficient models for geoid in the Cordillera Blanca. Proceedings of Geoscience and Remote
undulation determination using a spherical harmonic representa- Sensing Symposium, IEEE International 4, 26052607.
tion of the height anomaly/geoid undulation difference. Journal of Welch, R., Jordan, T., Lang, H., Murakami, H., 1998. ASTER as a
Geodesy 71 (5), 282289. source for topographic data in the late 1990's. IEEE Transactions
Raup, B.H., Kieffer, H.H., Hare, T.M., Kargel, J.S., 2000. Generation on Geoscience and Remote Sensing 36 (4), 12821289.
of data acquisition requests for the ASTER satellite instrument for Wilson, J.P., Gallant, J.C., 2000. Digital terrain analysis. In: Wilson, J.P.,
monitoring a globally distributed target: glaciers. IEEE Transac- Gallant, J.C. (Eds.), Terrain analysis - principles and applications.
tions on Geoscience and Remote Sensing 38, 11051112. Wiley, New York.
Rentsch, H., Welsch, W., Heipke, C., Miller, M., 1990. Digital terrain Wood, J.D., Fisher, P.F., 1993. Assessing interpolation accuracy in
models as a tool for glacier studies. Journal of Glaciology 36 (124), elevation models. IEEE Computer Graphics and Applications 13
273278. (2), 4856.
Rivera, A., Casassa, G., 1999. Volume changes on Pio XI glacier,
Patagonia: 19751995. Global and Planetary Change 22 (14),
233244.