Previous applications

The methods for glacier delineation with Landsat TM data used in previous investigations can be divided into three distinct groups (Reference PaulPaul, 2001): (1) segmentation of ratio images from various TM band combinations, (2) unsupervised and (3) supervised classification techniques. Method 1 is used, for instance, by Reference Bayr, Hall and KovalickBay rand others (1994) with digital numbers (DN) from TM4 and TM5 as input, or with the planetary reflectance at the satellite sensor of the same bands by Reference Hall, Chang and SiddalingaiahHall and others (1988) or Reference Jacobs, Simms and SimmsJacobs and others (1997). Reference RottRott (1994) created a glacier mask after the thresholding of a ratio image from TM3 and TM5 but used the atmospherically corrected spectral reflectance of each channel. Reference Aniya, Sato, Naruse, Skvarca and CasassaAniya and others (1996) used method 2 to classify the whole of Hielo Patagónico Sur (southern Patagonia icefield) (ISODATA clustering with TM1,TM4 and TM5 as input). Method 3 was used by Reference Gratton, Howarth and MarceauGratton and others (1990) and Reference Sidjak and WheateSidjak and Wheate (1999) (maximum-likelihood classification). The latter authors also investigated the use of principal-component analysis (PCA) and a normalized-difference snow index (NDSI). Reference Serandrei Barbero, Rabagliati, Binaghi and RampiniSerandrei Barbero and others (1999) created a glacier classification scheme using fuzzy-set theory and a DEM within a GIS framework. So far, however, most methods have been applied only to a smaller number of glaciers (550), and all have been unable to classify the debris-covered ice of a glacier.

Comparison of different glacier-mapping methods from Landsat TM

In this study we apply different glacier-mappingmethods (see below) to a subset of a Landsat TM scene (path 195, row 28) from 12 September 1985. The test region (15 km by 15 km) is located in the Weissmiesgroup in the Saas valley, Swiss Alps (Fig. 1). This region is typical for glaciated environments in Switzerland. It is characterized by steep relief (1500–4500 ma.s.l.), with cast shadow and debris cover on some glaciers. Together with abundant small snowfields, these three influences on glacier-mapping accuracy, known to be critical from previous studies, can be examined in this test region.

Fig. 1. Location of the test area, the Weissmies group of mountains in Switzerland (see inset), and as seen with Landsat TM (band 3) on 12 September 1985 (contrast enhanced). Black lines indicate the glacier outlines from the digitized glacier inventory of 1973. Size of imagery is about 15 km by 15 km. Landsat TM data © Eurimage.

In Figure 2a–d we present the results from different glacier-mapping methods, comparing two glacier maps at a time in each figure: (a) segmentation of a ratio image from TM3/TM5 vs TM4/TM5 using the DN (Fig. 2a); (b) as (a), but using the spectral reflectance instead of DN (Fig. 2b); (c) an unsupervised ISODATA clustering with 20 classes vs TM4/TM5 from DN (Fig. 2c); and (d) a supervised maximum-likelihood classification with eight classes vs TM4/TM5 from DN (Fig. 2d). Further methods were applied (NDSI, PCA, usage of atmospheric-corrected TM bands), but results are not shown because they were less accurate.

Fig. 2. Glacier masks comparing two methods at a time. Areas in light grey were identified by both methods as being a glacier, dark grey areas only by the first method, and black areas only by the second method. The compared methods are: (a) TM3/TM5 and TM 4/TM5 from DN; (b) as (a) but all TM bands from satellite planetary reflectance; (c) with an unsupervised ISODATA clustering (20 clusters) algorithm and TM4/TM 5 from DN; (d) with a supervised maximum-likelihood classification of training areas with eight classes and TM4/TM5 from DN.

A glacier map (black = glacier, white = other) is created by interactive thresholding of the ratio images for methods (a) and (b). The 20 classes of (c) were separated into glacier and other by visual interpretation. For the maximum-likelihood classification (d), training areas in eight classes were chosen: glacier 1 and snow 1 (in sunlight), glacier 2 and snow 2 (in shadow), forest, meadow, terrain and cast shadow. For the final glacier map the latter four were converted to other and the first four classes to glacier. For each of Figure 2a–d two glacier maps were combined with the following colour scheme: glacier onboth maps: light grey; glacier only on the first/second map: black/dark grey; and, other on both maps: white. To improve the quality of the classification, a 3 by 3 median filter was applied to all glacier maps before combination. A more detailed analysis of 32 glaciers reveals that the average change in glacier area by the median filter is –0.4%, if glaciers smaller than 0.1km² were not considered.

All methods other than TM4/TM5 with DN reveal problems with regions in cast shadow (indicated by arrows in Fig. 2a), where they map too much (Fig. 2a, b and d) or too little (Fig. 2c) glacier area. Additional regions within cast shadow are mapped from both methods displayed in Figure 2b. Small snowfields were mapped with the methods displayed in Figure 2c and d. The accuracy of all investigated methods could be improved partly by changing the relevant parameters (thresholds, training areas, number of clusters) but at the cost of more incorrect results in other places. All methods fail to detect debris-covered ice because of the spectral similarity to the surrounding terrain. The accuracy of the glacier classification with segmentation of a TM4/TM5 ratio image using the raw DN proved to be the best method with respect to glacier areas in cast shadow or assigning snowfields to other.

Accuracy of the best glacier-mapping method

To evaluate the accuracy of this best-suited classification method, the TM-derived glacier areas were compared with areas derived manually from a higher-resolution SPOT Pan scene (10 m). Unfortunately, this scene (path 55, row 256, acquired on 17 September 1992) does not cover the Weissmies test area, but it shares a small region with a TM scene (path 195, row 28), acquired only 2 days prior to the SPOT scene. Because of the good temporal coincidence, another test site (located to the south of the Nufenenpass) depicted on both the TM and SPOT scenes was selected. For 32 glaciers within this site, the automatically TM-derived areas were 2.3% smaller (on average) than on the manually analyzed SPOT image. This deviation is well within the accuracy of the manual glacier delineation regarding small snowpatches or the delineation of debris-covered areas. Thus, for debris-free ice, the accuracy of the glacier areas inferred from TM is better than about 3%.

The accuracy is illustrated in Figure 3 for a region 5 km by 7 km in size showing Cavagnoli (C) and Basodino glaciers (B). The outline as derived from the TM4/TM5 glacier-mapping method (black) is superimposed on the SPOTscene together with the glacier outline of 1973 (white) from the digitized Swiss glacier inventory. The depicted overlay suggests the following:

Fig. 3. Cavagnoli (C) and Basodino glaciers (B) with outline from TM (black) from 15 September 1992 and the Swiss glacier inventory from 1973 (white) on a SPOT Pan scene from 17 September 1992. This area (about 5 km by 7 km) is located 45 km northeast of the test area Weissmies. SPOTdata ß SPOT Image.

1. (1) the TM-derived glacier outlines fit quite well to the visible glaciers,

2. (2) small isolated ice fields are not classified as glaciers,

3. (3) the smallest glaciers shrank through disintegration into snowpatches,

4. (4) there is a differentiated retreat of the larger glaciers,

5. (5) the largest glacier (Basodino) was even larger in 1992 than in 1973.

Requirements and possibilities

A DEM has two main functions within the SGI 2000: the orthorectification of the satellite imagery, and the derivation of three-dimensional glacier parameters within a GIS. The orthorectification is mandatory for at least four different tasks:

1. (1) To eliminate the effects of perspective distortion, terrain elevation has to be considered during georectification of imagery of rugged terrain. For instance, a pixel with a height of 3000ma.s.l., located 90 km from the nadir point, is shifted by 370 mfrom its real position in the uncorrected image.

2. (2) The borders between individual glaciers were assigned to the classified TM image from a (georeferenced) vector layer containing digitized glacier basin outlines.

3. (3) Overlay of TM scenes from other years or with scenes from other sensors.

4. (4) Fusion with the DEM itself used to derive glacier parameters.

All scenes for SGI 2000 were orthorectified with a set of ground-control points to a residual rms error of about half a pixel, using a DEM with 25m spatial resolution (SGI 2000 DEM) from the Swiss Federal Office of Topography.

Three-dimensional glacier parameters, like minimum and maximum elevation, glacier length, the median and 2:1 altitude (equilibrium-line altitude), and a detailed hypsography, can be computed automatically with a DEM. Slope and aspect of each glacier can be obtained as an average for the entire glacier or as a percentage of selected zones (e.g. accumulation and ablation area). Moreover, average illumination or the percentage area in cast shadow during a day can be calculated. In this way it is possible to achieve a more thorough understanding of topographic influences on monitored changes in glacier area or length.

Comparison with a reference DEM

To analyze the vertical accuracy of the SGI 2000 DEM, a comparison with a reference DEM directly inferred from stereo-photogrammetry was performed (Reference KääbKääb, 2000). An illuminated version of the SGI 2000 DEM is shown in Figure 4a for a small area (5.7 km by 5.0 km) within the test region. Also indicated are the outlines of seven glaciers, analyzed in the discussion below. Artefacts from the interpolation process between the originally digitized contour lines are visible on the illuminated SGI 2000 DEM. They are notably pronounced in gradient products like slope, as illustrated in Figure 4b, which shows the difference in slope to the reference DEM.

Fig. 4. (a) Illuminated version of the DEM used for the SGI 2000 in a subsection of the test area Weissmies with outlines of seven glaciers (numbered) from the digitized inventory of 1973. DEM artefacts are clearly visible. (b) Differences in slope between the DEMin (a) and a reference DEM ranging from –53˚ (black) to +65˚ (white).Most of the artefacts are located outside of glacier areas. Elevation data: DEM 25 © Swiss Office ofTopography (BA013305).

