Large-Scale Debris Cover Glacier Mapping Using Multisource Object-Based Image Analysis Approach

: Large-scale debris cover glacier mapping can be efﬁciently conducted from high spatial resolution (HSR) remote sensing imagery using object-based image analysis (OBIA), which works on a group of pixels. This paper presents the spectral and spatial capabilities of OBIA to classify multiple glacier cover classes using a multisource approach by integrating multispectral, thermal, and slope information into one workﬂow. The novel contributions of this study are effective mapping of small yet important geomorphological features, classiﬁcation of shadow regions without manual corrections, discrimination of snow/ice, ice-mixed debris, and supraglacial debris without using shortwave infrared bands, and an adaptation of an area-weighted error matrix speciﬁcally built for assessing OBIA’s accuracy. The large-scale glacier cover map is produced with a high overall accuracy of ≈ 94% (area-weighted error matrix). The proposed OBIA approach also proved to be effective in mapping minor geomorphological features such as small glacial lakes, exposed ice faces, debris cones, rills, and crevasses with individual class accuracies in the range of 96.9–100%. We conﬁrm the portability of our proposed approach by comparing the results with reference glacier inventories and applying it to different sensor data and study areas.


Introduction
Detailed and accurate information about the size and spatial extent of debris-covered glaciers is vital for climate change and other glaciological applications such as observing glacier changes, conducting mass balance studies, predicting glacial lake outburst floods, and many more [1][2][3]. However, there appears to be a lack of large-scale (1:1000-1:5000) maps of debris-covered glaciers [4,5]. Recently, several remote sensing-based studies were carried out for glacier cover mapping using various pixel-based classification approaches [6][7][8]. However, object-based image analysis (OBIA), which works on the groups of homogeneous pixels called image objects, emerged as the alternative for the accurate extraction of fine details from high spatial resolution (HSR) images. It has significant advantages over traditional pixel-based classification approaches in utilising spatial information such as image texture, contextual information, pixel proximity, and geometric characteristics from HSR imagery [9].
Applications of OBIA in glaciology using remote sensing data received much attention in the past decade (e.g., [10][11][12][13][14][15][16][17][18]). Still, the literature is silent on mapping rugged Himalayan  The Bara Shigri Glacier, located at 32 • 10 12"N and 77 • 40 48"E, is also debris-covered in the Lahaul-Spiti district of Himachal Pradesh, northwestern Himalaya [28]. Figure 2 shows the lower ablation zone of the glacier, with the zoomed-in views showing the various glacier cover classes. It is the largest glacier in the Indian state of Himachal Pradesh, possessing 28 km in length. This glacier experiences hindrance by the Pir Panjal ranges from the south for the Indian summer monsoon and by the Great Himalayan ranges from the north for mid-latitude westerlies [28].

Datasets
The datasets used in this study are listed in Table 1. The datasets utilised for both glaciers include WorldView-2 imagery, Landsat TM thermal infrared (band 6) imagery, and ASTER Global DEM v2. The WorldView-2 multispectral imagery (2 m spatial resolution) covers different glacier cover classes. The expensive WorldView-2 data compelled us to procure a scene constituting only the main glacier trunk of Gangotri Glacier for experi-Remote Sens. 2022, 14, 3202 4 of 28 mental purposes. Additionally, we procured Resourcesat-2 Linear Imaging Self-Scanning System (LISS)-4 imagery of Gangotri Glacier to validate the robustness of the proposed multisource OBIA approach.
Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 31 The Bara Shigri Glacier, located at 32°10′12″N and 77°40′48″E, is also debris-covered in the Lahaul-Spiti district of Himachal Pradesh, northwestern Himalaya [28]. Figure 2 shows the lower ablation zone of the glacier, with the zoomed-in views showing the various glacier cover classes. It is the largest glacier in the Indian state of Himachal Pradesh, possessing 28 km in length. This glacier experiences hindrance by the Pir Panjal ranges from the south for the Indian summer monsoon and by the Great Himalayan ranges from the north for mid-latitude westerlies [28].

Datasets
The datasets used in this study are listed in Table 1. The datasets utilised for both glaciers include WorldView-2 imagery, Landsat TM thermal infrared (band 6) imagery, and ASTER Global DEM v2. The WorldView-2 multispectral imagery (2 m spatial resolution) covers different glacier cover classes. The expensive WorldView-2 data compelled us to procure a scene constituting only the main glacier trunk of Gangotri Glacier for experimental purposes. Additionally, we procured Resourcesat-2 Linear Imaging Self-Scanning System (LISS)-4 imagery of Gangotri Glacier to validate the robustness of the proposed multisource OBIA approach.
The images of the post-monsoon month were selected to ensure minimum cloud cover and snowfall over the study area. The ancillary information such as brightness, temperature, and slope was fulfilled in this study through the Landsat TM band 6 imagery (to integrate with WorldView-2 imagery for both Gangotri and Bara Shigri glaciers), Landsat 8 TIRS band 10 imagery (to integrate with LISS-4 imagery of Gangotri Glacier), and the ASTER GDEM v2 (30 m), respectively. These data formed the basis for multisource classification for mapping the glacier cover classes and geomorphological features using OBIA. The images of the post-monsoon month were selected to ensure minimum cloud cover and snowfall over the study area. The ancillary information such as brightness, temperature, and slope was fulfilled in this study through the Landsat TM band 6 imagery (to integrate with WorldView-2 imagery for both Gangotri and Bara Shigri glaciers), Landsat 8 TIRS band 10 imagery (to integrate with LISS-4 imagery of Gangotri Glacier), and the ASTER GDEM v2 (30 m), respectively. These data formed the basis for multisource classification for mapping the glacier cover classes and geomorphological features using OBIA.
The reference glacier outlines for both glaciers were generated by performing on-screen digitisation in the ArcGIS environment utilising HSR multispectral imagery (0.5 m) derived by pan-sharpening WorldView-2 multispectral (2 m) and panchromatic (0.5 m) images using the widely accepted Gram Schmidt method [29]. While performing digitisation, utmost care was taken to consider visual interpretation elements, different combinations of false-colour composites, and image spectra to digitise the glacier boundary with efficacy. We also selected the Randolph Glacier Inventory (RGI) 6.0 global inventory published on 20 July 2015 [30] to validate the glacier boundaries derived using the proposed multisource OBIA approach.
All the datasets were subjected to radiometric and atmospheric corrections as described in Figure 3. The brightness temperature and slope images were resampled to 2 m and 5 m and cropped to complement the spatial resolution and area of the WorldView-2 and LISS-4 images, respectively. The ASTER Global DEM v2 supported the 3D visualisation of the glacier, which assisted in the manual demarcation of the glacier boundaries in the regions covered with shadows, debris, and ice divides. We also corrected the RGI 6.0 global inventory for seasonal snow cover and ice cliffs for our study areas.  . OBIA conceptual workflow adopted for extracting glacier cover classes (BT: brightness temperature; GLCM: grey level co-occurrence matrix; NIR: near-infrared).  Figure 3. OBIA conceptual workflow adopted for extracting glacier cover classes (BT: brightness temperature; GLCM: grey level co-occurrence matrix; NIR: near-infrared).

