Next Article in Journal
Landsat-8 Sea Ice Classification Using Deep Neural Networks
Next Article in Special Issue
Long-Term Landsat-Based Monthly Burned Area Dataset for the Brazilian Biomes Using Deep Learning
Previous Article in Journal
Hyperspectral Anomaly Detection Using Deep Learning: A Review
Previous Article in Special Issue
Assessing the Wall-to-Wall Spatial and Qualitative Dynamics of the Brazilian Pasturelands 2010–2018, Based on the Analysis of the Landsat Data Archive
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Mapping Three Decades of Changes in the Tropical Andean Glaciers Using Landsat Data Processed in the Earth Engine

1
Instituto del Bien Común (IBC), Lima 15072, Peru
2
Programa de Doctorado en Recursos Hídricos (PDRH), Universidad Nacional Agraria La Molina, Lima 15024, Peru
3
Fundación EcoCiencia, Quito 170517, Ecuador
4
Fundación Amigos de la Naturaleza (FAN), Santa Cruz 2241, Bolivia
5
Earth Observation and Geoinformatics Division, National Institute for Space Research, São José dos Campos 12227-010, Brazil
6
Subdirección de Información y Análisis (SDIA), Instituto de Investigación en Glaciares y Ecosistemas de Montaña (INAIGEM), Ancash 02001, Peru
7
Instituto do Homem e Meio Ambiente da Amazônia (Imazon), Belem 66055-200, Brazil
*
Author to whom correspondence should be addressed.
Remote Sens. 2022, 14(9), 1974; https://doi.org/10.3390/rs14091974
Submission received: 16 March 2022 / Revised: 11 April 2022 / Accepted: 19 April 2022 / Published: 20 April 2022
(This article belongs to the Special Issue State-of-the-Art Remote Sensing in South America)

Abstract

:
The fast retreat of the tropical Andean glaciers (TAGs) is considered an important indicator of climate change impact on the tropics, since the TAGs provide resources to highly vulnerable mountain populations. This study aims to reconstruct the glacier coverage of the TAGs, using Landsat time-series images from 1985 to 2020, by digitally processing and classifying satellite images in the Google Earth Engine platform. We used annual reductions of the Normalized Difference Snow Index (NDSI) and spectral bands to capture the pixels with minimum snow cover. We also implemented temporal and spatial filters to have comparable maps at a multitemporal level and reduce noise and temporal inconsistencies. The results of the multitemporal analysis of this study confirm the recent and dramatic recession of the TAGs in the last three decades, in base to physical and statistical significance. The TAGs reduced from 2429.38 km2 to 1409.11 km2 between 1990 and 2020, representing a loss of 42% of the total glacier area. In addition, the time-series analysis showed more significant losses at altitudes below 5000 masl, and differentiated changes by slope, latitude, and longitude. We found a more significant percentage loss of glacier areas in countries with less coverage. The multiannual validation showed accuracy values of 92.81%, 96.32%, 90.32%, 97.56%, and 88.54% for the metrics F1 score, accuracy, kappa, precision, and recall, respectively. The results are an essential contribution to understanding the TAGs and guiding policies to mitigate climate change and the potential negative impact of freshwater shortage on the inhabitants and food production in the Andean region.

Graphical Abstract

1. Introduction

