Exploration of Glacial Landforms by Object-Based Image Analysis and Spectral Parameters of Digital Elevation Model

Glacial landforms are a significant element of landscape in many regions of Earth. The increasing availability of high-resolution digital elevation models (DEMs) provides an opportunity to develop automated methods of glacial landscape exploration and classification. In this study, we aimed to: 1) identify glacial landforms based on high-resolution DEM datasets; 2) determine relevant geomorphometric and spectral parameters and object-based features for the mapping of glacial landforms; and 3) develop an accurate workflow for glacial landform classification based on DEM. The developed methodology included the extraction of secondary features from DEM, feature selection with the Boruta algorithm, object-based image analysis, and random forest supervised classification. We applied the workflow for three study sites: one in Svalbard and two in Poland. It allowed the identification of six categories of glacial landforms: till plains, end moraines, hummocky moraines, outwash/glaciolacustrine plains, valleys, and kettle holes. The majority of relevant secondary features represented DEM spectral parameters calculated from 2-D Fourier analysis. The supervised classification models with the highest performance exhibited up to 96% overall accuracy with regard to a groundtruth dataset. This study showed that glacial landforms can be identified using novel image-processing methodology and spectral parameters of high-resolution DEM. The complete classification workflow developed herein provides a solution for the transparent generation of thematic maps of glacial landforms that may be reproducible and transferrable to various glacial regions worldwide.

technological improvements in data collection and processing in recent years allow an accurate description of landforms, providing a basis for the study of the origin of different types of terrain. High-resolution satellite-derived or Light Detection And Ranging (LiDAR) digital elevation models (DEMs) form the fundamentals for detailed inspections of terrain attributes to quantitatively and qualitatively describe the Earth's surface. Novel numerical modeling and machine learning techniques provide reliable and meaningful processing tools for remote sensing datasets, allowing the generation of thematic maps and drawing valuable geomorphological conclusions at various scales [1], [2]. In particular, the development of geomorphometry and feature extraction methodology is setting the current course in geomorphology [3]. Baker [4] included supervised classification, predictive modeling, and LiDAR mapping as some other trends in future geomorphology development. This article deals with all of the abovementioned disciplines, proposing innovative solutions for glacial geomorphology mapping that may be used worldwide.

A. State of the Art
The mapping of geomorphological structures has undergone substantial changes in recent years. It evolved from the manual digitization (vectorization) of the Earth's surface landforms [2] to the automatic or semi-automatic characterization of DEM using machine learning approaches [5]. Glacial landforms were manually determined based on multibeam bathymetry [6], [7], radar [8], [9], satellite imagery [10], or a combination of all these sources [11]. Whereas traditional results provide sufficient information for landform recognition and identification, the process of vectorization is very timeconsuming [12]. By comparison, automatic or semi-automatic approaches allow the rapid delineation of linear or polygonal characteristics, especially using object-based approaches. The advantages of automatic digitization are the lack of subjectivity, increased processing speed, repeatability, and consistency in interpretation [13]. The introduction of the object-based image analysis (OBIA) approach to geosciences has improved the objectivity of thematic mapping by changing the unit of analysis from pixels to objects [14]. Image objects have advantages, such as their own statistics, shape characteristics, and relations to other objects, and they are closely related to real-world objects [15]. When reference ground-truth samples are available, they are often used to classify and assess the accuracy of generated thematic maps [16]. The recent availability of high-resolution photogrammetric or LiDAR DEM for terrestrial terrains creates an opportunity for the automated classification and extraction of landform features from the landscape.
Terrestrial remote sensing studies usually benefit from the extraction of many features collected from different sensors, such as satellite-derived spectral signatures, spectral indices (e.g., Vegetation Indices and Geology Indices), and textural analysis [17]. However, it should be clarified that, in this research, we used the term "feature" as typical in remote sensing literature, describing the predictor variable for image classification [18]. Unlike in many other remote sensing literature, we use the term "spectral" as a descriptor of the DEM [19].

B. Rationale and Motivation
The recent advances in geomorphological mapping are largely being driven by a change from analog to digital techniques [20], significantly increasing the amount and size of available datasets [21]. It is therefore necessary to focus on the identification of new secondary features of DEM and automatic classification protocols providing precise thematic maps of landforms. Secondary features mean predictor derivatives that may be extracted from DEM (which is the primary feature), such as geomorphometric attributes, spectral parameters, or object-based variables [18]. In this study, we propose a comprehensive, repeatable workflow that allows the creation of reliable geomorphologic maps from highresolution LiDAR and satellite-derived photogrammetric DEM surfaces.
Feature selection methods can successfully adjust a feature set, improving thematic mapping accuracy in remote sensing studies. Unfortunately, feature selection methods have rarely been applied in geomorphological research with LiDAR or photogrammetric DEM [22]. Whereas some recent studies have applied the Boruta feature selection algorithm for the land cover mapping of forests [23]- [25], in this study, we propose that the Boruta technique is also suitable for geomorphological mapping based on LiDAR or photogrammetric datasets.
The identification and prediction of glacial landforms based on LiDAR datasets have rarely been studied using a classic pixel-based approach [26], [27]. Whereas OBIA has been occasionally applied for the determination of landslides [22], [28] and change detection [29], it has been rarely applied for the exploration of glacial geomorphology [23], [30] from high-resolution DEMs. Therefore, this research is one of the first to involve OBIA for this purpose.
Although the use of semiautomatic supervised classifiers in remote sensing studies is still growing [31], the random forest (RF) and support vector machine (SVM) techniques are currently among the most powerful methods for the accurate classification of geospatial datasets [32]. In geomorphological applications with high-resolution DEMs, RF has been applied for the investigation of, for example, beach geomorphology [33], fluvial geomorphology [34], and the modeling of barrier-island habitats [35]. However, to our knowledge, ours is the first study to employ an object-based RF classification methodology for the exploration of glacial landforms.