Methodology
This section explains the methodology workflow developed ( Figure 3) to achieve the research objectives of this study. The segmentation and classification were performed using multiresolution segmentation and ruleset-based classification, all incorporated in the most popular software for OBIA, i.e., eCognition. The advantages of ruleset-based OBIA classification are the incorporation of expert knowledge in the classification, formulation of complex class descriptions, and transferability. The sequential steps involved in the OBIA methodology adopted here for executing the multisource classification of glacier cover mapping include (1) identification of glacier cover classes using visual and spectral analysis; (2) image segmentation and attributes selection; (3) multisource OBIA classification with refinement, and (4) accuracy assessment.
We first experimented with the WorldView-2 + ancillary dataset covering the Gangotri Glacier. To determine the transferability of our approach, we then applied it to the LISS-4 + ancillary dataset of Gangotri Glacier and another WorldView-2 + ancillary dataset covering the Bara Shigri Glacier with few modifications in the rulesets.

Descriptive Analysis of Glacier Cover Classes
Visual and spectral identification of glacier cover classes over the study area is required to define their characteristics precisely to formulate rulesets in OBIA for each class.

Visual Interpretation
The classes that could be mapped in this study were identified on the false-colour composite (R: Near-Infrared (NIR-1), G: Red, B: Green) of WorldView-2 imagery, including snow/ice, IMD, SGD, PGD, valley rock, glacial lakes, exposed ice-faces (EIFs), crevasses, rills, debris cones, and meltwater streams (refer to Figure 1). Here, the snow/ice class refers to clean-ice. The mountain shadows, too, were considered a class. We could also identify outwash deposits as marked in Figure 1; however, we considered it in the PGD class.
The SGD is the debris on the glacier's surface [23], and the PGD is the debris found along the glacier's margins [23]. Valley rock, which appears red, is a source for both the SGD and PGD; therefore, the SGD and PGD are spectrally indistinguishable. Snow/ice forms the clean-ice part of the glacier and appears brightest. IMD, which appears in cyan, is a mixture of the SGD and ice and forms when the ice starts giving way to the SGD in the ablation zone [31]. Glacial lakes appear in blue, and mapping of glacial lakes formed at the glacier's terminus is essential because they may breach in a glacial lake outburst flood and can cause immense harm to the life and property downstream [32]. EIFs are the steep ice slopes formed by slumping or undercutting by surface water and are crucial in studying the mass balance of debris-covered glaciers [33]. They appear in blue but with some debris mixed with it and also little more protruded than glacial lakes. A crevasse is a deep crack or fracture in the ice, which may either be open or hidden by snow. The crevasses, which appear as black elongated features in the imagery, pose difficulties and dangers for mountaineers and glacier travellers. A debris cone is conical-shaped sedimentation formed at the base of valley-side slopes [34] and plays a crucial role in the development and evolution of debris-covered glaciers [35]. Debris cones appear brighter than the PGD. Rills that act as actual glacier drainage systems are shallow channels formed in debris cover formed by flowing water [36]. They are found in the PGD and appear as black and elongated like crevasses. The meltwater stream is a source of water emerging out of the snout and can help track a glacier snout's position [37]. It appears as a dark blue elongated feature. Mountain shadows that appear black pose a challenge in debris-covered glacier mapping because snow/ice, IMD, SGD, PGD, and valley rock information seem to be hidden under them.

Spectral Analysis
An in-depth analysis of the spectral properties of the glacier cover classes to be mapped was performed by selecting several representative samples for each class. The average spectral response plots of all the glacier cover classes in the optical and thermal infrared bands are shown in Figure 4. The reflectance and brightness temperature values were normalised to percentages for straightforward representation. gated feature. Mountain shadows that appear black pose a challenge in debris-covered glacier mapping because snow/ice, IMD, SGD, PGD, and valley rock information seem to be hidden under them.

Spectral Analysis
An in-depth analysis of the spectral properties of the glacier cover classes to be mapped was performed by selecting several representative samples for each class. The average spectral response plots of all the glacier cover classes in the optical and thermal infrared bands are shown in Figure 4. The reflectance and brightness temperature values were normalised to percentages for straightforward representation. It can be seen in Figure 4 that the main confusion arises in the spectral separation of the SGD, PGD, and valley rock from each other. The SGD shows a higher spectral response in the yellow and red-edge wavelength regions (WorldView-2 bands 4 and 6, respectively), which is similar to the PGD and valley rock, as these are compositionally the same. Therefore, the brightness temperature and slope profiles ( Figure 5) of these spectrally similar classes help segregate them.
It can be seen from Figure 5A that the spectral response of the SGD deviates remarkably from the PGD and valley rock in the thermal infrared region (Landsat TM band 6). This is because the SGD has a lower temperature than the PGD and valley rock due to the glacier ice beneath it [23,38]. Similarly, it can be observed from Figure 5B that the slope information also aids in separating the glacier surface (SGD) from the non-glacier surface (PGD and valley rock) as the former can be found at lower slopes than the latter [2,39]. It can be seen in Figure 4 that the main confusion arises in the spectral separation of the SGD, PGD, and valley rock from each other. The SGD shows a higher spectral response in the yellow and red-edge wavelength regions (WorldView-2 bands 4 and 6, respectively), which is similar to the PGD and valley rock, as these are compositionally the same. Therefore, the brightness temperature and slope profiles ( Figure 5) of these spectrally similar classes help segregate them.
It can be seen from Figure 5A that the spectral response of the SGD deviates remarkably from the PGD and valley rock in the thermal infrared region (Landsat TM band 6). This is because the SGD has a lower temperature than the PGD and valley rock due to the glacier ice beneath it [23,38]. Similarly, it can be observed from Figure 5B that the slope information also aids in separating the glacier surface (SGD) from the non-glacier surface (PGD and valley rock) as the former can be found at lower slopes than the latter [2,39].
Thus, analysis of brightness temperature and slope profiles shown in Figure 5A,B, respectively, clearly depicts a considerable difference among the SGD, PGD, and valley rock. Thus, to compensate for spectral similarities, brightness temperature and slope data were included as ancillary layers with the spectral data to improve the delineation between the SGD, PGD, and valley rock.

