A Comparative Analysis of Machine Learning Techniques for National Glacier Mapping: Evaluating Performance Through Spatial Cross-Validation in Perú.

1. Introduction

Tropical glaciers are sensitive indicators of climate change [1] and vital dry season sources for drinking water, agriculture, and livelihood of many people that depends upon them [2].

Tropical Andean glaciers are among the fastest shrinking and largest contributors to sea level rise on Earth. Over the period from 1975 to 2020 Southern Peruvian Andes have receded by ∼32 % and are now at less than half their original size [3,4]. Most recent studies reported a glacier recession in the Cordillera Blanca of −46 % between 1930 and 2016 [5]. This has bought effects in spatiotemporal alterations in both quantity and quality of mountain water resources [2,6,7] specially this reliance increases sharply during drought conditions [8].

Accurate extraction of glacier areas and continuous monitoring of glacier morphology are fundamental prerequisites for glacier research, as they provide crucial information for assessing hydrologic risks and facilitating climate change mitigation efforts [9].

For example, recent studies focusing on hydrologic modeling of glacier contribution to watersheds in the tropical Andes have utilized multi-temporal estimates of glacier area derived from remotely sensed data to describe the effect that glacier change has had and will have on water resources for this region [10,11].

Due to the arduous, time-consuming, and subjective nature of manually delineating glacier boundaries, numerous techniques have been devised to automatically delineate glacier outlines using primarily multispectral imagery [1], this technique has been proven to generate equivalent accuracy compared to manual digitization techniques when a large sample of glaciers is analyzed [12]. Several studies have been carried out to map glaciated areas on a local and regional scale in Peru using remote sensing techniques [1,6,11,13] as well as in continental scale [9].

The prevailing and currently operational technique employed for delineating debris-free glaciers in the Peruvian Andes relies on the utilization of a Normalized Difference Snow Index (NDSI) threshold [12,14]. However, it is important to note that different threshold values can yield varying results [15]. Although the threshold limit may need to be adjusted for specific regions [10,13], it is commonly kept constant at approximately 0.35 to 0.55 in the literature [1].

Machine learning (ML) is commonly used to build geospatial prediction models [16,17,18,19,20] and has been universally employed in geoscientific research such as global soil properties mapping [21,22,23], landslide susceptibility [24], climate change studies [25], and wildfire risk mapping. Specially several machine learning approaches have been suggested for glacier mapping studies [26,27,28,29]. To be sure that observed glacier changes are related to real changes rather than caused by imprecise determination of the outline, the accuracy of the outlines must be known [30]. This measure is also appropriate to assess the significance of any derived relative changes in glacier size.

Geospatial prediction of environmental variables often relies on map accuracy assessment [31,32], classical map accuracy assessment is rooted in sampling theory [33,34] in which an unbiased estimate of map accuracy is obtained.

Several approaches can be chosen in order to accomplish accuracy assessment [35]. A classical procedure is to evaluate model performance and associated errors by randomly selecting a number of test observations that are set aside at the model calibration stage and only used to quantify model prediction error [36]. While the best method is simply using a completely independent test sample, this is not always feasible [33,34].

A possible solution is K-fold cross validation (K-CV)[37,38,39]. K-CV is a resampling-based technique for the estimation of a model's predictive performance [38]. The basic idea behind K-CV is to split an existing dataset into training and test sets using a user-defined number of partitions. First, the dataset is divided into k partitions or folds. The training set consists of k − 1 partitions and the test set of the remaining partition. The model is trained on the training set and evaluated on the test partition. A repetition consists of k iterations for which every time a model is trained on the training set and evaluated on the test set.

Each partition serves as a test set once. The procedure is repeated k-times and finishes when each unique subset has been used for testing once. A further improvement of this approach is the repetition of the K-CV: the whole procedure is repeated, producing a set of random samples each time. This approach is called the repeated k-fold CV method (RK-CV).

In environmental sciences, observations are often spatially dependent [24,40]. Subsequently, they are affected by underlying spatial autocorrelation by a varying magnitude. Although cross-validation is a particularly well-established approach for model assessment of supervised machine-learning models [39,41,42,43] it is generally agreed that most ML methods in spatial applications do not consider relative location and neighborhood features and that they analyze pixels regardless of their surroundings [44], and hence completely ignoring the spatial dependent and heterogeneity of spatial processes [16]. Therefore, the direct application of ML to geospatial data without accounting for the potential spatial autocorrelation could lead to biased outcomes. It is clear from many studies [40,44,45] that unaddressed spatial autocorrelation generates problems, such as overoptimistic fit of models, omitted information and/or biased (suboptimal) and therefore model performance estimates from conventional cross-validation leads to optimistically biased estimates of map accuracy [20,36,43] due to the similarity of training and test data in a non-spatial partitioning setup when using any kind of cross-validation for tuning or validation [36,43,46].

Several papers including [36,47] have investigated the over-optimistic accuracy estimates by promoting SP-CV methods. These methods start from the premise that spatial proximity of data points in the calibration and test data folds is to be avoided. This is commonly achieved by spatial blocking in K-CV [39,47]. The general overview from the literature [44] is that visible progress has been made recently in the development of spatial machine learning modelling. The current standard for modeling is to use a classic ML algorithm, and some sort of spatial cross-validation method (SP-CV). Therefore, cross-validation approaches that adapt to this problem should be used in any kind of performance evaluation when spatial data is involved [48,49].

A common approach is block-CV (BLK-CV) which divides all samples into contiguous blocks, and then avoids the selection of samples within the same block as both training and validation samples [47]. As an example [36,43] used K-Means clustering to split samples into five folds based on sample data’s locations. The validity of these procedure was empirically tested by [36,40,43]. They found that BLK-CV better approximate the error obtained when predicting species distribution and above ground forest biomass (AGB).

