The Impact of Digital Elevation Model Preprocessing and Detection Methods on Karst Depression Mapping in Densely Forested Dinaric Mountains

: Karst landscapes have an abundance of enclosed depressions. Many studies have detected depressions and have calculated geomorphometric characteristics with computer techniques. These outcomes are somewhat determined by the methods and data used. We aim to highlight the applica-bility of high-resolution relief laser scanning data in geomorphological studies of karst depressions. We set two goals: geomorphometrically to characterize depressions in different karst plateaus and to examine the inﬂuence of data preprocessing and detection methods on the results. The study was performed in three areas of the Slovene Dinaric Karst using the following steps: preprocessing digital elevation models (DEMs), enclosed depression detection, calculating geomorphometric characteristics, and comparing the characteristics of selected areas. We discovered that different combinations of methods inﬂuenced the number and geomorphometric characteristics of depressions. The range of detected depressions in the three areas were 442–491, 364–403, and 366–504, and the share of the depressions’ area conﬁrmed with all the approaches was 23%, 29%, and 47%, which resulted in different geomorphometric properties. Comparisons between the study areas were also inﬂuenced by the methods, which was conﬁrmed by the Mann–Whitney test. We concluded that preprocessing of high-resolution relief data and the detection methods in karst environments signiﬁcantly impact analyses and must be taken into account when interpreting geomorphometric results.


Introduction
Karst landscapes cover an extensive area of the Earth. According to Goldscheider et al. [1], about 15.2% of carbonate rock on land (outside glaciated areas) is at or near the surface and on which karst can develop. The best-known karst areas in the world include the Dinaric Karst, which stretches along the western coast of the Adriatic Sea from Italy to Albania. The karst of the southern part of Slovenia is part of the Dinaric Karst. Karst covers 9530 km 2 of Slovenia [2], which is just under half of the country. Karst boasts an abundance of landforms, most often dolines and other enclosed karst depressions, mostly located in forest-covered areas [3]. Karst depressions are often used as one of the diagnostic indicators to determine a karst area [4]. Different terms are used for terrain depressions (karst depressions, dolines, sinkholes, etc.); however, in view of the diversity of Dinaric karst depressions, we decided to use a general term-enclosed karst depressions (hereinafter karst depressions).
In the Dinaric Karst, there are differences between karst depressions-as well as natural dolines and collapse dolines, dolines anthropogenically altered for agricultural purposes are also common near settlements. The last is the so-called cultivated dolines, and they have anthropogenically levelled bottoms with a deeper soil layer, which make them more suitable for farming than the surrounding areas [5]. In general, the natural rocky and rugged karst landscape is often a limiting factor for agriculture and other economic activities, and it is often difficult to cross and map natural phenomena in the field. Due to the diversity of the bedrock and the aforementioned human influence, karst depressions have different shapes [6][7][8][9]. A review of the recent literature shows a number of karst depression detecting analyses, mostly focusing on the most common types of depressionround and ellipsoid [10][11][12][13][14][15][16][17][18][19]. In addition, round-shaped depressions have also been detected in non-karstic environments [20].
Most studies are based on high-resolution topography data, such as digital elevation models (DEMs), with a cell resolution of a few meters (e.g., [10,12]) to LiDAR elevation models, with a cell resolution of up to one meter (e.g., [13][14][15][16][17][18]). Studies have primarily focused on high-resolution DEMs, especially LiDAR-derived, which are becoming more and more detailed and fine-scaled, and enable geomorphometric analyses at different spatial scales, including the micro-scale for studying micro-features, such as solution features on rocks.
In our study, we focused on medium-sized karst features-enclosed karst depressionsas the most widespread and well-known karst surface landforms and, lately, also the most studied karst landforms [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In this decade, we have witnessed a switch toward the application of high-resolution LiDAR DEMs and the fast development of detection methods and geomorphometric analyses. In the future, technological development in remote sensing and laser scanning technology will enable even more detailed and fine-scaled DEMs. This quality has recently been achieved by using terrestrial laser scanners (TLS), which are mostly used for local studies covering small areas. There are some obstacles in forested areas due to dense and high vegetation cover. Removing the vegetation cover is also a methodological issue when using aerial laser scanning data because we consider that fine-scale laser DEM is influenced by vegetation noise. Such noise can be decreased by preprocessing. However, it can be assumed that with the progress in DEM data acquisition, preprocessing methods will become even more necessary.
Secondly, a number of computer methods can be utilized for detecting depressions; with the advancement of technology and improved resolution of spatial data, authors are increasingly choosing semi-automatic and automatic methods instead of manual digitization [19]. However, even new machine learning techniques need an input of learning samples in order to train their algorithm [7]. A certain determination of karst depressions is therefore still needed. Bauer [15] divided methods into three categories: (1) delineation of depressions based on the outermost contour line (e.g., [8]), (2) watershed-based doline delineation [19,21,22], and (3) filled-DEM based doline delineation [11][12][13][14]16]. One of the most recent automatic studies detected karst depressions for the entire surface area of Slovenia [9,23] by using the U-Net machine learning method.
The results of analyses of complex data can also differ significantly [24]. In general, in a set of detected landforms, the calculation of their geomorphometric or other characteristics, and their categorization, are influenced by the methods used (for example, when classifying landscapes [25,26]). From the geomorphometric perspective, the most problematic is the upper rim definition, which has already been discussed by several authors [27][28][29][30][31][32].
Input data from the same data source can also be preprocessed in many different ways. Telbisz et al. [16] compared different data sources (topographical map, contour line interpolation, laser scanning point cloud) and observed differences among them. Doctor et al. [10] also used preprocessing of DEM in order to remove some artificial influences (e.g., roads). When preparing data (e.g., removing noise), the method of DEM smoothing or small depression filling could also be used on karst terrains [12,16,18,31]. DEM smoothing on rocky karst includes some methodological issues since this type of landscape is characterized by underground water flow, creating concentric depressions on the surface. This is even more noticeable when analyzing forested areas, where the density of laser scanned points is lower, and data noise is present due to the classification of cloud points and vegetation removal, so preprocessing with smoothing or filling of the DEM before analysis is necessary (e.g., [16]). At the same time, smoothing obscures certain Remote Sens. 2022, 14,2416 3 of 20 depressions, so a comparison with the basic, non-smoothed, or non-filled DEM is also necessary. Consequently, this is a very complex process since the seemingly numerous flaws (noise) on the surface can actually be depressions of various sizes (even under 1 m). The problem is that accidental flaws (noise) are difficult to separate from actual depressions on the surface, especially in forests. A possible solution to improve the results is to smooth the acquired polygons of depressions (e.g., [8]) or sub-select detected polygons (depressions) based on pre-defined criteria of depth and area [10], which constitutes an additional intervention into the acquired data.
Since computer approaches require some kind of input (e.g., learning samples), a researcher cannot avoid influencing the geomorphometric results by selecting the methods. In other words, a computer can detect certain relief forms, but learning samples still have to be defined by a researcher [7]. The influence on methods can be related to the preparation of the data stage (e.g., preprocessing in order to avoid noise) and/or to the object detection stage (e.g., enclosed depression detection). In our study, we tried to combine different methods in these two stages and also include different study areas (with different geological characteristics).
We aimed to identify, characterize, and compare karst depressions of three densely forested study areas in the NW Dinaric Karst (Slovenia): Kras (with termophilous deciduous forest), Logaško-begunjski ravnik (Dinaric fir-beech forest) and Matarsko podolje (Mediterranean montane mixed forest). Accordingly, we compared various methods used for enclosed karst depression detection while also taking into account the multitude of ways laser scanning data of relief can be preprocessed. In other words, the selection of data preprocessing and detection methods might influence the results. Our goal was, therefore, to test and compare the results of different existing methods of karst depression detection, taking into account different preprocessing options of the laser scanning data. The main object of our study was karst depressions with diameters above 10 m and a depth of over 2 m. The contribution of our study is threefold: comparison and evaluation of differences in geomorphometric results based on preprocessing of high-resolution DEM, based on detection methods and based on study area properties.
The remainder of the paper is structured in four sections. In Section 2, the overall research design is presented, with the research area, data, and methods. Section 3 is dedicated to the results, which are presented in several subsections and provide information on the number of detected depressions, their detailed geomorphometric characteristics for each study area, and the overlapping of depressions. In Section 4, the discussion is focused on geomorphometric differences between study areas and specificities in the computerbased geomorphometric research of a karst environment. Section 5 concludes the paper with its main points.

