*3.3. Land Cover Classification*

The supervised land cover classification has been carried out in two time steps, 1992 and 2019, while using Landsat 5 TM for the earlier date and for the latter Landsat 8 OLI. To account for annual variation in the snow and glacier coverage, data from August and September, when snow and glacier coverage have their annual minimum, have been used. A cloud mask was applied to remove cloud contamination.

The Random Forest Classifier (RFC) method [39] was applied using Google Earth Engine (GEE) for the classification [40]. For constructing the study wide cloud free mosaic, the median function of GEE has been used, which takes the median value of each pixel in available image temporal stack. In order to achieve higher classification accuracy, we followed the method of [41]: Gray-level co-occurrence matrix (GLCM) texture features [42,43] and spectral indices were produced in order to serve as collective variable predictors for the classification algorithm. The texture characterizes the variance of the pixel DN value over space, so it needs to be measured in a multiple pixel neighborhood. Within this neighborhood of pixels can be found the following three elements: tonal (DN) difference between pixels, the distance over which this difference is measured, and directionality [42]. These neighboring pixels are considered a window or kernel, and usually have a square area with an odd number of pixels for practical reasons. It is important to define the window size, because the larger window size can include edges or patches with different textures; this is particularly applicable for larger window sizes. For the first time in 1973, Haralick et al. [42] proposed GLCM textures, which are co-occurring or second-order texture measures. The calculation is based on tonal (DN) differences in a spatially defined relationship between pairs of pixels, taking into account all pixel pairs within the neighborhood. Hall-Beyer explains that second-order measurements can distinguish two pixels wide vertical stripes from one pixel wide stripes, given uniform DN values in each stripe; first-order texture measurements are not able to perform this [42]. The GLCM can account for all three elements of texture and that is one of its advantages. GLCM can be calculated while using single input layer and defined window size (i.e., 7 × 7), selected by the user, and can deliver to one or more output layers based on the selected measurements (i.e., variance, homogeneity, entropy, etc.). Based on the empirical result, a window size of 7 × 7 yielded a better result for generation of GLCM textures and the following textures features were generated: Variance, Inverse Difference Moment, which measures the homogeneity, Contrast, Dissimilarity, Entropy, Correlation, and Angular Second Moment. which

measures the number of repeated pairs [42]. The GLCM of band 3 and Band 4 of Landsat 8 were used for 2019 land cover classification, and GLCM band 4 of Landsat 5 TM was used for 1992 to generate the texture features. The spectral indices used include: Modified Normalized Difference Water Index (MNDWI) [44], Enhanced Vegetation Index (EVI) [45], Normalized Different Moisture Index (NDMI) [46], Green Optimized Soil Adjusted Vegetation Index (GOSAVI) [47], Built-up Area Extraction Index (BAEI) [48], and the Normalized Difference Bareness Index (NDBai) [49].

The Smile RFC method was applied using 200 decision trees and eight variables per split, which accounts for two-third of all variables. The number of input variables for both years have been filtered according to the variable importance function of the RFC and only the variables that contributed most have been selected in order to produce the final study area land cover for the list of input variables.

Training data were collected from annual land cover data by ESA [50]. Stratified random sampling techniques were used to collect 500 points per class with a total 10,000 points. The overall classification accuracy reached over 80% for all time steps.