Overall, the intention of this work is to demonstrate the need for spatial cross validation when using machine learning in geospatial classification of glacier land cover when the aim is to estimate a biased-reduced predictive performance. The following objectives (and hypotheses) are addressed:

• The study compares the predictive performance of machine learning algorithms: k-nearest neighbors, Random Forest, Gradient Boosting Machines, and classical statistic models: logistic regression in glacier mapping. It is expected that the machine learning algorithms will demonstrate significantly higher predictive performance.
• The study also examines the predictive performance of classification algorithms when spatial and non-spatial cross-validation methods are used. It is hypothesized that non-spatial partitioning methods will yield over-optimistic results due to the presence of spatial autocorrelation.
• Additionally, the study investigates the effects of spatial clustering on the distribution of errors in the analyzed algorithms for glacier mapping. It is anticipated that non-spatial models will exhibit spatial autocorrelation in their errors.

2. Materials and Methods

2.1. Study Area

Peru encompasses the largest concentration of tropical glaciers worldwide. About 70 % of all tropical glaciers, covering an area of 1602.96 km2 are located there [50]. According to [51] the glacierized areas of the Peruvian Andes can be divided into three subregions based on their temperature, precipitation, and humidity characteristics:

• R1 is the northern wet outer tropics, with a high mean annual humidity of 71 %, nearly no seasonality of the temperature and a total annual precipitation of 815 mm. R1 includes the Cordillera Blanca, Central, Chonta, Huagoruncho, Huallanca, Huayhuash, Huaytapallana, La viuda and Raura.
• R2 is the southern wet outer tropics, with moderate mean annual humidity of 59 %, an annual variability of the mean monthly temperature of about 4 ◦C and a total annual precipitation of 723 mm. R2 includes the Cordillera Vilcabamba, Urubamba, Vilcanota, Carabaya and Apolobamba.
• R3 is the dry outer tropics, with low mean annual humidity of 50 %, a mean annual temperature of −4.0 ◦C and low total annual precipitation of 287 mm. R3 includes the Cordillera Ampato and Coropuna.

The three subregions experience a dry season from May to September, coinciding with the austral winter, and a wet season from October to April, corresponding to the austral summer [51]. During the wet season, the glaciers predominantly accumulate mass, while the lower portions of the glaciers experience ablation consistently throughout the year.

This study focusses in R1 and R2 glacier subregions (Figure 1), excluding glaciers form R3 subregion, hence reducing the number of Cordilleras analyzed. From the entire National Glacier Inventory (hereafter NGI) [52], we identified 9 glacier Cordilleras as show in Table 1. Nevertheless, in order to ensure the broad applicability of our findings, we carefully select Cordilleras that provide a comprehensive representation of the glacier distribution in Peru, particularly all selected study areas span proximally 90% of the glaciated surface of Cordilleras of Perú [52]. Blanca (40.20%), Central (3.80%), Huallanca (0.47%), Huayhuash (4.75%), Huaytapallana (1.92%), Raura (2.29%), Urubamba (2.11%), Vilcabamba (9.05%) and Vilcanota (22.96%).

Additionally, the majority of selected glaciers have been studied in term of hydrology, ablation and dynamics, especially the Cordillera Blanca and Cordillera Vilcanota have been intensely studied the in recent years [3,11,50,53,54].

2.2. Data Acquisition

2.2.1. Glacier Inventory

The analysis was conducted using the National Glacier Inventory (NGI) dataset provided by [52] as the reference data. This dataset consisted of polygonal vector data representing glacier areas delineated using mostly Sentinel 2 and Landsat data and high-resolution Google Earth images with field survey studies.

NGI is a shape file of the 2017-2018 glacier outlines for the main glaciated mountain regions in the Peruvian Andes, covering the entirety of Peru with information concerning their code name, date of acquisition and remote sensing sensor used for delineation. To prepare the data for classification, the glacier outlines from the NGI were used to generate a binary raster mask I(x) of 30 m grid size with values 0/1 corresponding to non-glacier/glacier pixels which is achieved after Equation (1):

$$I\left(\mathit{x}\right)=\left\{\begin{array}{c}1, if pixel belongs to glacier class \\ 0, otherwise.\end{array}\right.$$

(1)

Finally, the labels were projected to Universal Transverse Mercator (UTM) projection (UTM Zone 18S for subregion R1 and UTM Zone 19S for subregion R2). This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

2.2.2. Landsat Data and Processing

Landsat data corresponding to the period 2017–2018 were obtained from Landsat Collection 2 Level 2 and Tier 1 surface reflectance (SR) products [55] available online: https://www.usgs.gov/landsat-missions/landsat-collection-2-surface-reflectance (accessed on 1 July 2023). The Landsat image collections used were accessed and processed in the Google Earth Engine (GEE) platform [56]. For each selected glacier sub-region, a monthly composite was created using atmospherically corrected and topographic calibrated Landsat 8 OLI reflectance data imagery from the Tier 1 LANDSAT/LC08/C02/T1_L2 collection. The Landsat image collections used were accessed and processed in the Google Earth Engine (GEE) platform [56].

Table 1 provides an overview of the spectral characteristics of the OLI bands of Landsat-8. A more in-depth explanation of the Landsat Data processing steps is given below.

All available images from 2017 with a cloud cover less than 80% during the dry season (June and July) were filtered and selected. This period corresponds to the end of the ablation period, minimizing the potential confusion caused by transient snow cover, which could degrade the quality of the data and hinder accurate glacier discrimination [11,57]. In our study case, the implementation gets least cloudy scenes during the desired month in the dry season and nearby dates.

For the cloud masking we used the Band Quality Assessment (QA_PIXEL) information available in the Landsat Collection 2 Level 2 and Tier 1 surface reflectance (SR) images. Then we applied a spectral transformation by scale and offset parameters based on [58]. We created annual image composites for 2017 year using GEE’s image-reducing methods, having as an input all the cloud-masked Landsat scenes available in GEE’s Landsat collection that met our filter criteria. Finally, we got an annual Landsat mosaic composite conformed of Blue, Green, Red, NIR, SWIR1, SWIR2 reflectance bands for each AOI.

2.2.3. Normalized Difference Snow Index (NDSI)

The Normalized Difference Snow Index (NDSI) is a widely used index to separate snow from other land coverages, and therefore it is also useful to identify glacier coverage. Subsequently, NDSI has been extensively used for glacier mapping [6,9,10,11,59]. We computed NDSI values from each reflectance composite using Equation (2):

$$NDSI=({\mathsf{\rho }}_{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}-{\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}})/ ({\rho }_{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}+{\rho }_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}),$$