Glacier meltwater plays a vital role in providing water resources to populations downstream of tropical glacier mountain ranges. However, glaciers are decreasing in extent and mass during the dry season at an accelerated rate [1,2,3], due to the increase in global temperature caused by anthropogenic and natural causes [4,5,6]. As a result, glacier retreat is an essential indicator of climate change [2] and climate variability [7,8]. The tropical Andean glaciers (TAGs) are undergoing a rapid coverage reduction, which can economically impact populations in the tropical Andes, with effects on agriculture, drinking water, electricity generation, ecosystem integrity [9,10,11,12], cultural assets, and tourism activities [13,14,15], due to their cultural and tourism value. In addition, the accelerated retreat of glaciers can generate an increase in glacial lake outburst flood (GLOF) risks [7,16].
Most high mountains with glaciers have experienced glacier shrinkage [17]. Understanding the current status of glacier retreat is vital for better management of water resources [18] and adaptation to climate change. Many tropical glaciers, such as those in Peru and Bolivia, are critical buffers against reduced precipitation during the dry season. Ninety-nine percent of tropical glaciers are located in the Andes of South America, extending over Venezuela, Colombia, Ecuador, Peru, Bolivia, Chile, and Argentina [19].
The direct mapping of glaciers is difficult and dangerous due to the remoteness and inaccessibility of the terrain and challenges for conducting extensive fieldwork [20]. The application of remote sensing in glaciers began in the 1970s, initiating a new phase in glacier monitoring and solving to some extent the limitations and difficulties of glacier mapping [19]. The data obtained by remote sensors provide valuable, constantly updated, and almost real-time information on glaciers and associated geographic features [19,20,21].
At global and regional scales, the extent of glaciers can be derived through automated classification algorithms applied to multispectral satellite data [20], such as those from Landsat missions with more than three decades of freely available global information (http://landsat.usgs.gov/ (accessed on 1 March 2021)). Automated methods for glacier classification from multispectral data have been developed and tested in the last decade, and have proven to be suitable for detecting clean or slightly dirty glacier ice and fresh snow [20,22,23]. However, the efforts of mapping glacier coverage do not have a unique record of glacier area. They contain information from different dates on each glacier [24], making their temporal comparison difficult, as in the Global Land Ice Measurements from Space (GLIMS) and World Glacier Monitoring Service (WGMS) initiatives [19,20]. Therefore, it is necessary to have a comparable temporal record of glacier areas to observe their response to climate change and variability [21].
In recent years, the availability, accessibility, and temporal and spatial continuity of satellite data (such as Landsat and Sentinel) [25,26,27,28] and the development of cloud computing platforms, such as Google Earth Engine (GEE) [29], have created an opportunity to map the historical evolution of land cover change, such as water bodies and glaciers [16,18,21,26,30,31,32,33,34,35,36]. Despite the progress made in recent years in the assessment of glacier areas, there are still many cases in which the automatic delimitation of glacier boundaries in satellite images is complex, error prone, or time consuming [13,20,37]. In most cases, glacier cover change estimates cover 5- or 10-year periods [10,35,38].
This research estimates the annual extensions of the entire tropical Andean glaciers in the last three decades, based on a semi-automated classification method and the complete Landsat satellite images historical archives using the GEE platform. Unlike previous studies, this work’s spatial scope covers the entire tropical Andes, temporally comprising a 36-year time series by using an innovative methodological approach to map exposed glaciers based on cloud computing technologies to retrieve and process large quantities of data and the application of post-classification filters. Based on its yearly results, a detailed analysis of all the annual changes of the TAGs is presented, considering site morphometric characteristics (altitude, aspect, and slope) and geographic location.

2. Materials and Methods

2.1. Study Area

The study area refers to the tropical Andean glaciers (TAGs), geographically located between the Tropic of Cancer and the Tropic of Capricorn (between latitudes 23°26′13.3″N and 23°26′13.3″S) within the Intertropical Convergence Zone (ITCZ) [2,19]. According to the Randolph Glaciers Inventory (RGI), the TAGs have a surface area of 2337.55 km2 representing approximately 99% of all tropical glaciers [19,39]. According to the classification of glaciers from a climatic approach [40,41], the TAGs contain four climatic zones: Northern Wet Outer Tropics, Southern Wet Outer Tropics, Dry Outer Tropics, and Inner Tropics (Figure 1).
For the delimitation of the TAGs study area, we used the RGI (Version 6.0) [39], records to which we applied a catchment area or buffer of 1.5 km around each glacier to include the surrounding areas. We also visually identified and included glacier areas that were not mapped by the RGI in southern Peru, western Bolivia, and northern Chile, totaling an area addition of approximately 15%. The scope of our research precisely concerns the clean or marginally contaminated glacier surfaces. According to Herreid and Pellicciotti (2020) [42], debris-covered glaciers have a very low representativeness, accounting for approximately 5.4% of the total tropical glacial areas. Hence, the debris-covered glacial areas have not been dealt with in this re-search, for they also require specific methods for their mapping, as reported by Shukla et al. (2017) [43]. Finally, we ought to mention as well that glacier surges are not approached in this work. Contrary to other Andean glaciers [44], where these events are more commonly observed, glacier gains in the TAGs, usually resulting from surges, are sporadic [8] and occur to a minor extent.

2.2. Landsat Data

As the present study is an analysis of glacier coverage change at a multiannual level, the following data corresponding to the period 1985–2020 [26] were obtained from the image collections of the Landsat sensors: Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM +), and Operational Land Imager (OLI), on board the Landsat 5, Landsat 7, and Landsat 8 satellites, respectively [45]. We processed Landsat Collection 1 Level 1 and Tier 1 surface reflectance (SR) products [46], already orthorectified and with a spatial resolution of 30 m in the multispectral bands. The Landsat image collections used were accessed and processed in the Google Earth Engine (GEE) platform [29] (https://earthengine.google.com/ (accessed on 20 June 2021)).

2.3. Landsat Data Processing Workflow

Figure 2 summarizes the methodological workflow of our approach for the TAGs’ change detection. It comprises six stages: selection and masking of clouds and cloud shadows in the satellite images, calculation of the Normalized Difference Snow Index (NDSI), creation of annual Landsat mosaics, annual glacier classification, post-classification, and finally, accuracy assessment.

2.4. Images Selection and Cloud Mask

To select multiple images for each analyzed year, we carried out a sequence of procedures proposed by [20,47]. All Landsat images with cloud cover equal to or less than 70% acquired during the three-decade analyzed period were filtered from the Landsat archive and selected (Table S1 and Figure S1). After selecting the images based on the procedures above, we masked the clouds and cloud shadows. For cloud masking, we used the CloudScore technique [48,49]; in the case of cloud shadow, we used the Temporal Dark Outlier Mask (TDOM) algorithm [49,50] and the Band Quality Assessment (BQA) information available in the Landsat Collection. The BQA was not used for cloud masking, since clouds are overestimated in the glacier edge areas. We implemented these techniques in the GEE platform, which have been thoroughly described by the MapBiomas Project [51].

2.5. Normalized Difference Snow Index (NDSI)

The Normalized Difference Snow Index (NDSI) [52] is a widely used index to separate snow from other coverages [49], and therefore it is also useful to identify glacier coverage. We applied the NDSI index to each selected and masked image in the previous steps, based on Equation (1).
NDSI = ( Green SWIR Green + SWIR )

2.6. The Compositing and Mosaicking of the Annual Images

We created annual image composites for each year of the time series (36 years) using GEE’s image-reducing techniques, having as an input all the cloud-masked Landsat scenes available in GEE’s Landsat collection with less than 70% cloud cover. We applied the image reducers to the NDSI, RED, and NIR bands, as these variables are most suitable for glacier classification [53]. The NDSI values in each pixel of each year were reduced to a band of N D S I m i n (annual minimum NDSI), and the annual reduction in the RED and NIR bands were based on the 25th NDSI percentile. This percentile value was defined on a trial and error basis, aiming to reduce cloud cover and eliminate noise at the same time. In brief, the time-series data dimensionality was reduced by being condensed into a three-band artificial image for each year with the minimum snow cover, composed of manifold pixels from different dates throughout each of the analyzed years [51]. These three bands are the N D S I m i n , R E D m i n s n o w , and N I R m i n s n o w , and the two latter ones were calculated using Equations (2) and (3), respectively.
R E D m i n s n o w = median ( RED   with   the   25 th   NDSI   percentile )
N I R m i n s n o w = median ( NIR   with   the   25 th   NDSI   percentile )
The R E D m i n s n o w reduction is interpreted as the RED band with minimum snow in the year, and the N I R m i n s n o w as the NIR band with minimum snow. We used these criteria to reduce glacier confusion with temporary snow [20]. A similar work about exploiting a time-series stack of images and its corresponding reduction was reported by Winsvold et al. (2016) [54].

2.7. Annual Glacier Classification

Annual glacier classification for the 1985–2020 period was obtained using an empirical decision tree (EDT) constructed using samples of N D S I m i n , R E D m i n s n o w , and N I R m i n s n o w variables. The sampling strategy was performed for both Landsat 5 and 7 in 2000, and for Landsat 8 in 2020, respectively, and a visual sample collection was performed for the classes: glacier, water, and other coverages. Samples were collected in the entire study area considering the classes spectral variation (e.g., pixels in different slopes) and prioritizing pixels where there is greater confusion among the three classes (e.g., shaded pixels). In total, 9800 samples were collected for each sensor. Then, we conducted a statistical analysis of the sample set for each year using a boxplot violin diagram (Figure 3 and Figure 4) in order to find the optimal thresholds.
The decision tree used is shown in Figure 5. We empirically defined the decision tree nodes (1–1, 2–2, and 3–2) with thresholds α, β, and γ defined for each TAG classification region, classification year, and sensor type—TM, ETM+, and OLI (see Figure S2 for the classification regions details). EDT nodes were constructed with the objective of detecting glacier ( variables   N D S I m i n ) and eliminating confusion with water bodies and other coverages ( R E D m i n s n o w and N I R m i n s n o w   variables) [53].
The average value and standard deviation of the threshold α are 0.20 and 0.082, respectively, while the β and γ thresholds were set as constant at 0.11 and 0.10, respectively, which are very similar values to those used by Macander et al. (2015) [53]. The thresholds were adjusted using visual inspection in each classification region and year.

2.8. Post-Classification

To remove noise and information gaps due to cloud persistence, we designed and applied the following sequence of five post-classification filters:
  • Gap-fill filter to replace information gaps due to clouds with consistent temporal classification values;
  • Temporal filter to identify spurious class transitions (e.g., change in short periods from glacial to non-glacial);
  • Frequency filter to identify stable pixels [51];
  • Temporal permanence filter;
  • Spatial filter to remove salt-and-pepper classification noise without changing the boundary of the classified patches.
These filters were applied to reconstruct land cover and land use classes over large time-series, as in this study [51]. Below, we provide the details of each filter’s utility.
The filter sequence begins with the gap-fill technique (Figure 6). Severely cloud-affected regions commonly present information gaps in the annual TAGs classifications in long time-series, as a result of the cloud masked regions. The gap-fill tool reduces this effect by replacing ‘no data’ pixels with the temporally closest class value starting from the near past and moving forward in time (example Figure 6). The use of this filter allows for the achievement of comparable glacier areas in the time series [46].
Next, a temporal filter was applied [51]. The temporal filter uses sequential rankings in a one-way moving window of 3, 4, or 5 years’ kernel sizes to identify the temporally inconsistent transitions (e.g., glacier/non-glacier expansion or contraction when the historical record shows the continuity of either class) based on the intermediate-, first-, and last-year rules (Figure 7). The first- and last year rules were initially applied to respectively eliminate time inconsistencies in the years 1985 and 2020 (Figure 7a,b), followed by the intermediate years’ rules (i.e., 1986 to 2019) in the following order: the three-year rule (Figure 7c), the four-year rule (Figure 7d), and finally the five-year rule (Figure 7e).
Next, the frequency filter was applied. This filter improves the classification according to the historical occurrence of the class over the analyzed 36-year period. For this purpose, it considers the frequency percentage of the glacier class in the time series, from which it updates the ‘non glacier’ class occurrences, which can be safely considered classification errors in most cases. This study used a 70% frequency threshold; in other words, occurrences smaller than or equal to 30% of ‘non-glacier’ values in the time series were replaced by glacier values. This filter was applied only to the first and intermediate years; years at the end of the time series were excluded to prevent the loss of information of transitions occurring at the end of the time series (e.g., glacier loss in recent years) (Figure 8).
Based on the GLIMS definition, glaciers are defined as permanent ice or snow-mass [20]. We proposed applying a temporal permanence filter to separate temporary snow from glacier coverage (Figure 9).
Finally, to eliminate salt-and-pepper classification errors [37], a spatial filter based on GEE’s ‘connectedPixelCount’ function was applied. The spatial filter identifies sets of pixels that share the same pixel value based on their neighborhood. Consequently, only pixels that do not meet a predefined connection criterion (i.e., minimum number of identical-value connected pixels) are defined as isolated pixels and reclassified [51]. We considered a minimum of 5 grouped pixels for glaciers, representing a minimum area of approximately 0.5 ha.

2.9. Accuracy Assessment

Thematic accuracy analysis is the primary way to evaluate the quality of maps [55]. The accuracy assessment of annual maps in this study used a random sampling of 7299 points. This sample size guarantees a confidence level of 95%, with a normal distribution (Z = 1.96), a level of significance of 5%, and a maximum tolerable error below 1% (e = 0.5%). A minimum of 200 sample points per class and country was assured. Each point was interpreted and assigned a class using the QGIS Earth Engine plugin (qgis-earthengine-plugin) [56] and the thematic map accuracy evaluation plugin AcATaMa [57]. The reference dataset was obtained by means of visual inspection of the sample points, conducted by independent expert image interpreters, who relied on an RGB image composite consisting of the following bands: SWIR with the 25th NDSI percentile, NIR with the 25th NDSI percentile, and RED with the 25th NDSI percentile, for each analyzed year. This allowed us to calculate the error matrix and other accuracy metrics. We also conducted accuracy analyses according to the work of Olofsson et al. (2014), on good practices for estimating area and assessing accuracy of land change [58].
The accuracy of the classifications was tested using a set of binary metrics. The possible outcomes are true positive (TP), true negative (TN), false positive (FP), and false-negative (FN). The metrics tested are given in the following Equations (4)–(7) [59,60]:
Precision :   TP TP + FP  
Recall :   TP TP + FN
Accuracy :   TP + TN TP + TN + FP +   FN
F 1 score = 2 Precision Recall Precision + Recall = 2 TP 2 TP + FP + FN

3. Results

3.1. The Area Statistics of the TAGs

We estimated the annual area of the TAGs from 1985 and 2020, comprising 36 years of annual mapping. This annual TAGs time series is the largest one currently available, to the best of our knowledge. In order to have a good quality base year, we analyzed the results from the year 1990 onwards, using ten-year intervals until 2020 (Table 1). The smallest TAGs (minimum) for the period 2011–2020 had a total surface of 1409.12 km2 (Table 1).
We analyzed the results of the TAGs coverage at the country level using the Global Administrative Unit Layers (GAUL) country boundary available in GEE (https://developers.google.com/earth-engine/datasets/catalog/FAO_GAUL_2015_level2 (accessed on 14 January 2021)). We found that the most extensive glacier area in 2020 was found in Peru and Bolivia, with 72.76% and 20.35% of the total TAGs surface, respectively. Ecuador, Colombia, Chile, and Argentina combined represented 6.89% of the TAGs coverage in the Andes Cordillera (3.89%, 2.18%, 0.78%, and 0.04%, respectively). Venezuela has less TAGs coverage, with a percentage of less than 0.01%, accounting for nearly 0.03 km2 (Table 2 and Figure 10a).
The TAGs occur at different altitudes (Figure 10b). We verified that, in 2020, most of the glaciers are distributed between 4801 and 5600 masl, representing 82.63% of the total TAGs, while approximately only 1.72% at altitudes less than or equal to 4800 masl, and 15.64% are found above 5600 masl.
The results show that, in 2020, 19.47% of the TAGs face north, 21.56% east, 31.30% south, and 27.67% face west (Figure 10c). Most TAGs (73.31%) occur in areas with slope ranges from 11° to 40° (Figure 10d). Flat and low slopes terrains (i.e., from 0° to 10°) contain 11.57% of the TAGs, while areas with pronounced slopes (above 41°) concentrate 15.11% of the TAGs (Figure 10d).
The results show that, in the year 2020, TAGs are located between latitudes −24° and 11° (Figure 11a). The largest glacier areas occur in the range between −14° to −10°, which contains 44.75% of the total TAGs; the Peruvian Cordilleras (mountain ranges), Blanca, Vilcanota, Vilcabamba, and Urubamba, are located within these latitudes. The longitude range of the TAGs varies from −78° to −67° (Figure 11b).

3.2. The Spatial and Temporal Patterns of the TAGs

The total area of the TAGs considerably reduced between 1990 and 2020, following a linear decline pattern. We found an average glacier area reduction of 28.42 km2/year with an R2 of 0.98. In 1990, the mapped area of the TAGs reached 2429.38 km2, reducing by 42% by 2020 to an area of 1409.11 km2 (Figure 12a). Of the last three decades, the decade of 2011–2020 witnessed the greatest loss of glacier areas (associated with the decade’s minimum surface value) (Figure 12b). Figure 13 presents a map for the Vilcanota Cordillera (mountain range) in Peru. The retreat buffer width varies along the TAGs edges and over the years. The smallest TAGs (i.e., those equal to or less than 10 ha) lost more than 90% of their area in the period from 1990 to 2020 (Figure 13).
For the analysis of glacier loss or retreat, 1990 was adopted as the base year, and annual maps of losses were created (Table S2). The area losses of the TAGs by country are shown in percentages in Figure 14, where, in the year 2020, Peru, Bolivia, Ecuador, Colombia, Chile, Argentina, and Venezuela presented retreats of 41.19%, 42.61%, 36.37%, 60.19%, 47.24%, 45.47%, and 96.93%, respectively, with respect to the base year of 1990.
There are three distinct patterns of TAGs retreat when area loss is disaggregated by country (Figure 14, Table S3). First, Venezuela, Argentina, and Chile showed a rapid decline in the TAGs between 1990 and 1995, with approximately 20% of glacier area loss during that specific period (Figure 14), although this can be due to the low satellite image quality in those initial years. By 2020, Venezuela’s TAGs losses exceeded 90%, but in the case of Chile and Argentina, the losses reached 45.47% and 47.24%, respectively, since the year 1990. A second pattern is found in Colombia’s TAGs, which showed an acceleration from 1995 onwards, reaching a loss of 60.19% by 2020. The third pattern of gradual linear TAGs loss occurred in Peru, Bolivia, and Ecuador (Figure 14). Ecuador showed a lower rate of retreat.
We analyzed the glacier area loss by different latitudes and longitudes, with respect to the average loss of all TAGs (Tables S4 and S5). We found that the latitudes farther south (−24° to −20°) have higher losses with regard to the total mean loss, with more than 62.54% of the glacier loss observed during the analyzed time series (Figure S3), however these glaciers represent less than 1% of the total area of the TAGs in 2020. The TAGs in the southern latitudes, from −9° to 0°, had a reduction of 31.61%, lower than the total average (42.00%); these latitudes account for 33.6% of the total area of the TAGs in 2020. At the same time, we verified that glaciers with latitudes ranging from 1° to 11° north had an accelerated glacier retreat in recent years (since 1996), with a 54.27% loss from 1990 to 2020 (Table S4). We also observed differentiated change patterns when analyzing glaciers at different geographic longitudes. The glaciers between longitudes −74° and −73° experienced changes of 59.94%, well above the total average (i.e., mean total retreat TAGs over the period 1990–2020), while longitudes −78° and −77° had losses of 31.68%, below the total average (Table S5 and Figure S3).
The area loss of the TAGs with regard to aspect (i.e., different glacier slope orientations) was lower in the south and east (i.e., 41.59% and 40.63%, respectively) and higher in the north and west orientations (42.86%) (Table S6 and Figure S4). Our analyses also estimated that glaciers with a larger slope than 80° rapidly reduced by more than 30% in 1996, with regard to the 1990 TAGs state. The retreat then slowly decelerated in the following years, until it reached a total loss of 41.56% by 2020. On the other hand, the TAGs with slopes ranging from 61° to 70° had a greater acceleration of loss from 2015 onwards, representing a loss of 46.37% by 2020. Nevertheless, in other cases, we could not find a slope range that controls the retreat of the TAGs area (Table S7 and Figure S5). In summary, we also observed that the TAGs areas below 4600 masl lost an average of 92.41% of their extension in 1990, while at altitudes above 6000 masl, the losses have been lower, 9.73% on average (Table S8 and Figure S6).
Due to the importance of the Amazon basin in South America, we verified that the glaciers in this region have an approximate area of 869.59 km2 in 2020, representing 61.71% of the TAGs total area. We calculated and compared the evolution of the TAGs area loss in the period from 1990 to 2020 in the Amazon basin, and concluded that it has lost 43.07% of its glacier area, unlike other non-Amazonian basins, which have lost 40.18% along the analyzed time series. We found a more significant acceleration in recent years, approximately from 1995 onwards, where the loss of the Amazon basin exceeds that of other basins (Table S9 and Figure S7).

3.3. Accuracy of the TAGs Mapping

We conducted the accuracy assessment of the TAGs annual maps over the entire time series of the reference dataset generated for the study area based on 7299 random points and visual inspection, as previously described. The average annual accuracy of the original classification was 95.29% (Table S10). After the application of the filters, we observed a gradual improvement in map accuracy to 95.30%, 96.31%, 96.37%, 96.27%, and 96.32% for the classifications relying on the cumulative sequential use of the gap-fill, temporal, frequency, temporal permanence, and spatial filters, respectively (Figure 10, Table S11). Although there was a slight decrease in the average accuracy, the standard deviation of the time-series classification accuracy was reduced from 1.24 to 0.44 with the filtering techniques, hence exemplifying the utility of the post-classification filtering routines. The marginal gain of the filtering techniques (i.e., 1.03 ± 0.44) is justified with the better spatial and temporal consistency in the classification of glacier coverage for multiannual comparisons. The average classification accuracy of the final filtered product attained 92.81%, 96.32%, 90.32%, 97.56%, and 88.54% for the validation metrics F1 score, accuracy, kappa, precision, and recall, respectively. The kappa index and the accuracy showed the highest values, while recall was the lowest (Figure S8).
This indicates an omission error of 11.46% and a commission error of 2.44%. The omission errors mostly occurred in shaded areas and in small glaciers, while the commission errors took place at the glaciers’ borders. The highest omission and commission errors were found in the first and last years (1985 and 2020), caused by the low amount of Landsat images in the early years (Figure S3), since the filters still have limitations in the early and late years (Figure 15a). The accuracy reduction valleys in given moments throughout the time series (Figure 15a) are related to the mosaic quality in those years (e.g., 2002), due to the limited availability of cloud-free images, which had an important effect on the accuracy of those years’. The first five years of the time series were also impacted by the limited image availability.
The accuracy analyses considering the area proportions according to [58] are presented in Table S12, which also shows the annual accuracy metrics disregarding these proportions. The difference between these two types of mean accuracy is only 0.42%. This minor discrepancy is probably due to the fact that we are dealing with a binary classification.

4. Discussion

The thematic accuracy of the obtained results present acceptable values for the multitemporal maps, with an average accuracy of 96.32%. Similar studies in the literature reported inferior accuracies in general. Albert (2002) [22] conducted several semi-automatic experiments of ice-area classification in the Quelccaya Ice Cap, Peru, and compared their results to a hand-digitization. The employed methods comprised pixel-based supervised and unsupervised techniques, including non-parametric approaches, such as fuzzy logic and empirical decision trees (thresholds). Although the hand digitization reached an accuracy of 99.8%, the semi-automated results obtained accuracies varying from 73.7% (fuzzy technique) to 96.7% (Spectral Angle Mapper—SAM). The mean accuracy was 84.8%. Specifically regarding the empirical results for the decision trees, the accuracies ranged from 93.0% to 93.9%. In their turn, Rittger et al. (2013) [60] compared the accuracies of three MODIS global products, namely, the binary and fractional snow product (MOD10A1), based on the NDSI, and another fractional snow product (MODSCAG) based on spectral mixture analysis. The study areas relate to the mountain ranges located in the USA (Sierra Nevada, Rocky Mountains, the Upper Rio Grande) and in the Nepal Himalaya, southern Asia, which are not included in tropical regions. The attained precision values ranged from 71.2% to 99.9%, and the mean precision values for the four study areas were 95,8% for the binary MOD10A1, 98.4% for the fractional MOD10A1, and 99.2% for the fractional MODSCAG. We obtained a mean precision (considering all TAGs and all years of the time series) of nearly 97.6%.
Post-classification filters specially adapted for the detection of glacial cover according to the GLIMS definition were also applied [20]. The application of filters helped filling the gap of non-observed areas and correct temporal and spatial inconsistencies, improving the classification accuracy by up to 6.91% in some years. Nevertheless, the gap-filling filters did not improve the map accuracy in the first and last years of the time series due to the smaller size of the filter’s temporal kernel and the limited number of images in the first years of the time series [37]. The first years of the time series were not consistent with the remaining time series due to the limited number of images. As a result, we used the year 1990 as the baseline for the gain–loss analysis and removed the mapping between the years 1985 and 1989. The accuracy assessment, particularly in the intermediate years, indicated that the proposed methodology generates comparable maps, presenting spatial and temporal coherence, useful for glacier change analysis.
As mentioned in [8] and other studies about the mass balance analysis of the tropical Andean glaciers, although there have been some sporadic gains in the TAGs, the general trend is negative, i.e., glaciers are retreating in the Andes in the past decades. Our study corroborates these findings. Detecting clouds and seasonal snow is critical to obtain consistent estimates of glaciers and measure their change over time with satellite imagery [42]. To overcome this obstacle, we refined cloud masking methods (i.e., CloudScore and TDOM) and applied temporal and spatial filters to fill the unobserved gaps due to cloud (using scenes with cloud cover <70%).
The annual reduction in the NIR and RED bands based on the 25th percentile of the NDSI values allowed seasonal snow detection. Our algorithm excluded it as soon as the non-glacial ground became exposed. This procedure enabled us to select the pixels with less snow coverage in the image collection, improving the Minimum NDSI composite used by Huang [49]. The method proposed by Huang [49], based on the minimum reducer in the NIR and RED bands, generates noise when there are pixels where clouds and cloud shadow have not been properly masked out, leading to underestimation of the glacier cover in the classification. Moreover, an appropriate selection of images and spectral bands [8] was also crucial for mapping glaciers to avoid confusion with temporary and seasonal snow cover. Finally, using regional thresholds in the empirical decision tree led to better classification results than a similar study for this region that exclusively relied on a single threshold-based decision tree [7].
When referring to the distribution of tropical glaciers, a worldwide reference is the RGI of the GLIMS initiative [19,39]. However, despite being widely used in many studies, this reference has a significant limitation to calculating the distribution of tropical glaciers and making comparisons of multitemporal glacier loss, since the RGI data are from different dates in each place as it is targeted at producing a glacier inventory. Therefore, in our study, we reconstructed the distribution of the tropical Andean glaciers in the last three decades and we updated the distribution of glaciers to the year 2020. Our new data set allows comparing glacier change across different countries. Peru and Bolivia comprise 72.76% and 20.35% of the total TAGs, while Ecuador, Colombia, Chile, and Argentina account for 3.89%, 2.18%, 0.78%, and 0.04%, respectively. In turn, Venezuela is the country with the smallest glacier extent. The countries with less glacier coverage in 1990 are the ones that have had the most significant glaciers losses in percentage by 2020.
Overall, the TAGs are experiencing a steady retreating trend, and their extent reduced by an average of 42% between 1990 to 2020, which indicates a loss of nearly half of the tropical Andean glaciers extension in the last three decades. This fast loss rate of the TAGs is an unprecedented change in tropical glaciers [8], potentially linked with the ongoing warming of the planet as a consequence of climate change and climate variability [19,38]. The increase in global temperature [61] as well as the variability of sea surface temperatures (SSTs) of the Pacific Ocean can be the potential drivers of the TAGs retreat [8] we presented in this study. Therefore, the retreat of the TAGs constitutes an early case of the need to adapt to climate change to reduce the economic and social impacts associated with the loss of glacier cover and the freshwater they produce [11].
TAGs dynamics are also influenced by altitude, latitude, longitude, slope, aspect, and country. We realized that the tropical Andean glaciers below 5000 masl are the ones that changed the most in the last three decades. This result is consistent with previous studies, such as a glacier in Bolivia and Peru [10,35]. The TAGs located below 5000 masl lost almost 80.25% of their initial area in the year 1990. These lower-height glaciers retreated faster than higher ones, which can be explained by the high negative impact of the net short-wave radiation linked with cloud cover and surface albedo [8,35]. The albedo feedback is suggested as the primary cause of accelerated retreat [8]. It is worth mentioning that the slope variable also influences the glacier loss, since areas with a slope above 80° have lost more glacier cover than average. Therefore, if the current environmental conditions remain, it is very likely that in the future the glacier surface will disappear at low elevations below 5000 masl in the TAGs, accounting for approximately 171.26 km2, representing 12.15% of the TAGs in the year 2020.
We also found an accelerated retreat of a glacier in the hillsides belonging to the Amazon basin. Climatic and non-climatic factors could explain this significant change in the Amazon in recent years. One of the non-climatic factors could be the increase in forest burning in recent years in the Amazon, which generates black carbon [62,63] that can accelerate glacial retreat when entering the glacier surface [19,62]. However, it can also be explained by solar radiation in the tropical zone, which is more or less constant throughout the year [8] and lately intensified by climatic changes. Even though the hydrological contribution of glaciers to the Amazon basin is minor, this may influence the rivers fed by the Andean glaciers in the Amazon basin [35], as is the case of the Madeira River in Brazil [64], fed by the glaciers in Peru and Bolivia.
The countries with less TAG coverage, such as Venezuela, Argentina, Colombia, and Chile, are the ones that have lost their glacier coverage more rapidly in the last three decades. Furthermore, this may have a possible relationship with the geographic latitude and longitude of these TAGs. We observed that glaciers with latitudes farther from the equator are the ones that have lost more glacier coverage with respect to the total TAGs mean. In addition, we found that glaciers with latitudes ranging from −9° to 0° lost less area than the total average, and we noted that they are glaciers with more extensive areas in the tropical zone. The glaciers with latitudes closer to the Pacific Ocean, where the Blanca Cordillera (mountain range) in Peru is located, had lower losses than the total average, which could be explained by a possible association of glacier cover loss in the TAGs with the climatic behavior of the Pacific Ocean [8].
Several studies have shown the high (negative) correlation between annual air temperature and glacier mass balance in tropical regions [65]. Notably, the freezing level heights (FLHs), i.e., the elevation where the 0 °C air temperature is found, has been used as a proxy for the TAGs mass balance [15,66,67]. Higher FLH values correspond to glacier mass loss (or more negative mass balance). The continuous tropical FLH increase is basically linked with the Pacific ocean’s sea surface warming trend [67,68]. In some years during the ENSO (El Niño and La Niña events, with recurrence cycles from 2 to 7 years) occurrence, the Pacific ocean’s warming conditions could produce significant changes in the glaciers´ surface extension [67,68], but with different impacts on each country and region. Those sudden and considerable changes of the TAGs surface detected in the previous analyses can be due, in several cases, to ENSO effects, but this is out of the scope of the present research.
As limitations of the proposed methodology, we ought to mention that the data dimensionality reduction depends on the quantity and quality of the input images. The adopted filters do not add information into the classification process, but merely correct the temporal and spatial inconsistencies. The assessment of the glaciers´ retreat was conducted on a 2D basis, being restricted to the observable surface in the images. It would be highly recommendable to advance this research into a volumetric estimate of glacier loss, and this initiative would probably rely on sensors dedicated to measure the ice-sheet thickness (e.g., IceSat-2/Atlas). Finally, the classification methodology is not fully automatic, depending on expert image interpreters and heuristic routines for setting the decision tree thresholds.

5. Conclusions

In this research, we reconstructed the annual extensions of the TAGs for the last three decades, based on a semi-automated method and the historical Landsat series in the GEE platform. This approach built a new and updated dataset of the TAGs based on remote-sensing time-series data. The change detection analysis presented in this study presented a new understanding of annual glacier coverage changes across different zones of the tropical Andean Mountain Ranges. As a result, we provided the first multiannual and spatially detailed analysis covering the total extent of the TAGs. The observed changes highlight the worrisome retreat of glaciers throughout the tropical Andes, which will likely lead to considerable environmental, economic, and cultural problems in this region. In addition, the results represent an important source of information for developing water management strategies capable of meeting the challenges of climate change. The results shown in the present work are a source of data for understanding the relationship between glacier recession and different factors that determine either the progressive loss or permanence of tropical glaciers. Future steps in our research project are to provide the TAG dataset through a dashboard to support better decision making and water resource planning and management, climate change adaptation, and improving glacier monitoring.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs14091974/s1, Table S1: Number of images per sensor and year in the study area; Table S2: Access to multiannual classification products and glacier-loss maps; Table S3: Results of glacier surface coverage in the TAGs by country; Table S4: Results of glacier surface coverage (km2) in the TAGs at different geographic latitudes; Table S5: Results of glacier surface coverage (km2) in the TAGs at different geographic longitudes; Table S6: Results of glacier surface coverage (km2) in the TAGs by aspect; Table S7: Results of glacier surface coverage (km2) in the TAGs on different slopes; Table S8: Results of glacier surface coverage (km2) in the TAGs at different altitudes; Table S9: Glacier surface coverage results in the TAGs by basin; Table S10: Validation results of the glacier cover classification without filters; Table S11: Validation results of the glacier cover classification with filters; Table S12: Validation metrics disregarding and considering area proportions; Figure S1: Number of images per sensor and year in the study area; and Figure S2: Benchmark glaciers and classification regions. The division of regions is based on the criteria of similar geoclimatic conditions; Figure S3: Evolution of annual glacier area coverage loss in the TAGs from 1987 to 2019: (a) loss at different latitudes and (b) loss at different longitudes; Figure S4: Evolution of annual glacier area coverage loss in the TAGs from 1990 to 2020 at different glacier slope orientations: (a) polar plot of glacier loss, and (b) change in annual glacier area loss; Figure S5: Loss of the annual glacier area coverage in the TAGs from 1990 to 2020 in different slope ranges; Figure S6: Annual glacier loss in the TAGs from 1990 to 2020 at different altitude ranges; Figure S7: Evolution of the annual loss of glacier area in the TAGs from 1990 to 2020 by basins: (a) percentage loss of glacier area and (b) cumulative loss of glacier area in km2; and Figure S8: Multitemporal validation metrics for the final multiyear glacier cover maps.

Author Contributions

Conceptualization, E.Y.T.C., M.O.B., R.E.-V., N.M., R.C., C.A., K.H., C.Y. and C.M.S.J.; methodology, E.Y.T.C., R.E.-V., C.Y. and C.A.; software, E.Y.T.C.; validation, E.Y.T.C., M.O.B., N.M., R.C., C.A. and K.H.; formal analysis, E.Y.T.C., R.E.-V., M.O.B., C.A. and C.M.S.J.; data curation, E.Y.T.C., M.O.B., R.E.-V., N.M., R.C., C.A. and K.H.; visualization, E.Y.T.C., R.E.-V., C.A. and C.M.S.J.; writing—original draft, E.Y.T.C., R.E.-V., C.A. and C.M.S.J.; writing—review and editing, E.Y.T.C., M.O.B., R.E.-V., N.M., R.C., C.A., K.H., C.Y. and C.M.S.J.; supervision, R.E.-V., C.A., C.Y. and C.M.S.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research is part of the MapBiomas project implemented by the Amazon Geo-referenced Socio-Environmental Information Network (RAISG) with the support of the Gordon & Betty Moore Foundation, Norway’s Inter-national Climate and Forest Initiative (NICFI), Quadrature Climate Foundation, Good Energies Foundation, Regnskoqfondet, and Instituto Clima e Sociedade (ICS); the study was also conducted as part of La Lagunas de Origen Glaciar en el Perú: Evolución, Peligros e Impacto del Cambio Climático (Proyecto GLOP), which is a project jointly funded by the UK Natural Environment Research Council (NERC; Grant: NE/S01330X/1) and the Consejo Nacional de Ciencia, Tecnología e Innovación (CONCYTEC; Esquema financiero E031-2018-01-NERC, contrato No. 007-2019-FONDECYT) and UNAM research group (Geomática Ambiental). This research also has the support of the Brazilian National Council for Scientific and Technological Development—CNPq (Grant: 303523/2018-2; 311324/2021-5).

Data Availability Statement

Not applicable.

Acknowledgments

We would like to thank the team involved in the MapBiomas Amazonia Collections. We gratefully acknowledge USGS and NASA for the courtesy of the Landsat images. We are also grateful to Google for providing access to the Google Earth Engine platform that enables the processing of all data. We would like to thank the institutional partners of the Amazonian Network of Georeferenced Socio-Environmental Information (RAISG). We also thank Danilo Iturralde who helped with the photo interpretation of Ecuador’s data validation process. Finally, we would like to express our gratitude to Sandra Rios, Saul Cuellar, Cicero Cardoso Augusto, Marcos Rosa, João Victor N. Siqueira, Julia Shimbo, Carmen Josse, Richard Chase Smith and Tasso Azevedo for their collaboration, support, and direction of the Amazon MapBiomas Project.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Huh, K.; Baraër, M.; Mark, B.; Ahn, Y. Evaluating Glacier Volume Changes since the Little Ice Age Maximum and Consequences for Stream Flow by Integrating Models of Glacier Flow and Hydrology in the Cordillera Blanca, Peruvian Andes. Water 2018, 10, 1732. [Google Scholar] [CrossRef] [Green Version]
  2. Kaser, G.; Osmaston, H. Tropical Glaciers; Cambridge University Press: Cambridge, UK, 2002; ISBN 0-521-63333-8. [Google Scholar]
  3. Calizaya, E.; Mejía, A.; Barboza, E.; Calizaya, F.; Corroto, F.; Salas, R.; Vásquez, H.; Turpo, E. Modelling Snowmelt Runoff from Tropical Andean Glaciers under Climate Change Scenarios in the Santa River Sub-Basin (Peru). Water 2021, 13, 3535. [Google Scholar] [CrossRef]
  4. Barr, I.D.; Lynch, C.M.; Mullan, D.; De Siena, L.; Spagnolo, M. Volcanic Impacts on Modern Glaciers: A Global Synthesis. Earth-Sci. Rev. 2018, 182, 186–203. [Google Scholar] [CrossRef] [Green Version]
  5. Marzeion, B.; Cogley, J.G.; Richter, K.; Parkes, D. Attribution of Global Glacier Mass Loss to Anthropogenic and Natural Causes. Science 2014, 345, 919–921. [Google Scholar] [CrossRef]
  6. Eyring, V.; Gillett, N.P.; Achutarao, K.; Barimalala, R.; Barreiro Parrillo, M.; Bellouin, N.; Cassou, C.; Durack, P.; Kosaka, Y.; McGregor, S. Human Influence on the Climate System: Contribution of Working Group I to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. In IPCC Sixth Assessment Report; Cambridge University Press: Cambridge, UK, 2021. [Google Scholar]
  7. Seehaus, T.; Malz, P.; Sommer, C.; Lippl, S.; Cochachin, A.; Braun, M. Changes of the Tropical Glaciers throughout Peru between 2000 and 2016—Mass Balance and Area Fluctuations. Cryosphere 2019, 13, 2537–2556. [Google Scholar] [CrossRef] [Green Version]
  8. Rabatel, A.; Francou, B.; Soruco, A.; Gomez, J.; Cáceres, B.; Ceballos, J.L.; Basantes, R.; Vuille, M.; Sicart, J.-E.; Huggel, C.; et al. Current State of Glaciers in the Tropical Andes: A Multi-Century Perspective on Glacier Evolution and Climate Change. Cryosphere 2013, 7, 81–102. [Google Scholar] [CrossRef] [Green Version]
  9. Bradley, R.S.; Vuille, M.; Diaz, H.F.; Vergara, W. Threats to Water Supplies in the Tropical Andes. Science 2006, 312, 1755–1756. [Google Scholar] [CrossRef]
  10. Veettil, B.K.; Wang, S.; Simões, J.C.; Pereira, S.F.R. Glacier Monitoring in the Eastern Mountain Ranges of Bolivia from 1975 to 2016 Using Landsat and Sentinel-2 Data. Environ. Earth Sci. 2018, 77, 452. [Google Scholar] [CrossRef]
  11. Vergara, W.; Deeb, A.; Valencia, A.; Bradley, R.; Francou, B.; Zarzar, A.; Grünwaldt, A.; Haeussling, S. Economic Impacts of Rapid Glacier Retreat in the Andes. Eos Trans. AGU 2007, 88, 261–264. [Google Scholar] [CrossRef]
  12. Escanilla-Minchel, R.; Alcayaga, H.; Soto-Alvarez, M.; Kinnard, C.; Urrutia, R. Evaluation of the Impact of Climate Change on Runoff Generation in an Andean Glacier Watershed. Water 2020, 12, 3547. [Google Scholar] [CrossRef]
  13. Ramírez, N.; Melfo, A.; Resler, L.M.; Llambí, L.D. The End of the Eternal Snows: Integrative Mapping of 100 Years of Glacier Retreat in the Venezuelan Andes. Arct. Antarct. Alp. Res. 2020, 52, 563–581. [Google Scholar] [CrossRef]
  14. Wang, S.-J.; Zhou, L.-Y. Integrated Impacts of Climate Change on Glacier Tourism. Adv. Clim. Chang. Res. 2019, 10, 71–79. [Google Scholar] [CrossRef]
  15. Vuille, M.; Carey, M.; Huggel, C.; Buytaert, W.; Rabatel, A.; Jacobsen, D.; Soruco, A.; Villacis, M.; Yarleque, C.; Elison Timm, O.; et al. Rapid Decline of Snow and Ice in the Tropical Andes—Impacts, Uncertainties and Challenges Ahead. Earth-Sci. Rev. 2018, 176, 195–213. [Google Scholar] [CrossRef] [Green Version]
  16. Wood, J.L.; Harrison, S.; Wilson, R.; Emmer, A.; Yarleque, C.; Glasser, N.F.; Torres, J.C.; Caballero, A.; Araujo, J.; Bennett, G.L.; et al. Contemporary Glacial Lakes in the Peruvian Andes. Glob. Planet. Chang. 2021, 204, 103574. [Google Scholar] [CrossRef]
  17. Baťka, J.; Vilímek, V.; Štefanová, E.; Cook, S.J.; Emmer, A. Glacial Lake Outburst Floods (GLOFs) in the Cordillera Huayhuash, Peru: Historic Events and Current Susceptibility. Water 2020, 12, 2664. [Google Scholar] [CrossRef]
  18. Funaki, S.; Asaoka, Y. Long-Term Change in Ablation Area of Tropical Glaciers by Landsat Data. Procedia Eng. 2016, 154, 168–175. [Google Scholar] [CrossRef] [Green Version]
  19. Veettil, B.K.; Kamp, U. Global Disappearance of Tropical Mountain Glaciers: Observations, Causes, and Challenges. Geosciences 2019, 9, 196. [Google Scholar] [CrossRef] [Green Version]
  20. Racoviteanu, A.E.; Paul, F.; Raup, B.; Khalsa, S.J.S.; Armstrong, R. Challenges and Recommendations in Mapping of Glacier Parameters from Space: Results of the 2008 Global Land Ice Measurements from Space (GLIMS) Workshop, Boulder, CO, USA. Ann. Glaciol. 2009, 50, 53–69. [Google Scholar] [CrossRef] [Green Version]
  21. Taloor, A.K.; Kothyari, G.C.; Manhas, D.S.; Bisht, H.; Mehta, P.; Sharma, M.; Mahajan, S.; Roy, S.; Singh, A.K.; Ali, S. Spatio-Temporal Changes in the Machoi Glacier Zanskar Himalaya India Using Geospatial Technology. Quat. Sci. Adv. 2021, 4, 100031. [Google Scholar] [CrossRef]
  22. Albert, T.H. Evaluation of Remote Sensing Techniques for Ice-Area Classification Applied to the Tropical Quelccaya Ice Cap, Peru. Polar Geogr. 2002, 26, 210–226. [Google Scholar] [CrossRef]
  23. Paul, F.; Kääb, A.; Haeberli, W. Recent Glacier Changes in the Alps Observed by Satellite: Consequences for Future Monitoring Strategies. Glob. Planet. Chang. 2007, 56, 111–122. [Google Scholar] [CrossRef]
  24. Marta, S.; Azzoni, R.S.; Fugazza, D.; Tielidze, L.; Chand, P.; Sieron, K.; Almond, P.; Ambrosini, R.; Anthelme, F.; Alviz Gazitúa, P.; et al. The Retreat of Mountain Glaciers since the Little Ice Age: A Spatially Explicit Database. Data 2021, 6, 107. [Google Scholar] [CrossRef]
  25. Hansen, M.C.; Loveland, T.R. A Review of Large Area Monitoring of Land Cover Change Using Landsat Data. Remote Sens. Environ. 2012, 122, 66–74. [Google Scholar] [CrossRef]
  26. Veettil, B.K. Glacier Mapping in the Cordillera Blanca, Peru, Tropical Andes, Using Sentinel-2 and Landsat Data. Singap. J. Trop. Geogr. 2018, 39, 351–363. [Google Scholar] [CrossRef]
  27. Veetil, B. Un Análisis Comparativo Del Retroceso Glaciar En Los Andes Tropicales Usando Teledetección. Santiago Investig. Geogr. Chile 2016, 52, 3–36. [Google Scholar] [CrossRef] [Green Version]
  28. Zhu, Z. Change Detection Using Landsat Time Series: A Review of Frequencies, Preprocessing, Algorithms, and Applications. ISPRS J. Photogramm. Remote Sens. 2017, 130, 370–384. [Google Scholar] [CrossRef]
  29. Gorelick, N.; Hancher, M.; Dixon, M.; Ilyushchenko, S.; Thau, D.; Moore, R. Google Earth Engine: Planetary-Scale Geospatial Analysis for Everyone. Remote Sens. Environ. 2017, 202, 18–27. [Google Scholar] [CrossRef]
  30. Li, K.; Xu, E. High-Accuracy Continuous Mapping of Surface Water Dynamics Using Automatic Update of Training Samples and Temporal Consistency Modification Based on Google Earth Engine: A Case Study from Huizhou, China. ISPRS J. Photogramm. Remote Sens. 2021, 179, 66–80. [Google Scholar] [CrossRef]
  31. Souza, C.; Kirchhoff, F.; Oliveira, B.; Ribeiro, J.; Sales, M. Long-Term Annual Surface Water Change in the Brazilian Amazon Biome: Potential Links with Deforestation, Infrastructure Development and Climate Change. Water 2019, 11, 566. [Google Scholar] [CrossRef] [Green Version]
  32. Kamp, U.; Pan, C.G. Inventory of Glaciers in Mongolia, Derived from Landsat Imagery from 1989 to 2011. Geogr. Ann. Ser. A Phys. Geogr. 2015, 97, 653–669. [Google Scholar] [CrossRef]
  33. Paul, F.; Winsvold, S.; Kääb, A.; Nagler, T.; Schwaizer, G. Glacier Remote Sensing Using Sentinel-2. Part II: Mapping Glacier Extents and Surface Facies, and Comparison to Landsat 8. Remote Sens. 2016, 8, 575. [Google Scholar] [CrossRef] [Green Version]
  34. Burns, P.; Nolin, A. Using Atmospherically-Corrected Landsat Imagery to Measure Glacier Area Change in the Cordillera Blanca, Peru from 1987 to 2010. Remote Sens. Environ. 2014, 140, 165–178. [Google Scholar] [CrossRef] [Green Version]
  35. Kozhikkodan Veettil, B.; de Souza, S.F. Study of 40-Year Glacier Retreat in the Northern Region of the Cordillera Vilcanota, Peru, Using Satellite Images: Preliminary Results. Remote Sens. Lett. 2017, 8, 78–85. [Google Scholar] [CrossRef]
  36. Borja, M.O.; Camargo, R.; Moreno, N.; Turpo, E.; Villacis, S. A Long-Term Land Cover And Land Use Mapping Methodology For The Andean Amazon. In Proceedings of the 2020 IEEE Latin American GRSS & ISPRS Remote Sensing Conference (LAGIRS), Santiago, Chile, 21–26 March 2020; IEEE: Santiago, Chile, 2020; pp. 386–391. [Google Scholar]
  37. Crawford, C.J.; Manson, S.M.; Bauer, M.E.; Hall, D.K. Multitemporal Snow Cover Mapping in Mountainous Terrain for Landsat Climate Data Record Development. Remote Sens. Environ. 2013, 135, 224–233. [Google Scholar] [CrossRef] [Green Version]
  38. Francou, B.; Ramirez, E.; Cáceres, B.; Mendoza, J. Glacier Evolution in the Tropical Andes during the Last Decades of the 20th Century: Chacaltaya, Bolivia, and Antizana, Ecuador. AMBIO A J. Hum. Environ. 2000, 29, 416–422. [Google Scholar] [CrossRef]
  39. RGI Consortium. Randolph Glacier Inventory—A Dataset of Global Glacier Outlines: Version 6.0; National Snow and Ice Data Center: Boulder, CO, USA, 2017. [Google Scholar] [CrossRef]
  40. Sagredo, E.A.; Lowell, T.V. Climatology of Andean Glaciers: A Framework to Understand Glacier Response to Climate Change. Glob. Planet. Chang. 2012, 86–87, 101–109. [Google Scholar] [CrossRef]
  41. Veettil, B.K.; Wang, S.; Florêncio de Souza, S.; Bremer, U.F.; Simões, J.C. Glacier Monitoring and Glacier-Climate Interactions in the Tropical Andes: A Review. J. S. Am. Earth Sci. 2017, 77, 218–246. [Google Scholar] [CrossRef]
  42. Herreid, S.; Pellicciotti, F. The State of Rock Debris Covering Earth’s Glaciers. Nat. Geosci. 2020, 13, 621–627. [Google Scholar] [CrossRef]
  43. Shukla, A.; Gupta, R.P.; Arora, M.K. Estimation of Debris Cover and Its Temporal Variation Using Optical Satellite Sensor Data: A Case Study in Chenab Basin, Himalaya. J. Glaciol. 2009, 55, 444–452. [Google Scholar] [CrossRef] [Green Version]
  44. Falaschi, D.; Bolch, T.; Lenzano, M.G.; Tadono, T.; Lo Vecchio, A.; Lenzano, L. New Evidence of Glacier Surges in the Central Andes of Argentina and Chile. Prog. Phys. Geogr. Earth Environ. 2018, 42, 792–825. [Google Scholar] [CrossRef]
  45. Wulder, M.A.; White, J.C.; Loveland, T.R.; Woodcock, C.E.; Belward, A.S.; Cohen, W.B.; Fosnight, E.A.; Shaw, J.; Masek, J.G.; Roy, D.P. The Global Landsat Archive: Status, Consolidation, and Direction. Remote Sens. Environ. 2016, 185, 271–283. [Google Scholar] [CrossRef] [Green Version]
  46. United States Geological Survey. Landsat Collection 1 Surface Reflectance. 2021. Available online: https://www.usgs.gov/landsat-missions/landsat-collection-1-surface-reflectance (accessed on 1 March 2021).
  47. Paul, F.; Bolch, T.; Kääb, A.; Nagler, T.; Nuth, C.; Scharrer, K.; Shepherd, A.; Strozzi, T.; Ticconi, F.; Bhambri, R.; et al. The Glaciers Climate Change Initiative: Methods for Creating Glacier Area, Elevation Change and Velocity Products. Remote Sens. Environ. 2015, 162, 408–426. [Google Scholar] [CrossRef] [Green Version]
  48. NASA. Crisis Mapping Toolkit (CMT) V1; Environmental Science Earth Air Space Exoplanet; NASA: Washington, DC, USA, 2021.
  49. Huang, L.; Li, Z.; Zhou, J.M.; Zhang, P. An Automatic Method for Clean Glacier and Nonseasonal Snow Area Change Estimation in High Mountain Asia from 1990 to 2018. Remote Sens. Environ. 2021, 258, 112376. [Google Scholar] [CrossRef]
  50. Housman, I.; Chastain, R.; Finco, M. An Evaluation of Forest Health Insect and Disease Survey Data and Satellite-Based Remote Sensing Forest Change Detection Methods: Case Studies in the United States. Remote Sens. 2018, 10, 1184. [Google Scholar] [CrossRef] [Green Version]
  51. Souza, C.M.Z.; Shimbo, J.; Rosa, M.R.; Parente, L.L.; Alencar, A.A.; Rudorff, B.F.T.; Hasenack, H.; Matsumoto, M.; Ferreira, L.G.; Souza-Filho, P.W.M.; et al. Reconstructing Three Decades of Land Use and Land Cover Changes in Brazilian Biomes with Landsat Archive and Earth Engine. Remote Sens. 2020, 12, 2735. [Google Scholar] [CrossRef]
  52. Hall, D.K.; Riggs, G.A. Normalized-Difference Snow Index (NDSI). In Encyclopedia of Snow, Ice and Glaciers; Singh, V.P., Singh, P., Haritashya, U.K., Eds.; Encyclopedia of Earth Sciences Series; Springer: Dordrecht, The Netherlands, 2011; pp. 779–780. ISBN 978-90-481-2641-5. [Google Scholar]
  53. Macander, M.J.; Swingley, C.S.; Joly, K.; Raynolds, M.K. Landsat-Based Snow Persistence Map for Northwest Alaska. Remote Sens. Environ. 2015, 163, 23–31. [Google Scholar] [CrossRef]
  54. Winsvold, S.H.; Kaab, A.; Nuth, C. Regional Glacier Mapping Using Optical Satellite Data Time Series. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2016, 9, 3698–3711. [Google Scholar] [CrossRef] [Green Version]
  55. Congalton, R.G.; Green, K. Assessing the Accuracy of Remotely Sensed Data: Principles and Practices; CRC Press: Boca Raton, FL, USA, 2009; ISBN 978-1-4200-5512-2. [Google Scholar]
  56. Donchyts, G.; Llano; Baart, F.; Hereñú, D.; Braaten, J. Gee-Community/Qgis-Earthengine-Plugin: V0.0.4 (0.0.4-Alpha-Doi); GitHub: San Francisco, CA, USA, 2022. [Google Scholar] [CrossRef]
  57. Llano, X.C. AcATaMa—QGIS Plugin for Accuracy Assessment of Thematic Maps. 2019. Available online: https://github.com/SMByC/AcATaMa (accessed on 1 March 2021).
  58. Olofsson, P.; Foody, G.M.; Herold, M.; Stehman, S.V.; Woodcock, C.E.; Wulder, M.A. Good Practices for Estimating Area and Assessing Accuracy of Land Change. Remote Sens. Environ. 2014, 148, 42–57. [Google Scholar] [CrossRef]
  59. Olson, D.L.; Delen, D. Advanced Data Mining Techniques; Springer: Berlin/Heidelberg, Germany, 2008; ISBN 978-3-540-76916-3. [Google Scholar]
  60. Rittger, K.; Painter, T.H.; Dozier, J. Assessment of Methods for Mapping Snow Cover from MODIS. Adv. Water Resour. 2013, 51, 367–380. [Google Scholar] [CrossRef]
  61. Randall, D.A.; Wood, R.A.; Bony, S.; Colman, R.; Fichefet, T.; Fyfe, J.; Kattsov, V.; Pitman, A.; Shukla, J.; Srinivasan, J.; et al. Climate Models and Their Evaluation. In Climate Change 2007: The Physical Science Basis. Contribution of Working Group I to the Fourth Assessment Report of the IPCC (FAR); Cambridge University Press: Cambridge, UK, 2007; pp. 589–662. [Google Scholar]
  62. de Magalhães, N.; Evangelista, H.; Condom, T.; Rabatel, A.; Ginot, P. Amazonian Biomass Burning Enhances Tropical Andean Glaciers Melting. Sci. Rep. 2019, 9, 16914. [Google Scholar] [CrossRef]
  63. Barboza Castillo, E.; Turpo Cayo, E.Y.; de Almeida, C.M.; Salas López, R.; Rojas Briceño, N.B.; Silva López, J.O.; Barrena Gurbillón, M.Á.; Oliva, M.; Espinoza-Villar, R. Monitoring Wildfires in the Northeastern Peruvian Amazon Using Landsat-8 and Sentinel-2 Imagery in the GEE Platform. IJGI 2020, 9, 564. [Google Scholar] [CrossRef]
  64. da Rocha, N.S.; Veettil, B.K.; Grondona, A.; Rolim, S. The Influence of ENSO and PDO on Tropical Andean Glaciers and Their Impact on the Hydrology of the Amazon Basin. Singap. J. Trop. Geogr. 2019, 40, 346–360. [Google Scholar] [CrossRef]
  65. Vuille, M.; Francou, B.; Wagnon, P.; Juen, I.; Kaser, G.; Mark, B.G.; Bradley, R.S. Climate Change and Tropical Andean Glaciers: Past, Present and Future. Earth-Sci. Rev. 2008, 89, 79–96. [Google Scholar] [CrossRef] [Green Version]
  66. Schauwecker, S.; Rohrer, M.; Huggel, C.; Endries, J.; Montoya, N.; Neukom, R.; Perry, B.; Salzmann, N.; Schwarb, M.; Suarez, W. The Freezing Level in the Tropical Andes, Peru: An Indicator for Present and Future Glacier Extents: The Freezing Level in the Tropical Andes. J. Geophys. Res. Atmos. 2017, 122, 5172–5189. [Google Scholar] [CrossRef]
  67. Bradley, R.S.; Keimig, F.T.; Diaz, H.F.; Hardy, D.R. Recent Changes in Freezing Level Heights in the Tropics with Implications for the Deglacierization of High Mountain Regions. Geophys. Res. Lett. 2009, 36, L17701. [Google Scholar] [CrossRef] [Green Version]
  68. Yarleque, C.; Vuille, M.; Hardy, D.R.; Timm, O.E.; De la Cruz, J.; Ramos, H.; Rabatel, A. Projections of the Future Disappearance of the Quelccaya Ice Cap in the Central Andes. Sci. Rep. 2018, 8, 15564. [Google Scholar] [CrossRef]
Figure 1. Location of tropical Andean glaciers based on the RGI (Version 6.0) [39] and their climatic classification based on Sagredo and Lowell [27,40,41].
Figure 1. Location of tropical Andean glaciers based on the RGI (Version 6.0) [39] and their climatic classification based on Sagredo and Lowell [27,40,41].
Remotesensing 14 01974 g001
Figure 2. Methodological steps for the TAGs classification implemented in Google Earth Engine.
Figure 2. Methodological steps for the TAGs classification implemented in Google Earth Engine.
Remotesensing 14 01974 g002
Figure 3. Example of a violin box plot diagram for the Landsat TM and ETM+ sensors in 2000.
Figure 3. Example of a violin box plot diagram for the Landsat TM and ETM+ sensors in 2000.
Remotesensing 14 01974 g003
Figure 4. Example of a violin box plot diagram for the Landsat 8 OLI sensor in 2020.
Figure 4. Example of a violin box plot diagram for the Landsat 8 OLI sensor in 2020.
Remotesensing 14 01974 g004
Figure 5. Empirical decision tree used for the classification of the TAGs.
Figure 5. Empirical decision tree used for the classification of the TAGs.
Remotesensing 14 01974 g005
Figure 6. Effect of the gap-fill filter application: (a) time series of the original classification (before the filter application), and (b) time series of the filtered classification (gap-filled).
Figure 6. Effect of the gap-fill filter application: (a) time series of the original classification (before the filter application), and (b) time series of the filtered classification (gap-filled).
Remotesensing 14 01974 g006
Figure 7. Effect of the temporal filter application for each rule: (a) first-year rule, applied only to the year 1985; (b) last-year rule, applied only to the year 2020; intermediate-year rules with (c) a three-year moving window; (d) a four-year moving window; and (e) a five-year moving window. The first- and last-year rules maintain class continuity at the beginning and end of the series, while the intermediate rules identify the inconsistent transitions and correct them, such as the sudden appearance of glacier values in a non-glacier area, or vice versa. Given that glaciers tend to be slow-changing entities, sudden changes in a time series are most likely classification inconsistencies.
Figure 7. Effect of the temporal filter application for each rule: (a) first-year rule, applied only to the year 1985; (b) last-year rule, applied only to the year 2020; intermediate-year rules with (c) a three-year moving window; (d) a four-year moving window; and (e) a five-year moving window. The first- and last-year rules maintain class continuity at the beginning and end of the series, while the intermediate rules identify the inconsistent transitions and correct them, such as the sudden appearance of glacier values in a non-glacier area, or vice versa. Given that glaciers tend to be slow-changing entities, sudden changes in a time series are most likely classification inconsistencies.
Remotesensing 14 01974 g007
Figure 8. Effect of the frequency filter application: (a) time series of the original classification (before the filter application), and (b) time series of the filtered classification.
Figure 8. Effect of the frequency filter application: (a) time series of the original classification (before the filter application), and (b) time series of the filtered classification.
Remotesensing 14 01974 g008
Figure 9. Example of the temporal permanence filter application: (a) time series of the original classification, and (b) time series of the filtered classification.
Figure 9. Example of the temporal permanence filter application: (a) time series of the original classification, and (b) time series of the filtered classification.
Remotesensing 14 01974 g009
Figure 10. TAGs results for the year 2020: (a) percentage of the spatial distribution of glaciers by country; (b) occurrence of glaciers by altitude ranges; (c) distribution of glaciers by terrain aspect; and (d) percentage of spatial distribution of glaciers by terrain slope.
Figure 10. TAGs results for the year 2020: (a) percentage of the spatial distribution of glaciers by country; (b) occurrence of glaciers by altitude ranges; (c) distribution of glaciers by terrain aspect; and (d) percentage of spatial distribution of glaciers by terrain slope.
Remotesensing 14 01974 g010
Figure 11. TAGs results for the year 2020: (a) spatial distribution of glaciers by latitude, and (b) longitude.
Figure 11. TAGs results for the year 2020: (a) spatial distribution of glaciers by latitude, and (b) longitude.
Remotesensing 14 01974 g011
Figure 12. TAGs retreat between 1990 and 2020: (a) linear regression of the glacier area throughout the time series, (b) boxplot of the glacier area loss by decades.
Figure 12. TAGs retreat between 1990 and 2020: (a) linear regression of the glacier area throughout the time series, (b) boxplot of the glacier area loss by decades.
Remotesensing 14 01974 g012
Figure 13. Map of glacier cover loss by decade from 1990 to 2020 in the Vilcanota Cordillera (mountain range) in southern Peru.
Figure 13. Map of glacier cover loss by decade from 1990 to 2020 in the Vilcanota Cordillera (mountain range) in southern Peru.
Remotesensing 14 01974 g013
Figure 14. Evolution of the TAGs’ annual coverage loss by country in the period from 1990 to 2020.
Figure 14. Evolution of the TAGs’ annual coverage loss by country in the period from 1990 to 2020.
Remotesensing 14 01974 g014
Figure 15. Accuracy after filter application (a) multiannual accuracy plot, and (b) accuracy box plot for the cumulative sequential use of temporal filters.
Figure 15. Accuracy after filter application (a) multiannual accuracy plot, and (b) accuracy box plot for the cumulative sequential use of temporal filters.
Remotesensing 14 01974 g015
Table 1. Glacier-mapping summary statistics estimated for decennial periods.
Table 1. Glacier-mapping summary statistics estimated for decennial periods.
Period (# Year)Minimum (km2)Maximum (km2)MeanStandard
Deviation (s)
(km2/Year)
1991–20002089.202364.802220.3595.58
2001–20101815.572069.411928.3789.03
2011–20201409.121807.831670.35123.14
Table 2. Glacier-mapping summary statistics estimated for decennial periods per country.
Table 2. Glacier-mapping summary statistics estimated for decennial periods per country.
CountryGlacier Area (km2)
1990200020102020
Peru1743.451503.671320.851025.29
Bolivia499.58434.75364.02286.69
Ecuador86.1073.7068.5254.78
Colombia77.2160.1849.2130.73
Chile20.7715.7312.1510.96
Argentina1.150.730.680.63
Venezuela1.120.440.140.03
Total2429.382089.201815.571409.12
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Turpo Cayo, E.Y.; Borja, M.O.; Espinoza-Villar, R.; Moreno, N.; Camargo, R.; Almeida, C.; Hopfgartner, K.; Yarleque, C.; Souza, C.M., Jr. Mapping Three Decades of Changes in the Tropical Andean Glaciers Using Landsat Data Processed in the Earth Engine. Remote Sens. 2022, 14, 1974. https://doi.org/10.3390/rs14091974

AMA Style

Turpo Cayo EY, Borja MO, Espinoza-Villar R, Moreno N, Camargo R, Almeida C, Hopfgartner K, Yarleque C, Souza CM Jr. Mapping Three Decades of Changes in the Tropical Andean Glaciers Using Landsat Data Processed in the Earth Engine. Remote Sensing. 2022; 14(9):1974. https://doi.org/10.3390/rs14091974

Chicago/Turabian Style

Turpo Cayo, Efrain Yury, Maria Olga Borja, Raul Espinoza-Villar, Nicole Moreno, Rodney Camargo, Claudia Almeida, Kathrin Hopfgartner, Christian Yarleque, and Carlos M. Souza, Jr. 2022. "Mapping Three Decades of Changes in the Tropical Andean Glaciers Using Landsat Data Processed in the Earth Engine" Remote Sensing 14, no. 9: 1974. https://doi.org/10.3390/rs14091974

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop