Glacier mapping from multispectral satellite data

Two decades of glacier mapping from multispectral satellite data has resulted in a wealth of methods and data. Although only cloud-free imagery from the end of the ablation period can be used for efficient glacier mapping, appropriate scenes from nearly all the world’s glaciers are available today, due to the long time period covered (starting with TM in 1982). Given that many glacierized regions are not covered by appropriate ASTER/ETM scenes, the archived dataset should also be considered and analyzed systematically for the GLIMS project.

Apart from manual glacier delineation by on-screen cursor tracking, most methods utilize the low reflectance of ice and snow in the middle infrared part of the spectrum for glacier classification (e.g. Reference Bayr, Hall and KovalickBayr and others, 1994; Reference PaulPaul, 2002b; Reference KääbKääb and others, 2003a). This spectral region is covered by several space-borne sensors (e.g. ASTER, IRS-1C/ D, Landsat TM/ETM, SPOT 4/5), and the methods are thus widely applicable. Latest reviews and comparisons of the applied methods were given by Reference Sidjak and WheateSidjak and Wheate (1999), Reference Gao and LiuGao and Liu (2001), Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002) and Reference BishopBishop and others (2004). In particular, thresholded ratio images from various TM band combinations have proven to be a very accurate, fast and robust method, at least for clean glacier ice (Reference AlbertAlbert, 2002; Reference Paul, Huggel, Kääb and KellenbergerPaul and others, 2003). Debris-covered glacier ice is not detectable by multispectral classification alone, but the combination with DEM information (slope) and neighbourhood analysis has shown promising results (Reference Paul, Huggel and KääbPaul and others, 2004 and references therein).

From a remote-sensing perspective, the special challenges at our test site are deep shadows, lakes of varying turbidity adjacent to glacier ice, and debris-covered glaciers. In this section we focus on mapping of glaciers in cast shadow using the ETM scene. In view of the accurate glacier maps obtained from simple band ratios and the absence of an appropriate DEM, we have compared three ratios for glacier mapping: ETM3/ETM5, ETM4/ETM5 and (ETM2-ETM5)/(ETM2 + ETM5), where the latter is also known as the normalized-difference snow index (NDSI). We have also investigated the effect of subtracting path radiance previous to the classification and applied a median filter (noise reduction) to all glacier maps shown in the following comparisons.

In Figure 2a our test region is shown in a contrast-stretched ETM band 3 image (cf. Fig. 1 for location). The most critical regions (shadow, debris, snowfields) are indicated by arrows. The comparisons presented in Figure 2b-d depict the glacier maps from two methods at a time: Figure 2b compares ETM3/ETM5 with ETM4/ETM5, Figure 2c the NDSI with and without subtraction of path radiance in ETM3, and Figure 2d compares ETM3/ETM5 with the NDSI with subtracted path radiance in ETM3. The comparison shows that ETM4/ETM5 is not capable of mapping snow or ice in cast shadow, while the NDSI without path radiance subtraction maps all regions in cast shadow (including bare rock). The visual inspection reveals two ‘winners’ which are compared in Figure 2d, revealing only small differences. The NDSI performs somewhat better in regions with ice in cast shadow, while ETM3/ETM5 is more accurate in regions with snow in cast shadow. Lastly, we use ETM3/ETM5 for glacier mapping, as final manual editing (e.g. for debris cover) has to be applied for many glaciers anyway. Glacier mapping in cast shadow remains the most critical point. Thus, we suggest that efficient glacier mapping benefits most from a careful selection of the threshold parameter in such regions. This will markedly reduce time-consuming manual corrections in these regions.