Image Segmentation and Attributes Selection
The fundamental step in OBIA is determining the optimal segmentation parameters for the multiresolution segmentation process [40], as optimal segmentation reduces the under-and over-segmentation errors. The inputs to the segmentation process, that is, reflectance, brightness temperature, and slope, were weighted as 1, 2, and 2, respectively. High weights to brightness temperature and slope were assigned to discriminate the SGD, PGD, and valley rock accurately. The eight bands of WorldView-2 multispectral imagery were weighted equally to consider all the spectral information made available by each band.
Thus, analysis of brightness temperature and slope profiles shown in Figure 5A,B, respectively, clearly depicts a considerable difference among the SGD, PGD, and valley rock. Thus, to compensate for spectral similarities, brightness temperature and slope data were included as ancillary layers with the spectral data to improve the delineation between the SGD, PGD, and valley rock.

Image Segmentation and Attributes Selection
The fundamental step in OBIA is determining the optimal segmentation parameters for the multiresolution segmentation process [40], as optimal segmentation reduces the under-and over-segmentation errors. The inputs to the segmentation process, that is, reflectance, brightness temperature, and slope, were weighted as 1, 2, and 2, respectively. High weights to brightness temperature and slope were assigned to discriminate the SGD, PGD, and valley rock accurately. The eight bands of WorldView-2 multispectral imagery were weighted equally to consider all the spectral information made available by each band.

Selecting Optimal Segmentation Parameters
The multiresolution segmentation process based on a mutual-best-fitting approach begins with one pixel and iteratively merges them in pairs to larger sizes until the upper threshold of the homogeneity criterion does not exceed locally. The homogeneity criterion is a combination of spectral and shape homogeneity and is defined by the three segmentation parameters: scale, shape, and colour [41]. This study followed the iterative evaluation procedure [42] to derive the optimal segmentation parameters.
Further, one major characteristic of multiresolution segmentation is its ability to segment an image into segments of similar scales at the same level and objects of different scales across levels. This concept is referred to as multiscale segmentation [43]. Multiscale segmentation is based on the inheritance concept, i.e., the attributes of a parent (or super) segment are inherited by a child (or sub) segment. The technique offers the possibility of creating hierarchy levels that facilitate accurate extraction by allowing topology relationships to be built between segments produced at finer scales and those generated at coarser scales. Therefore, we extracted objects using multiscale segmentation, and the scale parameter influences the maximum and minimum size of individual objects. A larger value of the scale parameter results in larger image segments and vice-versa [41,44]. At a low scale value, image segments contain a few pixels and thus are too small to represent any class of interest in the imagery effectively. At a high scale value, image segments become so large that they may represent more than one glacier class. After moving over the complete range of scale values, a scale parameter of 200 is found to be appropriate for extracting glacier cover classes of larger size, such as snow/ice, IMD, SGD, PGD, and valley rock, whereas a scale parameter of 60 was determined to be suitable for extracting glacier cover classes of smaller size such as glacial lakes, EIFs, crevasses, rills, and meltwater streams for this study. Figure 6 explains the segmentation results. Eventually, values ranging from compactness = 0.1/smoothness 0.9 to compactness = 0.9/smoothness 0.1, were tested iteratively in steps of 0.1. Lower compactness (higher smoothness) values allow objects to increase in length and become narrower. On the contrary, raising the value of the compactness parameter minimises the perimeter/area ratio of image objects and results in compact and regular objects. A compactness value of 0.4 (and smoothness of 0.6) is optimal here. Iterative evaluation of the segmentation parameters thus resulted in the most optimum set of parameters such as scale = 200 (for large size glacier cover classes)/60 (small size glacier cover classes), shape = 0.2, and compactness =  and EIFs) at a scale parameter of 60, respectively. On the other hand, Figure 6B shows the under-segmentation of the large-size glacier cover classes (SGD, PGD, and valley rock) at a smaller scale than 200. Figure 6C shows the over-segmentation of the small-size glacier cover classes (glacial lakes and EIFs) at a scale larger than 60. As seen in Figure 6D, the segments fit well with the glacier cover classes to be extracted, and the overlapping with the surrounding glacier cover classes was reduced.
The shape parameter defines the textural homogeneity of segments and is a composite of smoothness and compactness. Colour (spectral value) is the segment's pixel value or reflectance. The shape and colour parameters are related in a way that their total value is equal to one. Similarly, smoothness and compactness parameters are also related so that their total value equals one [41,44]. Values for the shape/color parameter ranging from shape = 0.1/color = 0.9 to shape = 0.9/color = 0.1 in increments of 0.1 were tested for the study area. At the shape value of 0.1, the colour parameter having a 0.9 value allows the resulting image segments to grow so that it maximises the homogeneity of image spectral values. The value of a colour parameter is decreased and increased, maintaining the uniformity of the shape parameter, which in turn is regulated by the compactness/smoothness parameter, as discussed below. A relatively low value of 0.2 for the shape parameter (colour 0.8) was selected, as glacier cover classes usually exhibit irregular shapes due to the complex landforms.
Eventually, values ranging from compactness = 0.1/smoothness 0.9 to compactness = 0.9/smoothness 0.1, were tested iteratively in steps of 0.1. Lower compactness (higher smoothness) values allow objects to increase in length and become narrower. On the contrary, raising the value of the compactness parameter minimises the perimeter/area ratio of image objects and results in compact and regular objects. A compactness value of 0.4 (and smoothness of 0.6) is optimal here. Iterative evaluation of the segmentation parameters thus resulted in the most optimum set of parameters such as scale = 200 (for large size glacier cover classes)/60 (small size glacier cover classes), shape = 0.2, and compactness = 0.4. The segmentation levels based on the optimum scale parameters were named L200 and L60 for subsequent processing (Figure 3).

Selection of Attributes
Once an optimum segmentation is achieved, a class hierarchy is prepared ( Figure 3) to establish hierarchical and contextual relationships among the objects of interest. The class SGD was divided into subclasses (child classes) containing more detailed glacier cover classes such as glacial lakes, EIFs, and crevasses. Similarly, the PGD class was broken into debris cones, rills, and meltwater streams. Thus, the classification was performed in two levels of segmentation as described in the methodology (refer to Figure 3): (1) L200 level, to classify large-size segments (i.e., snow/ice, IMD, SGD, PGD, valley rock, and shadows); (2) L60 level, to classify small size segments (i.e., glacial lakes, EIFs, crevasses, debris cones, rills, and meltwater streams).
The beauty of OBIA is that one can compute object attributes per layer such as the mean value, spectral ratio, standard deviation, and so on; besides, there are geometry attributes such as shape, area, and the like and texture attributes such as the standard deviation and grey level co-occurrence matrix (GLCM) statistics. As described in Section 3.1, "identification of glacier cover classes", the brightness temperature, and slope data are stacked with reflectance data to deal with the spectral similarity of the SGD, PGD, and valley rock. Therefore, attributes such as the mean reflectance value, standard deviation, brightness temperature, and slope form the inputs to OBIA to produce the glacier cover map. Other derived attributes include the NIR-2/Yellow band ratio, shadow detection index, normalised difference water index (NDWI), GLCM mean, Max. Diff., length/area, and their combinations that are found to be suitable in this study to map other glacier cover classes depending on their characteristics. All the attributes (with their descriptions) defined in this study to map the glacier cover classes are listed in Table 2. Table 2. Spectral and spatial attributes selected and their usage.