(2)

wheree ${\mathsf{\rho }}_{\mathrm{G}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$is the surface reflectance in the green band and ${\mathsf{\rho }}_{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}}$ is the surface reflectance in the SWIR band. The use of NDSI as only mean of glacier detection is a common approach in glacier mapping generally with acceptable results [57,60,61]. In this work, we utilized the Normalized Difference Snow Index (NDSI) both as a covariate and as an independent method for delineating glacier outlines.

2.2.4. Digital Elevation Model (DEM)

The significance of topographic characteristics in glacier classification has been demonstrated, as the distinction between debris-covered and non-covered glaciers based solely on reflectance properties is practically impossible [12,30,59]. Accordingly, 10 topographic parameters including elevation, slope, aspect, profile curvature, plan curvature, longitudinal curvature, cross-sectional curvature, maximum curvature and minimum curvature were generated using SAGA GIS [62] using the 30-meter resolution JAXA’S ALOS WORLD 3D DEM downloaded from https://opentopography.org/ through the R package elevatr (https://github.com/jhollist/elevatr) for each glacier region under study.

2.2.5. Pixel Sampling and Classification Matrix Assembly

We build a label-predictors classification matrix by combining the NGI raster mask, Landsat surface reflectance derived products and topographic parameters as is explained in the following section.

Although it is common approach in glacier mapping to use all available pixel in a given reflectance scene [9,57], processing each pixel might be redundant for regional or national scale mapping specially with high spatial resolution (i.e.,30 meters) images and specially in an experiment-benchmark based studies. As such, computing resources were prioritized to process only a sample of pixels in each Cordillera.

Therefore, a stratified random sampling approach (SRS) [40,63] was employed to select a balanced sampling of pixels for each glacier region under study, hence the training samples were a subset of the available pixels within each scene, this approach is common in glacier mapping via machine learning [29,64]. Finally, all covariates and label rasters were harmonized to UTM Zone 18S and UTM Zone 19S glacier Cordilleras within subregion R1 an R2 respectively. This was done before cropping, resampling to 30 meters of spatial resolution and stacking. After preparing the data-cubes, the complete stacks of label-covariate rasters were discretized to 5,000 training pixels over all studied glacier Cordilleras (Table 2). With this harmonized dataset, each glacier area was processed one by one in the main workflow, allowing the implementation of modeling methods for individual study areas. Overall, all Landsat mosaic composite images cover a total area of 42 135 km².

2.3. Machine Learning Classifiers

A large variety of machine learning (ML) methods have been used for spatial prediction and more specially in glacier mapping [28,42,64,65]. However, the focus of this paper is on the evaluation of CV methods. Thus, similar to previous cross validation comparison studies [36,43,45], we selected just a small amount of ML methods for our experiments: Random Forest (RF), Gradient Boosted Machines (GMB) and Weighted k-Nearest Neighbors (KKN), additionally a statistical model, Logistic Regression (LG). A detailed explanation of each model is outside of the scope of this study, but a brief summary is given for each model with relevant references for each.

2.3.1. Logistic Regression

The binary logistic regression takes the dependent variable in the form of binary data as presence (1) and absence (0) and linearly relates to the independent variable(s) [66]. Using an exponential function (sigmoid function), it calculates the probability of each input sample for classification. The multinomial binary logistic regression can be written as follows:

$$Z=\mathrm{a}+\sum _{i=1}^n{b}_i{x}_i$$

(3)

$$P\left(Z\right)=\frac{1}{1+e^{-Z}},$$

(4)

where P(Z) is the probability; Z is a parameter; $\mathrm{a}$ is the intercept and $b_i$’s are the coefficients for independent variables $x_i$’s; the i index is for each covariate. Usually, the probability value of 0.5 is used to categorize data of either presence or absence class, and thereby compute the classification metrics (i.e., accuracy).

2.3.2. k-Nearest Neighbors

The Weighted k-Nearest Neighbors (KNN) algorithm is a supervised learning method which uses kernel functions to weight the neighbors according to their distances [67]. KNN classifies data points based on the closest k training samples in the feature space. It considered as non -parametric because KNN makes no statistical assumption about the data. KNN is a common classification algorithm used in remote sensing data mining applications and it has been widely used for mapping glacier surfaces [64,68].