Fig. 2. (a) Test region for glacier mapping in a contrast-stretched ETM band 3 image with critical regions indicated by numbered arrows (image size is 10.5 by 9km): 1a. cast shadow on ice; 1b. cast shadow on snow; 1c. cast shadow on rock; 2. debris-covered ice; 3. snowpatches. (b–d) Glacier maps comparing two methods at a time. Light-grey areas were identified by both methods as being a glacier, black areas only by the first method, and dark-grey areas only by the second method. The compared methods are: (b) ETM3/ETM5 vs ETM4/ ETM5 from digital number; (c) NDSI with and without correction for path radiance; and (d) NDSI (with path radiance correction) vs ETM3/ETM5.

GIS-based glacier basin delineation

A GIS acts at the interface between the satellite-derived glacier map and the data obtained for individual glaciers. In particular, the digitizing of glacier basins (ice divides) and the automatic data extraction for each glacier, including 3-D parameters obtained from a DEM, are greatly facilitated (Reference Kääb, Paul, Maisch, Hoelzle and HaeberliKääb and others, 2002; Reference PaulPaul, 2002b; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). Here we also use the digitizing capabilities of a GIS to delineate debris-covered glaciers and LIA trimlines or lateral and terminal moraines.

While for the SGI 2000 the digitized Swiss glacier inventory from 1973 was used to define the ice divides (Reference PaulPaul, 2004), such information was not available for Baffin Island. Moreover, the definition of ice divides is hampered, as many small ice caps or larger glaciers are connected in their accumulation area without visible separation. When separation into individual glacier parts was feasible, we used topographic shading to guide the delineation of ice divides. In other regions we used preliminary ice divides (Fig. 3a) as a first step for defining a glacier basin. For a later glacier inventory, these divides should be adjusted according to the former Canadian inventory, which is straightforward with the digital data in a GIS. The divides were extended to closed polygons, roughly surrounding the actual glacier and its possible LIA extent. Contrast-stretched RGB (red, green, blue) composites (bands 3, 2, 1) and the classified glacier outlines are used as background to assist in the digitizing process. With the ETM3/ETM5 method, nearly all lakes are classified as glacier, independent of their turbidity. They are detached from the glaciers in the course of the basin delineation, as they are clearly visible in the background image. Often a glacier lake fills the depression enclosed by the LIA moraine. In this case, the glacier basin polygon encloses the lake, and an additional line is digitized for later lake separation (cf. Fig. 3b). Obvious snowfields within a glacier basin are also enclosed from a polygon and assigned with a label point for later exclusion from glacier area calculation.

Fig. 3. Digitized glacier basins (black) with automatically classified glacier outlines (white) on a pan-sharpened ETM scene. (a) The preliminary locations of ice divides in regions with unclear separations are indicated by arrows. (b) Splitting of misclassified lakes with glacier basin outlines (top arrow) or additional glacier outlines (bottom arrow). The digitized LIA maximum extent of some glaciers is given by a thick grey line. The lake at the top has been excluded from further calculations, while the outline from the lake at the bottom has partly been used to map the LIA glacier extent.

Little Ice Age moraines

As a result of the harsh climatic conditions and the recent deposition, the lateral and terminal moraines of many glaciers appear very fresh and are free of any vegetation or even lichen cover. Moreover, most ice caps are surrounded by a lichen-free seam (trimline) of variable extent which is quite obvious even in coarser-resolution satellite imagery (Reference Andrews, Davis and WrightAndrews and others, 1976; Reference Locke and LockeLocke and Locke, 1977). However, when these lichen-free areas are not connected to a recent glacier, they may be attributed to persistent snowfields (Reference KoernerKoerner, 1980). A general recession from the LIA maximum extent did not start until the 1920s, as reported by Reference BirdBird (1967) and Reference DowdeswellDowdeswell (1995). Investigations from lichen dating have shown that the LIA maximum glacier extent (about 130–200 years ago) has been quite similar to mid-Holocene or in some places even late Wisconsinan glacier extent several thousand years ago, and lateral moraines from both periods can be found side by side (Reference MillerMiller, 1973; Reference DykeDyke, 1979). With respect to the sensor resolution of ETM and ASTER (15 m), we are not able to discern these two moraines from space and there might be some interferences in interpretation.