C. Research Objectives
The main objectives of this study included: 1) the identification of glacial landforms based on LiDAR and satellite-derived DEM; 2) the determination of relevant secondary features from geomorphometric features, spectral parameters calculated on the basis of 2-D Fourier analysis, and object-based shape features; 3) the development of an accurate workflow for glacial-landform classification based on high-resolution DEM; 4) the mapping of glacial landforms that is reproducible and may be transferrable to different study sites. We investigated the abovementioned aspects using the methodology described in Section II.

A. Study Sites
We analyzed three study sites with a well-developed glacial landscape: one located in the presently glaciated Arctic region (Svalbard) and two located in the area of northern Poland with glacial relief formed during the Late Pleistocene ( Fig. 1). We chose these study sites to check our workflow with various glacial landforms formed by a contemporary glacier and palaeo-ice sheet, based on high-resolution DEMs available for these regions.
The study site within the presently glaciated area is located in the Kaffiøyra region-a coastal plain in Oscar II Land at western Spitsbergen. There are several glaciers located in the Kaffiøyra, among which the biggest land-terminating glacier is the Elise glacier which contains a particularly well-developed marginal zone [ Fig. 1(A)]. Conspicuous glacial landforms located on the foreland of Elise glacier were formed as a result of glacier retreat from its maximum Little Ice Age (LIA) extent, which ended around Svalbard at the beginning of the 1920s [45]. Thus, the study site represents a very fresh glacial landscape, comprising landforms that were formed about 100 years ago. To date, the foreland of the Elise glacier was investigated by classical geomorphological mapping [37]- [39] and sedimentological investigations of glacial landforms in the field [40]. The study site consists of a conspicuous arc of end moraines, indicating the position of Elise glacier's margin during the LIA, extensive and relatively flat outwash plain on the foreland of the end moraines, as well as till plains, outwash plains, and hummocky moraines located between end moraines and present position of Elise glacier's margin [ Fig. 1(A)]. The end moraine belt is 180-500 m wide and up to 30-35 m high above the surrounding outwash plains. Till plains have undulated surface with elongated hills (flutes) mostly 100-200 m long and usually a few meters high. Hummocky moraines have undulating, irregular topography with hills and depressions. Hills are tens of meters wide and long, up to 5 m high, and are characterized by various degrees of elongation [see Fig. 1(A) and Table I]. The study sites with glacial relief formed in the late Pleistocene are located in northern Poland, within regions occupied by the last Scandinavian Ice Sheet (SIS) during the Last Glacial Maximum and exposed as a result of the subsequent ice-sheet recession. The first area [ Fig. 1(B)] is situated close to the Baltic Sea coast within the Gardno-Leba Plain, and the second area [ Fig. 1(C)] is located close to the maximum extent of the last SIS on the northern slope of the Lubawa Upland. The study site in the Gardno-Leba Plain consists of a conspicuous arc of end moraines formed during the recession of the last SIS as a consequence of local ice readvance ca. 15-14-ka BP [41], [42]. Other glacial landforms typical for marginal zones of the ice sheet, such as: till plains, outwash/glaciolacustrine plains, hummocky moraines, or valleys, also occur there [see Fig. 1(B) and Table I]. The most prominent landform in this study site is a 0.6-1.8 km wide and up to the 60-m high belt of end moraines.
Valleys are up to 500-600 m wide and up to 10-15 m deep. Hummocky moraines consist of irregular hills and depressions up to 8-12 m high/deep. Till plains are flat or slightly undulated, with relative heights which do not exceed a few meters. Outwash/glaciolacustrine plains are also flat or undulated (if secondary landforms such as paleochannels or young erosional incisions occur within them) with maximum local denivelation in order of 5-10 m (Table I). This glacial landforms assemblage is one of the most conspicuous geomorphological records of the palaeo-ice sheet marginal zone in northern Poland. So far, it was identified and analyzed mainly based on traditional field-based geological and geomorphological mapping [43] investigations of the structural sedimentology of sections exposed within the end moraines [44]. The study site in the Lubawa Upland consists of hummocky moraines, till plains, kettle holes, and valleys [see Fig. 1(C) and Table I]. Local denivelations within till plains do not exceed a few  Table I). The final shaping of these landforms is associated with the recession of the last SIS from its LGM position at ca. 22-18-ka BP [45]. This study site was previously mapped using traditional methods of delineation and interpretation of particular landforms based on geological data (boreholes and hand drillings) and topographic maps [46]. High-resolution DEM data were previously used only for a visual inspection and interpretation of the landforms visible in the study sites. However, the knowledge on the glacial landforms in the area is relatively solid due to the availability of geological and geomorphological maps or aerial and satellite images. Our approach using high-resolution DEM and semiautomated classification was used for the first time in these selected areas representing glacial relief formed by the contemporary glacier around a hundred years ago (Spitsbergen), and by palaeo-ice sheet ca. 22-14-ka BP (northern Poland).

