Glaciers are sensitive indicators of climate change, growing and wasting in response to changes in temperature and precipitation (Reference Fitzharris, Hay and JonesFitzharris and others, 1992; Reference HaeberliHaeberli, 1995; Reference Oerlemans and ReichertOerlemans and Reichert, 2000). Regular glacier inventories are needed for understanding global and local climate change patterns, as well as monitoring water resources (Reference ChinnChinn, 2001; Reference KargelKargel and others, 2005). These inventories should be repeated at intervals of a few decades (Reference Haeberli, Hoelzle, Paul and ZempHaeberli and others, 2007). Regular glacier inventories can be made by the use of satellite remote sensing, capable of acquiring comprehensive, uniform and frequent global observations of glaciers and ice sheets. In New Zealand the first and only glacier inventory including the two main islands’ glaciers was made from aerial photographs recorded in 1978 (Chinn, unpublished). The mapping showed that the Southern Alps hosted 3144 glaciers exceeding 0.01 km2, totalling an area of 1158 km2 and an estimated ice volume of 53.3 km3 (Reference ChinnChinn, 2001). Only 0.5% of these glaciers were >10 km2, and 78.2% were <0.2 km2. In agreement with the repeat interval suggested for national glacier inventories, an update of the 1978 inventory was recommended.
The goals of international glacier monitoring have evolved during the past century. At the initiation of worldwide collection of information about ongoing glacier changes, at the 6th International Geological Congress in Zürich, Switzerland, in 1894, it was hoped that long-term glacier observations would give insight into processes of climatic changes. In 1986 the World Glacier Monitoring Service (WGMS) took over the maintenance and collection of information on ongoing glacier changes from the World Glacier Inventory (WGI). WGMS collects standardized observations on changes in mass, volume, area and fluctuations as well as statistical information on glacier inventories. Even though satellite imagery offers an unprecedentedly powerful and efficient medium with which to study glaciers that are usually located in remote, inaccessible and inhospitable environments (Reference Gao and LiuGao and Liu, 2001), Reference AlbertAlbert (2002) underlines that satellite sensors are still limited in their measurement resolution, and some variables can be difficult to observe from these platforms. For snow and ice purposes, these could include proper discrimination between glacier ice and the surrounding glacial moraines and proper recognition of glacier ice in shaded steep mountain walls. Reference AlbertAlbert (2002) therefore stated that ground-based observations and field research would continue to play a crucial role in the coming years. The Global Land Ice Measurements from Space (GLIMS) program was initiated in 1999 by the US Geological Survey with the aim of mapping world glaciers and extracting glaciological information through the use of satellite imagery. Over 60 institutions across the globe are involved in GLIMS, including institutions in New Zealand (Reference KargelKargel and others, 2005; Reference RaupRaup and others, 2007). GLIMS primarily uses data from optical satellite instruments, such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) sensor on board the Terra spacecraft. This instrument includes shortwave infrared (SWIR) spectral bands and has a spatial resolution up to 15 m in the visible, making it well suited to meet the accuracy requirements for New Zealand conditions in which many glaciers are small and located on steep slopes (Reference KargelKargel and others, 2005).
Multispectral imaging in the visible and infrared is unsurpassed for systematic global study of glaciers (e.g. Reference AlbertAlbert, 2002; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002; Reference RaupRaup and others, 2007). Several approaches for glacier classification of multispectral imagery have been tested (e.g. manual digitizing, band ratios, normalized-difference snow index (NDSI), supervised classification) (Reference KääbKääb, 2005). From a wealth of experience and a comprehensive review, Kääb concluded that human interpretation remains the best tool for extracting higher-level information from satellite imagery for many glacier types. The analyst is able to include experience, knowledge and complex logical rules in the decision and digitizing process, also relying on nonspectral data or knowledge. Manual delineation is also often needed to complement and correct automatic classifications, such as the extraction of debris-covered glaciers (Reference KääbKääb, 2005). Automated glacier extraction using band ratios (e.g. near-infrared (NIR)/SWIR) with certain thresholds or normalized band differences like the NDSI (Green–SWIR)/(Green + SWIR) has proven to be an accurate, fast and robust method for detecting clean glacier ice (Reference AlbertAlbert, 2002; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002; Reference KääbKääb and others, 2003; Reference ShangguanShangguan and others, 2006; Reference Ye, Kang, Chen and WangYe and others, 2006; Reference Bolch, Menounos and WheateBolch and others, 2010). These methods are based on the distinct low reflectance of ice and snow in the SWIR part of the spectrum and its high reflectance in the visible part. In supervised classification, training areas of the categories to be mapped are operator-selected and the spectral signatures of these classes are used to automatically segment the entire imagery (Reference KääbKääb, 2005). Yet, like all spectral classifications, it fails to efficiently separate debris-covered ice from the surrounding non-vegetated glacier forefield (Reference Paul, Huggel and KääbPaul and others, 2004a; Reference ShangguanShangguan and others, 2006).
The overall aim of this work was to quantify the glacier area changes in the central Southern Alps of New Zealand between 1978 and 2002 and to test various automated classification methods for glacier extraction from multi-spectral ASTER data. The glacier changes were obtained by comparing the New Zealand glacier inventory produced from aerial photographs recorded in 1978 (Chinn, unpublished) with the updated inventory for the central parts of the New Zealand Alps produced in this study from an ASTER image recorded in February 2002.
2. Study Area
The central Southern Alps of New Zealand (43.36–44.03° S, 169.79–170.75° E; Fig. 1) were targeted in this study as they include by far the most heavily glacierized area in the country, as well as New Zealand’s largest glaciers. In addition, the study area is the most accessible glacierized area of New Zealand, although the mountains of New Zealand are generally very remote. Fieldwork often requires costly helicopter logistic support. All New Zealand glaciers are distributed along the Southern Alps, from 42° S to 45° S, except for some glaciers on the volcano Mount Ruhapehu on the North Island (Chinn, unpublished). The number of individual glaciers is relatively high, but the average size is small because of the steep topography, which favours individual small cirque glaciers and ice patches on separate peaks (Reference ChinnChinn, 2001). The average summit height across the range varies from 1850 m a.s.l. in Fiordland in the south to 3000 m a.s.l. in the central Southern Alps (including the country’s highest summit, Mount Cook, at 3754 m a.s.l.) and 2000 m a.s.l. north of the central Southern Alps.
Precipitation is generally high, with a strong east–west gradient induced by the north–south orientation of the range across the path of prevailing westerly winds. Precipitation increases rapidly from ∼3 m a−1 on the western coast to >10 m a−1 just west of the Main Divide, which separates western and eastern slopes, and drops to ∼1 m a−1 40 km east of the Main Divide (Reference Griffiths and McSaveneyGriffiths and McSaveney, 1983). This strong west–east orographic precipitation gradient creates a steep eastward rise of glacier equilibrium-line altitudes (ELAs). New Zealand glaciers are mainly high-activity maritime types with precipitation at or well above 3 m a−1. Extreme maritime glaciers occur west of the Main Divide, with ‘dry’ balance glaciers and rock glaciers prevailing east of the divide (Reference ChinnChinn, 1999).
Recent work by Chinn and others (unpublished information) suggests that ice volume in New Zealand glaciers decreased 15% from 54.5 km3 in 1976 to 46.1 km3 by 2008. Annual changes are believed to be due mainly to mass-balance changes driven by climate (precipitation and temperature), but over the entire period much of this decrease was attributed to loss of ice volume for larger glaciers due to calving into growing glacial lakes and trunk downwasting (Chinn and others, unpublished information). Over a much longer period, Reference Hoelzle, Chinn, Stumm, Paul and HaeberliHoelzle and others (2007) reported a general pattern of ice retreat in New Zealand since the Little Ice Age (LIA). They used a parameterization scheme with four variables – surface area, total length, maximum and minimum altitude – to estimate specific mean mass balance and glacier volumes between the ‘1850 extent’ and the volume in the mid-1970s. They calculated that area change in this period was a decrease of 49%, with a corresponding volume loss of 61%. The calculated loss in net mass balance varied between 0.67 m w.e. a−1 in the most maritime area and 0.54 m w.e. a−1 in the most continental area. However, some of the glaciers in New Zealand also experienced advances, commencing in the early 1990s and ceasing around 2000, which were more extensive than any other since the end of the LIA. The positive glacier balances were associated with an increase in the strength of westerly atmospheric circulation which brought increased precipitation (Reference Chinn, Heydenrych and SalingerChinn and others, 2005). Evidence of glacier advance despite the general picture of glacier retreat was shown by Reference Allen, Owens and SirgueyAllen and others (2008). In that study, multitemporal glacial mapping revealed some small steep glaciers significantly advancing between 2002 and 2006, up to 100–140 m for Stocking and Eugenie Glaciers, east of the Main Divide.
A complete and unique inventory of New Zealand glaciers was carried out by T. Chinn using aerial photographs from 1978 (Chinn, unpublished). Glacier outlines were photo-interpreted from the oblique aerial photographs recorded in 1978 and were digitized onto topographic base maps at a scale of 1 : 63 360 (imperial system where 1 in = 1 mile). The contour interval of these maps was 30.5 m (100 ft). In addition to mapping all the glaciers, 48 index glaciers, well distributed along the Southern Alps, were chosen for annual snowline survey. Oblique aerial photographs of the index glaciers have been taken every year since 1978, except in 1979, 1990 and 1991.
For updating the New Zealand glacier inventory in the main glacierized area, the central Southern Alps, an image from the ASTER sensor was used. ASTER captures high spatial resolution data in 14 bands, operating in three subsystems: visible and near-infrared (VNIR) consisting of three spectral bands with a spatial resolution of 15 m; SWIR consisting of six spectral bands with a resolution of 30 m; and thermal infrared (TIR) consisting of five spectral bands (Table 1). The best available image was selected considering the following criteria: (1) cloud-free, (2) end of summer, and (3) year with a negative mass balance. The image chosen was recorded on 14 February 2002 and was a processed level-1B image (Fig. 1, image ID: SC:AST_L1B.003:2019534558). The year 2002 was a very suitable year for compiling the glacier inventory since it had the most negative mass balance between 2001 and 2005 (Reference Chinn, Heydenrych and SalingerChinn and others, 2005). This meant that most of the surrounding snowpatches remaining in late summer in other years were gone and could not be misinterpreted as glaciers (Reference ShangguanShangguan and others, 2006). The image contains ∼40% of the total glacier area mapped in the 1978 glacier inventory. Six of the forty-eight index glaciers used for the annual snowline monitoring are also situated within the image (Fig. 1). Landsat images might have permitted us to investigate a longer time series, but were not considered here due to their lower spatial resolution (30 m). ASTER was believed to have the appropriate spatial resolution to map New Zealand glaciers. Changes in glacier area of the small, alpine glaciers, which contribute by far the largest numbers of individual glaciers in New Zealand, would have been more difficult to detect.
The image was orthorectified (PCI Geomatica OrthoEngine software) to ensure an appropriate mapping accuracy. A total of 31 survey points were collected in the field from late February until late April 2005 using a Trimble Pro-XR GPS instrument. These features were relatively easily detectable on the satellite image and were chosen as ground control points (GCPs). These were, for instance, centers of bridges, sharp intersections between scree slopes and vegetation, or characteristic rock features in the fringe of glaciers. Only 22 of these GCPs were later found to be of high enough quality to be used for further processing, i.e. easy to separate from the surrounding pixels in the image and easy to recognize in the field. The GPS data were post-processed using the carrier phase technique (sub-metre theoretical accuracy) with reference data collected from the Trimble 40 GPS Base Station located at the University of Otago. The field-collected survey points were divided into GCPs and independent check points (CPs). The first orthorectification performed with a digital elevation model (DEM) from the Land Information New Zealand (LINZ) topographic database, with a pixel size of 15 m, showed large shifts in pixel position when overlaying geographical features from LINZ. This was likely due to the planimetric and altimetric extrapolation of the geometric sensor model resulting from the spatially clustered and uneven three-dimensional (3-D) distribution of the field-collected control points (Reference ToutinToutin, 2004). To improve the sensor model and orthorectification result, height points from the New Zealand topographic database were added. Finally 28 GCPs and 15 CPs were used for the orthorectification. The final orthorectification error (root-mean-square (RMS)) for the GCPs (X and Y) was 15.3 and 14.5 m, respectively, and for the CPs (X and Y) was 17.1 and 14.7 m, respectively.
3.2. Extracting land ice area using classification methods
The automatic extraction of debris-free glaciers from multi-spectral satellite imagery relies on the spectral uniqueness of ice and snow compared with surrounding targets (soil, vegetation, rock): high reflectance in the visible (0.4–0.7 μm), and medium and low reflectance in the NIR and SWIR (0.7–2.5 μm), respectively. Glacier boundaries have been mapped successfully in various places in the world using different techniques: (1) photo interpretation and manual delineation on colour composites (Reference Williams, Hall and ChienWilliams and others, 1997; Reference Khromova, Osipova, Tsvetkov, Dyurgerov and BarryKhromova and others, 2006; Reference BeedleBeedle and others, 2008); (2) image segmentation (or thresholding) of principal components (Reference Sidjak and WheateSidjak and Wheate, 1999), single bands (e.g. Thematic Mapper 3 (TM3); Reference Bronge and BrongeBronge and Bronge, 1999), single band ratios (e.g. ASTER 3/4; Reference KääbKääb and others, 2003) or spectral indices (e.g. NDSI; Reference Hall, Bayr, Bindschadler, Shöner, Hardy and FrankensteinHall and others, 2001); and (3) unsupervised and supervised classifications (e.g. ISODATA, maximum likelihood algorithm; Reference Aniya, Sato, Naruse, Skvarca and CasassaAniya and others, 1996). Thresholding of band ratios or spectral indices is generally the preferred approach because of its simplicity, low requirements in terms of time and processing power, and ease of implementation in an operational context. Only one study reports on the extraction of land ice in New Zealand and it used the NDSI approach (Reference Mathieu, Chinn and FitzharrisMathieu and others, 2009). However, the delineation of ice/non-ice targets was not assessed in great detail as the main aim was to assess the suitability of ASTER data to extract ELAs. Advanced image-processing techniques such as spectral unmixing have been used efficiently for retrieving subpixel fraction cover either for snow or ice (Reference Nolin, Dozier and MertesNolin and others, 1993; Reference Klein and IsacksKlein and Isacks, 1999). However, a per pixel approach was considered more appropriate in this work because subpixel fraction cover derived from a 15 m spatial resolution ASTER image would not provide any additional accuracy benefits when comparing with the 1978 national glacier inventory. Indeed, the glacier boundaries in the earlier inventory are believed to be at best ±30 m accurate, corresponding to the accuracy of the topographic base maps (1 : 63 000) used to report the photo-interpreted features.
In this study, we assess the suitability of three common classification techniques to extract land ice in the New Zealand context: (1) band ratio (ASTER 3/4), (2) NDSI (ASTER 1 − ASTER 4)/(ASTER 1 + ASTER 4), and (3) supervised classification (maximum likelihood algorithm). The ratio of TM4/5 (ASTER 3/4 equivalent) was found to be the most appropriate method for the new Swiss glacier inventory (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002). As we were only working with one image, no atmospheric correction was performed. Further, ratios were performed directly on raw digital numbers as this was found to be the most robust approach in an operational context (e.g. Reference AlbertAlbert, 2002; Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002, Reference Paul, Huggel and Kääb2004a).
In the case of the ASTER 3/4 ratio, thresholds of 1.8, 2.0, 2.2 and 2.4 were tested to separate glacier and non-glacier. For the NDSI, thresholds of 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65 and 0.7 were tested. Nearly all lakes in the area were classified as glaciers, independent of their turbidity in the NDSI-based glacier products. Turbid water, like the Tasman Lake, was still classified as glacier even with the highest threshold tested. However, this was not a significant problem for the ASTER 3/4 ratio. All lakes were clearly visible in the colour composite and were therefore manually digitized and removed from the NDSI results. The supervised classification was only tested for one small test area of the image. Supervised classification demands substantial work for selecting the best possible training sites.
Significant misclassifications on the glaciers were expected because of the presence of heavily debris-covered glacier tongues in many of the largest glaciers in the area, in addition to the extreme vertical relief causing large shadows. Consequently, it was decided to manually digitize the glacier areas of the entire image. This option was considered the most reliable approach to quantify glacier changes between 1978 and 2002. The manual delineation was done by photo-interpreting the image displayed as colour composite using the VNIR channels (NIR, red and green). It is important to note that individual glacier boundaries were not mapped in this study, but only the land ice area as a whole. Glacier area changes were calculated from glacier area polygons as measured in 1978 and 2002, using simple change overlay techniques. In addition, more detailed maps were produced to analyse how the changes were distributed in space.
3.3. Validation of the classification methods
The classification methods were validated by (1) qualitatively comparing GPS outlines for selected glaciers with the manual and automated classifications and (2) quantitatively comparing the manual classification (i.e. digitized product) with the automated classifications (indices and supervised classification). GPS outlines of selected glaciers were surveyed in the field for the validation of the classification methods. To reduce the potential negative effects of the time difference between the fieldwork (2005) and the image collection (2002), index glaciers were preferred since aerial photographs of these glaciers existed for each year between 2002 and 2005. The index glaciers are monitored annually during the snowline survey, and possible changes could therefore be carefully considered. The surveyed glaciers are well distributed over the study area, covering both eastern and western glaciers (Fig. 1). Four of the six index glaciers located in the area were of convenient size and accessibility for GPS ground mapping: Jalf/Bauman, Chancellor Dome, Langdale, and Ridge glaciers (Fig. 2). In addition, Tasman (which is also an index glacier), Hooker, Fox, Franz Josef and Stocking glaciers were also investigated in the field (Fig. 3). GPS outlines of the individual glaciers were obtained by walking along the boundary between snow and rock. The GPS receiver was set to record a new position every 5 s, calculated from at least five satellites. This criterion was decided in order to ensure a high positional accuracy of the GPS data. The glacier outlines were post-processed using the carrier phase technique with reference data collected from the Trimble 40 GPS Base Station. Both the carrier phase processing technique and the constant lock of the GPS receiver with at least four satellites ensured systematic sub-metre accuracy, even with a limited number of position records for a single point (i.e. dynamic GPS data collection). Higher accuracy is achieved by using the higher frequency of the L1/L2 carrier signal as opposed to the pseudo-random code. Hence, the accuracy of the GPS glacier outlines was much more dependent on the field interpretation rather than on inaccuracies in the GPS recordings. For some glaciers and because of safety issues, the GPS positions could not be recorded at the exact boundaries (e.g. frontal part of Hooker and Tasman Glaciers (Fig. 3)). In these cases, the average distance to the exact boundaries was noted and the positions were later adjusted in Geographical Information System (GIS) software. Once the manually digitized map was assessed with the field-surveyed GPS outlines, this map was used as an efficient way to select the best thresholds for the index-based classifications and to validate all automated classifications.
4.1. Performance of manual interpretation
In this section we first present the assessment of the manual interpretation because it served as a validation set for the automated classification methods. The manually digitized glacier outlines were compared with the GPS-collected glacier outlines in order to assess the accuracy of the interpretation. For the smaller glaciers, there is a good correlation between the manual digitizing and the field outlines (except from Stocking Glacier, which has been retreating between 2002 and 2005). This indicates that the manual digitizing was fairly accurate for these small glaciers, and that they have changed little during the 3 years (Fig. 2). The validation data as well as the operator experience suggest that the manual digitizing is accurate within 2 pixels (30 m) (±15 m). This would ignore the fact that glacier/nonglacier boundaries are not hard but tend to be fuzzy in nature, and that small discrepancies might still be attributable to time difference between the imagery and the fieldwork. Although difficult to evaluate, this error is likely to increase in the case of steeper glaciers, in shadow where the snow/ice and rock contrast will not be as good, and in the case of debris-covered areas. The termini of the larger glaciers had changed too much in the same period to be used for validation purposes. Nevertheless, the lateral margins of these glaciers also seemed to match relatively well (within 2–3 pixels), despite the difficulty in separating rock and debris-covered ice (Fig. 3). Considering the previous observations the overall accuracy of the digitizing of the 2002 glacier outlines is estimated to be within 3 pixels (±22.5 m).
4.2. Performance of automated classification techniques
The automated classifications were compared visually with the GPS outlines and the manually digitized outlines. Then confusion matrices were computed between the manually digitized glacier maps (acting as a reference dataset) and the automated classifications (Tables 2, 4 and 6). Producer’s, user’s and overall accuracy were calculated from the confusion matrices (Tables 3, 5 and 7). The overall accuracy is the total proportion of pixels classified correctly. User’s and producer’s accuracy are equivalent to type I and type II errors (for details see Skidmore, 1999). For simplicity, accuracy measures are only presented for the threshold values selected as optimum (2.0 for ASTER 3/4 ratio and 0.5 for NDSI).
In the confusion matrices the ideal threshold is where both the omission error (non-glacier detection) and commission error (false glacier inclusion) are as low as possible. However, if choosing the lowest value of one of the errors, the other will be very high. Therefore a choice has to be made when choosing thresholds based on what one wants as an outcome. If glacier mapping over a large and heterogeneous area is the main target (as for this case), the ideal threshold should not exclude too much (what is glacier) nor include too much (what is not glacier). If only the small and very steep glaciers were of interest it could be wise to choose a threshold that included more of the glaciers in shaded areas (and a higher commission error would have to be accepted). In the New Zealand context, with a high degree of debris-covered ice on many long valley glaciers, manual correction was required, as most of these areas were left out of the classifications. Thus, for the ASTER 3/4 ratio and the NDSI, the task was to find the threshold where a minimum manual adjusting of land ice boundaries was needed. By inspection of the ASTER 3/4 ratio image, the best threshold was found to be ∼2.0 (Fig. 4), this result being also confirmed with the confusion matrix (Table 2). However, major changes in glacier areas will be detected even though the threshold used is not optimal.
In Tables 2 and 4 (ASTER 3/4 and NDSI results for the entire image) we note that >20% of the glacier area defined in the photo interpretation is not detected in the classification (Omission error). Most of this error by far can be attributed to debris-cover glacierized areas of long valley glaciers and others which are not detected as glaciers. We also note that applying a 3 × 3 median filter to this classification slightly increases the accuracy (Table 3). This filtering was efficient for removing noise in areas with dense vegetation, but outside these areas it did not influence the classification results significantly. Areas defined as non-glacier in the photo interpretation, but classified as glaciers in both the ASTER 3/4 and the NDSI maps (Commission error), were typically very bright rocks or relatively small snow/icy patches, often located on terrain features where it was not likely to be a glacier (e.g. small hills separated from the surrounding topography). The initial test of supervised classification did not give very promising results. We compared the results obtained for the supervised classification over a limited area around Franz Josef Glacier (Tables 6 and 7; Fig. 5) with the results obtained over the same area for the best thresholds applied to the ASTER 3/4 ratio and the NDSI (results not shown). The supervised classification both excluded more area classified as glacier in the digitizing and included more area not classified as glacier in the other two methods. Table 6 shows this clearly, where both the omission and commission error are very high. The overall accuracy given in Table 7 is also very low. The strips clearly visible in Figure 5a are due to residual errors in the radiometric correction of the SWIR band, and are in the order of 1 digital number.
4.3. Glacier area change from 1978 to 2002
In the 2002 ASTER image the area of the manually digitized glacier outlines totalled 428 km2. This represented almost 11% of the entire image. Of this 428 km2, 242 km2 were located on the eastern side of the Main Divide and 186 km2 on the western side of the Main Divide (Table 8). The glacier polygons of the 1978 inventory were clipped to match the ASTER image. Within the limits of the ASTER image, a total of 513 km2 of glacier area was calculated, which accounted for ∼41% of the glacierized area in New Zealand in 1978. Of this total, 296 km2 were located east of the Main Divide and 217 km2 west of the Main Divide (Table 8).
Changes during the 24 year period in the central Southern Alps indicate that the total glacier area decreased from 513 km2 to 428 km2 (a loss of 16.6%). We estimated the uncertainty for the calculated areas to be ±1.1 % (or ±5 km2) for the 1978 glacier outlines and ±0.6% (or ±2 km2) for the 2002 glacier outlines. The uncertainty in the total glacier area change was ±6.8% (or ±6 km2) (see section 5.2 for details of uncertainty estimations). The glacier retreat was smaller on the western side of the Main Divide (14.3%) than on the eastern side (18.3%). The two main reasons for this are that (1) the two largest glaciers on the western side (Fox and Franz Josef Glaciers) advanced during the study period, and (2) the development or enlargement of proglacial lakes on the large eastern glacier rapidly shrank the glaciers by lake calving. On the western side, Fox and Franz Josef Glaciers advanced 300–650 m in the time period studied, while La Perouse and Victoria Glaciers both retreated 1900 m (Fig. 6). On the eastern side all the largest glaciers retreated 1–2.1 km, most significantly Tasman Glacier (Tasman Glacier 2100 m, Murchison Glacier 1000 m, Hooker and Classen Glaciers 1700 m; Fig. 7). The enormous loss of ice mass due to the recent growth of these glaciers’ proglacial lakes contributes to such large area losses on the eastern side that the large retreat on most of the other valley glaciers on the western side of the divide (e.g. Victoria and La Perouse Glaciers (Fig. 6)) is small by comparison. If the eastern glaciers with proglacial lakes were discounted, the difference in relative decrease of glacier areas between the eastern and the western side of the Main Divide would be much smaller or even might be reversed.
In contrast, the smaller, high-elevation glaciers showed only minimal retreats compared with the large valley glaciers. Figure 8 illustrates the 1978 and 2002 glacier outlines for the five smallest glaciers studied in the field. It is not easy to detect significant glacier area changes. Some of these glaciers show small retreats, some have approximately the same size, and some even seem to have increased in some parts between 1978 and 2002. The latter observation is consistent with a trend to positive mass balances over the past two decades (Reference Chinn, Heydenrych and SalingerChinn and others, 2005).
5.1. Mapping accuracy of automated classification techniques
Mapping debris-covered glacier ice in a multispectral image is a recurring challenge because the spectral signature of the debris often prevails over the ice or snow signature if a certain percentage of the glacier pixels are debris-covered. Often confused with surrounding rocks, debris-covered glacier tongues are commonly manually added to classifications as a post-classification step (e.g. Reference Williams, Hall and ChienWilliams and others, 1997). Reference Paul, Huggel and KääbPaul and others (2004a) describe a semi-automatic method for delineating debris-covered glaciers which combines the advantages of automated multispectral classification for extracting clean glacier ice with slope information derived from a DEM. One important question raised by our research is how debris-covered glacier tongues can be mapped accurately from satellite imagery, even by using manual delineation, when it is nearly impossible to recognize their exact boundaries in the field. However, even in the worst cases misinterpretation in the field will not deviate more than a few pixels from the true boundary, as long as we assume that the ice beneath the debris is part of the glacier.
Snow in shadow tends to have reflectance properties similar to those of bright rock, and is another challenge for the automatic extraction of glaciers. Manual detection of these surfaces was relatively straightforward as long as the interpretation was careful, and topographic maps and other terrain information (e.g. DEMs and aerial photographs) were used. However, these areas were not always easy to extract with automatic methods, depending on the type of snow or ice (clean snow versus dirty ice) and steepness of the slope (snow in steeper slopes tended to be harder to extract). Most of the shaded clean snow could be accurately delineated using the lowest threshold tested for the ASTER 3/4 ratio. However, this threshold would also include large areas which obviously were not glaciers, so a trade-off was required to balance omission and commission errors. Many dirty snow/ice patches in shadow were left out. The NDSI was also shown to detect shaded clean snow relatively well, but not as well as the ASTER 3/4 ratio.
When we applied the NDSI, nearly all lakes in the area were classified as glaciers (Table 4). Turbid water, like the Tasman Lake, was still classified as glacier ice even with the highest threshold tested. The reason for this is the similar reflectance of very turbid lakes (with high load of rock flour in the water) and glacier ice in the green spectral domain. Increasing the NDSI threshold increased the omission error significantly (Table 4). However, this was not found to be a problem with the ASTER 3/4 ratio. In New Zealand the lakes are clearly visible in the VNIR composite image, and the water bodies were therefore manually digitized and removed from the NDSI results. After the water mask was applied to the NDSI results, this method revealed results similar to the ASTER 3/4 ratio (Fig. 4). In some areas the NDSI could be even more efficient. Similar problems of misclassifications of lakes as glaciers were reported by Reference Paul and KääbPaul and Kääb (2005) in Arctic Canada when using the Landsat Enhanced Thematic Mapper (ETM) 3/4 method. The best threshold for the NDSI in this study was ∼0.5 (Table 4), in accordance with other studies (e.g. Reference Mathieu, Chinn and FitzharrisMathieu and others, 2009). Band ratios have been reported to be an especially simple, robust and fast method for glacier mapping over large areas (e.g. in Switzerland (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002) and Peru (Reference AlbertAlbert, 2002)), and were confirmed as good performers in this study in New Zealand. The ASTER 3/4 ratio with a threshold around 2.0 was the most efficient (Table 2). This is also similar to other studies (e.g. Reference Paul, Kääb and HaeberliPaul and others, 2007).
5.2. Accuracy of glacier change detection between 1978 and 2002
The accuracy of glacier change detection between 1978 and 2002 was constrained by the difficulty of photo-interpreting glacier boundaries in both datasets. Figure 8 shows a comparison of the 1978 glacier outlines drawn from oblique aerial photographs and the manual digitizing of glaciers from the 2002 ASTER image for a selected number of glaciers, also delineated in the field. Drawing glacier outlines from oblique photographs is time-consuming and the accuracy depends on the local knowledge of the interpreter as well as the availability of good cartographic information for the area, like topographic maps and DEMs. This is apparent in Figure 8 where more details are evident in the digitizing of the satellite image compared with the 1978 interpretation, even though the ASTER image has a lower spatial resolution than the oblique photographs. The 1978 mapping was done at a coarser scale, mostly because details were practically impossible to transfer to the paper topographic maps. The 1978 digitizing was done on topographic base maps at a scale of 1 : 63 360 (imperial system where 1 in = 1 mile), where 1 mm corresponds to ∼63 m, a distance that includes four ASTER pixels. For the 1978 inventory we believe that in most areas the accuracy of the digitizing onto the base maps is within 1 mm. However, some corners are slightly cut off and some of the outlines are slightly generalized (e.g. Fig. 8). We conservatively estimated the accuracy of the 1978 glacier boundaries to be within 2 mm (±1 mm), giving us a ground-equivalent distance of ±63 m.
The accuracy of the 2002 manual digitizing firstly depends on the quality of the orthorectification. In the final orthorectified product, we achieved a root-mean-square error (RMSE) of 17.1 and 14.7 m for the CPs, in the X and Y directions respectively. In practical terms, these results are slightly more than a pixel (15 m), a standard level of accuracy when orthorectifying satellite imagery. The consistency of these results was verified for the entire image by visual inspection of the orthorectification overlaid by rivers and other vector information from the New Zealand topographic database. The other source of errors is related to the digitizing of the glacier boundaries. We estimated the uncertainties of the 2002 digitized glacier outlines for the most part to be within two pixels (30 m) (i.e. ±15 m). No corner was cut and every pixel was accounted for (Fig. 2). In areas of doubt (e.g. steep areas with glaciers in shadows), aerial photographs, 3-D model and field observations were used to support our interpretation. However, in order to accommodate difficulties of interpreting debris-covered tongues and the problems of interpreting steep mountain faces (e.g. east face of Mount Cook; Figs 9 and 10), we estimated the total uncertainty in the manual digitizing to be within 3 pixels (±22.5 m). By combining the uncertainty in the manual digitizing with the uncertainty of the orthorectification, we obtained an overall uncertainty estimation for the 2002 glacier outlines of ±28.5 m.
In order to calculate the resulting uncertainty in area, we took each single glacier polygon (for both inventories) and generalized them into circles (>600 digitized polygons within the study area in 2002, many containing several glaciers). This is a simplification, since the polygons certainly have more complex shapes; however, it should give us a relatively good error estimate. The uncertainties in area were calculated for each circle, both for the 1978 and the 2002 data. The uncertainty in the total glacier area change was ±6.8% (±6 km2) (Table 8).
Snow pixels show up as bright white in the visible imagery because of their high reflectance. However, whether this snow is part of a glacier or belongs to an occasional snowpatch is impossible to determine from the spectral signature itself. The interpretation is subjective, depending on what the interpreter sees as a glacier or not. In positive mass-balance years, with ample snow and a cold ablation season, many snowpatches will survive the summer and be extracted as glaciers from a satellite image by automated classification. If these areas were separated from the glacier, most could easily be removed by setting a minimum size for glacier polygons (e.g. 0.01 km2). However, if these snowpatches are connected to the glaciers, this will cause a larger area to be classified as glaciers compared with the true glacier area. To minimize these potential confusions, it is important to select an image from a negative mass-balance year.
In addition to the difficulty of distinguishing debris-covered glacier parts from rocks, it was also problematic to photo-interpret the upper extent of glaciers near steep rock faces. For the 1978 glacier inventory, the definition was that if these steep faces contributed ice by flow, they were part of the glacier, but if they were disconnected and only contributed by winter snow avalanches, the faces were not part of the glacier (Chinn, unpublished). However, this distinction was almost impossible to make only from the ASTER image because of the coarse resolution of the image relative to the features. Steep rock faces were not excluded from the total glacier area unless it clearly seemed to be bare rocks. The comparison of the 1978 and 2002 inventory maps shows many red areas around the highest peaks along the Main Divide (see example around Mount Cook in Fig. 9a), suggesting that glaciers would have grown in these areas during the study period, something which is not likely. Most of these red areas were classified as glaciers in 2002 because they were snow-covered at the time of image acquisition, whereas in reality they were just snow-covered rocks (cf. Fig. 9). The arrow in Figure 9 points to an area that was most likely wrongly classified as a glacier in the 1978 inventory. After closer investigation we found that this area is in fact a corridor of active rock avalanche. This is another illustration of the complexity of tracing glacier outlines from aerial photographs or satellite images. However, considering the two inventories as a whole, circumstances like active rock avalanches (leading to classifying part of glaciers as non-glaciers) and snow-covered mountain faces (leading to classifying rocks as glaciers) are minor and hardly influence the total outcome of the two inventories. When estimating the uncertainties in the two inventories, we have accounted for all possible error sources, and our estimates were based on conservative assessments. These estimates point to small uncertainties in the glacier area calculations and we can argue that the overall pattern and magnitude of glacier retreat would not change significantly after a further refinement of these outlines.
5.3. Change of glacier area in New Zealand from 1978 to 2002
The results show a measured retreat of the glaciers in the central Southern Alps of ∼17% during the period 1978–2002. Similar results have been reported in other parts of the world. The new Swiss glacier inventory (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2004b) measured the relative loss in glacier size from 1973 to 1998/99 as ∼20%. However, in Switzerland, glaciers smaller than 1 km2 (100 ha) contributed to ∼40% of the total loss, while they cover only 15% of the area. About 18% of this change occurred from 1985 to 1999, illustrating a recent acceleration of the glacier retreat. This strongly contrasts with the New Zealand situation where the largest glaciers have contributed to the largest ice loss. Reference Quincey and GlasserQuincey and Glasser (2009) reported a retreat of Tasman Glacier by 3.5 km between 1990 and 2007. They also found that the area of the proglacial lake doubled between 2000 and 2007 alone. In Norway, the maritime glaciers along the western coast with a large annual mass turnover had a large mass surplus between 1962 and 2000. Conversely, the continental glaciers with small summer and winter balances had a mass deficit over the same period (Reference Andreassen, Elverøy, Kjøllmoen, Engeset and HaakensenAndreassen and others, 2005). Since 2001, all monitored glaciers in Norway have shown a considerable mass deficit. Overall the Norwegian glaciers have retreated during the 20th century. Continental glaciers have, with few exceptions, retreated throughout the century, while many maritime glaciers have been through periods of advance and recession, although recession has been the main process (Reference Andreassen, Elverøy, Kjøllmoen, Engeset and HaakensenAndreassen and others, 2005). Reference Andreassen, Paul, Kääb and HausbergAndreassen and others (2008) reported glacier area shrinkage at investigated glaciers in Jotunheimen National Park, Norway, by 23% since the 1930s for 38 glaciers. The largest reduction occurred in the first part of the period, between 1930 and 1965. Since 1965 the area reduction has been 12% for 164 glaciers. Reference Paul and KääbPaul and Kääb (2005) calculated a decrease in area for 225 glaciers in Cumberland Peninsula, Baffin Island, Arctic Canada. They compared Landsat ETM+ and ASTER scenes from August 2000 with manual delineation of LIA trimlines and moraines. They found an average area loss of 11% over the period. The relative loss of area increased towards smaller glaciers (<10 km2) and was more or less constant for larger glaciers. Reference Khromova, Osipova, Tsvetkov, Dyurgerov and BarryKhromova and others (2006) assessed glacier recession in the east Pamirs, central Asia, using historical data and space images from Russian satellites and ASTER. They found a 7.8% decrease in glacier area during 1978–90, and 11.6% during 1990–2001. In addition to the decrease in glacier area and retreat of glacier fronts, they also found increased debris-covered area and the development of new lakes. All these reports of glacier retreat are similar to what is found in this study, but the reasons behind the recessions might be very different.
Qualitative analysis of our results suggests that the smallest glaciers in New Zealand do not contribute to the largest ice loss (Figs 6–8). The large heavily debris-covered valley glaciers have lost enormous areas over the past few decades, whereas some of the smaller glaciers may even have shown a small advance. This last observation was corroborated by Reference Allen, Owens and SirgueyAllen and others (2008) who reported that Stocking and Eugenie Glaciers (east of the Main Divide) advanced up to 100–140 m between 2002 and 2006.
The clustering of these large valley glaciers in the central Southern Alps covered by the 2002 ASTER image indicates that the area change for the study period is likely to be larger here than for the rest of the glacierized areas in New Zealand. This suggests that if the 2002 glacier inventory had been made for the whole of New Zealand, the total area changes from 1978 would probably have been lower. The response time for alpine glaciers is ∼30 years (Reference Cuffey and PatersonCuffey and Paterson, 2010), while response times for high-velocity glaciers like Franz, Fox and Stocking are in the order of 5–8 years (Reference Fitzharris, Lawson and OwensFitzharris and others, 1999). Some of the large valley glaciers maintained their LIA areas (1850–90) up until the 1970s, giving a response time of ∼100 years. The relatively thick ice of the lower reaches of these glaciers is well insulated beneath thick protective mantles of debris and became largely relict. Ice loss was shown by downwasting (surface lowering) alone, with little to no terminus retreat. These glaciers are in a state of extreme disequilibrium with the present climate. However, during the past few decades, many have reached irreversible tipping points that have triggered proglacial lake expansion into terminus areas (Reference Quincey and GlasserQuincey and Glasser, 2009). These lakes, with their rapid expansion by both melt and glacier calving, are destroying the lower trunks of the large valley glaciers, making the glacier behaviour not directly linked to the climate (Reference Hochstein, Watson, Malengrenau, Nobes and OwensHochstein and others, 1998; Reference Kirkbride and WarrenKirkbride and Warren, 1999; Reference Purdie and FitzharrisPurdie and Fitzharris, 1999; Reference RöhlRöhl, 2006; Reference Winkler, Chinn, Gärtner-Roer, Nussbaumer, Zemp and ZumbühlWinkler and others, 2010). Reference ChinnChinn (1999) reported a reversal of the past century glacier recession trend. Recent interannual fluctuations of ELAs in New Zealand indicate that alternate positive and negative mass-balance years (1998–2008) have followed two decades of positive mass balance (1977–98) (Reference Willsman, Salinger and ChinnWillsman and others, 2008). Figures 6 and 7 show the changes in glacier area for the western and eastern glaciers in the 24 year period. Only the two largest of the west coast glaciers have advanced, whereas all the other large glaciers have retreated. It may be inferred that the western climate is more favourable for glacier advance, compared with east of the Main Divide. However, Reference ChinnChinn (2001) argues that the notion that the ‘western’ glaciers behave differently from the ‘eastern’ glaciers, which is frequently mentioned in the literature, is unfounded. This false impression arises from the many simplistic comparisons made between only two groups of glaciers that are of different types and have very different responses to climate: ‘Frontal changes at the steep and exceptionally responsive Franz Josef and Fox Glaciers, which provide enhanced responses to climate within 5 to 8 years (Fizharris and others, 1999), have repeatedly been compared with the large, low-gradient valley glaciers of the Mount Cook area, which have slow, subdued response times in the order of 100 years’ (Reference ChinnChinn, 2001). Fox and Franz Josef Glaciers have very steep mass-balance gradients and enormous accumulation basins compared with the other valley glaciers presented in Figures 6 and 7. An illustration of their enormous accumulation basins is given in Figure 11 and a comparison is made with the other main valley glaciers west of the Main Divide. The two glaciers with the largest accumulation basins (Fox and Franz Josef Glaciers) have advanced, those with reasonably large basins (e.g. Spencer Glacier, Fritz Glacier, Salisbury Snowfield) seem to be stable, while the valley glaciers with the smallest accumulation basins (e.g. La Perouse and Victoria Glaciers) have shown a marked retreat. It is important to highlight in this context that the hypsometry (i.e. the distribution of glacier area over its altitudinal range) for Fox and Franz Josef Glaciers is also very different from that of the other western valley glaciers. The main part of the accumulation basins for these two glaciers is located above 2000 m a.s.l., while only very small sections of the other valley glaciers are above this elevation. This would make Fox and Franz Josef Glaciers very responsive to small changes in climate, much more so than the other western valley glaciers. Because of their fast reaction time, scientists have tried to correlate the retreat and advance of these two glaciers with climatic parameters. Reference Hooker and FitzharrisHooker and Fitzharris (1999) found a strong link between atmospheric circulation changes, climate variables and the behaviour of Fox and Franz Josef Glaciers. However, causality effects between climate and glacier behaviour are complex, and hypsometry is only one factor in this interplay (Reference Dyurgerov, Meier and BahrDyurgerov and others, 2009). The processes leading to glacier changes were not the subject of this study, but the reader can consult Reference Chinn, Heydenrych and SalingerChinn and others (2005) for a comprehensive summary of annual mass balances in New Zealand and their possible causes during the past few decades.
6. Conclusion and Perspectives
This study has stressed the large potential of glacier extraction from ASTER imagery, but at the same time underlined the problems involved. Thresholding of the ASTER 3/4 ratio proved to be the most efficient automatic classification method, and when combined with the manual extraction of debris-covered tongues it yielded reasonably accurate results over a large area. The optimum threshold was ∼2.0, which is in agreement with other studies (e.g. Reference Paul, Kääb and HaeberliPaul and others, 2007). For the NDSI, there was a large misclassification of water bodies, especially for turbid water. After manually masking out all the water bodies, the NDSI gave results very similar to those of the ASTER 3/4 method. The optimum threshold for the NDSI was ∼0.5, similar to values reported in other studies (e.g. Reference Mathieu, Chinn and FitzharrisMathieu and others, 2009). The initial test of supervised classification did not provide very promising results. Shaded vegetation was classified as glacier ice, and shaded snow was poorly detected by the supervised classification. The method for semi-automatic extraction of glacier debris described by Reference Paul, Huggel and KääbPaul and others (2004a) should be tested for the New Zealand environment, as debris was found to be one of the most significant limitations for remote-sensing based automatic glacier mapping in the study area.
The glaciers in the study area show an overall area loss of ∼17% during the 24 year period 1978–2002. Glaciers west of the Main Divide of the Southern Alps lost ∼14% of their 1978 area, whereas glaciers east of the Main Divide lost ∼18%. The overall loss in the study area is similar to that previously reported from other parts of the world. Despite the large climatological precipitation difference between glaciers west and east of the Main Divide, they show relatively similar area losses. Large glaciers show the greatest area losses, whereas some of the small, alpine, high-altitude glaciers show no change or in some cases signs of advance, which is in agreement with the reported positive trend in mass balance from the New Zealand annual snowline survey over the past three decades (Reference Chinn, Heydenrych and SalingerChinn and others, 2005). It is important to state that if all the large valley glaciers that have developed proglacial lakes during recent decades were removed from the analysis, the changes in glacier area between 1978 and 2002 would have been minimal. As underlined by Reference Khromova, Osipova, Tsvetkov, Dyurgerov and BarryKhromova and others (2006), the importance of in situ studies and historical glaciological and climatic data in order to best assess glacier changes from satellite data cannot be underestimated. Thus, more fieldwork should be undertaken on the New Zealand glaciers to obtain a better understanding of the current processes, for instance the rapid retreat of glaciers with proglacial lakes.
This work was carried out in collaboration between the Department of Geoscience, University of Oslo, and the School of Surveying, University of Otago. We are grateful to the School of Surveying for the use of their facilities and field support during the stays of the lead author in New Zealand. Some field support was provided by the Department of Geoscience. We thank P. Sinden for field assistance. We are grateful to Glacier Explorers in Mount Cook Village for boat transportation on Tasman Lake and for acting as back-up security during the fieldwork. We thank The University Centre in Svalbard (UNIS) and A. Hormes for providing support for finalizing this work. Finally, we thank B. Raup, an anonymous reviewer and scientific editor T. Scambos for constructive, thorough comments and suggestions that appreciably improved the manuscript.