Research Area
Slovenia boasts great landscape diversity since it is located at the contact of different European landscape units [33]. Karst relief covers 43% of Slovenia; different types of karst are consequently formed, with numerous surface and underground geomorphological features [3,34]. Geographically and geologically speaking, the Dinaric Karst is part of the Dinaric Mountains.
Considering the geological and geomorphological diversity of karst landscapes, consequent differences in characteristics of enclosed karst depressions (such as different types of depressions based on their genesis, size, shape, presence, or absence of anthropogenous intervention, etc.) are known [34] and were also expected in our DEM analysis. The research, therefore, encompassed three geologically (different types of carbonate rocks), geographically (location, elevation, and climate) and biogeographically (vegetation cover) different karst areas in the Slovene Dinaric Karst ( Figure 1): Kras, Logaško-begunjski ravnik, and Matarsko podolje. Four square kilometers of the forested area were randomly selected in each region. All three selected areas are geomorphologically classified as carbonate karst analysis. The research, therefore, encompassed three geologically (different types of carbonate rocks), geographically (location, elevation, and climate) and biogeographically (vegetation cover) different karst areas in the Slovene Dinaric Karst ( Figure 1): Kras, Logaško-begunjski ravnik, and Matarsko podolje. Four square kilometers of the forested area were randomly selected in each region. All three selected areas are geomorphologically classified as carbonate karst plains in the midst of larger, predominantly levelled karst areas marked by karst depressions. They are characterized by tectonic stability, thick bedrock layers, rock outcrops, and a deep groundwater table (over 100 m) [34]. Kras is a low Dinaric plateau, composed of Cretaceous limestones. It lies at an elevation of between 150 m to the NW and 450 m to the SE. In terms of diversity of karst depressions, solution dolines and collapse dolines dominate, with an average density of 60 dolines per km 2 , taking up approximately 12% of the surface area of Kras. They typically exhibit steep, rocky slopes with large and levelled bottoms, a consequence of anthropogenic intervention; they are also called cultivated dolines [5,9,34]. The Kras Plateau was traditionally cultivated, and before World War II (WWII), the region was characterized almost as karst-rocky desert, where the natural forest had been completely cut down centuries ago. After WWII, the agricultural land was gradually abandoned, and, since then, natural afforestation has been in progress. Most of the dolines have recently been covered by young (up to 70 years old) termophilous deciduous forests, dominated by Quercus pubescens [5].
Logaško-begunjski ravnik is a wide karst plain at an elevation of between 480 and 610 m. It has formed on Cretaceous limestones, Jurassic dolomites alternating with Jurassic limestones, and Triassic dolomites that follow in belts from west to east. The network of dolines, which is mostly solution and collapse by genesis, is very dense with Kras is a low Dinaric plateau, composed of Cretaceous limestones. It lies at an elevation of between 150 m to the NW and 450 m to the SE. In terms of diversity of karst depressions, solution dolines and collapse dolines dominate, with an average density of 60 dolines per km 2 , taking up approximately 12% of the surface area of Kras. They typically exhibit steep, rocky slopes with large and levelled bottoms, a consequence of anthropogenic intervention; they are also called cultivated dolines [5,9,34]. The Kras Plateau was traditionally cultivated, and before World War II (WWII), the region was characterized almost as karst-rocky desert, where the natural forest had been completely cut down centuries ago. After WWII, the agricultural land was gradually abandoned, and, since then, natural afforestation has been in progress. Most of the dolines have recently been covered by young (up to 70 years old) termophilous deciduous forests, dominated by Quercus pubescens [5].
Logaško-begunjski ravnik is a wide karst plain at an elevation of between 480 and 610 m. It has formed on Cretaceous limestones, Jurassic dolomites alternating with Jurassic limestones, and Triassic dolomites that follow in belts from west to east. The network of dolines, which is mostly solution and collapse by genesis, is very dense with up to 400 dolines per km 2 . The dolines constitute about 35% of the plain surface area [9]. The area is traditionally covered by mature and well-developed Dinaric fir-beech forests (Omphalodofagetum) [34]. The most widespread anthropogenic activity is traditional forestry.
Matarsko podolje is a karst plain, levelled out at an elevation between 500 and 680 m. It is formed on Cretaceous and Paleogene limestones [8]. The network of dolines, which are almost exclusively solution by genesis, is less dense, with only a few tens of dolines per km 2 . The dolines cover only 7% of the entire plain. As on Kras, the area of Matarsko podolje also features characteristic anthropogenically reshaped dolines with a levelled bottom. Compared to the dolines on Kras, the dolines on Matarsko podolje are shallower and have smaller dimensions [9,19]. Traditionally, the area was used as pastures for centuries, which were abandoned in the first half of the 20th century and overgrown by Mediterraneanmontane forests dominated by Fagus sylvatica and Corylus avellana in karst depressions and Pinus negra on shallower soils outside of depressions [5].
The selected study areas are strongly dominated by dense forests; the relief is very rugged, with rock outcropping mostly covered by thin soil or even without soil, which creates some areas that are hard to cross and access. Methods for automatic detection of karst depressions and visualizations based on LiDAR DEM are, therefore, the most promising and possible ways that allow us to explore karst depressions in these three areas, even better than with field survey.