Attributes
Usage Description

Mean and standard deviation
Preliminary classification of all glacier cover classes The mean value and standard deviation attributes describe the spectral properties of image objects [45].

Mean Brightness Temperature (K)
To differentiate between SGD, PGD, and valley rock SGD has a lower temperature as compared to PGD and valley rock due to the underlying glacier ice ( Figure 5A).
To map snow/ice and IMD In the thermal infrared region, IMD shows a higher spectral response than snow/ice because IMD contains debris that raises its temperature ( Figure 4).
To map shadows The temperature of shadows is remarkably lower than all the other classes ( Figure 4).

Mean Slope (Degrees)
To differentiate SGD from PGD and valley rock The glacier surface (SGD) is formed at lower slopes than the non-glacier surface (PGD and valley rock).

NIR-2/Yellow
To differentiate IMD from snow/ice Snow/ice and IMD have a high spectral response in the yellow wavelength region, whereas it is low in the NIR-2 region. Further, the spectral response of snow/ice in both these regions is relatively higher than that of IMD ( Figure 4). Therefore, band ratio NIR-2/Yellow was developed to distinguish IMD from snow/ice.

Shadow Detection Index To map shadows and EIFs
Shadow detection index developed by Shahi et al. [46] using WorldView-2 multispectral data could effectively map shadows. The formula for the Shadow detection index is defined in Equation (1).

Max. Diff. To map glacial lakes
Max. Diff. attribute [47] is the result of the difference between the maximum value of an object and its minimum mean value. The mean values of all layers (WorldView-2 reflectance imagery, brightness temperature, and slope) belonging to an object are compared to get the maximum and minimum values. Subsequently, the result is divided by the brightness. Glacial lakes show the highest values of Max. Diff. attribute.

NDWI
To map EIFs NDWI used in this study is a spectral index developed by Wolf [48] using WorldView-2 multispectral data. The formulation of NDWI is as given in Equation (2). EIFs were successfully mapped using NDWI.

GLCM mean
To map crevasses GLCM mean is the average expressed in terms of the GLCM. The pixel data are weighted by the frequency of its occurrence in combination with a certain neighbor pixel data [49]. Crevasses were mapped using GLCM mean.

NIR-1 To map debris cones
The spectral response of debris cones is higher than that of PGD in the NIR-1 region ( Figure 4). Therefore, debris cones are classified using mean values of NIR-1.

Length/Area
To map rills and meltwater streams The length/area is a geometric attribute, which defines elongation [47]. Since rills and meltwater streams (that originate from the snout) are elongated features, length/area attribute was used to classify these classes.

Multisource OBIA Classification
The following two-step classification is performed to obtain the large-scale glacier cover map. The knowledge base on the attributes as defined in Table 2 was utilised in defining the rulesets. The thresholds of the attributes were determined through trial and error (for threshold values, refer to Figure 3) using the feature selection tool of eCognition.

Mapping of Objects at L200 Level
Objects of interest such as the SGD, PGD, valley rock, snow/ice, IMD, and shadows were delineated at the L200 segmentation level since they constitute the basis for extracting other classes. At first, the SGD and PGD were delineated using the brightness temperature and slope information. Valley rock found at a higher slope than the SGD and PGD was classified using the slope information alone. To classify snow/ice, the brightness temperature was found to be valid, whereas to classify IMD, which is a mixture of snow/ice and debris, the NIR-2/Yellow band ratio was found to be applicable. Mountain shadows that exhibited lower temperatures than other classes were classified successfully using a combination of the shadow detection index (defined in Equation (1)) and brightness temperature.

Mapping of Objects at L60 Level
At the L60 segmentation level, objects of interest such as glacial lakes, EIFs, crevasses, debris cones, rills, and meltwater streams were classified. Glacial lakes were mapped using the Max. Diff. attribute. Since EIFs exhibit similar properties to glacial lakes, a combination of NDWI (defined in Equation (2)) and the shadow detection index was used to delineate EIFs from glacial lakes. To map crevasses, the GLCM mean was used. Debris cones that show spectral similarity with the PGD were differentiated from it using the mean values of the NIR-1 band. The length/area attribute was used to extract rills and meltwater streams.

Refinement
The results obtained using the multisource OBIA were refined using the merge region and region growing algorithms of eCognition software. Some objects on the SGD were misclassified as the PGD, and a few PGD objects were mislabeled as the SGD. Both misclassifications may occur owing to the slight overlap of brightness temperature values in the SGD and PGD. To correct the SGD misclassified objects, we used a ruleset as "relative area of SGD ≥ 0.75" in the region growing algorithm. It says that all the objects classified as the PGD and having a "relative area of SGD ≥ 0.75" be grown into the SGD. Likewise, to rectify the PGD misclassified objects, we applied the ruleset as "relative border to PGD > 0.5" in the region growing algorithm. It says that all the objects classified as the SGD and having a "relative border to PGD > 0.5" be grown into the PGD.
Some finer meltwater stream objects were misclassified as the PGD, specifically, outwash deposits (as described in Section 3.1.1 "visual interpretation"). Meltwater streams at the snout are usually loaded with debris and possess high velocities. However, outside the glacier boundary, the water loses its velocity, causing some of the debris to settle down, allowing the water stream to split into many channels enclosed by sediments. This process forms braided stream networks, commonly called outwash (here, the PGD) [37]. Therefore, the thinner meltwater stream objects are highly likely to be wrongly classified as outwash. We developed a ruleset "0.1 < NDWI < 0.2" in the region growing algorithm to correct it. It means that all the meltwater stream objects with "0.1 < NDWI < 0.2" be produced as the PGD.
Eventually, the objects of each glacier cover class are merged into a single region using the merge region algorithm. Thus, the final glacier cover map containing a total of twelve glacier cover classes is shown in Figure 7.