Clustering with the k-means algorithm has the significant advantages of ease of interpretation, a high degree of flexibility and computational efficiency; however, its main disadvantage lies in the need to specify a priori the number of k clusters (Kopczewska, 2022). Moreover, three hyper-parameters have to be set: the number of nearest neighbors (k), the distance and the kernel are needed to set before training a KNN. In this research we implemented KNN using the knnn package for R (https://github.com/KlausVigo/kknn)

2.3.3. Random Forest

Random forest (RF) belongs to the group of ensemble learners, which builds upon a large number of basic model structures called Decision Trees (DT) that are trained in parallel. RF and merges the results of these DT together to get a more accurate and stable prediction compared to a single DT [69]. For each DT, a random subset of the training sample and a random subset of the classification features are used [69], hence each DT is trained independently.

Currently RF is the most common machine learning method used for geospatial predictions and is widely acknowledge for its general accuracy and successful applications in diverse geoscientific problems [17,19,22,49,70,71,72,73,74].

To train an RF model, some hyper-parameters should be selected, among others, the number of individual decision trees (ntree) and the number of features selected at each split of the trees (mtry). In this study the number of trees is kept at a moderate size of 100 to 500 trees for a balance between computational efficiency and predictive stability [75].

2.3.4. Gradient Boosting Machines

Gradient Boosting Machines (GBM) were first presented in [76], is another type of model based of ensemble learners. To generate the final prediction results, GBM could use weak learners in a sequential learning process, in the form of an ensemble of weak predictions such as DT. Unlike RF, in GBM Decision Tree learners (DT) are trained sequentially. GBM is a useful practical tool for different predictive tasks, and it can consistently provide more accurate results than the conventional single machine learning models [77,78,79]. A downside of GBM compared to RF, is that GBM has more hyperparameters that need to be calibrated.

2.4. Implementation

All modeling analyses were conducted using the R programming language, an open-source statistical software [80]. The algorithm implementations from several packages were employed, including gbm (https://github.com/harrysouthworth/gbm) for boosted regression trees [76], kknn [81] for weighted k-nearest neighbor classification (KNN), ranger [82] for random forest (RF) modeling based on [69].

2.5. Model Evaluation Metrics

2.5.1. Matthews Correlation Coefficient (MCC)

Evaluating binary classifications is a pivotal task in statistics and machine learning, because it can influence decisions in multiple areas.

To evaluate the accuracy of the models in binary classification problems, we employ the Matthews correlation coefficient (MCC), Equation 3, as an established and robust error measurement metric [83]. The MCC is derived from Cramér's V and is applied to a 2 × 2 traditional confusion matrix, consisting of true positives (TP), true negatives (TN), false negatives (FN), and false positives (FP):

$$\mathrm{M}\mathrm{C}\mathrm{C}=\frac{TP\times TN-FP\times FN}{\sqrt{\left(TP+FN\right)\times \left(TP+FN\right)\times \left(TN+FP\right)\times \left(TN+FN\right)}}$$

(5)

Given the number of experiments we ran in each glacier region, we focused solely on MCC for sequent comparative analysis because this was proved to be a balanced, more reliable summary of binary classification error than accuracy, F1 score and Kappa[83,84]. We choose not to use Kappa coefficient as is has proven to not to be a reliable measure of accurate classification and is difficult to interpret [85].

The Matthews correlation coefficient instead, generates a high score only if the classifier correctly predicted most of the positive data instances and most of the negative data instances, and if most of its positive predictions and most of its negative predictions are correct. A high Matthews correlation coefficient (close to +1) means always high values for all the four basic rates of the confusion matrix: true positive rate (TPR), true negative rate (TNR), positive predictive value (PPV), and negative predictive value (NPV).

It should be emphasized that in the classification results presented in this study, no manual modification of glacier outlines/areas was conducted, contrary to common practice. Although such modifications are typically carried out to improve classification accuracy, evaluating the effectiveness of these procedures is challenging due to their irreproducibility [12]. Moreover, attributing the accuracy results solely to manual corrections would undermine the primary objective of this study, which is to assess the accuracy of the algorithm and cross-validation method itself. This was important in our work as we needed to automatically analyze and compare different algorithm-dataset combinations.

2.5.2. Moran’s I

Another way to assess the quality of a model applied on geographical data is that degree of spatial autocorrelation of errors [20]. Preferably the spatial autocorrelation of errors should be minimized or even eliminated, which would imply that the model performs similarly in space and that there are minimal (or no) subregions with strong patterns of over- or underestimated values [45,49]. The Moran’s I metric (MI) can be used for detect and quantify global spatial autocorrelation [86]. Therefore, Moran’s I were estimate from the absolute error of each model averaged for all repetitions for each glacier area.

MI varies from −1 to +1. A positive value means similar values happen in close proximity, while a negative value shows that dissimilar values are spatially grouped close to each other. A value close to zero shows no spatial autocorrelation for the variable, indicative of a spatially random process; in the latter case, the assumption of independence, essential for many statistical methods, is met [73].

In all cases Moran’s I was tested under Monte Carlo simulation and were derived from nearest neighbor spatial weights matrices based on k = 5 nearest neighbors using the R library spdep [87].

2.5.3. Spatial Autocorrelation of Glacier Class

Although cross-validation is a particularly well-established approach for model assessment of supervised machine-learning models [39,41,42,43] a particular issue is that spatial autocorrelation in the data can invalidate model validation approaches because an observation point cannot serve as spatially independent validation of nearby training data points [45]. Even so when the data points are spatially clustered, conventional cross-validation leads to optimistically biased estimates of map accuracy [45,49]. Therefore, when machine learning models are directly applied to spatial data we need to considering if spatial autocorrelation is present and to what extent [36].

To better understand the effect of spatial autocorrelation in our validation approaches, we test for spatial autocorrelation on class labels using indicator variogram analysis (Equation 6). The indicator semivariogram ${\gamma }_I$ can be inferred from the indicator class data using the estimator [32,88,89,90,91]:

$${\gamma }_I\left(U_k;\mathit{h}\right)=\frac{1}{2n\left(h\right)}\sum _{i=i}^{n\left(\mathit{h}\right)}{\left[i\left(U_k;\mathcal{x}\right)-i\left(U_k;\mathcal{x}+\mathit{h}\pm ∆\mathit{h}\right)\right]}^2;k=1,\dots ,m-1$$

(6)

where $\left\{U_1,\dots ,U_m\right\}$possible clases exist (in our case two classes for glacier and non- glacier surface, i.e., m = 2), $n\left(\mathit{h}\right)$ is the number of pairs of indicator data that are a distance $\mathit{h}\pm ∆\mathit{h}$ apart, finally the indicator random function or indicator transform $i\left(U_k;\mathcal{x}\right)$ of N ($\mathcal{x}$) is defined using Equation 7:

$$I\left(U_k;\mathit{h}\right)=\left\{\begin{array}{c}1 if N\left(\mathcal{x}\right)\in \left\{U_1,\dots ,U_k\right\} \\ 0 if N\left(\mathcal{x}\right)\in \left\{U_{k+1},\dots ,U_m\right\}\end{array}\right.;k=1,\dots ,m-1.$$

(7)

where $\mathrm{N}\left(\mathcal{x}\right)$ denotes the value of the random variable at pixel $\mathcal{x}$. Since we already have binary class labels, we can readily apply equation 3 without apply the indicator transform in equation 2.

In order to model the indicator semivariograms for each glacier dataset, we employed various model specifications including the exponential, spherical, Gaussian, and Matern functions. Variograms were fitted by weighted least squares using $N_{\mathrm{j}}/h_j^2$as weights, where $N_{\mathrm{j}}$denotes number of points pairs in the j-th lag and $h_j^2$ is the corresponding lag distance. The model that yielded the smallest residual sum of squares when compared to the sample variogram was selected as the final model. The variogram models were implemented using the gstat package for R [92,93].

2.6. Spatial Vs. Non-Spatial CV Comparison Experiment

Our experimental benchmark (illustrated in Figure 2) is mainly composed of 4 steps. Step 1 deals with the construction of training data and prediction locations from the entire dataset. In step 2 the prediction locations are applied to calculate the reference prediction error ${MCC}_{ref}$. This reference error is used to check which CV method can provide more accurate and unbiased predictions. In step 3 we implement both NSP-SP and SP-CV methods on the training dataset. This implementation provides estimated prediction errors ${MCC}_{CV}$ of each CV method. In step 4 we calculate the absolute difference ($d_{CV}$) of the two error values, ${MCC}_{ref}$ and ${MCC}_{CV}.$

2.6.1. Step 1: Construction of Sample Data and Prediction Locations

In order to assess the accuracy of different CV methods, it is crucial to incorporate an unbiased reference prediction error is an indispensable element to include in the benchmarks. Therefore, CV methods cannot be implemented directly in the entire dataset, the dataset should be divided into two parts, one containing training data to implement CV methods, and the other containing prediction locations (test data) to provide a standard quasi real prediction error i.e., (${MCC}_{ref}$). Therefore, the initial stage of the experiments involves the creation of sample datasets and prediction datasets for each Cordillera [46].

For the test data, we employed simple random sampling to generate 1000 independent sampled pixels for each glacier region being investigated, which served as the prediction locations. In contrast, the training samples were selected using a designed stratified sampling approach to ensure a balanced representation of the two classes [33,94]. This was done because it has been proven that class imbalanced in machine learning model can produce highly overoptimistic biased result [95].

2.6.2. Step 2: Calculate Reference Prediction Error

Since randomly independent prediction locations are directly used as the test set, the calculated value of reference prediction error ${MCC}_{ref}$ can completely depict all model’s unbiased performances [34,40]. First, all the available training samples are used to build the benchmark models. Then, these models are used to predict in the independent test locations to obtain the reference error ${MCC}_{ref}$based on the prediction of the machine learning models at the prediction locations and the true class values. As mentioned earlier, the incorporation of the reference prediction error enables us to assess and compare the potential biases of machine learning models in the context of glacier classification.

2.6.3. Step 3: Calculate the Prediction Error of Each CV Method

In this stage we employed two cross-validation strategies to assess the predictive power of our models and generate the prediction error (${MCC}_{CV}$). The first strategy involved a standard K-fold cross-validation approach (NSP-CV), where observations were randomly divided into V subsets, ignoring any structure of spatial dependence in the data. The models were then trained iteratively on V-1 sets, with each subset V serving as a test set in a different iteration. The predicted vector of glacier and non-glacier land cover classes was used to calculate cross-validation error statistics, we focused on Matthew correlation coefficient (MCC). In our study, we used K = 5, representing five-fold cross-validation.

The second strategy, known as spatial K-fold CV (SP-CV), differs from random K-fold CV (NSP-CV) by considering the spatial structure of the data when partitioning observations into subsets. The objective is to group observations into spatially homogeneous clusters that are larger in size than the range of spatial autocorrelation, thereby achieving independence between cross-validation folds. In this work, we adopted the spatial cross-validation approach proposed by [47], and utilized by [43] which utilizes k-means clustering to mitigate the influence of spatial autocorrelation. In contrast to non-spatial cross-validation, spatial cross-validation partitions the data into spatially disjoint subsets, reducing the impact of spatial autocorrelation.

For the second CV scenario (i.e., SP-CV) we developed novel code to accommodate the specific characteristics of our dataset and research objectives. Since k-fold cross-validation is commonly used and well-known among machine learning practitioners [41,96], we provided a detailed explanation of the spatial cluster-based implementation instead in Figure 3:

In all experiments, k was set to 5. Both 10 and 5 are the most commonly used values in CV [36]. Then, each CV method was implemented 50 times (i.e., 50 repeated five-fold partitioning setting was chosen for performance estimation of ${MCC}_{CV}$ for each model). Thus, 250 models were fitted and tested at each glacier region and the average ${MCC}_{CV}$ was derived to account for random errors [36].

2.6.4. Step 4: Calculate CV Method’s Prediction Error Difference

After step 2 and step 3, we obtain the reference prediction error ${MCC}_{ref}$ and the evaluation result of every CV method (${MCC}_{CV}$). By comparing them, we can find out which CV method produces more realistic results. For this purpose, we use the CV method’s prediction error difference ($d_{CV}$) as a quantitative metric. We calculate by subtracting from and getting the absolute value ($d_{CV}=\left|{MCC}_{CV}-{MCC}_{ref}\right|$). When is closer to zero, the corresponding CV method’s performance is considered better.

2.7. Satistical Comparison of Model Results

In order to determine the influence of the two different validation scenarios on model performance, a paired t-test was carried out to see if there were significant differences in the results between models. The analysis was based on modified paired t-test [97] considering as factor the model and validation strategy used and as dependent variable the classification performance measure (i.e., MCC). The null hypothesis establishes that average differences between models’ results (i.e., MCC) are negligible, then both cross-validation approaches behave equally. For analyzing the results, we examined the p-values. Prior to applying the t-test, a Shapiro-Wilk test was conducted to confirm the normal distribution of the experimental results, as required by the test as in [98].

3. Results and Discussion

3.1. Spatial Autocorrelation

Semivariograms were calculated for the glacier class indicator variable for each glacier area under study to confirm the presence of spatial autocorrelation in the datasets.

Figure 4 shows the experimental and fitted indicator variograms for each glacier data set for 1000 test pixels, the semivariograms were rescaled by variance of each indicator variable to facilitate comparison across datasets. The parameters of the indicator variogram models of the 9 data sets are given in Table 3.

In general, it is clear that the glacier class sampled at 30 m resolution presents a significant spatial correlation, the smallest range of spatial dependence was found for cordilleras Central, Urubamba and Huallanca with ranges from 0.37, 0.68 and 0.87 Km respectively, and the largest range was found for cordillera Blanca (5.42 Km).

For some glacier regions the indicator semivariograms display a very large nugget effect and a short-range structure (i.e., Central, Huaytapallana and Urubamba), as well as bigger kappa parameters, indicating a fairly constant spatial process in those cases. But in general, relative larger ranges and smaller nugget effects occur for the majority of other regions. For instance, in the cordillera Blanca case, the indicator variogram shows that at a 30 m resolution km spatial resolution, glacier class presents a significant spatial correlation up to 5.42 km (Figure 4). This spatial autocorrelation can notably be observed in almost all cases, where clusters of homogeneous glacier class values are present. The cordillera Blanca region exhibited a larger range due to its considerably larger area (13,963 Km2). The absence of similar range behavior in other areas can be attributed to the substantial differences in the overall size and extent of each individual region because glaciers assume a size and flow rate that are in balance with local climate [15].

Given the relatively high sampling intensity (and resulting proximity) of glacier class pixels selected for this work, and the long range of spatial autocorrelation in the data (i.e., 5 km), it is obvious that any given randomly selected test pixel will not be independent from a large number of neighboring pixels, thereby violating the core hypothesis of model validation (i.e., the independence between training and test sets). This result probably doesn’t hold for the regions for which variograms revealed poor spatial structure like Central, Huaytapallana and Urubamba Cordilleras. This suggest that a spatial cross validation method could possibly be useful in this context [36].

It is worth to note that no stationarity assumptions are needed since the indicator variograms are simply used as a means to describe the spatial structure of the data, rather than to perform a model based inferential spatial prediction of a spatial process.

3.2. Model Specification

For the benchmark, all models were trained using either default hyperparameter values and/or recommendations from the literature that were specifically tailored to our case data. A detailed summary of the models and hyperparameters is shown in the Table 4.

3.3. Spatial Vs. Non-Spatial CV Comparison Experiment

Figure 5 shows the results of our experimental benchmark. It is possible to identify the models that exhibit superior and inferior performance in both scenarios. In NSP-CV settings, K-nearest neighbors (KNN) was the overall best model across regions consistently demonstrating the highest performance followed by logistic regression (LR) and random forest (RF). The SP-CV settings, generally shown lower MCC values, K-nearest neighbors (KNN) remained as the overall best model in almost all regions, followed by gradient boosting machines (GBM) and logistic regression (LR), in Cordillera Blanca and Central respectively. Some studies have found that KNN shows the best predictive performance in spatial settings [81,99] although this is generally not the case [43,100,101,102].

Although SP-CV generally shows lower MCC results than NSP-CV, this is not always the case, especially in Cordillera Central, Raura, Urubamba, and Vilcanota. We hypothesize that these results are connected with the degree of spatial autocorrelation/clustering of the glacier class in these specific Cordilleras. For example, Cordilleras Central and Urubamba present almost a pure nugget effect, which could have generated the mixed response of MCC to both validation approaches. Hence, it is impossible to distinguish the effect of SP-CV from NSP-CV in the presence of weak spatial autocorrelation.

Surprisingly logistic regression (LR) demonstrates good performance in some cases in this quasi-linear problem, without clear evidence of being surpassed by models such as random forest (RF) and boosted trees (GBM) in some cases, as commonly acknowledged in other geospatial contexts [17,43,74]. This result shows the importance of traditional parametric approaches in spatial modeling which make this algorithm a valid choice, especially if the differences in predictive accuracy compared to black-box models are small

Concerning the poor performance of random forest (RF), the literature generally agrees upon its general applicability in geospatial contexts as a "go-to" model [17,43,60,74,103]. RF uses "bagging" (bootstrap aggregation). As spatial data is correlated, this resampling violates the assumption of independent and identically distributed (i.i.d.) data units, which is fundamental to bootstrapping. [104] provided evidence that these limitations could lead to inferior prediction performance of RF under spatial dependence, and this could be the reason for the observed performance of RF.

There were notable differences in the distribution of the estimated MCC between the spatial and nonspatial settings (i.e., ${MCC}_{SP-CV}$ and ${MCC}_{NSP-CV}$). Figure 6 shows the results for both CV methods for each glacier region under study along the reference error ${MCC}_{ref}$ in red doted horizontal lines.

First and foremost, it can be seen than almost in all cases the SP-CV and NSP-CV methods produced quite different results in terms of MCC, it is also remarkable that the proposed SP-CV, which consider the clustered nature of the data, produce closer evaluation results to the reference prediction MCC than NPS-CV in almost all cases. This suggests that SP-CV may produce bias reduced spatial predictions when using ML models for glacier mapping, especially when training locations area scare and highly clustered as is usual in glacier monitoring and mapping. As expected, the SP-CV results shown high variance for all models.

Upon careful examination of Figure 6, Cordillera Central, Urubamba, and Vilcabamba, it can be found quite high differences between CV-MCCs and the test MCC (dotted red lines in Figure 6), regardless of the validation approach. It seems that, in those particular cases, both CV methods are incapable of estimating the true error of the models. This indicates the presence of biased results, even when employing the spatial cross-validation approach suggested in this study. Although there could be multiple reasons to explain these results, we hypothesize that they are due to possible overfitting of the models in those specific areas.

A more detailed analysis can be found in Table 5 that shows the mean MCC grouped by Cordillera and CV method as well the difference regard to the test reference MCC (in parenthesis). For example, Cordillera Blanca mean MCC for NSP-CV models is 0.949, while mean MCC for SP-CV models is 0.928, and their difference with respect to the reference independent test (i.e., ${MCC}_{ref}$ 0.844) are 0.1054 and 0.0845 respectively.

Overall, the spatial 5- fold CV led to a slightly decline in model’s MCC (i.e., about - 4%) respect to the non-spatial k fold method. Likewise, the SP-CV yields a sharp decline in biases $(d_{CV}$). compared to the reference MCC evaluated at the test-sets (i.e., about - 18%).

Furthermore, SP-CV’s evaluation results were much closer to the reference prediction error than the results of NSP-CV for all glacier regions under study. This shows that the proposed method, which considers the spatial structure of the data in the evaluation of the model could indeed provide a reasonable unbiased result. To further strengthen this observation, it is important to analyze the individual model’s discrepancy between the SP-CV and NSP-CV MCC estimates obtained in the benchmark’s settings, which will be explored in the subsequent section.

3.4. Satistical Comparison of Model Results

We carried out a comparative statistical analysis of the results using non-parametric statistical tests to compare the performances of all the considered models. Table 6 shows the t-statistics and the p-values (in parenthesis) of paired t-test comparison for the MCC results obtained after the 50-repeteated 5-fold CV using the two types of cross-validation for each glacier region (p-values < 0.05 means that there is a significance statistical difference between cross-validated MCC results and p-value > 0.05 means that both cross-validation approaches yield equal results).

Based on the statistical analysis, considering the results of the paired t-test, it is not definitively conclusive to assert that spatial cross-validation (SP-CV) yields significantly different Matthew’s correlation coefficient (MCC) estimates compared to non-spatial cross-validation (NPS-CV) in the majority of cases with exception of remarkable cases that need further clarification:

For Cordillera Blanca only GBM and KNN models showed significant differences (p-value < 0.05). For Cordillera Huallanca, Huayhuash and Vilcabamba RF was the only model that showed significant differences. For Cordillera Huaytapallana only GBM showed significant differences. For Cordillera Urubamba only LR showed significant differences.

Some algorithms showed no sensitivity to the cross-validation method (i.e., KNN and LR). Since KNN is the best model overall in all cases regardless the CV method, this suggests that KNN can produce more reliably estimates of error, and at the same time the best performing predictions. On the other hand, some algorithms were quite sensible to the cross-validation method. RF showed, especially RF showed significant differences between SP-CV MCC’s and NP-CV MCC’s. In all those cases, the null hypothesis establishes that average differences between error’s models (i.e., MCC) are non-negligible. Therefore, this difference can be attributed to an overoptimistic bias in nonspatial performance estimates in the presence of spatial autocorrelation [36,43,46,72,74].

However, it should be noted that these results may slightly vary depending on the data distribution within each fold and the randomization in each repetition, which can impact the comparison tests [43]. Nonetheless, given the considerable number of repetitions in the experiments, these results are generally robust.

While our results align with previous studies, such as those by [43,105] and [45] suggesting that non-spatial performance estimates tend to be significantly “superior´’ to spatial performance estimates, it is important to note that our experiments do not allow for such a definitive statement in all cases. The observed differences in performance between spatial and non-spatial estimates in our study may not be as pronounced or statistically significant regardless of the classification algorithm.

3.5. Spatial Autocorrelation of Errors

Is traditionally acknowledge that when applying standard ML models in a spatial prediction setting, errors could remain depended in space due to failing to account for spatial dependency or heterogeneity in important explanatory variables or in model architecture [72,106,107]

Correspondences between predicted and label class data were indicator coded. If the class of the test sampled pixel matched with that of the prediction class, an indicator code 1 was assigned to that sample pixel. Conversely, an indicator code 0 was given to sites where the predicted class differed from the test class. These values were identified as classification errors and obtained through the 50 repeated spatial 5-fold either spatial and non-spatial cross-validation (i.e., SP-CV and NSP-CV) and then aggregated through a simple majority vote approach for all repetitions. Next, we analyzed spatial autocorrelation of the indicator-coded data using Moran's I.

Figure 7 shows that regardless of the type of modeling algorithm, there were low MI values with statistically non-significant p-values, (i.e., p values greater to 0.05). At first it seems that SP-CV could possibly reduce the spatial dependence of errors compared to NSP-CV which shows generally higher MI values (i.e., spatial correlated errors). Thus, suggesting that in our SP-CV validation method, spatial structure of errors is completely eliminated in all glacier regions but not in NSP-CV. But a closer inspection of p-values (not shown in figure) suggest that the Null hypothesis should not be rejected and the spatial error patterns are almost random in both cross-validation scenarios.

Urubamba was the only exception in terms of significant MI p-values, especially when KNN model was applied, indicating the occurrence of spatial correlated errors. Overall, after either SP-CV and NSP-CV methods, the errors’ spatial structure was completely absorbed by almost all the models in all the glacier regions (statistically non-significant p-values) regardless the CV approach.

These results suggest that our expectation that the inclusion of spatial information in SP-CV method is supposed to capture the spatial autocorrelation of errors better than the NSP-CV doesn’t hold for glacier classification errors in the Tropical Andes using machine learning algorithms.

3.6. Spatial Predictions

Finally, we generated glacier outline maps over all studied regions using trained models to visually inspect and compare the predictions. For predict glacier class we took the approach of [108] described in [45]. We used all available pixels of each area to fit the final prediction models. As the error estimates from such a model are invalid, error estimates derived from the blocked cross-validation methods should still be used. This approach favors final prediction quality over perfect accuracy of error estimates. It has the advantage of using all the data and thus likely being the best predictor, particularly for smaller datasets. It has the disadvantage that the error estimates from the cross-validation no longer apply perfectly to the predictions, as they were made with slightly different models.

To analyze the predictions in greater detail we show a set of predictions on four small areas within Cordillera Blanca (Figure 8). Overall, the spatial distribution of the model predictions is quite similar. Here, KNN is not clearly the best predictor all around, as all the outlined maps exhibit a strong agreement with the NGI outline. There are, however, certain areas showing overestimation and underestimation of glacier surfaces. Notably, overestimation is noticeable, particularly near the glacier periphery of NGI in inset 3. In inept 2 we observe some areas of underestimation, where glaciers areas were misclassified as non-glacier, especially in the case of KNN and RF algorithms.

These classification errors could be attributed, to some extent, to the compositing and cloud filtering of Landsat 8 images in GEE. This filtering process might lead to the improper detection of transient snow areas, which are then mistaken for glaciers. Importantly, this phenomenon is not observed in other prediction areas. This highlights the significance of the cloud filtering algorithm, as snow pixels can seriously undermine the performance of machine learning models of glacier classification. Additionally, the spatial resolution of Landsat imagery used in this study can introduce a substantial number of mixed pixels—pixels that encompass both the glacier and the surrounding terrain. To address the former issue, a series of recommended algorithms has been previously suggested [57].

The debris-free outline map generated using the NDSI approach exhibits a high degree of similarity with the NGI outlines, although not perfect, probably due to NDSI being the current operational method for glacier mapping in Perú [11,52]. Although high similar, differences exist because our approach uses Landsat 8 imagery, whereas the NGI uses Sentinel 2 and Landsat 8 as well.

Prior studies have already concluded that debris-free ice can be accurately mapped using simple methods, such as the band ratio of the Red/NIR bands of Landsat or S2 data [26,42,53,64,109]. While our primary objective was to compare the predictive classification errors of machine learning algorithms using spatially distributed glacier class data, it's important to note that NDSI and band ratio approaches remain robust and widely used methods in remote sensing-based glacier monitoring. As this research demonstrates, these methods should not be dismissed. Nevertheless, evaluating the classification errors of glacier mapping using these approaches requires a proper assessment, typically involving a comparison with high-resolution satellite images [30,109], a procedure which is highly manual and inefficient. Therefore, we hope that a statistically sound approach could enhance current national glacier mapping efforts.

3.7. Methodological Limitations of the Study

Although machine learning (ML) has been effectively used to map glaciers in recent years [110], for example, [64] used approaches such as K-nearest neighbors (KNN), while [111,112] employed Random Forest (RF). The use of deep learning segmentation algorithms has been recently suggested for glacier mapping [42,110]. However, the computational cost can be high, limiting nationwide mapping efforts at the moment. Additionally, deep learning approaches have disadvantages compared to conventional machine learning algorithms, as spatial validation methods cannot be readily applied to them.

Some limitations detected in this work concern primarily the training/testing pixel sampling schema. When using multi-spectral optical satellite data for glacier mapping, an independent validation sample is usually not used, and maps are validated either manually [5,53,57] or through a comparison with high-resolution satellite images [30,109]. On the other hand, some studies used an independent sample [9] as a means of validation. We would argue that the size of the training and validation samples would inevitably affect the classification and validation results.

We acknowledge that hyperparameter tuning would likely yield better results or slightly alter the best model's outcomes [107]. A few previous studies have tested the effects of hyperparameter tuning and spatial cross-validation. For example, [43] implemented hyperparameter tuning in both spatial and non-spatial settings with mixed results. Therefore, here we fixed the hyperparameters for ease of comparison between the proposed validation methods, a similar approach that has been used by [36].

We also hypothesize that the size of the validation blocks needs calibration, as previous results suggested [36]. The size of the blocks should be substantially larger than the range of spatial autocorrelation in the model residuals to provide a reliable error estimate.

In general, we suggest that more work needs to be done to (1) incorporate different pixel sampling schemes and sample sizes into consideration, (2) develop experimental benchmarks that incorporate hyperparameter tuning, and (3) understand the effect of the spatial dimensions of validation blocks on the possible reduction of error estimates.

4. Conclusions

Here we have experimentally explored the effects of machine learning algorithm and cross validation in terms of glacier classification performance. If spatial dependence in the validation approach is not taking into account, our study provides that with the exception of Cordillera Blanca, Huayhuash and Vilcanota where Random Forest, GBM and logistic regression respectively outperforms in terms of Matthew correlation coefficient, KNN outperformed other models. In the SP-CV settings, the results change sightly, K-nearest neighbors (KNN) remained as the overall best model in almost all regions, followed by gradient boosting machines (GBM) and logistic regression (LR) in Cordillera Blanca and Central respectively. More importantly, we found that non-spatial cross validation leads to overoptimistic performance results (up to 18%), especially in the presence of spatial autocorrelation. Although these differences are not always statistically significant, we would recommend to use spatial CV instead of non-spatial CV for estimating the prediction performance of machine models when using spatial data, as only this ensures the assessment of bias-reduced predictive performance results, and as the error estimates are always conservative, this is especially important when the corresponding results form the basis of conservation and policy making.