The well-preserved indicators of a former glacier extent form a promising base for manual delineation, in particular as such data have not been obtained thus far for a larger sample of glaciers in this region. However, there are at least three critical points for the delineation:

1. 1. the boundary of the trimlines is quite diffuse (ice caps),

2. 2. there are two or more possible locations for the boundary (lateral/terminal moraines), or

3. 3. they are not visible due to shadow.

In the first case, doubtful parts of the former extent were neglected (Fig. 4a); in the second case, the outline is placed at the inner rim (Fig. 4b); and in the last case, the outline stopped or followed the overall glacier shape (Fig. 4c). Thus, the glacier change calculated here is a minimum value. To delineate former glacier extent, we used pansharpened ETM band 4, 3, 2 RGB composites, obtained by the intensity–hue–saturation (IHS) fusion technique. Here the average intensity from the three 30 m bands 2–4 is replaced with the 15m pan band and than transformed back to the RGB colour space. This composite shows much finer spatial details than without fusion, and the near-infrared (NIR) composite has much better contrast in regions with vegetation than the true colour band 3, 2, 1 composite. The digitizing process includes three steps: the initial digitizing, removal of former misclassified boundaries and correction of dislocated connecting points (node errors) or dangling segments.

Fig. 4. Examples for manual delineation of LIA glacier outlines (thick white lines), also showing satellite-derived outlines (thick black lines) and manually digitized debris-covered glacier parts (thin white lines). (a) The arrows indicate questionable trimlines along an ice cap that has not been used for delineation. (b) While the white arrows indicate probable trimlines, the black arrows hint at the outer rim of the lateral moraine. We have used the crest of this moraine for delineation of the LIA outline. (c) Our glacier boundary delineation in regions of cast shadow (arrows) is speculative. The delineation follows the general shape of a glacier and the topography.

General remarks

Some of the major improvements of the ASTER sensor compared to TM or ETM are: (1) the along-track back-looking NIR band (3B) for generation of DEMs, (2) the improved sensor resolution in the visible/near-infrared (VNIR) bands (15 m), (3) the pointing capability of the VNIR sensor (±24°), and (4) the possibility of changing the gain settings of individual bands (Reference KääbKääb and others, 2003a). Among the disadvantages are the smaller ground area covered (factor of nine compared to Landsat) and the missing band in the blue part of the spectrum (useful for lake detection). Apart from DEM creation, the enhanced capabilities of the ASTER sensor have been utilized in numerous glaciological studies (for an overview see Reference KääbKääb and others, 2003a; Reference BishopBishop and others, 2004): measurement of glacier velocities from sequential images (Reference Dowdeswell and BenhamDowdeswell and Benham, 2003; Reference Kääb, Lefauconnier and MelvoldKääb and others, 2005), assessment of glacier hazards (Reference Huggel, Kääb, Haeberli, Teysseire and PaulHuggel and others, 2002; Reference Kääb, Wessels, Haeberli, Huggel, Kargel and KhalsaKääb and others, 2003b), analysis of supraglacial lakes (Reference Wessels, Kargel and KiefferWessels and others, 2002) and mapping of debris-covered glaciers (Reference Paul, Huggel and KääbPaul and others, 2004). The generation and accuracy assessment of ASTER-derived DEMs has been discussed in detail before (e.g. Reference KääbKääb, 2002; Reference ToutinToutin, 2002a, b) and is thus only briefly summarized below.

We apply ASTER data in this study for two major objectives: generation of a DEM and comparison of the glacier mapping results with ETM. While glacier mapping within the ASTER scene is somewhat hampered by snowfields, the DEM creation benefits from the low gain setting in the VNIR bands. The DEM has several potential applications, including orthorectification of satellite imagery, 3-D glacier parameters, perspective visualizations, topographic normalization, potential solar radiation, identification of ice divides, and mapping of debris-covered glaciers. Here we apply the DEM for orthorectification and calculation of 3-D glacier parameters for 340 glaciers with the method described in Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002). The semi-automatic method for delineating debris-covered glaciers (cf. Reference Paul, Huggel and KääbPaul and others, 2004) did not work well in this region and has thus not been applied. For the ETM scene (no DEM available) we applied manual delineation of debris cover for 85 glaciers using a pan-sharpened band 4, 3, 2 composite as background.