Accuracy Assessment
The essence of ground truth data for validating classification results was highlighted by Foody [50]. However, due to the high cost of field expenditure and inaccessibility of the Gangotri Glacier, in-situ data could not be collected for validation in this study. To deal with this, a simple stratified random sampling scheme [51] was followed to collect the testing samples on the same multisource dataset used for classification. The collection of testing samples was aided by visual interpretation, spectral signature plots, and the knowledge of experts.
The critical difference in evaluating pixel-based and object-based classification is the reference data type which can be either a point (pixel) or an object. The use of point data tends to underestimate the accuracy of the classification result obtained from OBIA, and there is uncertainty about whether that class is consistent across the entire object [52]. Therefore, we used objects as references for assessing the efficacy of the proposed OBIA approach.
Further, we adapted the area-weighted error matrix proposed explicitly for the accuracy assessment of OBIA by Radoux et al. [24] for our case study and compared it with the traditional pixel-based error matrix. The area-weighted error matrix accounts for the variable object sizes, unlike the conventional pixel-based error matrix, and its implementation is shown by MacLean and Congalton [53].
A total of 469 testing sample objects were collected on the WorldView-2 + ancillary dataset in the case of Gangotri Glacier. The error matrices for defining the accuracies are shown in Table 3.

Transferability of Proposed Multisource OBIA Approach
To ascertain the transferability of the proposed multisource OBIA approach, we applied the same rulesets but with slight alterations in the threshold values on (1) same study area but different sensor data and (2) same sensor data but different study area. Slight modifications in the thresholds and/or rulesets are required to be performed because local temperature regime and topography are bound to fluctuate with a change in the study area (glacier), as supported by Shukla and Ali [2]; Bolch et al. [22]; and Keshri et al. [31]. The rulesets are also constrained by the spatial and spectral characteristics of the satellite imagery used. For example, fine details can be extracted from high spatial resolution imagery (e.g., WorldView-2 having 2 m spatial resolution), which is not the case with comparatively lower spatial resolution imagery (e.g., LISS-4 having 5m spatial resolution).

Same Study Area, Different Sensor Data
In this scenario, we used LISS-4 imagery covering the Gangotri Glacier to implement the proposed rulesets. The ancillary information of brightness temperature and the slope was collected from Landsat 8 TIRS band 10 and ASTER GDEM v2, respectively. The spatial resolution of LISS-4 data limited the classification to snow/ice, IMD, SGD, PGD, valley rock, glacial lakes, debris cones, and mountain shadows as other small features such as crevasses, EIFs, and rills could not be identified on the LISS-4 imagery ( Figure 8A). Regarding image segmentation, it also means that the scale of analysis varies with the spatial resolution of the imagery because the scale parameter depends on the spatial resolution. Objects may form at different scales of analysis in many scenarios, even in the same image [54]. Hence, the LISS-4 + ancillary dataset analysis had to be modified accordingly.  Here, we used WorldView-2 imagery covering the Bara Shigri Glacier to implement

LISS-4 + Ancillary Dataset Analysis
In the segmentation process, iterative testing facilitated the eventual selection of scale = 190 (for large size glacier cover classes such as snow/ice, IMD, SGD, PGD and valley rock)/40 (small size glacier cover classes such as glacial lakes/ponds), shape = 0.2, and compactness = 0.4. Given that the LISS-4 imagery does not contain coastal blue and yellow spectral bands, the rulesets for the LISS-4 + ancillary dataset had to be revised. This also implied that the spectral index, such as the shadow detection index and NIR-2/Yellow band ratio, had to be replaced. The rest of the OBIA approach stays similar to the WorldView-2 + ancillary dataset. The rulesets altered to classify the LISS-4 + ancillary dataset are listed in Table 4. To refine the SGD and PGD classes using the region growing algorithm, the same rulesets were utilised as that for the WorldView-2 + ancillary dataset of Gangotri Glacier with certain adjustments in the threshold values as tabulated in Table 4. The objects of every class were merged into a single region using the merge region algorithm. The final large-scale thematic map of Gangotri Glacier using the LISS-4 + ancillary dataset is shown in Figure 8B.
To assess the accuracy of the classified image, a total of 309 testing objects were collected using stratified random sampling. The area-weighted error matrix shown in Table 5 demonstrates the individual accuracies of the glacier cover classes. Table 5. User's and producer's accuracies of individual glacier cover classes of Gangotri Glacier using LISS-4 + ancillary dataset with an overall accuracy of 89.9%.

Same Sensor Data, Different Study Area
Here, we used WorldView-2 imagery covering the Bara Shigri Glacier to implement the proposed rulesets. The ancillary information of brightness temperature and the slope was collected from Landsat TM band 6 and ASTER GDEM v2, respectively. We could map all the glacier cover classes in Bara Shigri Glacier as extracted in Gangotri Glacier, except crevasses as they are missing in the scene of Bara Shigri Glacier. Here also, we could observe outwash deposits (refer to Figure 2) which were included in the PGD class. Further, as can be observed from Figure 2, the snow/ice and IMD are found outside the glacier boundary in the scene. Therefore, the rulesets used to classify snow/ice and IMD were modified suitably, as described in Table 6.
An additional attribute of brightness with brightness temperature aided in classifying the snow/ice class due to its high brightness values. On the contrary, IMD exhibits less brightness and a lower brightness temperature than snow/ice. Another adjustment was made to the glacial lakes' ruleset. Unlike Gangotri Glacier, which is full of supraglacial lakes, we also observed few periglacial lakes in Bara Shigri Glacier other than the supraglacial lakes (refer to Figure 2). However, the ruleset defined for classifying supraglacial lakes in Gangotri Glacier alone served the purpose of classifying both the supraglacial and periglacial lakes of Bara Shigri Glacier. The threshold range of the Max. Diff attribute was revised consequently, as referred to in Table 6 below. Later on, we simply merged them to signify glacial lakes in general.
In the refinement process, the classification of the SGD, PGD, and meltwater streams was improved by implementing the same rulesets in the region growing algorithm as that of the WorldView-2 + ancillary dataset of Gangotri Glacier. Only specific modifications in the threshold values were performed, as presented in Table 6. The objects of every class were merged into a single region using the merge region algorithm. The final thematic map of Bara Shigri Glacier is presented in Figure 9.

Results
In this section, we start with the evaluation of the spectral patterns and spatial properties of the glacier cover classes. This is then accompanied by the qualitative and quantitative analysis of the glacier cover classes mapped using the proposed spectral and spatial   To perform accuracy assessment, the area-weighted error matrix was defined in Table 7 to evaluate the individual accuracies of the classes by collecting a total of 178 testing sample objects.