Data
The entire country was scanned with laser scanning (LiDAR) technology in 2011, 2014, and 2015 in the framework of a national project financed by the Ministry of the Environment and Spatial Planning. Data layers for surface (relief with buildings and vegetation) and relief only were prepared with the same methodology [35] and are freely available and distributed by the Slovenian Environment Agency in the form of a point cloud and readyto-use relief raster tiles with a resolution of 1 m. The estimated density of ground points varies across the country and can be 0.5 per m 2 or more [35]. For our areas, the average densities were between 2.4 and 3.0 per m 2 . We used the available LiDAR point cloud already classified as ground (relief) by the provider and created a high-resolution DEM with a cell resolution of 1 × 1 m for each study area. The raster DEM was used for further preprocessing to detect enclosed karst depressions and to calculate their geomorphometric characteristics.

Methods
The methodological part of the study can be divided into the following sections ( Figure 2

Preprocessing the Input Laser Scanning DEMs
The basic raster data layer of the DEM was preprocessed in two different ways for each area in order to remove any potential noise in the data. The following tools were used: (1) Focal statistics-Smoothing using a five-cell circular radius filter; (2) Fill-Filling the depressions that were shallower than 1 m.
As suggested by Telbisz et al. [16], we used the Focal statistics tool for smoothing the DEM. This tool makes a calculation for each cell and ascribes to it an average value of all the surrounding cells in a certain radius. Following Telbisz et al. [16] work with LiDAR data, we also decided to use a five-cell radius for relief smoothing.
The second tool fills the depressions. This process detects cells that are encircled by cells with a higher elevation on all sides and therefore represents a theoretical sink. This kind of depression is filled up to the depth of the z-limit parameter. In our case, the depressions were filled further up to 1 m.
The two approaches are two of the different methods used to prepare a DEM prior to geomorphological analyses, such as detecting the connectivity of surface water flow [36] or (karst) depressions [16,18,32], and which are sensitive to any potential noise in the data [16].

Preprocessing the Input Laser Scanning DEMs
The basic raster data layer of the DEM was preprocessed in two different ways for each area in order to remove any potential noise in the data. The following tools were used: (1) Focal statistics-Smoothing using a five-cell circular radius filter; (2) Fill-Filling the depressions that were shallower than 1 m.
As suggested by Telbisz et al. [16], we used the Focal statistics tool for smoothing the DEM. This tool makes a calculation for each cell and ascribes to it an average value of all the surrounding cells in a certain radius. Following Telbisz et al. [16] work with LiDAR data, we also decided to use a five-cell radius for relief smoothing.
The second tool fills the depressions. This process detects cells that are encircled by cells with a higher elevation on all sides and therefore represents a theoretical sink. This kind of depression is filled up to the depth of the z-limit parameter. In our case, the depressions were filled further up to 1 m.
The two approaches are two of the different methods used to prepare a DEM prior to geomorphological analyses, such as detecting the connectivity of surface water flow [36] or (karst) depressions [16,18,32], and which are sensitive to any potential noise in the data [16].
For the sake of better result systematization, code FS5 marks the set of enclosed karst depressions that were detected based on the smoothed DEM and code F1 marks the set of enclosed karst depressions that were detected based on the filled DEM.

Enclosed Karst Depression Detection
The karst depressions were detected in two ways. The first method, Filled-DEM based doline delineation, detects karst depressions with theoretical filling and is carried out in several phases; it detects depressions of the first, second, and third or higher level [10,11]. The authors who used this method subsequently modified and adapted it or used only part of the method, since their studies were based mostly on detecting depressions of the first level, i.e., dolines (e.g., [14][15][16]). For the purposes of this article, we also used only part of the method and detected karst depressions of the first level. Using the input data, we used the hydrology toolset to calculate the flow direction, identify sinks or areas with no surface flow, determine watersheds, and theoretically filled them to the point at which the water would spill over the depression rims. Based on the difference in subtracting the theoretically filled depressions from the input DEM, we delineated the depressions. After preparing the final layer, we retained the robustness of the polygons in line with the resolution of the input DEM (1 m).
The contour line method is based on searching for completely closed contour lines. The method involves finding the highest closed contour line, which in theory represents the rim of a karst depression. A similar approach was applied by Bauer [15] and Verbovšek and Gabor [8]. Our research pre-set the criterion that depressions had a round shape. After a trial-and-error examination, we set the circularity index range at 0.75-1.25. The value for a square is 1.27 and the value for a circle is 1.00; in practice, the most realistically detected depressions had values up to 1.25, which also led us to set the bottom limit at 0.75. The circularity index is provided in Table 1 (code CIRCB). This helped us to exclude larger depressions of irregular shapes. To avoid detecting smaller depressions, we set an additional condition that the polygon must be at least 10 m 2 . As with the method of filled-DEM depression delineation, we also prepared the final layer with the robustness of the polygons according to the resolution of the input DEM (1 m). Both methods returned a set of depressions. We then eliminated from the sets any depressions that were too small in relation to the diameter or were too shallow, pointing to the likelihood that they represented noise in the data and were not truly karst depressions: The criteria were set based on the literature, which determined that enclosed karst depressions and dolines are those with a depth greater than or equal to 2 m and a diameter greater than or equal to 10 m [10,34]; Kobal et al. [14], for example, selected polygons of karst depressions to be analyzed based on the same criterion. There was no top limit in terms of the size (area) or depth of depressions in this study.
The final selection of depressions included only those polygons that were completely within the area of the 2 km × 2 km studied area.
For further analysis and a better systematization of the results, we distinguished the results (the set of karst depressions) that were detected with the filled-DEM method by designating them with the code S1 and the results (set of karst depressions) that were detected using the contour line method with the code R1.

Calculating the Geomorphometric Characteristics of Karst Depressions
We ascribed to the detected karst depressions the selected geomorphometric variables based on data regarding the elevation, slope, and other geomorphometric characteristics (e.g., ruggedness, shape), such as area, perimeter, diameter, depth, elongation, length, width, orientation, the lowest, highest, and average elevation, slope, and others. The geomorphological variables listed in Table 1 were used for further analyses. For these metrics, mean, median, and coefficient of variation were calculated in order to provide basic statistical parameters for each study area, with a specific combination of methods (settings). This allowed us to geomorphometrically characterize each study area. The values of the geomorphometric variables were calculated using geoinformation tools in ArcGIS Pro 2.8.2.

Overlapping the Detected Karst Depressions
An overlapping analysis was carried out on the layers of the detected karst depressions for each area using the Count overlapping features tool. We added all the cells that were recognized in each individual approach and calculated the surface area in m 2 . The final result of the overlapping analysis is the surface area of the karst depressions as detected by all six detection approaches.
Additionally, karst depressions detected with different approaches were overlaid and compared with reference karst depressions that had been digitized based on topographic maps (1:5000). These maps are the most detailed topographic maps available for the study areas.

Analyzing the Differences between the Areas
The set of detected karst depressions and their calculated geomorphometric characteristics (see Section 2.3.3) enabled us to compare the geomorphological characteristics between the study areas and examine whether there is an influence of different approaches for detecting depressions: that is, (1) the differences in the results of different detection methods with the same input data (DEM), and (2) the differences in the results of different input data (DEM) using the same detection method.
The areas were compared based on the set of karst depressions with their calculated geomorphometric characteristics. In order to determine whether statistically significant differences exist between two study areas according to a certain variable, the Mann-Whitney test was applied. We also calculated the basic statistical characteristics for each individual area, such as the average, median, and variation coefficient, to compare how the data varies.

The Number of Detected Karst Depressions within the Study Areas
Enclosed karst depressions with diameters above 10 m and depths of over 2 m were the main object of the study. We calculated their geomorphometric characteristics and used them further to compare three different areas and to recognize whether data preprocessing and two different depression detection methods significantly affect the acquired geomor-Remote Sens. 2022, 14, 2416 9 of 20 phometric comparisons. We prepared several layer sets of detected karst depressions in three selected karst areas. The layers representing the results for three different areas were compared statistically to determine whether we could reach more or less the same conclusions regarding the features of the depressions regardless of the approach (data preprocessing method and detecting method).
According to Table 2 (see also Supplementary Figures S1-S18), the greatest number of detected depressions by both methods was recorded in the area of Logaško-begunjski ravnik. Method S1 detected between 475 and 488 karst depressions, and method R1 detected between 442 and 491 depressions. This area, therefore, also has the greatest depression density. The lowest number of depressions by both methods was detected in Kras, where the depression density was between 91 and 100.8 depressions per km 2 . The situation in Matarsko podolje was exactly the opposite, since approaches S1, FS5, and R1, FS5 detected the same number of depressions, 504, while approaches S1, F1, and S1, OS detected significantly fewer depressions (about 130) than approaches R1, F1, and R1, OS ( Table 2). Depression detecting approach S1, F1 S1, FS5 S1, OS S1, F1 S1, FS5 S1, OS S1, F1 S1, FS5 S1, OS In the area of Logaško-begunjski ravnik and Kras, approaches S1, F1, and R1, F1, and S1, OS and R1, OS were the least divergent in the number of detected depressions. The differences in the detected depressions of these approaches were no more than 13. A greater deviation was detected between approaches S1, FS5, and R1, FS5, with which the difference in the detected depressions in the area of Logaško-begunjski ravnik was 33, and 37 in the area of Kras.

Geomorphometric Characteristics of Karst Depressions within Individual Areas
The basic geomorphometric characteristics (average, median, variation coefficient) of karst depressions for the test area in Logaško-begunjski ravnik (

Logaško-Begunjski Ravnik
Based on the results, the smoothing of the input DEM causes the formation of milder and larger depressions, while the concave landform is mostly preserved. As can be seen from Table S1, some statistical characteristics are very diverse. For example, the mean surface area (AREA) in method S1 with the basic DEM (OS) is 598 m 2 , and 794 m 2 in the R1 method and the smoothed DEM (FS5); the same case is with mean volume (VOL) being the greatest for smoothed DEM (FS5) (1738 m 3 ) and exceeded the values for non-processed and filled DEM by almost 400 m 3 .
Far fewer differences can be observed, for example, in the mean elongation (ELONG), which ranged between 1.22 and 1.27 in all the results, which is much more comparable, so the preprocessing had a smaller impact on this variable. A difference in the utilized methods R1 and S1 is noticeable, for example, in the mean circularity (CIRCB). Since the R1 method is based on determining round contour lines, the mean circularity values (CIRCB) are expectedly lower in the R1 than in the S1 method, which points to rounder depressions.
The mean depth also has relatively equal values everywhere (just over 4 m). There are very small differences when comparing the results with the basic and the filled DEM.
Here it must be emphasized that karst depressions with a depth of under 2 m had already been excluded at the beginning, but we still anticipated that greater differences would be detected.
The mean shares of leveled depression bottoms (SH_3), which can potentially indicate differences between natural and cultivated dolines, have similar values-about 1%. Due to traditional forest land cover, cultivated dolines are not characteristic of this study area.

Kras
Similar to the case of Logaško-begunjski ravnik, some differences between the geomorphometric characteristics were also detected for the area of Kras (Table S2). The mean surface area (AREA) was between 671 m 2 (S1, FS5) and 1325 m 2 (R1, FS5), while the mean volume (VOL) ranged from 1648.63 m 3 for S1, OS and up to 5101.63 m 3 for R1, FS5. The coefficients of variation for surface area and volume are the highest (up to 1.51 and 3.29), showing a much higher diversity of depressions than in the other two study areas. In general, depressions on the Kras Plateau are much bigger, and almost all were anthropogenically modified throughout the centuries and used for cultivation. This can also be marked in shares of leveled depression bottoms (SH_3), which are, on average, much higher (up to 5%) than in the other two study areas (about 1%).
Similarly, the mean depth values (DEPTH) ranged from 5.23 (S1, OS) up to 6.26 m (R1, FS5). There are small differences when comparing the results with the basic and the filled DEM. However, the mean elongation (ELONG) is once again more similar between areas, with values being between 1.20 and 1.33.
DEMs apparently do not have a lot of noise, and there are also not so many smaller real depressions.

Matarsko Podolje
The mean surface area (AREA) in Matarsko podlje is between 727 m 2 (S1, F1) and 1162 m 2 (R1, FS5), while the mean volume (VOL) ranges from 1537 m 3 for S1, F1 up to 2667 m 3 for R1, FS5. While the trend for surface area and volume is pretty much similar to Kras Plateau, the coefficients of variance for surface area and volume are lower, showing lower diversity of depressions than on Kras. In general, depressions in this part of Matarsko podolje had not been anthropogenically modified since they were used as pastures. The mean shares of leveled depression bottoms (SH_3) are about 1-3%.
The depths on Matarsko podolje are very similar to those on Logaško-begunjski ravnik, i.e., between 4.39 and 4.54 m, while the mean surface areas are more diverse (between 726 and 1162 m 2 ). As expected, the influence of DEM filling with 1 m also had less influence on the end result (Table S3).

Overlapping of Karst Depressions within Individual Areas
Some differences in karst depressions within individual areas can also be visualized by overlapping the resulting layers. The overlapping exposed firstly differences in the number of detected landforms and, secondly, in their surfaces. Analysis of karst depression overlapping (Table 3, Figures 3-5) determined how much surface was covered by the depressions that were detected with each detection approach (each combination of preprocessing and detection methods). The numbers (1-6) denote how many times the area was detected as a depression. In two of the selected areas, Kras and Matarsko podolje, the most surfaces were detected with four out of six approaches, while the least surfaces were detected with five out of six approaches. The majority of the surfaces on Logaško-begunjski ravnik were detected in all six ways and the least surfaces in five out of six approaches. The shares of the depressions' area confirmed with all the approaches were 23%, 29%, and 47%.   The numbers (1-6) denote how many times the area was detected as a depression (e.g., number 5 denotes that an area was marked as a depression with five approaches). Figure 3. Overlapping of detected karst depressions for the area of Logaško-begunjski ravnik. The numbers (1-6) denote how many times the area was detected as a depression (e.g., number 5 denotes that an area was marked as a depression with five approaches).

Figure 4.
Overlapping of detected karst depressions for the area of Kras. The numbers (1-6) denote how many times the area was detected as a depression (e.g., number 5 denotes that an area was marked as a depression with five approaches). (1-6) denote how many times the area was detected as a depression (e.g., number 5 denotes that an area was marked as a depression with five approaches).
When comparing detected karst depressions with the reference layer of karst depression based on topographic data (1:5,000), it can be observed that differences in the recall are diverse-between 59% and 95% (Table 4). In general, the detection method R1 provided a higher rate of matching with topographically defined karst depressions. There are only two cases in which S1 had a higher match-when preprocessing DEM with FS5 for Logaško-begunjski ravnik and Kras. We noticed that, in general, preprocessing DEM with FS5 resulted in better matching when using detection method S1. In the case of Matarsko podolje and Kras, the detection method R1 provided the best results when combined with preprocessing method FS5; in the case of Logaško-begunjski ravnik, this (1-6) denote how many times the area was detected as a depression (e.g., number 5 denotes that an area was marked as a depression with five approaches).
When comparing detected karst depressions with the reference layer of karst depression based on topographic data (1:5,000), it can be observed that differences in the recall are diverse-between 59% and 95% (Table 4). In general, the detection method R1 pro-vided a higher rate of matching with topographically defined karst depressions. There are only two cases in which S1 had a higher match-when preprocessing DEM with FS5 for Logaško-begunjski ravnik and Kras. We noticed that, in general, preprocessing DEM with FS5 resulted in better matching when using detection method S1. In the case of Matarsko podolje and Kras, the detection method R1 provided the best results when combined with preprocessing method FS5; in the case of Logaško-begunjski ravnik, this combination is the worst. Depression detecting approach S1, F1 S1, FS5 S1, OS S1, F1 S1, FS5 S1, OS S1, F1 S1, FS5 S1, OS The precisions are above 75% for Kras and Matarsko podolje, but mostly less than 30% for Logaško-begunjski ravnik. The reason might be in a low-quality production of topographic maps for that area, which is the most densely vegetated.
The differences in matching confirm that the selection of preprocessing and detection methods are crucial and have a great influence on the results. After comparing the results with different preprocessing DEMs, we noted that in most cases (except for one), recalls and precisions with FS5 preprocessing are better than those that were made based on the original DEM.

Discussion on Geomorphometric Differences between the Study Areas
When comparing the differences between areas (e.g., answering the question: 'Do karst depressions in Logaško-begunjski ravnik differ from depressions in Kras according to the geomorphometric characteristics?'), we observed, based on the Mann-Whitney test (Table 5), that the results were partially coordinated. Table 5. Analysis of the differences between areas according to the selected geomorphometric variables. Statistical significance of the Mann-Whitney test: the grey cells indicate that the areas differ significantly from one another (with p < 0.05). After the Mann-Whitney test, analysis of the differences between areas according to individual geomorphometric variables showed that the karst depressions of the study areas differ significantly among one another in most cases, i.e., according to most geomorphometric characteristics. There are some noticeable differences between the analyses of different sets of karst depressions. Some results can therefore be designated as different for some pairs of comparisons if a different method of karst depression detection and/or a differently preprocessed input data (DEM) ( Table 4) is used. Some discrepancies should be highlighted. The most noticeable inconsistency was observed for the elongation (length of the major axis/length of the minor axis) (ELONG) variable, in which eleven comparisons determined that there were no differences, but seven comparisons confirmed the said differences. A similar scenario occurred with the circularity (CIRCB) variable, with which a similar number of comparisons (eight vs. ten) were confirmed or rejected. In the other variables, the difference between the areas was confirmed, but there is at least one example in each that opposes that. The exception is the variable ratio between the horizontal distance to the centroid (2D) from the lowest point and the length of the major axis (RAT_DD), in which all the areas differ regardless of the approaches.

The Issues Related to Forest Cover and Vegetation Filtering
Many studies were referenced in this manuscript in relation to the detection and characterisation of enclosed karst depressions, mostly dolines, but only rarely have methods been tested in the presence of forest-tree canopies. As Kobal et al. [14] already highlighted, tree canopies can decrease the homogeneity of the spatial distribution obtained from ground hits. This effect depends on the distribution of the canopy density, especially where high-resolution surveys are conducted because higher emitter frequencies correspond with lower laser pulse energies and lower penetration capabilities [14].
Based on Tables S1-S3, it can be seen that the highest values of characteristics, such as surface area and volume of karst depressions for all three study areas, were calculated when detecting karst depressions based on smoothed (FS5) DEM. When detection methods were done based on filled (F1) or original (OS) DEM, the values of these characteristics were lower. It can therefore be confirmed that different preprocessing of DEMs has a noticeable impact on the final results of the general characteristics of depressions. Characteristics, such as volume and area, also had higher average values in cases in which karst depressions were detected based on the contour line method (R1). These results should be further considered in future geomorphological analyses of karst depressions in similar forested karst areas with complex geomorphology, a high density of karst depressions, and other karst phenomena.
When preparing data (e.g., removing noise), the methods of DEM smoothing or small depression filling can be used [16,18,32], as was done for karst depression (mostly dolines) detection in the area of the Slovak Karst National Park (Slovakia) by Hofierka et al. [18], in the Aggtelek Karst (Hungary) by Telbisz et al. [16], and in the Pipestem watershed (North Dakota, USA) by Wu et al. [32]. The preprocessing of the high-resolution LiDAR DEM is needed to guarantee that it is "hydrologically correct" for successive analyses [14]. In our case studies, when comparing areas, some inconsistencies between comparisons of geomorphometric characteristics based on original DEM and preprocessed DEM were observed. More inconsistencies occurred if DEM smoothing was used rather than small depression filling in the data preprocessing ( Table 4). The influence of smoothing, which was noticeable in our study, was also analyzed also by Erdbrügger et al. [39] for different relief types. They discovered an important impact of DEM smoothing on the results of flow direction analysis. Since our methods also used some steps of hydrological analysis (Section 2.3.1), an influence of smoothing was expected and also confirmed. It is also important to highlight that different values of filling or selecting different radii when smoothing the DEM, can impact the final results [16]. Telbisz et al. [16] compared different ways of filling in the input data of a laser scanning and elevation model based on topographic maps. They determined that the height of the fill affects the final number of detected depressions. They also concluded that the characteristics related to size (length, width, perimeter length) are similar between different sets but less so in terms of the circularity. Our research can confirm that comparisons of circularity (CIRCB) between study areas were notably influenced by the preprocessing and detection method.
The effect of the set of initial settings on the final result has, for example, already been confirmed in an analysis of the subjectivity of a DEM-based data segmentation method [40]. Lidberg et al. [36] analyzed the effect of different resolutions and data preprocessing to determine surface water flow modelling. They determined that a higher data resolution was better for modelling, whereas the breaching method turned out to be the most suitable data preprocessing method and depression filling the least. Gostinčar and Ciglič [41] compared the detection of relief features with different resolutions and determined that, according to the Peucker and Douglas method [42], the highest resolution is not necessarily the best option, since noise was present.

The Issues Related to Human Impacts on Terrain "Smoothness" Due to Agricultural Land-Use
When studying and geomorphometrically analysing karst areas, especially karst depressions, anthropogenic impacts need to be considered when selecting the preprocessing of data, as well as in the selection of the detection method. Hard carbonate rocks (limestone and dolostone) in which the karst depressions are shaped are extremely resistant to weath-ering, so anthropogenic remains (e.g., the formation of patterns of stone walls, vernicular architecture, etc.) can be preserved for millennia (e.g., the UNESCO world-heritage site Starogradsko polje on Hvar Island). Different agricultural practices, such as rock outcropping "clearing", are known in different karst landscapes in Europe, e.g., in the south of France [43] and also on the Kras Plateau in Slovenia [5,34], where the removal of rock outcroppings was carried out inside plots, and, years later, the effects of this activity are visible on a LiDAR image showing the "smoothness" of these plots [5].
For this study, we intentionally selected three karst landscapes with different carbonate rocks, characterized by karst depressions and covered by dense forest vegetation. In the study, the influences of preprocessing of DEM were less frequent when comparing Matarsko podolje and Logaško-begunjski ravnik and more noticeable when comparing these two areas with Kras (Table 4), where the surface has been traditionally more affected in the past due to agriculture (e.g., cultivated dolines, clearing of stones from grasslands). This difference can be confirmed by ruggedness measures for the whole study area (Table 6). In addition, the higher "smoothness" of the Kras area can be recognized through a measurably higher share of flat depression bottoms (share of cells with a slope ≤3 • ), characteristic of cultivated dolines (Tables S1-S3).

Remote Detection of Karst Depresions as the Best Alternative
The study found that the application of different preprocessing data and different detection methods can significantly affect the results of the comparison of depression characteristics in different areas. The influence of the method on the outcome of the analyses was already illustrated by Bauer [15], who used three different detection methods to extrapolate different sets of detected depressions (dolines), which also impeded the interpretation of the results. As with Bauer, in this research as well, all six sets of karst depressions (per each study area) showed visible and measurable (quantitative) differences in size (e.g., surface area) as well as shape (e.g., elongation) and volume (Figures 3-5). The difference in the utilized methods R1 and S1 is noticeable, for example, in the circularity (CIRCB). Since the R1 approach is based on determining round contour lines, the circularity values (CIRCB) were expectedly lower in R1 than in the S1 approach, which points to rounder depressions. The influence of the methodology selection was also observable in the difference in the number of detected depressions for each area. It should be mentioned that methods for DEM preprocessing and depression detecting, require certain settings. For example, smoothing requires the radius, shape, and type of filtering (e.g., mean, modus, etc.), and the filling requires a depth limit setting. Similarly, the detection method requires the definition of roundness, i.e., the selection of an appropriate measurement that defines circularity (see [31] for a list of geomorphometric parameters). These options provide numerous sub-settings-approaches-which cannot all be tested. Each geomorphometric approach can therefore involve some kind of bias, and the solution might be the use of different methods for a multiple-analysis approach.
Despite the high accuracy of LiDAR DEM, additional remote sensing data layers (e.g., satellite images, aerial images) can be included to increase the accuracy and success rate of detecting depressions. In order to detect denuded caves, which are more elongated depressions than round dolines, Grlj and Grigillo [22] combined the DEM and satellite data. Hoai et al. [44] used thermal depression imaging to detect dolines, which were further segmented and classified using machine learning. Unlike methods that focus on detecting individual shapes, so-called comprehensive classification methods focus on classifying the landscape as a whole. Analysts try to determine certain spatial patterns for larger areas of a certain landscape. For example, the visibility network method for a comprehensive classification of different landscape types (where the focus is not on an individual shape but a comprehensive feature of a selected landscape) has also been tested in synthetic DEMs [45]. The methods of comprehensive relief classification were developed by Drăguţ and Eisank [46], who classified the surface of the entire planet, and by Wieczorek and Migon [47], who classified SW Poland. The karst depressions that were detected using geoinformation technologies were either further classified by various authors into different types, or their various statistical characteristics were calculated (e.g., [8,9,16]). However, in addition to all the referenced calculation methods, different visualizations of high-resolution DEMs for remote areas are developing [48], which will give additional perspective to landform features in a study on high-resolution digital elevation models.
In any case, methods for karst depressions detection are suitable for identification, general spatial analysis, and calculation of their morphometric characteristics. Such methods allow the easy and inexpensive application to other areas, including over large scales, and thus provides several advantages over manual mapping, especially when study areas are difficult to access and forested, as well as when time investment and field costs are considered [49]. A combination of remote sensing with extensive field surveys will always provide the best and most reliable results. Unfortunately, there will always exist some remote and difficult of access regions on the Earth's surface and/or even on other planets (e.g., Mars) where, so far, remote sensing is the only existing and the safest method for the geomorphometric study.
Nevertheless, methods for karst depression detection and the use of remote sensing data have considerable potential for application in large-scale analyses in many fields in addition to geomorphology and karstology, including geology, botany, zoology, tourism, spatial planning, and nature conservation [5,[49][50][51][52].

Conclusions
The study highlights the influence of data preprocessing and various methods for karst depression detection. We confirmed that both have a significant effect on the results, in some cases to the extent that the results are statistically different when comparing different study areas. We observed that different combinations of methods influenced the number and geomorphometric characteristics of depressions. The range of detected depressions in the three study areas was 442-491, 364-403, and 366-504, and the shares of depressions' areas confirmed with all the approaches were 23%, 29%, and 47%. This also resulted in the calculation of different geomorphometric properties. We noticed that the preprocessing and detection methods used had an impact on most geomorphometric variables. Higher inconsistencies occurred in the case of smoothing DEM than if smaller depressions were filled in the data preprocessing. In light of the results, we determined that it is crucial for each study to clearly present the entire process, including data preprocessing, since this guarantees the most comprehensive transparency.
Nevertheless, the main questions remain partially open: Can we reach some general conclusions considering which method should be preferably used for individual karst type area? Which methodological approach should be used for different successive analyses? In theory, it would be easier to select the most suitable methodological process by comparing results with detailed geomorphometric measurements carried out in the field. This is the only way to determine which method comes closest to the real state in nature. However, it has to be highlighted that many areas in the Dinaric Mountains and other densely forested karst landscapes are not easily accessible, and remote sensing might be the only option to detect certain relief features. It can be expected that future studies will use even more precise and high-resolution DEMs, which will enable fine-scale geomorphometric analyses. With fast progress in DEM data acquisition, preprocessing methods will become even more necessary.
In addition to its methods, the results of this study also have great potential to be used in further environmental studies. Our future study will focus on the connectivity between karst depressions on the surface and underground karst caves. Specifically, the resulting layers of detected depressions will be compared (overlaid) with the layer of underground cave systems. In such a way, the locations of karst depressions above the cave systems will be determined and further studied in greater detail from (vertical) hydrological connectivity. Such potential connections can only be found by (semi) automatic determination of depressions for wide areas.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14102416/s1, Figure S1: Detected karst depressions in Kras area based on R1 method and DTM F1.; Figure S2: Detected karst depressions in Kras area based on R1 method and DTM FS5.; Figure S3: Detected karst depressions in Kras area based on R1 method and DTM OS.; Figure S4: Detected karst depressions in Kras area based on S1 method and DTM F1.; Figure S5: Detected karst depressions in Kras area based on S1 method and DTM FS5.; Figure  S6: Detected karst depressions in Kras area based on S1 method and DTM OS.; Figure S7: Detected karst depressions in Matarsko podolje area based on R1 method and DTM F1.; Figure S8: Detected karst depressions in Matarsko podolje area based on R1 method and DTM FS5.; Figure S9: Detected karst depressions in Matarsko podolje area based on R1 method and DTM OS.; Figure S10: Detected karst depressions in Matarsko podolje area based on S1 method and DTM F1.; Figure S11: Detected karst depressions in Matarsko podolje area based on S1 method and DTM FS5.; Figure S12: Detected karst depressions in Matarsko podolje area based on S1 method and DTM OS.; Figure S13: Detected karst depressions in Logaško-begunjski ravnik area based on R1 method and DTM F1.; Figure S14: Detected karst depressions in Logaško-begunjski ravnik area based on R1 method and DTM FS5.; Figure S15: Detected karst depressions in Logaško-begunjski ravnik area based on R1 method and DTM OS.; Figure S16: Detected karst depressions in Logaško-begunjski ravnik area based on S1 method and DTM F1.; Figure S17: Detected karst depressions in Logaško-begunjski ravnik area based on S1 method and DTM FS5.; Figure S18: Detected karst depressions in Logaško-begunjski ravnik area based on S1 method and DTM OS.; Table S1: The geomorphometric characteristics of karst depressions for the test area in Logaško-begunjski ravnik (according to two detecting methods and input data preprocessing). CV-coefficient of variation.; Table S2: The geomorphometric characteristics of karst depressions for the test area in Kras (according to two detecting methods and input data preprocessing). CV-coefficient of variation.; Table S3: The geomorphometric characteristics of karst depressions for the test area in Matarsko podolje (according to two detecting methods and input data reprocessing). CV-coefficient of variation.
Author Contributions: All authors discussed the article structure and main goal. The methodological description and geoinformatic analyses were conducted by Š.Č. and R.C. The introduction, results, discussion, and conclusion were written by all authors. M.B.V. undertook the final proofreading of the article. All authors have read and agreed to the published version of the manuscript.

Funding:
The study was financed by the Slovenian Research Agency in project J6-2592 (Connectivity concept in karst: the doline-cave coupling in the context of human impacts) and program P6-0101 (Geography of Slovenia).

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.