The minimum (maximum) elevation differences are –96 (+74)m for the entire area, with a standard deviation of 8.6 m. The corresponding values for the slope differences are –53 (+60)˚ and 6.6°. Those deviations are not acceptable for the SGI 2000, but the large differences were mostly found at isolated locations or crests, usually not related to glacier coverage. To estimate the influence of the artefacts on the derived glacier parameters, some of them were calculated for the seven glaciers indicated in Figure 4a. The average differences of minimum (maximum) elevation are –2.7 (–1.3)m, with a standard deviation of 7.6 (7.4)m. The corresponding values for slope are 1.2 (–2.5)˚ and 2.4 (3.3)°. These deviations are acceptable for the glacier parameters in the SGI 2000.

Data preparationsFootnote 1

By using an orthorectified glacier map derived from TM and a suitable DEM, three-dimensional glacier parameters can be obtained automatically within a GIS. Before the GIS-based processing was started, all data products were converted into Arc/Info (ESRI, 1999) specific formats, and three GIS-related tasks were prepared for the SGI 2000:

1. (1) digitizing of the glacier outlines from the inventory of 1973 into a vector layer,

2. (2) creation of a vector layer with glacier basin boundaries (see below),

3. (3) calculation of DEM products (e.g. slope, aspect) for obtaining three-dimensional glacier parameters.

The glacier outlines were digitized from the original maps (scale 1:25 000) as individual arcs, with an average rectification error of each map of about 5 m (rms). The central flowlines and the reconstructed outlines from about 1850 were also digitized.

To separate connected glaciers in the classified TM image into individual glaciers, the ice divides between them must be defined. This is done by on-screen digitizing of a new vector layer (coverage) using the digitized Swiss glacier inventory and the classified TM image as background in arcedit. Ice divides were taken without modification from the digitized inventory and extended to a closed polygon roughly surrounding the glacier. All other glaciers were also surrounded by closed polygons, which are large enough to include possible future variations of glacier area. The thick black lines in Figure 5 represent these polygons. They are shown together with the digitized glacier areas (in grey). With these predefined glacier basins it is also possible to assign a unique identity (ID) to glacier groups (see below). A group of glaciers can originate from a single glacier through disintegration over time, already be established in a former inventory or consist of an entire glacier comprised of different streams.

Fig. 5. Glacier basins (bold black lines) and the digitized Swiss glacier inventory of 1973 (grey areas). Contiguous ice masses were separated into individual glaciers according to the 1973 inventory, and roughly surrounded to obtain closed polygons. These polygonal glacier basins can also be used to eliminate gross classification errors (e.g. proglacial lakes) and for glacier identification, especially to assign a unique ID to glaciers.

DEM products such as slope or aspect were calculated within the digital image-processing software. Some DEM products are further converted with short FORTRAN programs (e.g. the aspects are classified into eight sectors). The elevation products (e.g. median elevation or hypsography) are computed within Arc/Info.

The location of the glacier ID is taken from the revised database of Reference Maisch, Wipf, Denneler, Battaglia and BenzMaisch and others (1999). Because many glaciers have split up during recent years, the location of their ID (assigned in 1973) often lies outside their present outline. Therefore, and to handle groups of glaciers, the ID is assigned to the entire basin. This is done by converting the database table (ID, x and y-coordinates) to a point coverage with generate and intersecting this coverage with the glacier basin coverage. Thus, each glacier basin holds the glacier ID in the attribute table.

Data flow

The data processing can be separated into a general work flow between three modules and a more specific data flow within each module. The modules are (Fig. 6):

Fig. 6. Work flow between sequential modules and principal data flow within each module. In each module different programs are available to calculate glacier parameters and their changes. Most of them can be combined into a complete digital processing chain.

1. (1) the GIS module for calculating quantities within the GIS (e.g. areas of individual glaciers),

2. (2) the CONV module for conversion and recalculation of Arc/Info output tables (e.g. glacier areas of two different years into relative changes in area),

3. (3) the VIS module for creation of graphic files from data files (e.g. with XMGR, GMTor IDL).

Each module holds individual programs for individual tasks, each of which consists of an input, calculation and output part. The different modules can be combined into a complete digital chain, with the output from a program in one module as the input for a program in the next module (cf. Fig. 6).

The GIS module is shown in Figure 7 in more detail. Only three input layers are needed for the GIS module:

Fig. 7. The principal data flow within the GIS module using Arc/Info and including a DEM product for three-dimensional glacier parameters. Areas of individual glaciers are available after step 2.

The calculation section consists of the following steps (numbers refer to Fig. 7): The first step (1) is a raster-vector conversion of the glacier map with imagegrid and gridpoly into a coverage. This glacier coverage is then (2) combined with the glacier basin coverage with intersect to obtain the individual glaciers. Together with intersect the glacier basin coverage cuts each glacier out of the glacier coverage in correspondence with its basin. The coverage with the individual glaciers (or a selection of them) is then (3) converted with polygrid to a zonal grid where each zone corresponds to a glacier. The DEM (or a product of it) is (4) converted with asciigrid to a value grid, where each cell holds the value of the DEM at that location. The last step (5) is the combination of the value grid and the zonal grid with zonalstats. This command gives for each zone (glacier) statistic parameters (e.g. minimum, maximum, range, mean, standard deviation) according to the underlying value grid (elevation in case of the DEM). The resulting output table can be processed further with the CONV and VIS modules.