Results
In this section, we start with the evaluation of the spectral patterns and spatial properties of the glacier cover classes. This is then accompanied by the qualitative and quantitative analysis of the glacier cover classes mapped using the proposed spectral and spatial attributes developed in the OBIA environment.
The quantitative evaluation uses the error matrix-based statistical measures such as overall accuracy, user's accuracy, and producer's accuracy. Segmentation, the pivotal step in OBIA is also discussed in this section to assess the quality of the segmentation results and to elaborate on the essence of integrating ancillary layers such as temperature and slope. Figure 4 portrays the spectral response curves for each class in the WorldView-2 multispectral imagery and Landsat TM thermal infrared imagery. Snow/ice and IMD can be easily separated from other classes due to their relatively high spectral response in all wavelength regions except the thermal infrared region. Further, it was observed that the snow/ice and IMD classes have a high spectral response in the yellow region, whereas it is low in the NIR-2 region. Moreover, the spectral response of snow/ice in both these regions is relatively higher than that of IMD. Therefore, the band ratio NIR-2/Yellow was developed in this study to distinguish IMD from snow/ice. Similarly, the debris cones follow the spectral pattern of the PGD. However, in the NIR-1 region, the spectral response of debris cones is higher than that of the PGD. Therefore, the spectral values of NIR-1 are found to be adequate to classify debris cones. Further, glacial lakes and EIFs, both of which are formed on the SGD, depict spectral overlap in the NIR-1 region. However, lakes show a higher spectral response than EIFs in all the other wavelength regions. Therefore, we used NDWI defined in Equation (2) and developed by Wolf [48] to map them. However, applying a global threshold misclassified the glacial lakes as EIFs.

Analysis of Spectral and Spatial Properties
It is not surprising to reveal that crevasses, rills, and shadows resemble each other in their spectral responses. Still, shadows differ from the other two in the thermal infrared region, and therefore, brightness temperature is derived from the thermal infrared layer aided in discriminating the shadows. A temperature difference can also be seen between rills and crevasses; however, both are so small that the global brightness temperature threshold would not give accurate results. Therefore, mapping of glacier cover classes in such cases of ambiguities would depend on the spatial attributes, which allow a degree of separability among them. For example, glacial lakes show the highest values of the Max. Diff. attribute, and therefore they were distinguished successfully from EIFs using the Max. Diff. To deal with the misclassification between glacial lakes and shadows, shadows were refined using the shadow detection index developed by Shahi et al. [46] defined in Equation (1). Furthermore, since rills are elongated features, the length/area attribute that defines elongation [47] helped segregate them from crevasses, and crevasses were refined using the GLCM mean. Similarly, meltwater streams (that originate from the snout), which are also elongated features, were discriminated from glacial lakes using the length/area attribute.

Comparison between the Segmentation Results without and with Ancillary Data
The purpose of this section is to explain the impact of ancillary data used in this study. For this, the subsets of segmentation results using only multispectral (single-source) data and using an integrated dataset consisting of multispectral and ancillary (multisource) information of the brightness temperature and slope are shown in Figure 6E,F, respectively. It can be seen from Figure 6E that the objects of the SGD and PGD cannot be discriminated in some places, whereas in Figure 6F, after adding ancillary data, they can be demarcated. A similar observation was made in the case of the differentiation of the PGD and valley rock. The PGD and valley rock can be seen as a single object in Figure 6E, whereas these are shown as two different objects in Figure 6F after using the ancillary data. This demonstrates the use of brightness temperature and slope data for accurately discriminating the SGD, PGD, and valley rock from each other.
Further, the multiscale segmentation approach facilitated the delineation of individual classes having between-class spectral variation (such as the SGD, PGD, and valley rock) and within-class spectral variation (such as between glacial lakes and EIFs within the SGD parent class). The critical point is the temperature difference between the glacier and non-glacier regions and the slope difference between the flat glacier surface and the steep non-glacier surface. Their low spectral contrast is the main issue in differentiating the glacier surface (SGD) and non-glacier surface (PGD, valley rock). Assigning high weight to the brightness temperature and slope at the initial segmentation process assisted in differentiating glacier and non-glacier parts.

Assessment of Individual Class Accuracies
In glaciology, assessing classification quality using error matrix-based statistical accuracy measures is very common and widely used [7,55]. Foody [56] argues that comparing reference and classified labels formed from differently sized sampling units may lead to dissimilar accuracy estimates. In the case of OBIA, where the reference sample units are polygons, i.e., segments (and not pixels), the size of each reference sample unit varies [24]. Therefore, the traditional pixel-based error matrix that deals with the pixels (which are of the same size) cannot tackle the variable sizes of the polygons in the accuracy assessment. Therefore, Radoux et al. [24] exclusively designed the area-weighted error matrix for assessing the accuracy of OBIA in which each cell of the traditional pixel-based error matrix is weighted by the area of the reference units [53]. Thus, this study inclined towards OBIA discusses the results of classification based on the area-weighted error matrix. Furthermore, high individual accuracies of user's and producer's accuracy with a negligible or minimal difference between them suggest accurate delineation of the class. Therefore, in this study, we present the results by evaluating the individual accuracies of each class.

Gangotri Glacier
The glacier surface covers snow/ice, IMD, and the SGD; the non-glacier surface includes the PGD and valley rock. The features that define the geomorphology of debriscovered glaciers include glacial lakes, EIFs, debris cones, rills, meltwater streams, and shadows. The analysis of the Gangotri Glacier is discussed as follows using both the WorldView-2 and LISS-4 images.

•
WorldView-2 + Ancillary Dataset The developed multisource OBIA approach achieved high user's accuracy in the range 91.1-99.4% for mapping of snow/ice and IMD. Distinguishing the SGD is difficult because of its spectral overlap with the PGD and valley rock. However, the OBIA approach commendably mapped the SGD with a very high user's accuracy of 99.0% and producer's accuracy of 96.8%. The 2.2% difference between the user's accuracy and producer's accuracy may be linked to the miscategorisation of some pixels as the SGD owing to their lower temperatures and slope gradients. The fine supraglacial features such as glacial lakes, EIFs, and crevasses, which are part of glacial geomorphology, were also mapped with high accuracies ranging between 98.4% and 100%.
The non-glacier surface features also showed good performance, with the PGD exhibiting high values of both accuracies (user's accuracy = 94.5% and producer's accuracy = 85.4%). The approach also precisely mapped valley rock with individual accuracies of user's accuracy = 95.2%, producer's accuracy = 91.9%. The geomorphological features such as debris cones and rills were classified with accuracies varying between 88.1% and 100%. The meltwater stream that gushes out from the snout was mapped with the accuracy of the user's and producer's, respectively, of 92.1% and 100%. •