DEM generation

The DEM was obtained from the ASTER 3N and 3B channels at 30 m spatial resolution using 17 ground-control points (GCPs) and 25 tie points within standard image-processing software (PCIs Orthoengine). Elevations and locations of the GCPs are taken from the 1 : 250 000 scale Canadian topographic map (sheet 16L) for several mountain peaks and some points at sea level. The average difference of calculated and mapped elevations is –8.6m, with a standard deviation of 155 m. Visual inspection of the shaded DEM reveals only few locations with outliers, mostly located at steep, shadowed slopes and not on glaciers. Small gaps in the DEM are filled by interpolation in the course of geocoding. High-frequency noise in the resulting DEM has been reduced by a Gaussian low-pass filter which changed the original elevations by only a few metres.

Glacier mapping with ASTER

The glacier map from ASTER has also been obtained from a simple band ratio using a thresholded band 3/4 image and an additional threshold in band 1. The latter removes nearly all of the wrongly classified shadow pixels that are located over bare rock, due to the markedly higher reflection of shadowed ice and snow pixels. For a small section shared with the ETM scene, we compared the glacier mapping results of both sensors. Visual inspection reveals very good agreement of the two outlines (Fig. 5a), apart from small shifts due to the missing orthorectification of the ETM scene. A more quantitative evaluation has been performed for 62 glaciers >0.1 km², yielding an average relative area difference of -2.4%, with a standard deviation of 5.6% (Fig. 5b). Thus, ASTER-derived glacier areas are in general somewhat larger, but the deviations can also be explained by differing pixel sizes (Reference Paul, Huggel, Kääb and KellenbergerPaul and others, 2003). For glaciers larger than 1 km² the differences are within ±5%. For clean glacier ice, previous studies have revealed area mapping errors with TM of <5% (e.g. Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others 2002, Reference Paul, Huggel, Kääb and Kellenberger2003). However, this error is largely theoretical, as nearly every glacier is partly covered by debris.

Fig. 5. Comparison of glacier mapping from ETM and ASTER. (a) Overlay of glacier outlines mapped from ETM (black) and ASTER (white) on a small subset (11.5 by 8.5 km) of the ASTER scene (cf. Fig. 1). (b) Scatter plot of relative area differences (ETM– ASTER) vs glacier size for 62 glaciers larger than 0.1 km². The glacier area mapped with ETM is in general a few per cent smaller, and the area difference exhibits increasing scatter towards smaller glaciers.

Size distribution

The distribution of glaciers by number and area for seven distinct area classes is depicted in Figure 6a for 770 glaciers as obtained from the ETM scene. While the percentage of glacier-covered area clearly increases towards larger glacier size classes, the percentage of glaciers by number is more-or-less similar for all size classes less and larger than 5 km², respectively. For Alpine glaciers, nearly the opposite is the case, with more-or-less similar percentages by area for all size classes (summing up all glaciers smaller than 1 km²) and a nearly exponential decrease of the percentage by number towards larger glaciers (Reference PaulPaul, 2004). Thus, we suggest that calculation of global glacier distribution from universal scaling paradigms (Reference BahrBahr, 1997) may give erroneous results. There is not much change in this general pattern if size classes are grouped in powers of two as in the World Glacier Inventory.

Fig. 6. (a) Distribution of glacier size by area (white bars) and number (black bars) summarized for seven area classes. The percentages are quite different from Alpine glacier distribution. (b) Relative change of glacier size with glacier area for 225 individual glaciers. In contrast to Alpine glacier change, there is no increase of scatter towards smaller glaciers and thus a clear dependence of relative area change on glacier size.