B. Data Acquisition and Processing
This research benefits from elevation data from an open access photogrammetric ArcticDEM digital surface model of the Arctic [47] and LiDAR DEM available for the territory of Poland [48]. ArcticDEM was generated by applying stereo autocorrelation techniques to overlapping pairs of high-resolution optical satellite images. It was constructed from in-track and cross-track high-resolution (∼0.5 m) imagery acquired by the DigitalGlobe constellation of optical imaging satellites, mostly from the panchromatic bands of the WorldView-1, WorldView-2, and WorldView-3 satellites [47]. Stereo pair images were processed to DEMs using the Surface Extraction from TIN-based Searchspace Minimization software [49]- [51]. The resulting DEM data were downloaded from ArcticDEM explorer Website (https://livingatlas2.arcgis.com/arcticdemexplorer) as 2-m spatial resolution 32-bit GeoTIFF of mosaicked elevation data, blended and feathered to reduce void areas and edge-matching artifacts. Filtered IceSAT altimetry data have been applied to the raster files to improve absolute accuracy [47]. We applied ArcticDEM 2-m GeoTIFF elevation data to glacial landforms mapped on the foreland of the Elise glacier in Spitsbergen.
DEM data for the Gardno-Leba Plain and the Lubawa Upland study sites were downloaded from the database of the Head Office of Geodesy and Cartography in Poland (GUGiK) as a part of the ISOK flood protection project [52]. They were collected based on LiDAR scanning of the terrain surface conducted in Poland in 2011-2014 [48]. They represent bare surfaces without any objects, such as buildings, trees, or plants. The Airborne Laser Scanning (ALS) point cloud had a density of 4 points/m 2 , a vertical accuracy of ≤0.15 m and a horizontal accuracy of ≤0.50 m. Points within the ALS cloud were initially classified according to the type of objects they represented (e.g., ground, vegetation, and water). Then, based on ground and water points, the LiDAR DEM was created by the interpolation of point values to a density of 1 point/m 2 . Consequently, a 1-m DEM representing the elevation of the ground surface and water surface was downloaded as ASCII files (XYZ coordinates) from the GUGiK database. The original 1-m resolution DEM was used to map glacial landforms in the Lubawa Upland area (∼37 km 2 ), whereas a 5-m DEM based on resampling of original ASCII files was used for analyses in the Gardno-Leba Plain in order to speed up the calculation processes (especially spectral parameters) for this relatively extensive area (∼117 km 2 ) and to filter out some small-scale landforms originating from anthropogenic activity.
In this study, we followed the main assumptions of remote sensing image classification summarized by Lu and Weng [17]. A summary of the processing workflow explicitly developed in this study is provided in Fig. 2. Sections II-C to II-I contain methodological descriptions of all the steps of the workflow.

C. Ground-Truth Dataset
The ground-truth dataset was constructed based on visual inspections of DEMs, existing geological and geomorphological maps, orthophoto maps, and the identification and interpretation of distinct types of landform. Based on the  Table I, we sampled representative points indicating the location of a particular landform type within DEMs (Fig. 1).
For the foreland of Elise glacier, 167 samples were separated into four classes: till plains, outwash plains, end moraines, and hummocky moraines. The locations of all the ground-truth samples are provided in Fig. 1(A). In the Gardno-Leba Plain, the ground-truth dataset contained 246 point locations appropriate for the five landform categories: valleys, end moraines, till plains, outwash/glaciolacustrine plains, and hummocky moraines [ Fig. 1(B)]. Finally, for the Lubawa Upland, our ground-truth dataset contained 382 samples, divided into four categories: valleys, till plains, kettle holes, and hummocky moraines [ Fig. 1 Once determined, the ground-truth sample dataset was divided into training and validation (test) subsets, allowing the training of the supervised classifier and measurement of its performance in terms of accuracy. We separated the training and validation subsets with the ratio 70/30 with the Subset Features procedure in ArcGIS (Environmental System Research Institute (ESRI), 380 New York Street, Redlands, CA 92373, USA). The procedure allows for random separation of the point shapefile with an unbiased analysis.

D. Feature Extraction
In this study, we extracted 21 types of secondary features from DEMs. A list of all the extracted features is provided in Table II. In Sections II-E-II-I, we describe the methods for feature extraction with pixel-based and object-based separation.

E. Pixel-Based Features
All the pixel-based geomorphometric features were created using a 3 × 3 neighborhood window size in ArcGIS. The slope and aspect of DEMs were calculated using the methodology described by Burrough and McDonell [53]. The curvature was calculated using the algorithm proposed by Zevenbergen and Thorne [54]. Curvature was calculated perpendicularly to the slope and is expressed as the planar curvature, whereas the profile curvature was computed parallel to the slope [54]. The DEM roughness was calculated using the Surface Area to Planar Area (SAPA) algorithm proposed by Jenness [55].
Pixel-based spectral features were calculated using our own algorithms developed in MATLAB (MathWorks, 1 Apple Hill Drive, Natick, MA 01760-2098, USA). They were generally computed based on 2-D Fast Fourier Transform (2-D FFT). In order to avoid spectral leakage, 2-D spectra were multiplied by a function of Discrete Prolate Spheroidal (Slepian) Sequences (DPSS). The algorithm was executed using moving window sizes from 80 × 80 m to 500 × 500 m, with 90% of the window overlapping depending on the resolution and extent of the considered dataset. The 2-D FFT allowed the determination of cross sections producing 1-D spectra, every 5 • , in a range of 0 • -180 • , allowing the further determination of averaged spectral parameters according to Tegowski et al. [56].
For each of the 37 spectra, we calculated two spectral moments (m 0 and m 2 ) [57]. Whereas ω 0 represents the mean frequency of the spectrum, the concentration of the spectral energy density around the mean frequency is expressed as the spectral width (ν 2 ). The spectral skewness defined for central moments (γ s_centr ) is a parameter that may represent slight changes of the DEM surface [58]. The fractal dimension (Dfft) is the parameter determining the fractal characteristics of the LiDAR dataset [59]. Being a measure of spectrum peak "sharpness," the quality factor (Q) is a parameter representing combined characteristics of spectral moments [60]. Detailed descriptions of the spectral parameters of the DEM used in this study and the necessary formulas are provided by Trzcinska et al. [19]. In general, spectral parameters refer to the spatial frequencies that build up the total shape of rough terrain. Therefore, the spectral parameters characterize the scale and degree of surface roughness.

F. Object-Based Features
Object-based features were calculated using our rulesets developed in eCognition (Trimble, Inc., 935 Stewart Drive Sunnyvale, CA 94085, USA). They are provided as mean values per image object, created using a multiresolution segmentation algorithm with a scale factor of 25 or 40. The used segmentation methodology, with its parameters, is described in more detail in Section II-H.
All the extracted object-based features represented the shape properties of the image objects. Asymmetry described the dependence between the length of an image object relative to that of its regular polygon based on the ellipse surrounding the image object. A similar feature, Length/Width, represented the smallest of the relations of the eigenvalues of the covariance matrix or the Length/Width ratio calculated from the bounding box of the image object. The compactness measured the ratio between the magnitude of the image object's length and width to its area, expressed in pixels. It can take values from 0 to ∞, where 1 means identical compactness. The smoothness of an image object's border was described by the Shape Index feature; the lowest value it can take is 1, and a lower value means a smoother border. The comparison of an image object with an ellipse was performed according to the Roundness feature, taking values from 0 to ∞, where the lowest values mean ideal roundness. The similarity of an image object to an ellipse of similar size and shape was estimated by another measure, Elliptic Fit. The last object-based feature extracted in this study was Rectangular Fit, expressing the similarity of an image object's shape to a rectangle of similar size. A detailed description of the object-based features used in this study, with reference to specific formulas, is provided in eCognition's Reference Book [61].

G. Cross Correlation and Feature Selection
The initial data analysis includes cross correlation, which was the first step before feature selection. It was performed in advance to preclude highly correlated features from further processing so that they did not make a significant contribution in relation to other features. We performed cross correlation using the "caret" library in the R statistical software [62] (R Foundation for Statistical Computing, Vienna, Austria, Wirtschaftsuniversität Vienna, Welthandelsplatz 1, 1020 Vienna, Austria). Features with Pearson's correlations higher than 0.75 were considered highly correlated.
The relevance of all the extracted features for supervised classification was determined with the Boruta feature selection algorithm [63]. The algorithm, based on the RF classifier, allows the selection of the most important features, providing their importance or Z scores in relation to other irrelevant features generated as "shadow" variables [64], [65]. We used the R implementation of the algorithm provided in the "Boruta" library. A detailed description of the working principles of the used feature selection method is given by Kursa and Rudnicki [63]. We performed the Boruta algorithm based on sets of features before and after cross correlation.

H. Object-Based Image Segmentation and Classification
OBIA was the basis for the generation of object-based features and the determination of image objects with similar characteristics in terms of the DEM. As resulted from the processing workflow developed in this study, the properties of all the features were calculated based on the mean values of image objects belonging to ground-truth categories (see the process "export object statistics" in Fig. 2). Image objects were generated using the multiresolution segmentation algorithm proposed by Benz et al. [15] in eCognition. The algorithm uses an iterative, bottom-up merging technique to create image objects suitable for delineating specific features in remote sensing images (in our case, high-resolution DEM). Benz et al. [15] defined the parameters of multiresolution segmentation that affect the size and shape of the final image objects, such as the scale, shape, and compactness. The scale is a dimensionless parameter based on the relative sizes of image objects. The two other parameters represent dependences between image object properties such as color, shape, smoothness, and compactness. A detailed description of the multiresolution segmentation's working principles is provided by Benz et al. [15]. All further analyses were performed using image objects (instead of pixels) as the primary units of analysis.
Image classification was performed using a supervised approach related to decision trees and machine learning. Of several tested supervised classifiers (Classification and Regression Trees, SVM and K -Nearest Neighbor), SVM and RF are recommended, often providing the highest accuracy according to remote sensing literature [31], [66]. The principle of SVM is to separate the dataset after transformation into a multidimensional feature space where it seeks the best fitting hyperplane between the classes [67]. Being an ensemble classifier, RF generates multiple decision trees (often called forests) based on a randomly selected subset of ground-truth training samples and features [64]. In this study, the RF classifier was trained using a training subset of ground-truth samples and selected important features. We used the RF implementation in eCognition [61].

I. Accuracy Assessment
The accuracy of the supervised classification was calculated based on an error matrix and accuracy statistics using a validation subset of ground-truth samples [16]. The calculated accuracy statistics include the user's and producer's accuracy, omission error, commission error [68], [69], overall accuracy, and Cohen's Kappa coefficient [70]. The user's accuracy means the proportion of total ground-truth samples that were adequately classified in relation to a particular class. The producer's accuracy may be explained as the proportion of all classified pixels in the remote sensing image that were classified correctly. Whereas the overall accuracy is the sum of all the correctly classified ground-truth samples in relation to the total number of samples, Kappa is the other measure of overall accuracy that includes the possibility of misclassification by chance. Accuracy was assessed using eCognition [61]. Moreover, we calculated error matrices and accuracy assessment statistics for quantitative comparison with manual separations of all study sites. Our compliance with manual classification was calculated in ArcGIS based on cross-tabulation between automatic and manual classification areas. The results were expressed in percentage values. Manual classification of landforms in the foreland of Elise glacier was based on available geomorphological maps of the area [39] and visual interpretation of the DEM and orthophoto map (Fig. S1A). For the Gardno-Leba Plain and the Lubawa Upland, this interpretation was based on available maps of surface deposits (Fig. S1B, C) combined with a visual inspection of LiDAR DEMs.

A. Feature Extraction and Selection
Data analysis indicated that 10 of 15 extracted features were relevant for the foreland of the Elise glacier study site. Similarly, 10 of 12 features were selected as important for Gardno-Leba Plain. Features for these scenarios were selected after initial cross correlation. In comparison, 24 of 30 extracted features were relevant for the Lubawa Upland study site (Table III). These features were selected without cross correlation. We trained classifiers based on datasets before and after cross correlation, then we selected the results with the highest accuracy. The lists of all important features for the three considered study sites were provided in Table III. The importance of all selected features was presented as boxplots in Fig. 3. The selected features were used for further supervised classification.

B. Object-Based Image Segmentation and Classification
Object-based image segmentation was performed using the following parameters for multiresolution segmentation: shape, 0.1; compactness, 0.5. Whereas the multiresolution segmentation scale for the foreland of Elise glacier and the Gardno-Leba Plain was 40, the same parameter for the Lubawa Upland was 25. In all scenarios, the segmentation algorithm was applied based on a single DEM feature.
SVM and RF provided the best classification performances confirmed by accuracy assessment results from several tested supervised classifiers. Whereas the RF classifier was the most appropriate for the foreland of the Elise glacier scenario, SVM had the highest performance for the Gardno-Leba Plain and the Lubawa Upland case studies. Fig. 4 shows close-up results list of representative glacial landforms marked in areas provided in Fig. 1. In this outcome, we provided a comparison between the original DEM and the expected versus obtained results. Whereas the expected column in Fig. 4 refers to a manual classification of landforms based on existing satellite images, geological/geomorphological maps (Fig. S1), and visual interpretation of DEM, the obtained column is related to semiautomated results of the supervised classification performed according to our processing workflow.
We provided a comparison between hillshade reliefs, manual maps, and semiautomated maps for all study areas ( Fig. 5; A-the foreland of Elise glacier; B-the Gardno_Leba Plain; C-the Lubawa Upland). Maps presented in Fig. 5 include polygon data related to other landforms (nonglacial), water bodies, and glacier ice (see Fig. S1). Visual comparison of the classification results from Figs. 4 and 5 with spatial distributions of ground-truth samples (Fig. 1) suggests a high correlation for all scenarios.

C. Comparison With Relevant Spectral Parameters
From six relevant spectral features selected by the Boruta algorithm in the foreland of Elise glacier, spatial variability of spectral moment m 0 calculated in window sizes 160 m largely follows the spatial distribution of the main glacial landform types. High values of m 0 are correlated with diverse topography within the arc of end moraines and also partly within till plain with flutes, whereas low values of m 0 occur across much flatter areas of outwash plains (Fig. 1A, S2). Quality factor Q (calculated for 100-and 300-m windows) reveals high variability within all identified landform types, whereas spectral skewness defined for central moments γ s_centr (calculated for 100-, 200-, and 300-m windows) has high values for all landforms, with the highest variability and locally low values within flat areas of outwash plains (Fig. 1A, S3).
In the Gardno-Leba Plain, high values of spectral moment m 0 calculated for 500-m window sizes are also correlated with diverse topography within end moraine arcs and partly within hummocky moraines and along valley slopes, in contrast to spectral skewnessγ s calculated for a 250-m window, which is the lowest in areas of diverse topography. Flat areas of outwash/glaciolacustrine plains show the lowest  (Fig. 1B, S4). Mean frequency ω 0 and quality factor Q (calculated for 250-and 300-m window) reveal high variability within all glacial landforms identified in this study site.
The spectral moments m 0 and m 2 (calculated for window sizes 80 and 140 m) in the Lubawa Upland show high values within valleys and other local depressions such as kettle holes or parts of hummocky moraines (Fig. 1C, S5). The high value of the mean frequency parameter ω 0 means the presence of small forms, whereas its low value means large ones, with respect to the moving window size of 80 m. Therefore, high values for the 80-m window size mainly occur in valleys (being relatively small). Low values of spectral skewnessγ s and γ s_centr occur in valleys and kettle holes for both 80-and 140-m moving windows, high values were recorded mainly within till plains. The spectral width ν 2 and quality factor Q show high values mainly within valleys and till plains, whereas low values of these parameters occur mainly within hummocky moraines (Fig. 1C, S5).

D. Accuracy Assessment
The error matrices and accuracy assessment statistics for the SVM and RF models were presented in Tables IV-VI.  Tables VII-IX provide error matrices and accuracy statistics (Table IV). The other results showed overall accuracy and Kappa scores from 79% to 86% (Tables V and VI). Whereas the highest compliance with manual separation for the foreland of Elise glacier scenario reached 90% (Table VII), the overall accuracy statistic for two other areas indicated 63%-65% correctness (Tables VIII and IX). User's accuracy for individual landform types in Gardno-Leba Plain varies from 79% for end moraines and hummocky moraines to 100% for valleys (Table V), and in the Lubawa Upland from 79% for valleys to 100% for kettle holes (Table VI).
In general, these quantitative results show high to very high correspondence with ground-truth samples and lower compliance with manual geomorphological classifications.

A. Reference to the Main Hypothesis and Research Problems
This study confirmed that glacial landforms, such as till plains, end moraines, hummocky moraines, outwash/ glaciolacustrine plains, valleys, and kettle holes, can be identified using novel image-processing methodology and spectral parameters of high-resolution DEM. From a total of 122 extracted variables (62 for the foreland of Elise glacier, 30 for the Gardno-Leba Plain and the Lubawa Upland), we determined relevant types of spectral features of the DEM that may be diagnostic for the thematic mapping of glacial geomorphology. Besides the DEM, they include four types of geomorphometric features, eight types of spectral features, and four kinds of object-based features. The complete workflow for glacial landforms classification developed in this study allowed for the quick and transparent analysis of DEM. It provided a solution for the generation of thematic maps of glacial geomorphologic features that is reproducible and may be transferrable to other areas of interest.

B. Research Relevance in Relation to Other Works and Interpretations of the Main Findings
In comparison with other applications of the RF algorithm in remote sensing studies [33]- [35], the classifier's performance in this work was amongst the highest, suggesting high relevance of RF for investigations of glacial geomorphology based on DEM surfaces.
According to our best knowledge, spectral parameters were not extracted from LiDAR and satellite-derived DEM before. Object-based features are the other relevant group of features selected in this study, confirming their usefulness in light of the very few OBIA glacial geomorphology studies based on LiDAR datasets [30], [71]. Finally, the relevance of geomorphometric features, confirmed by their high importance, supplements current trends in the development of geomorphological science [3], including evaluation of geomorphometric attributes [72].
Moreover, the current trends also involve the technological development of automatic mapping methods presented for a small section of glacial landforms in this study [20], [21]. This method avoids the disadvantages of manual digitization, such as subjectivity and a lack of repeatability [71].
User's accuracy for individual landform types of the foreland of Elise glacier varies from 88% for hummocky moraines to 100% for outwash plains and end moraines (Table IV). The comparison with manual separation of particular landform types shows the slightly lower but still very high accuracy of our semiautomated classification-overall accuracy of 90% (Table VII). Robb et al. [71] used OBIA to the semiautomated classification of four types of glacial landforms (moraines, eskers, outwash, and till plains) on the foreland of Breiðamerkurjökull in Iceland, and they obtained an overall accuracy against the manual mapping of 77%. Our level of accuracy indicates that the proposed workflow for glaciallandform classification works particularly well in presently glaciated areas with young glacial landforms.
The obtained levels of accuracy in the Gardno-Leba Plain and Lubawa Upland relative to the ground-truth dataset are lower than on the foreland of Elise glacier in Spitsbergen but still high enough (≥80%) to consider our workflow as applicable for mapping the glacial geomorphology of formerly glaciated regions. The level of overall accuracy concerning reference datasets in the order of ca. 80% is often considered in remote sensing analyses as enough for reliable automatic maps [73]- [76]. The comparison with manual separation of landform types (manual maps) shows an overall accuracy of 63% in the Gardno-Leba Plain (Table VIII) and 65% in the Lubawa Upland (Table IX). Lower accuracy with respect to manual maps is understandable because traditional geomorphological mapping is based not only on an analysis of morphological patterns but also on the spatial distribution of surface sediments/rocks when designating particular landforms. While the approach we present here is a semiautomated classification based on morphological characteristics, secondary features, and ground-truth points but without additional input of surface sediments/rocks. However, analysis from other regions of Pleistocene glaciations such as North America shows the accuracy of supervised classification of glacial landforms with respect to traditional geological maps in order of 51%-61% [77].
In addition, the influence of three different DEMs resolutions on outputs of supervised classification is worth discussing. The highest accuracy was reached for the DEM, separated for four landform types with 2-m spatial resolution. Whereas an increase in spatial resolution often complicates the classification process, it is not surprising that the other DEM with a higher spatial resolution (Lubawa Upland) and the same number of classes had a lower accuracy. Besides, our results suggest that increasing classification classes may have a greater influence on the accuracy than a lower resolution of spatial datasets, which is visible especially in the Gardno-Leba plain. In this scenario, the DEM had the lowest spatial resolution while separated into five classes. It is also worth mentioning that all calculations were performed regarding the object-based approach, which by design should be independent of the spatial resolution of datasets [14], [15]. Therefore, our results seem to support the fact that classification complexity had a much more significant effect on the accuracy of the results than the spatial resolution of DEM datasets in the object-based approach.

C. Potential Sources of Errors
Accurate thematic maps may obviously be obtained when ground-truth samples have previously been determined in a precise manner. This is crucial because all the accuracy statistics were calculated in relation to the ground-truth-sample test dataset. Whereas ground-truth samples can be obtained through, for example, field work or expert knowledge, their determination should be as precise as possible. It should be noted that the interpretation of landform types for the ground-truth dataset has a degree of subjectivity. It is reduced by increasing the number of samples or increasing the use of high-resolution DEMs. The quality of the ground-truth dataset is higher when training points are inserted based on manual mapping using the most reliable remote sensing data in particular regions, e.g., orthophoto maps in unvegetated regions or high-resolution LiDAR DEMs representing a detailed picture of the terrain surface, even within areas with rich vegetation [78]. Ideally, ground-truth samples should be designed in a random, representative way and not spatially autocorrelated. At the same time, from a statistical point of view, they should be designed as large as possible and capture all occurring landform types. Whereas in the case of spatially extensive landforms, determination of such points does not usually cause any problems, it may be challenging for smaller, geomorphologically distinct forms. Examples of such differences in our study include, e.g., outwash/glaciolacustrine/till plains versus kettle holes and hummocky moraines. In order to ensure statistical significance, samples in classes occupying small areas had to be more densely clustered, thus exposing themselves to significant spatial autocorrelation. In this work, local spatial autocorrelation analysis (LISA, [79]) has shown that in the case of hummocky moraine class in the foreland of Elise glacier, samples of low spatial diversity were surrounded by other samples of low spatial diversity. In the case of the kettle hole class in the Lubawa Upland study site, the spatial autocorrelation of samples was statistically insignificant.
A closer look at the spatial distribution of kettle hole and valley classes in Lubawa Upland shows that some validation samples of kettle holes were found within valleys. On the other hand, no samples of other classes were found within the kettle holes classified areas (Table VI). This suggests that the kettle hole class should have a slightly larger surface area instead of the valley class. The other potential sources of errors may be related to the occurrence of ground-truth samples nearby the border of landforms. Considering the same example, the boundary of kettle hole-valley classes in the one site occurred in less than 0.2 m (comparing to 1-m raster resolution) from the one misclassified kettle hole sample. It is worth mentioning that mistakes at the border of two classes can always happen, but they are usually few. Some attempt to deal with them is to determine the locations of data samples away from potential boundaries of given areas, but ultimately we have no direct influence on how the boundaries will be determined in the process of image segmentation, especially when the border has gradational or fuzzy character.

D. Limitations of Our Research
In comparison with conventional geomorphometric studies, our thematic map did not delineate landforms as polylines along their crests or thalwegs [2], [80], but instead areas (polygons, see [10], [11], [81]) where geomorphometric features occur with high probability. This was caused by the working principles of the object-based methodology. The proposed approach is similar to land cover or habitat mapping and can be used as a diagnostic tool for the further (qualitative) delineation of landforms and decision-making.
Although the precision may affect the accuracy of classification [82], the details of how ground-truth samples should be defined appropriately were beyond the scope of this research. Another problem related to ground-truth samples is the training/test separation ratio used for image classification. Although the ratio of 70/30 used in this research is common in remote sensing studies [83], [84], some researchers emphasize the need for its further examination [82], [85], [86].
The proposed procedure produces maps with delimited areas of given geomorphometric features and identifies given glacial landforms. However, in terms of genetic interpretations of glacial landforms (e.g., linking their geomorphometric features with morphogenetic processes), it is essential to bear in mind that various geomorphological processes may produce morphometrically similar or even identical landforms [87]- [89]. For example, although different processes stay behind the formation of till plains and outwash/glaciolacustrine plains, their morphology (surface relief) may be very similar or even the same (Table I). Moreover, hills and hummocks produced in different glacial settings may also have similar morphology. In the foreland of Elise glacier, we identified landforms of comparable shape and dimensions belonging to end moraines, hummocky moraines, or till plains. Hills with similar length, width, and elongation ratio as well as height and cross-profiles were detected within end moraines belt and hummocky moraines or within end moraines belt and undulation of till plains [ Fig. 6(A)]. Also, undulations occurring within some parts of till plains and outwash/glaciolacustrine plains in the Gardno-Leba Plain are comparable, creating similar relief [ Fig. 6(B)]. Thus, except for pure geomorphometric parameters, the geological and geomorphological context of particular landform occurrences is always crucial for their adequate identification and interpretation. Information such as the spatial context of landforms distribution, their inner structure, or types of sediments/rocks may be critical for adequate interpretation. Geological investigations (e.g., classical landforms and sediments mapping in the field) may provide new details not provided by surface analysis. Therefore, the compliance of our semiautomated classification with manual landform mapping is in the order of 63%-65% for formerly glaciated regions. However, we have to bear in mind that classical geomorphological maps are also based on subjective interpretation and the detail of landform mapping depends on the scale of the mapping process, lists of landform types separated during mapping, author's experience, etc. Thus, manual separations are only an approximation of reality and should not be treated as an exact reference for automated glacial landforms extraction.
We observed a significant difference in accuracy level for the classification of glacial landforms in the presently glaciated arctic region (overall accuracy of 96%) and the region of Pleistocene glaciations (overall accuracy of 84%-86%). Lower accuracy in the region of past glaciations is probably related to a significant postglacial and anthropogenic transformation of relief captured on high-resolution DEMs [ Fig. 6(B)]. Therefore, our new methodology is particularly suitable for regions of present glaciation, where glacial landforms are young and clearly recognizable in the landscape.

E. Recommendations for Future Research
Although this research included the development of an RF and SVM supervised classifiers that performed with the highest accuracies, other classifiers were also tested. This corresponds with other authors suggesting the evaluation of different methods of classification [32], [90]. Therefore, we also recommend testing different methods of supervised classification, such as Decision Trees [91], Artificial Neural Networks [92], etc.
Whereas this research benefited from the 70/30 proportion for the training and test subsets of ground-truth samples, future research may consider tests of other separation ratios, e.g., 50/50 [93], [94], 60/40 [95], [96] or 80/20 [97]. The other approaches may include resampling methods, as proposed by Lyons et al. [85] and investigation on the impact of the spatial arrangement of ground-truth samples for classification accuracy [98].
The spectral parameters could potentially be extracted for other digital surfaces using different moving window sizes. We used nine different window sizes (from 80 to 500 m), and our feature selection results confirmed the almost all of them were relevant for classification. In any case, we recommend extraction from high-quality spatial datasets and testing different moving window sizes, considering tailoring them to the size of the expected glacial landforms, bearing in mind the significant increases in processing time and necessary computing power for smaller window sizes. However, once spectral parameters and other secondary features were developed, single running of key processes (segmentation, classification, and accuracy assessment) required several minutes for the study area. Possible further applications of the methodology developed in this study include other glacial landforms, such as drumlins, eskers, kames, tunnel valleys, marginal channels, that do occur in other areas in the European lowlands or beyond.
We appreciate the potential of this classification method for both land and sea applications [19]. The future potential application of the method (taking into account the increase in the development of information technology) includes a quick screening of DEM from different sources. It provides fast, accurate information about landforms supporting further interpretation and decision-making.

V. CONCLUSION
This study considered research and methodological approaches that may be used for the mapping of glacial landforms. The main assumptions can be further transferred and support other crucial remote sensing applications, like land-use decision making, that requires increasing processing time and computing power. An essential part of geomorphological studies, the automatic method of thematic mapping based on LiDAR and photogrammetric datasets is a response to the current needs of the research community. Due to the development of computer science, the technical possibilities of data processing have increased over time. Therefore, we believe that automatic methods of classification should be supported and developed.
The general achievement of this study was the development of a novel method for the mapping of glacial landforms. This unique method included the use of OBIA, RF, SVM, DEM spectral parameters, and Boruta algorithms for the accurate exploration of glacial geomorphology in the three study sites in Northern Poland and Svalbard. We identified and classified the following glacial landforms: till plains, end moraines, hummocky moraines, outwash/glaciolacustrine plains, valleys, and kettle holes. Our feature extraction and selection confirmed the relevance of 44 of 122 features, including spectral parameters, as primary diagnostics for the mapping of glacial landforms. For all investigated study sites, spectral parameters such as skewness, moment m 0 , and quality factor were relevant for the identification of analyzed glacial landforms.
The methodological solutions proposed in this article can be easily adapted and utilized for other landscapes, either present and from the past. The main advantages are reliability, transferability, and repeatability, giving transparent results.
Once the model is generated, the other advantage is the fast computation for a particular area. The principal requirement is a high-quality DEM surface and precise ground-truth-sample dataset.