LISS-4 + Ancillary Dataset
Implementing the proposed multisource OBIA approach on LISS-4 imagery extracted the glacier surface features, that is, snow/ice, IMD, and SGD, with great user's and producer's accuracies extending between 85.7% and 98.8%. Glacial lakes were also mapped with respective user's and producer's accuracies of 100% and 92.9%.
However, discrimination and precise mapping of the PGD were again the key conclusions in the case of the non-glacier surface. The PGD was extracted with the lowest user's accuracy of 78.1% relative to other non-glacier surface features. Though this reduced accuracy does not negatively affect the glacier boundary, this low accuracy is apparently due to the lower spatial resolution of LISS-4 imagery than the WorldView-2 imagery. It can also be due to the constant threshold of brightness temperature and slope applied over the entire scene for its differentiation. Even so, valley rock and debris cones achieved accuracies between 80.4% and 95.7%.

Bara Shigri Glacier
Here, the glacier surface covers only the SGD in the scene because snow/ice and IMD are seen surrounding the glacier surface area in the scene (Figure 2). The proposed multisource OBIA approach classified the SGD with an excellent user's accuracy of 98.5% and producer's accuracy of 100.0%. Snow/Ice and IMD were distinguished with high individual accuracies varying between 89.1% and 100%. EIFs, which are a part of the SGD, were mapped with 100% user's accuracy but low producer's accuracy of 69.0%. This lowest producer's accuracy of EIFs is due to the misclassification of an EIF as the SGD because EIFs are a mixture of snow/ice, water, and debris [33].
The non-glacier surface that comprises the PGD, valley rock, meltwater streams, and shadows demonstrated effective results with individual accuracies extending between 89.5 and 100%. The lowest accuracy of 89.5% possessed by the PGD is owing to the erroneous classification of a few PGD objects as valley rock which again does not disrupt the accuracy of the glacier surface outline. Valley rock was extracted with high corresponding user's and producer's accuracies of 96.3% and 90.4%. The geomorphological features, debris cones and rills, were extracted successfully with high individual accuracies lying between 81.4% and 100%. Even the meltwater streams were classified precisely, giving user's and producer's accuracies of 96.3% and 94.3%, respectively. Still, we found that some rills marked in Figure 2 could not be mapped. These are the objects which are, in reality, rills, but maybe the scale value at which the rills were classified was so coarse that it led to overestimation error and merged them with the PGD objects surrounding them. A finer scale could allow the optimal segmentation of rills to extract them precisely.
Glacial lakes, which include both supraglacial and periglacial lakes, attained user's and producer's accuracies of 87.3% and 100.0%, respectively. The 12.7% variation in the user's and producer's accuracy is attributable to the incorrect classification of some reference meltwater objects as glacial lakes, specifically periglacial lakes.
In conclusion, while glacier mapping becomes a difficult task due to mountain shadows, the shadow detection index derived using WorldView-2 could map the shadows with high accuracies for both the glaciers and without manual corrections. We could also map mountain shadows precisely using the brightness attribute in the LISS-4 + ancillary dataset case. Overall, the proposed multisource OBIA approach accurately mapped all glacier cover classes, especially for spectrally similar ones such as snow/ice, IMD, SGD, PGD, valley rock, and debris cones. Moreover, the proposed multisource OBIA approach enabled the accurate mapping of geomorphological features of glacial lakes, EIFs, rills, and crevasses using spectral and spatial attributes.

Appreciation of Glacier Boundaries
The classification accuracy was also assessed by comparing the derived glacier boundary with the manually digitised glacier boundary and the global inventory of glacier outlines provided by RGI. The HSR remote sensing data showed promising results in creating glacier boundaries [1]. However, due to the similar spectral property of the glacier surface with the surrounding area (PGD, valley rock), such boundaries are generally created using visual interpretation, which is collectively laborious and subjective. Here, we present an automatic method based on a combination of HSR multispectral data of WorldView-2, medium resolution Landsat thermal data, and ASTER Global DEM derived slope data using OBIA into one workflow.
It is common in glacier mapping studies to evaluate the classification accuracy by comparing the derived glacier boundary with the manually digitised outlines or global inventory of glacier outlines [18,21,57]. Therefore, the glacier boundaries obtained from the classification of different datasets and glaciers (denoted by yellow colour) were compared in this section against the reference datasets of the manually digitised (MD) glacier boundary (denoted by maroon colour) and the manually corrected vector boundary of RGI version 6.0 (denoted by blue colour) as shown in  Visual analysis of the boundaries of the glaciers and datasets depict the proximity of the multisource OBIA-generated glacier boundary with the MD reference dataset. However, a mismatch between the OBIA-generated glacier boundary and the RGI 6.0 was observed in some places. This deviation may be attributed to the positional errors in the RGI boundary, as described by Pfeffer et al. [58]. The comparison of the glacier boundaries computed from classification with the reference boundaries was also performed quantitatively using the error estimators of bias and percent error (Table 8) as explained below.
we present an automatic method based on a combination of HSR multispectral data of WorldView-2, medium resolution Landsat thermal data, and ASTER Global DEM derived slope data using OBIA into one workflow.
It is common in glacier mapping studies to evaluate the classification accuracy by comparing the derived glacier boundary with the manually digitised outlines or global inventory of glacier outlines [18,21,57]. Therefore, the glacier boundaries obtained from the classification of different datasets and glaciers (denoted by yellow colour) were compared in this section against the reference datasets of the manually digitised (MD) glacier boundary (denoted by maroon colour) and the manually corrected vector boundary of RGI version 6.0 (denoted by blue colour) as shown in Figures 10-12.
The total glacier area from the multisource OBIA approach is found to be 22.50 km 2 . This leads to an area mismatch of 0.41 km 2 and 0.25 km 2 compared to the RGI 6.0 and MD glacier boundary, respectively. It can be inferred from Table 8 that the percent error reduced drastically from 13.62% to 1.79% (in the case of RGI 6.0 boundary as the reference) and from 9.85% to 1.09% (in the case of MD boundary as the reference), which amounts to a match of 98.21% and 98.90%, respectively. This high accuracy depicts that the proposed multisource OBIA approach can be relied upon to map the glacier boundaries with high precision. •

LISS-4 + Ancillary Dataset
The total glacier area measured from the multisource OBIA approach is 18.71 km 2 . A comparison with the RGI 6.0 and MD glacier boundaries indicated the bias of 0.68 km 2 and 0.63 km 2 , respectively. This leads to a percent error of 19.39% in the case of RGI 6.0 boundary as a reference and 19.34% in the case of MD boundary as a reference, resulting in a corresponding match of 80.61% and 80.66%.