Area change since the LIA maximum extent

The relative change in glacier size with area is illustrated for 225 individual glaciers in Figure 6b. The relative loss of area increases towards smaller glaciers (if sizes are <10 km²) and is more-or-less constant for larger glaciers. Interestingly, there is no increase in scatter towards smaller glaciers as observed, for example, in the Swiss Alps (Reference PaulPaul, 2004). In total, the 225 glaciers have lost 217 km² of area (–11%) since their LIA maximum extent. As several authors have pointed out, general glacier recession in this area started in the 1920s (e.g. Reference DowdeswellDowdeswell, 1995), resulting in an average decadal area loss of 27km², or 1.4%. This value is somewhat lower than the 38 km² (–2.2%) area loss per decade found for Alpine glaciers for the period 1850–1973, but still indicates that glacier recession in this region is considerable.

3-D parameters

Several 3-D parameters were obtained for 340 glaciers from the ASTER-derived DEM and the glacier outlines obtained according to the method described in Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others (2002). In order to keep this overview short, the two scatter plots in Figure 7 are presented only as an example of this dataset. Future studies might explore the full potential of the DEM for further applications. From Figure 7a a clear dependence of minimum and maximum glacier elevation on glacier size is evident, although there is some scatter towards smaller glaciers. The scatter plot in Figure 7b illustrates about 300m higher mean elevations for south-oriented than for north-oriented glaciers. For all 340 glaciers, the average of the mean elevation is 1020 m (±217 m), the elevation range is 550 m (±385 m) and the mean slope is 23° (±7.7°).

Fig. 7. Scatter plots of various 3-D glacier parameters as obtained from the ASTER-derived glacier outlines and DEM. (a) Glacier size vs minimum (circles) and maximum (crosses) elevation, indicating a general increase (decrease) of maximum (minimum) elevation towards larger glaciers. (b) Mean glacier elevation vs aspect. A slight increase (about 300 m) of mean elevation towards southern aspects is visible.

Working hours

In view of the GLIMS goal to create a global glacier inventory from space, we have tracked the time required for glacier mapping, in particular glacier basin delineation and manual editing of glacier outlines (debris cover, shadow gaps). We have spent about 22 hours on basin delineation for 880 glaciers (1.5 min per glacier) and about 4 hours on debris mapping and final editing of 85 glaciers (3 min per glacier). This sums to an average 5 min per glacier, or 7 years for all estimated 160 000 glaciers worldwide (based on 240 8 hour days per year). Mapping of LIA outlines took an additional 20 hours for 230 glaciers (5 min per glacier). The mapping time required for a full scene with the ratio method is about 1 hour, including algorithm and threshold selection for one or two test regions. Thus, glacier mapping is very fast, but post-processing takes considerable time.

Additional considerations

While glacier basin delineation is straightforward for most cirques, as well as for valley glaciers which are separated by mountain crests (as in the Alps), there are problems defining the ice divides of many glaciers found in the Arctic. For assessment of glacier change, we have introduced preliminary divides which can be corrected in a later revision. LIA glacier outlines have been delineated because they are prominent in the satellite imagery and it is valuable to quantify glacier retreat in this area. For some LIA glacier outlines, we cannot exclude misinterpretation in this first assessment, but the overall pattern and magnitude of glacier retreat will not change significantly after a later refinement of glacier basins. The 15m ASTER VNIR bands show finer spatial details (cf. Fig. 5a) than the corresponding pansharpened 15m ETM bands. This permits improved delineation of LIA moraines and debris-covered glaciers. The numerous applications of the ASTER-derived DEM were not fully exploited in this overview. However, the examples shown illustrate the potential of such DEMs, even if the cartographic database is coarse (1 :250 000 scale map). From a GLIMS perspective, the ASTER DEM is a vital part of the project, and the combination with the SRTM-derived DEM will be most fruitful in other parts of the world.