Bara Shigri Glacier
The total area of Bara Shigri Glacier (lower ablation zone) appraised from the multisource OBIA approach was estimated to be 2.71 km 2 , and from the reference datasets of RGI and MD was measured correspondingly as 3.25 km 2 and 2.74 km 2 . This amounts to the individual bias of 0.54 km 2 and 0.03 km 2 , leading to respective percent errors of 16.62% and 1.09% tallying separately to a match of 83.38% and 98.91%.
We would like to point out that most of the available glacier boundaries provided by RGI have significant inaccuracies across the debris-covered part of the glacier, which means that the boundary between debris-free and debris-covered glacier is not precisely distinguished. The uncertainty reported in RGI data is ±5% [58]. Further, it is an apparent phenomenon in remote sensing-based studies of glacier mapping (especially debris-covered glaciers) that differences in acquisition dates and spatial resolution of remote sensing data bring new challenges. However, we presented an automated delineation of debris-covered glacier boundaries with improved accuracy.

Comparison with Other Debris-Covered Glacier Mapping OBIA Methods
The complexity of mapping fine details increases with the spatial resolution and also with the complex landforms such as debris-covered glaciers. This study is a step toward the detailed mapping of the debris-covered glacier from the multisource dataset, which includes not just the commonly classified classes such as snow/ice, IMD, SGD, glacial lakes, PGD and valley rock but also the supraglacial features such as crevasses and EIFs and the geomorphological features such as debris cones and rills.
Though OBIA was found suitable for mapping ground objects using the HSR images [9], some studies on the use of OBIA for the classification of the debris-covered glacier from lower spatial resolution data were also reported. Therefore, the real usage of OBIA to map minor but significant geomorphological features from HSR data might not have been explored. The examples are classification of glacier cover classes using ASTER and Landsat TM and Landsat Enhanced TM Plus data by Rastner et al. [12], Robson et al. [18], Sharda and Srivastava [15], and Ahmad and Baig [16].
A few studies demonstrated the use of OBIA in mapping glacier cover classes of debris-covered glaciers from HSR remote sensing data [13,14,19,33]. Since these studies were focussed on discriminating merely clean ice (i.e., snow/ice) and debris-covered ice (SGD) [14] or were limited to a few classes such as glacial lakes/ponds [13,19] or EIFs [33], spectral information was sufficient. Nonetheless, OBIA has the potential to discriminate classes using spatial attributes. Additionally, the mountain shadows and some valley rock areas needed manual corrections, making it a tedious task [21]. We could automatically mask shadowed regions from the WorldView-2 + ancillary dataset using the shadow detection index and brightness temperature rulesets for both the glaciers and from the LISS-4 + ancillary dataset using the brightness attribute with accuracies lying in between 96.0 and 100%. Furthermore, to the best of our knowledge, no work has been carried out for mapping smaller but necessary glacial units such as debris cones and rills from the visible and NIR spectrum. This study is unique in extracting debris cones and rills from HSR imagery using OBIA.
Recently, studies developed machine and deep learning approaches for debris-covered glacier mapping [57,59,60]. The most critical drawback of machine learning and deep learning approaches is that they highly depend on the training dataset collected for a given study area. Moreover, the deep learning approaches require sufficient training data to train the model.

Constraints and Potentials
This study admits the lack of ground truth for validation. However, in highly inaccessible regions of the world, where remote sensing is the most viable source of images for mapping, there is a lack of ground truth. Moreover, sometimes, it is expensive to conduct fieldwork. In such situations, researchers are forced to use less than ideal data or very fine spatial resolution imagery for reference purposes [50]. Nevertheless, the fairness with which the reference samples are collected in this study, such as taking care of the size of the reference sample units, accentuates the rigorous validation conducted. Further, this study followed the area-weighted error matrix which is particularly proposed for the accuracy assessment of the classification performed using OBIA.
The multiscale segmentation followed in this study effectively came up with the hierarchical classification, which could establish the contextual relationship among the objects. However, the multiscale segmentation depends on the scale parameter which depends on the spatial resolution of the imagery; hence, the scale of analysis and other segmentation parameters adopted for the WorldView-2 + ancillary dataset had to be revised with the LISS-4 + ancillary dataset. Moreover, it is often impossible to determine the correct scale of analysis in advance because different kinds of images require different scales of analysis [54].
Therefore, the transferability of the proposed multisource OBIA approach in other debris-covered glacier areas may require a thorough understanding of the basis behind each ruleset applied for sequential retrieval of various glacier cover classes. Similarly, the brightness temperature and slope thresholds defined for mapping the SGD, PGD, and valley rock in the Gangotri Glacier had to be altered with the changes in the local temperature regime and topography of the Bara Shigri Glacier [12,22]. However, with the high spatial and spectral capabilities of WorldView-2, this study re-emphasises the significance of deriving spectral and spatial attributes unique to the sensor to retrieve the utmost information from the HSR data. While the logic followed to extract the glacier cover classes would remain unchanged, the rulesets defined in this study could be portable to other debris-covered glaciers with scene-dependent, resolution-dependent, and seasondependent modifications, a widely accepted finding in glacier mapping [2,31].

Conclusions
This study attempted to explore the potential of OBIA to utilise the spectral and spatial characteristics of HSR data with ancillary information for the effective mapping of multiple glacier cover classes. The proposed multisource OBIA approach permitted the creation of segments of the desired shape and size by integrated evaluation of the information from different sources such as optical, thermal, and slope data. Allowing high weights to the brightness temperature and slope in the segmentation step facilitated the formation of optimal segments. Specifically, a contrast in the response of the SGD and PGD in the thermal infrared wavelength region and slope between the PGD and valley rock helped accurately discriminate between the SGD, PGD, and valley rock. This study made it plausible to use spectral properties and ancillary inputs effectively for distinguishing snow/ice, IMD, and the SGD without depending on shortwave infrared bands. Other significant highlights of this study are the classification of geomorphological features and the independence of manual corrections for extracting shadows. There is no full-proof method for classification using remote sensing data, not even LULC classification, which suffers from minor errors and needs manual rectification to present the classification more transparently. In this paper, no manual correction was adopted.
Moreover, the OBIA-generated boundaries follow the glacier boundary of the manual digitisation quite closely. Using the pan-sharpened band of WorldView-2 increased the accuracy of manually demarcated glacier boundaries. The importance of the proposed OBIA approach lies with the extraction of geomorphological features such as glacial lakes, EIFs, debris cones, rills, and crevasses using the spectral and spatial attributes. This contribution is novel in the context of the Indian Himalayas. To our knowledge, no study has contributed to exploiting the area-weighted error matrix for assessing the accuracy of glacier cover classes mapped using the OBIA approach. Though this study lacks field validation, the fair procedure adopted to collect the testing samples followed by visual analysis, evaluation of spectral signature plots, and appraisal of glacier boundary justifies the reliability of the proposed OBIA approach. The availability of good quality DEM and HSR thermal infrared data may further enhance the quality of results. Further improvements can be focussed on augmenting machine learning methods with OBIA [55] which can likely refine the identification of recurring patterns and textures of the debris-covered glacier features.