MAT: GIS-Based Morphometry Assessment Tools for Concave Landforms

Abstract

The rapid development of remote sensing technology for obtaining high-resolution digital elevation models (DEMs) in recent years has made them more and more widely available and has allowed them to be used for morphometric assessment of concave landforms, such as valleys, gullies, glacial cirques, sinkholes, craters, and others. The aim of this study was to develop a geographic information systems (GIS) toolbox for the automatic extraction of 26 morphometric characteristics, which include the geometry, hypsometry, and volume of concave landforms. The Morphometry Assessment Tools (MAT) toolbox in the ArcGIS software was developed. The required input data are a digital elevation model and the form boundary as a vector layer. The method was successfully tested on an example of 21 erosion-denudation valleys located in the young glacial area of northwest Poland. Calculations were based on elevation data collected in the field and LiDAR data. The results obtained with the tool showed differences in the assessment of the volume parameter at the average level of 12%, when comparing the field data and LiDAR data. The algorithm can also be applied to other types of concave forms, as well as being based on other DEM data sources, which makes it a universal tool for morphometric evaluation.

1. Introduction

Morphometric analysis is based on measurements and a mathematical approach to the shape and dimensions of landforms [1,2]. The possibility of using geographic information systems (GIS) and digital elevation models (DEMs) in calculations contributed to the development of geomorphometry as the science of quantitative land-surface analysis [3,4,5,6]. Morphometric analysis enables the assessment of the spatial diversity of landforms, as well as the interpretation of the relationships between them. The detailed analysis of morphometric parameters allows obtaining information about the genesis and evolution of forms and enables them to be classified into various types [7,8]. The analysis of morphometric characteristics is also helpful in the reconstruction of the palaeogeographic conditions of the areas where these forms occur [9,10,11].

One type of landform is concave forms, which include valleys, gullies, glacial cirques, sinkholes, kettle holes, pockmarks, craters, and other types of local depressions [12]. Generally, one can distinguish dry forms in which there are no permanent streams or watercourses and catchments, in which a river network occurs. In order to characterise them, their shape or geometric parameters are determined, such as length, width, area, and hypsometric parameters, such as relief [13,14]. To assess the characteristics of river catchments, linear parameters based on the drainage network are additionally used [15,16,17].

The measurement techniques used to calculate the morphometric parameters of landforms have evolved over the years. Techniques used in geomorphological mappings conducted in the field, such as profiles [18], total stations (TS) [19], differential GPS [11,20,21], or nowadays terrestrial laser scanning (TLS) [22,23,24,25,26] and Structure from Motion (SfM) photogrammetry [23,27,28,29,30], should be distinguished. Data for the analysis of the large areas can be obtained from sources such as topographic maps [1], aerial and satellite images [28,31,32,33,34], or airborne laser scanning (ALS) [22,26,30,35,36,37,38]. Researchers have compared measurement methods in terms of their accuracy, their time and cost requirements, and their efficiency [22,26,30,39,40]. The advantage of modern remote sensing methods is the acquisition of high-resolution data, such as point clouds and orthophotos of observed objects, which enable the generation of very high-quality digital elevation models (DEMs) for morphometric analysis [41].

The morphometric analysis performed in GIS allows for the calculation of many parameters. However, there is a need to use various geoinformation tools, often in several programs. In order to automate the computation of the morphometric characteristics of concave forms, the authors developed their algorithms in the form of models, scripts, plug-ins, or extensions to facilitate work in GIS software [38,42,43,44,45,46,47,48,49,50,51]. Based on the analysis of the functionalities of the mentioned tools, the following geometric parameters were considered as important for the morphometric evaluation of concave forms: area, perimeter, length, width, orientation, elongation ratio, circularity ratio, compactness coefficient, form factor, and the following hypsometric parameters: maximum, mean, and minimum height, elevation range, slope, aspect, Topographic Wetness Index (TWI), Terrain Ruggedness Index (TRI), hypsometric integral (HI), volume. The study assumed the use of ArcGIS software because of its broad capabilities for DEM data analysis. We examined the availability of the mentioned characteristics in ArcGIS tools provided as open source (

Table 1. Morphometric parameters implemented in the ArcGIS open-source tools for the concave forms analysis.

Tool	Morphometric Parameters
	Area	Perimeter	Length	Width	Orientation	Elongation Ratio	Circularity Ratio	Compactness Coefficient	Form Factor	Elevation	Slope	Aspect	TRI	TWI	HI	Volume
CalHypso [42]															×
ACME [43]	×	×	×	×	×		×			×	×	×			×
BaDAM [44]	×	×	×	×		×	×	×	×	×					×
rillStats [46]																×
SWPT [47]	×	×	×	×		×	×	×	×					×
Hypsometric Integral Toolbox [51]															×

).

The CalHypso script [42] and the Hypsometric Integral Toolbox [51] were designed to calculate hypsometric curves and their statistical moments. The CalHypso tool extracts also the hypsometric curve as a graph. The ACME toolbox allows determining geometric and hypsometric parameters of glacial cirques, including length, width, planar and 3D area, perimeter, circularity, elevation, slope, aspect, and hypsometric integral [43]. The BaDAM and SWPT toolboxes for the assessment of river catchments are based on indicators such as elongation ratio, circularity ratio, compactness coefficient, form factor, shape factor, drainage density, and stream frequency, among others [44,47]. SWPT also includes spatial indexes, such as Topographic Wetness Index (TWI). It should be noted that neither tool determines the value of the Terrain Ruggedness Index (TRI).

In the characteristics of concave forms, most of the parameters mentioned above can be clearly defined, but there are still parameters that are problematic due to the ambiguous definition of the boundaries of the analysis. Such parameters include the volume of concave forms. Only one of the tools presented in Table 1 allows the calculation of the volume parameter. The tool was developed as a toolbox, which, in addition to the volume, also determines the form cross-sections [46]. The standard approach used in calculating the volume changes of forms is the analysis of digital elevation models created for at least two time periods, called the DEMs of Difference (DoD) approach [52]. This method has been used in research on, among other things, gully erosion [26,30,40], rill erosion [23,53,54], mass movements on slopes [55,56,57], and geomorphological processes (erosion and accumulation) in river valleys [58]. However, if there is no digital elevation model from the period before a given form of concave originated, the calculations are more complicated. In this case, the original surface area should be reconstructed by developing a 3D cover plate [59]. D’Oleire-Oltmanns et al. [59] used for this purpose the digitisation of the boundary of the depression, the creation of a 3D polygon, and then the conversion of the polygon into a raster, thus obtaining an overlay elevation model. A similar methodology was developed by Báčová et al. [46] using the vertex of the form polygon boundary to create an overlay TIN model and then converting it to a raster.

Based on the analysis of the functionality of the toolboxes analysed above, it should be stated that there is still no tool that would allow for a comprehensive morphometric assessment including the geometry, hypsometry, and volume of any dry concave form in the ArcGIS software. The aim of this study was to develop an open-source geographic information systems (GIS) toolbox for the automatic extraction of 26 morphometric characteristics. To enable a complex evaluation of the concave forms, it was assumed that in addition to the parameter values, the tool would also calculate descriptive statistics, hypsometric curves as graphs, and raster layers of the spatial parameters (slope, aspect, TRI, TWI). The Morphometry Assessment Tools (MAT) were developed using the ArcGIS ModelBuilder visual programming language and the functions available in the ArcGIS software. The tools were implemented as a toolbox for use in the ArcGIS 10.7 version. The advantage of the tools is their speed of calculation, their ability to be used for many objects at the same time, the minimal set of input data required, and their user-friendly interface.

2. Study Area

The toolbox was tested on the example of small dry valleys formed within the morainic upland, outwash plains, and slopes of subglacial channels in the young glacial area located in the West Pomeranian Lakeland in northwest Poland (Figure 1). The relief of the studied area is the result of the functioning of the post-glacial denudation system. This system is based mainly on the drainage routes and directions formed on the foreland of the ice sheet, which withdrew from the range line of the Pomeranian phase of the Weichselian glaciation 16–17 ka, [60]. This area shows a well-developed relief, accompanied by numerous lakes of the Drawsko and Szczecinek Lake Districts. In the foreground of the marginal forms of the Pomeranian phase, there was a fluvioglacial drainage route which formed the vast outwash plains of the Drawa, Piława, and Gwda basins, and in the background of the marginal forms there is a trough-like, low-lying area of the Białogard Plain, drained by the upper Parsęta and its tributaries [61]. Outwash sands cover the moraine hills that were formed during the retreat of the Weichselian ice sheet before the Pomeranian phase [62]. As a result of the slope processes in the form of surface mass movements, slope wash, and gully erosion, strongly inclined surfaces were remodelled, resulting in dry or periodically dry valley forms cutting through them [63]. The morphometric evaluation was performed for 21 dry valleys located in the test areas, of which three were morainic upland within the catchment of the Dębnica River (Parsęta tributary), and one was located within the subglacial channel of Lake Rymierowo (Figure 1).

The dominant type of relief in the northern and central part of the study area is the undulating moraine upland. The area of Buślarki is located on the eastern slope of the Buślarskie Góry. From the top of the upland at 120 m a.s.l., the terrain descends to 60 m a.s.l. towards the Dębnica valley. To the east of the Dębnica valley, there are complexes of glacial crevasse forms of the kame plateau of Skowrończe Góry [64]. The altitude of the Piaski Pomorskie test area ranges from 70 to 145 m a.s.l. The monotonous relief of the moraine upland is diversified by the forms of crevasse accumulation and extensive kames [63]. The Koprzywno test area is located on the eastern edge of the subglacial channel of Lake Koprzywno. Its altitude ranges from 80 to about 160 m a.s.l. The Rymierowo test area is located on the western edge of the lake’s subglacial channel cut into the outwash surface. The altitude of the area is from 140 m to 165 m a.s.l. Basic hypsometric parameters were calculated for individual test areas (

Table 2. Hypsometric characteristics of the test areas.

Test Area	Hypsometric Parameters
	Minimum Elevation (m a.s.l.)	Maximum Elevation (m a.s.l.)	Elevation Range (m)	Maximum Slope (°)	Mean Slope (°)
Buślarki	56.2	121.9	65.7	51.0	7.1
Piaski Pomorskie	69.9	144.3	74.4	62.3	5.9
Koprzywno	80.3	159.8	79.5	56.7	6.8
Rymierowo	141.5	164.2	22.7	46.5	4.2

). The geological structure of the research area is dominated by glaciofluvial sandy sediments overlying glacial tills [64]. Currently, the test areas are forested. The afforestation of edge zones significantly limits the further development of valley forms, and their dry bottoms are shaped mainly by episodic runoffs in periods of snowmelt or intense rainfall [63].

3. Materials and Methods

3.1. Data Sources

Two DEM data sets were used to test the toolbox, one from field measurements and the other from the remote sensing system (LiDAR). This approach made it possible to compare the hypsometry results according to the data source.

3.1.1. Field Data

For 21 erosion-denudation valleys, precise total station measurements were made using the GPS System SR530-RTK Leica equipment (Leica Geosystems AG, Switzerland) and the Elta R55 (Zeiss AG, Germany). The errors in vertical and horizontal accuracy did not exceed ±0.05 m. The points were measured at the background of the landforms, on slopes, and at the bottom. The elevation point mapping density was, on average, 7 points per m². Data were located in the National Coordinate System, PUWG 1992 (EPSG 2180). The elevation data obtained as a result of the measurements were the basis for the creation of the digital elevation models with a spatial resolution of 0.5 m in the Surfer software using the Kriging interpolation method (Figure 2a,b). The DEM data were then imported into ArcGIS 10.7. The field data were also used as the reference data on which the boundaries of the landforms were determined.

3.1.2. LiDAR Data

The second data source for testing the developed tool was the ALS LiDAR point cloud. A point cloud was generated from aerial laser scans obtained from the Head Office of Geodesy and Cartography, Warsaw, Poland [65]. The LiDAR point clouds were acquired within the national ISOK (Informatyczny System Osłony Kraju) project [66]. The ISOK project was dedicated to the country’s protection against extraordinary threats [66]. Within the ISOK project, the data were classified according to the American Society for Photogrammetry and Remote Sensing standards [67] into the following classes: ground (2), low (3), medium (4) and high (5) vegetation, buildings (6), low points (7), model key-points (8), water (9) and others. Data classification in the ISOK project was performed automatically, but the classification process was followed by a visual inspection and manual correction of misclassified points [66]. The automatic classification of the ground points was performed using the progressive densification filtering algorithm developed by Axelsson [68]. Threshold values of 0.4 m and 2 m were chosen to distinguish low, medium, and high vegetation. The acceptable error for the ground point classification was 1%, whereas the threshold for the other classes was 6% [66].

Data were obtained in the *.las format. The point cloud density was 4–6 per m², which means there was an average point distance of 0.5 m. This density is sufficient to generate a precise model and map both small natural terrain forms and large anthropogenic forms, such as flood embankments, slopes, road embankments, and other similar structures [69]. The height accuracy (mean error) of the ALS point cloud after alignment was ≤0.15 m, while the situational accuracy (mean error) was ≤0.5 m. The coordinate system used was the National Coordinate System, PUWG 1992 (EPSG 2180).

The digital elevation model was generated from the classified point clouds using the LAS Dataset to Raster tool in the ArcGIS software. In the used tool there are two types of interpolation: binning and triangulation. In the binning method, the process of determining the value of a pixel relies on examining the points that fall within the defined pixel. The final value can be chosen as the minimum, maximum, or average of the points [70]. For the digital elevation model, the best setting is to use the minimum value because the algorithm chooses then the lowest point (bare ground) inside the pixel [71]. The triangulation method uses Delaunay triangulation to create a surface from a network of triangular facets defined by nodes and edges that cover the surface, which is then rasterised [70]. The binning method is used for high-density data, while the triangulation method is recommended for very low-density LiDAR data, when binning cannot be used to create an appealing surface [70]. The study conducted by Dawid et al. [72] examined the difference between these two interpolation methods in generating a ground surface. They stated the minor differences between the two types of interpolations using the LiDAR data with high density. On this basis they used the binning interpolation method. In this study the binning method was also applied. This method allows the voids to be filled when there are no points collected within the area represented by a pixel in the resultant raster [70]. The simple type of void filling was set. This option computes the average using up to eight neighbouring pixels (with values), and it is used for the filling of small voids [70]. The digital elevation model was generated with a resolution of 0.5 m.

3.2. Description of the Toolbox

The Morphometry Assessment Tools (MAT) toolbox consists of two modules: (a) Geometry and (b) Hypsometry. The modules were developed using the ArcGIS ModelBuilder visual programming language. The schematic workflow for the MAT toolbox presents the input data, main models, and output data in the standard symbology of the ModelBuilder application (Figure 3). The division into modules allows the user to run the computational processes separately depending on which metrics need to be determined. This also allows shortening the computation time if only basic parameters are needed. Individual modules use different extensions of ArcGIS and operate with different levels of detail. The Geometry module can be used for calculations based only on vector analysis tools, while the Hypsometry module requires a license for the Spatial Analyst extension. The toolbox was tested with the ArcGIS 10.7.1. Models are run from the ArcCatalog or ArcToolbox applications as dialogue boxes (Figure 4). The user can also run the models in edit mode; thus, all intermediate data will also be saved. The input data requirements, the algorithms implemented in the models, and the output data are described in detail below.

3.2.1. Input Data Requirements

The first input data set for the analysis is the boundary of the landform. The boundaries should be prepared as polygons and saved as an ESRI Shapefile (.shp) or geodatabase feature class. The boundary line should be situated on the edges of the concave form (Figure 5). The landform edge polygon can be created manually using vectorisation methods based on a digital elevation model, slope map, contour lines, or orthophotos. It is also possible to apply the methods for the semi-automatic delimitation of the edge polygon using, for example, landform classification based on the Topographic Position Index (TPI) [73] or the algorithms used in drainage network extraction, such as sink and depression evaluation and catchment delineation [74]. In this study, data collected from the field measurements were used to determine the form boundaries. Points measured on the highest ridges (so-called break lines) of the analysed forms were used to determine the vertexes of the boundary polygon of each valley. The data were imported to the ArcGIS software using the Add XY Data function. Then the data were converted to the polygon layer using the Points to Line and the Feature to Polygon tools.

The second data set is a digital elevation model such as DEM, DTM, or DSM, which can be prepared in one of the raster data formats accepted by the ArcGIS software (among others, GeoTIFF, ESRI GRID, PNG, or JPG). DEM should cover the form boundary over the entire area and preferably should also include an outside buffer zone. The polygon layer and DEM need to be in the same projected coordinate system with linear units (meters) to perform the correct calculations.

3.2.2. Geometry Module

The module of basic morphometric parameters is based on the geometry of the polygon features.

Table 3. Definitions of geometric parameters.

Geometric Parameters	Units	Formula	References
Area (A)	(m2)	The projected area enclosed by form boundary
Perimeter (P)	(m)	Length of the horizontal projection of form boundary
Length (L)	(m)	Length of the long axis
Width (W)	(m)	W = A/L	[76]
Orientation (O)	(°)	Orientation (azimuth) of the long axis
Elongation ratio (Re)	(-)	$\mathrm{Re}=1.128\sqrt{A}$/L	[16]
Circularity ratio (Rc)	(-)	Rc = 4πA/P2	[77]
Compactness coefficient (Cc)	(-)	$\mathrm{Cc}=0.282\mathrm{P}/\sqrt{A}$	[15]
Form factor (Ff)	(-)	Ff = A/L2	[76]
Shape Factor (Sf)	(-)	Sf = 1/Ff	[77]

presents the names of parameters and formulas used for the calculations as well as the references of the individual characteristics. The selected parameters are the ones that are most commonly used to characterise the geometry of dry concave landforms.

The area (A) and perimeter (P) of the form are calculated using the Python functions based on the coordinate system of the input layer. The length (L) of the form was defined as the length of the long axis of a rectangle that encircles the identified feature, known in the literature as the minimum bounding box [38,75]. The convex hull method was used for the calculation of the minimum bounding box geometry.

The Minimum Bounding Geometry function was used in the calculation of the length (L) and the orientation (O; azimuth of the major axis) of the polygons. In the literature, there are various approaches used for calculating the length of a form. It can be determined, for example, as a line along the main channel to the farthest point on the basin boundary [16]. The results based on the minimum bounding geometry are a simplified calculation method, but it is possible to apply to various forms and also to those in which the channel cannot be determined. The mean form width (W) is calculated by dividing the area by the perimeter of the form [76]. The elongation ratio (Re) is based on the formula presented by Schumm [16], defined as the ratio of the diameter of a circle of the same area as the form to the maximum form length. The value of Re varies from 0 (in a highly elongated shape) to 1 (in a circular shape). The circularity ratio (Rc) is based on the equation proposed by Miller [77] as the ratio of the area of the form to the area of the circle with the same perimeter. Index values below 0.5 indicate the elongation of the form, while values close to 1 indicate a high circularity. The next index is the form’s compactness coefficient (Cc), defined as the ratio of the perimeter of a form to the circumference of the circular area, which equals the area of the form [15]. The form factor (Ff) presented by Horton [76] determines the ratio of the form area to the square of its length. Form factor values are always below 0.754 (for a perfectly circular form); the smaller the values are, the more elongated the form is. The shape factor (Sf) is calculated as the reciprocal of the form factor [77]. The parameters W, Re, Rc, Cc, Ff, and Sf are computed using Python expressions.

All the geometric parameter values calculated for each polygon in the input layer are saved in the attribute table of the new output polygon layer created in the .shp format. Based on the parameter values derived for all features, the descriptive statistics—maximum, mean, minimum, range, and standard deviation—are calculated and saved in the output statistics table in the .csv format. Output files are created in the folder specified by the user. The names of the output files are created automatically starting from the name of the input layer and ending with the “Geometry” and “Geometry_Statistics” suffixes, respectively.

3.2.3. Hypsometry Module

The module of hypsometric parameters is based on the polygon features and the digital elevation model.

Table 4. Definition of hypsometric parameters.

Hypsometric Parameters	Units	Formula	References
Maximum elevation (Zmax)	(m)	Maximum elevation from DEM
Mean elevation (Zmean)	(m)	Mean elevation from DEM
Minimum elevation (Zmin)	(m)	Minimum elevation from DEM
Elevation range (Zrange)	(m)	Zrange = Zmax − Zmin
Maximum slope (Smax)	(°)	Maximum slope calculated from DEM	[78]
Mean slope (Smean)	(°)	Mean slope calculated from DEM	[78]
Minimum slope (Smin)	(°)	Minimum slope calculated from DEM	[78]
Mean aspect (Amean)	(°)	Mean aspect calculated from DEM	[78]
Maximum Terrain Ruggedness Index (TRImax)	(m)	Maximum Terrain Ruggedness Index calculated from DEM	[79]
Mean Terrain Ruggedness Index (TRImean)	(m)	Mean Terrain Ruggedness Index calculated from DEM	[79]
Minimum Terrain Ruggedness Index (TRImin)	(m)	Minimum Terrain Ruggedness Index calculated from DEM	[79]
Maximum Topographic Wetness Index (TWImax)	(-)	Maximum Topographic Wetness Index calculated from DEM	[80]
Mean Topographic Wetness Index (TWImean)	(-)	Mean Topographic Wetness Index calculated from DEM	[80]
Minimum Topographic Wetness Index (TWImin)	(-)	Minimum Topographic Wetness Index calculated from DEM	[80]
Hypsometric integral (HI)	(-)	$\mathrm{HI}=\frac{Zmean-Zmin}{Zmax-Zmin}$	[81,82]
Volume (V)	(m3)	V = $\sum \mathrm{A}p\times \left(\mathrm{Z}b- \mathrm{Z}a\right)$ Ap—cell area Zb—cell value of the overlay DEM (before surface) Za—cell value of the DEM (after surface)

presents the names of the parameters, formulas, and references of the individual indexes. Based on the literature, the selected parameters are commonly used to characterise the hypsometry of concave landforms.

The maximum (Zmax), mean (Zmean), and minimum (Zmin) height and the elevation range (Zrange) values are extracted directly from the DEM layer using the Zonal Statistics function. The Zrange parameter corresponds to the total relief metric. In order to determine the slope values, a slope raster is generated (in degrees; the values range from 0 to 90) [78]. Then, the Zonal Statistics function is applied to calculate the maximum (Smax), mean (Smean), and minimum (Smin) values in the area of each polygon feature. The next tool used in the Hypsometry model is the aspect, which identifies the direction of the downhill slope faces. The values of each cell in the output raster indicate the compass direction the surface faces at that location (clockwise in degrees from 0 due north to 360) [78]. The mean aspect (Amean) parameter is calculated by averaging all values across the entire surface of each form using the Zonal Statistics function. Mean aspect could be a preferred characteristic compared with the orientation (O) parameter when it comes to evaluating the total effect of solar radiation and wind direction on the form distribution [43].

Another spatial index in the Hypsometry module is the terrain ruggedness, which determines the variability or irregularity in elevation (highs and lows) within a sampled terrain unit. For the calculation, the Terrain Ruggedness Index developed by Riley et al. [79] was selected, as it is widely used, but it is not available in the ArcGIS Desktop system tools. TRI expresses the elevation difference between adjacent cells of a DEM. The metric measures the difference in elevation values from a centre cell and the eight cells directly surrounding it. Then, the eight elevation differences are squared and averaged. The square root of this average results is the TRI value for the centre cell. To calculate the TRI raster in the hypsometry model, an algorithm was developed that uses the Focal Statistics function, which is based on the weighted neighbourhood. The ASCII-weighted kernel files were prepared to sequentially calculate the differences between the centre cell and the eight surrounding it. To make the elevation differences positive, the Abs function that calculates the absolute value of the cells was implemented. For the output TRI raster, the Zonal Statistics function was used to determine the values of the maximum (TRImax), mean (TRImean), and minimum (TRImin) values in the area of each polygon feature.

Although the toolbox is intended for dry concave forms, one spatial index was introduced to assess the topographic control on hydrological processes within the form area. The Topographic Wetness Index (TWI) was developed by Beven and Kirkby [80] within TOPMODEL. Among other things, TWI is widely used to identify hydrologic flow paths, but it is still missing in the ArcGIS Desktop system tools. The index is a function of the slope and the upstream contributing area:

TWI = ln (a/tanb)

(1)

where:

• a—upstream contributing area;
• tanb—local slope in radians.

The index is a number expressed without units, for which higher values indicate the water accumulation areas while lower values represent crests and ridges. The workflow of the TWI raster calculation in the Hypsometry module includes the Slope tool, Hydrology tools, and Math tools. As for the previous spatial indexes, the Zonal Statistics function was used to determine the maximum (TWImax), mean (TWImean), and minimum (TWImin) values in the area of each landform.

The following parameter is the hypsometric integral (HI), which is defined as the area below the hypsometric curve, and it has been used in classifying areas according to their geological stage of development [81]. HI is an indicator of the erosion cycle. The erosion cycle is defined as the total time it takes to denudate a unit of terrain to the level of the base erosion. The erosion cycle can be divided into three stages: monadnock (old) (HI ≤ 0.3), equilibrium state or mature stage (0.3 ≤ HI ≤ 0.6), and inequilibrium or young stage (HI ≥ 0.6). The hypsometric curve used to quantify HI represents the relative proportion of the form area below a given height (Figure 6) [81,82]. The shape of the hypsometric curve is related to its stage of form development. Convex hypsometric curves are typical of a youthful stage, s-shaped curves are related to a mature stage, and concave curves are indicative of a peneplain stage [16,42,81]. The hypsometric curves can also show the dynamic of geomorphological processes [2,83].

In the workflow of the HI parameter calculation, we followed Matos and Dilts [51]. The hypsometric integral is calculated by slicing the area into elevation bands and plotting the cumulative area for each band. Additionally, an algorithm was developed to plot a graph for the hypsometric curve of each analysed landform. For this purpose, a graph template was prepared, which is used in the Make Graph tool in the Hypsometry model.

The last characteristic is the volume (V) of the concave landform. The first step of the analysis is the cover plate raster creation according to the upper edge of the form. The polygon boundary is converted to points based on vertexes. Then, for the points, the elevation values from the DEM are extracted. Using altitude values, the interpolation is performed using the Natural Neighbor algorithm. The choice of the method was preceded by testing other interpolation tools available in the ArcGIS software, such as IDW, Kriging, Spline, Trend, and also the TIN model. It was stated that the Natural Neighbor method gave the most appropriate results—i.e., the interpolated cover surface fitted well to the edge of the form and at the same time was smooth in the interior area (Figure 7). A comparison of the results according to the interpolation method used is presented in the results section of this study.

The second step of the analysis is the estimation of the material volume. The resulting raster is generated from two surfaces: the overlay DEM and the DEM of the form. The volume for each cell is calculated according to the formula:

V = Ap ∗ (Zb − Za)

(2)

where:

• V—volume [m³];
• Ap—cell area [m²];
• Zb—cell value of the overlay DEM (before surface);
• Za—cell value of the DEM (after surface).

The final volume (V) parameter is calculated by summing up all the cells in the area of each form.

All the output files of the Hypsometry module are saved in a folder specified by the user. Output files include an output polygon layer in the .shp format—the parameter values are saved in the attribute table of the layer; an output statistics table in the .csv format—descriptive statistics (maximum, mean, minimum, range, standard deviation) of parameter values calculated for all features in the polygon layer; output raster layers in the GeoTIFF (.tif) format—DEM, Slope, Aspect, TRI, and TWI—clipped to the area of the polygon features; output graphs data in the .jpg and .csv format—hypsometric curves determined for the areas of all polygon features (by FID field) (Figure 8). The names of the output files are created automatically, starting from the name of the polygon input layer and ending with the “Hypsometry”, “Hypsometry_Statistics”, “DEM’, “Slope”, “Aspect”, “TRI”, “TWI”, “Hypsometric_curve” suffixes.

4. Results

4.1. Landform Geometry Assessment

The results calculated with the use of the Geometry module of the MAT toolbox for 21 erosion-denudation valleys are presented in

Table 5. Geometric parameters of the erosion-denudation valleys.

Test Area	Valley	Geometric Parameters
		A Area (m2)	P Perimeter (m)	L Length (m)	W Width (m)	O Orientation (°)	Re Elongation Ratio (-)	Rc Circularity Ratio (-)	Cc Compactness Coefficient (-)	Ff Form Factor (-)	Sf Shape Factor (-)
Buślarki	1	3529.01	286.23	122.10	28.90	107.50	0.55	0.54	1.36	0.24	4.22
	2	457.07	125.34	56.85	8.04	95.30	0.42	0.37	1.65	0.14	7.07
	3	1019.36	222.39	104.78	9.73	103.10	0.34	0.26	1.96	0.09	10.77
	4	434.06	106.12	39.98	10.86	124.30	0.59	0.48	1.44	0.27	3.68
	5	2637.33	389.97	145.92	18.07	91.00	0.40	0.22	2.14	0.12	8.07
	6	681.23	136.02	57.66	11.81	85.60	0.51	0.46	1.47	0.20	4.88
	7	1588.49	266.74	110.49	14.38	106.80	0.41	0.28	1.89	0.13	7.69
	8	2372.08	251.17	112.16	21.15	104.10	0.49	0.47	1.45	0.19	5.30
Piaski Pomorskie	9	1368.46	208.62	94.32	14.51	170.90	0.44	0.40	1.59	0.15	6.50
	10	3565.03	260.54	95.02	37.52	140.30	0.71	0.66	1.23	0.39	2.53
	11	1476.07	207.78	85.48	17.27	25.70	0.51	0.43	1.53	0.20	4.95
	12	887.51	138.55	52.54	16.89	12.00	0.64	0.58	1.31	0.32	3.11
	13	971.33	212.99	74.71	13.00	129.30	0.47	0.27	1.93	0.17	5.75
	14	775.84	146.13	63.40	12.24	25.60	0.50	0.46	1.48	0.19	5.18
	15	2015.84	267.09	109.86	18.35	58.30	0.46	0.36	1.68	0.17	5.99
Koprzywno	16	2557.94	283.54	122.84	20.82	47.50	0.46	0.40	1.58	0.17	5.90
	17	4161.63	330.26	142.88	29.13	74.80	0.51	0.48	1.44	0.20	4.91
	18	1649.22	180.37	79.26	20.81	83.00	0.58	0.64	1.25	0.26	3.81
Rymierowo	19	3607.61	262.82	96.52	37.38	124.40	0.70	0.66	1.23	0.39	2.58
	20	1124.03	146.11	59.50	18.89	127.70	0.64	0.66	1.23	0.32	3.15
	21	847.35	111.42	39.89	21.24	155.70	0.82	0.86	1.08	0.53	1.88
Statistics	Max	4161.63	389.97	145.92	37.52	170.90	0.82	0.86	2.14	0.53	10.77
	Mean	1796.50	216.20	88.87	19.09	94.90	0.53	0.47	1.52	0.23	5.14
	Min	434.06	106.12	39.89	8.04	12.00	0.34	0.22	1.08	0.09	1.88
	Std	1149.36	76.89	31.77	8.23	42.82	0.12	0.16	0.28	0.11	2.13

.

Erosion-denudation valleys show a diversified area (A), ranging from 434.06 m² (valley 4—Buślarki) to 4161.63 m² (valley 17—Koprzywno) (mean, 1796.50 m²; std deviation, 1149.36). The length (L) of the forms, which is the result of linear erosion processes taking place at their bottoms [63], does not exceed 150 m. The maximum values were recorded for the valley of Piaski Pomorskie—145.92 m (mean, 88.87 m; std deviation, 31.77). They are generally short forms with a relatively narrow bottom. The mean widths (W) of the valleys range from 8.04 to 37.5 m (mean, 19.09 m; std deviation, 8.23). Compared with other forms of this type occurring in the young glacial area of the Polish Lowlands, these are small forms [84].

The form elongation ratio (Re) ranges from 0.34 to 0.82, with a mean of 0.53 (std deviation, 0.12). The values of the Re correspond very well with the circularity ratio (Rc), proving that most of the studied valleys have a more elongated bottom. The elongation of the form is the result of linear erosion processes occurring along the longitudinal axis [85,86]. In the case of valleys where the Rc value is above 0.5, it can be concluded that longitudinal processes, such as linear erosion, occurred with a similar intensity as transverse processes, such as mass movements, solifluction, and slope wash, responsible for, among other things, widening the bottom of the forms. Such a situation was recorded in the Rymierowo area, where the Rc value for valley 21 is 0.86. This form belongs to the shortest type, as confirmed by the fact that in its development, the longitudinal processes acted with the same intensity as the transverse processes, and the form is at the initial stage of development.

The dimensionless form factor (Ff) represents the ratio of the valley area to the square of its length. Low values of the form factor (Ff) confirm the elongation of forms. In the case of the studied valleys, this factor varies from 0.09 to 0.53 (mean, 0.23; std deviation, 0.11). The dimensionless shape factor (Sf), being the inverse of the form factor, takes the highest values for the most elongated forms. These factor values range from 1.88 to 10.77 (mean, 5.14; std deviation, 2.13).

4.2. Landform Hypsometry Assessment

The parameters calculated with the use of the Hypsometry module of the MAT toolbox for 21 erosion-denudation valleys are presented in

Table 6. Hypsometric parameters of the erosion-denudation valleys.

Test Area	Valley	Hypsometric Parameters
		Zmax Maximum Elevation (m)	Zmean Mean Elevation (m)	Zmin Minimum Elevation (m)	Zrange Elevation Range (m)	Smax Maximum Slope (°)	Smean Mean Slope (°)	Smin Minimum Slope (°)	Amean Mean Aspect (°)	TRImax Maximum TRI (m)	TRImean Mean TRI (m)	TRImin Minimum TRI (m)	TWImax Maximum TWI (-)	TWImean Mean TWI (-)	TWImin Minimum TWI (-)	HI Hypsometric Integral (-)
Buślarki	1	78.06	73.27	64.57	13.49	33.24	14.20	0.20	109.96	0.50	0.21	0.02	10.42	3.20	0.55	0.64
	2	78.53	74.56	69.95	8.58	30.95	17.32	2.52	112.92	0.48	0.27	0.10	8.41	2.77	0.62	0.53
	3	81.73	76.11	70.57	11.16	32.97	14.35	0.78	108.73	0.51	0.23	0.04	9.85	2.99	0.58	0.49
	4	81.40	78.67	75.87	5.53	27.85	13.97	2.95	139.81	0.41	0.22	0.08	7.99	2.77	0.71	0.50
	5	82.35	75.90	67.09	15.26	45.27	21.37	0.46	122.07	0.80	0.35	0.05	11.21	2.56	0.12	0.57
	6	80.96	77.36	73.40	7.56	37.90	18.36	3.69	112.71	0.58	0.29	0.08	8.90	2.60	0.34	0.52
	7	81.21	74.71	66.84	14.37	36.58	16.33	0.78	115.49	0.56	0.25	0.04	10.42	2.83	0.41	0.54
	8	80.99	75.25	69.41	11.58	42.12	16.80	0.82	103.15	0.72	0.26	0.03	13.51	2.84	0.16	0.50
Piaski Pomorskie	9	100.99	92.83	83.06	17.93	36.53	19.95	3.33	169.72	0.57	0.31	0.06	9.06	2.72	0.32	0.54
	10	100.62	92.18	84.97	15.65	39.17	17.42	0.20	166.60	0.64	0.25	0.02	11.94	3.07	0.33	0.46
	11	103.81	101.66	97.84	5.97	37.00	11.95	0.43	185.87	0.58	0.18	0.02	10.99	3.27	0.35	0.64
	12	100.60	97.51	93.99	6.61	43.38	22.31	1.47	181.89	0.75	0.36	0.06	9.54	2.39	0.06	0.54
	13	95.19	91.49	87.52	7.67	38.24	20.36	1.23	145.52	0.56	0.30	0.04	9.95	2.56	0.41	0.54
	14	98.48	95.34	94.22	4.26	40.56	16.93	0.10	202.81	0.63	0.23	0.01	17.19	3.26	0.23	0.37
	15	102.54	98.61	96.70	5.84	40.87	15.37	0.10	233.24	0.67	0.22	0.01	13.57	3.30	0.14	0.36
Koprzywno	16	129.08	117.25	105.16	23.92	52.04	34.62	2.49	226.44	0.95	0.56	0.06	10.47	2.11	−0.11	0.50
	17	128.74	116.56	104.22	24.52	50.95	33.67	2.02	222.43	0.96	0.53	0.12	10.91	2.27	−0.15	0.50
	18	94.91	90.20	85.23	9.68	44.27	22.84	0.62	222.71	0.75	0.35	0.02	12.87	2.54	0.09	0.51
Rymierowo	19	157.76	150.84	146.32	11.44	38.15	16.82	1.26	144.23	0.60	0.24	0.02	11.37	3.16	0.31	0.39
	20	159.41	156.51	151.22	8.19	43.19	24.44	0.10	134.17	0.73	0.39	0.01	11.85	2.47	0.12	0.64
	21	158.11	155.71	150.43	7.68	34.46	16.40	0.23	140.08	0.51	0.24	0.02	7.89	2.64	0.40	0.68
Statistics	Max	159.41	156.51	151.22	24.52	52.04	34.62	3.69	233.24	0.96	0.56	0.12	17.19	3.30	0.71	0.68
	Mean	103.59	98.22	92.31	11.28	39.32	19.32	1.23	157.17	0.64	0.30	0.04	10.87	2.78	0.29	0.52
	Min	78.06	73.27	64.57	4.26	27.85	11.95	0.10	103.15	0.41	0.18	0.01	7.89	2.11	−0.15	0.36
	Std	27.14	26.86	26.94	5.70	6.03	5.86	1.15	44.00	0.15	0.10	0.03	2.18	0.34	0.23	0.09

.

The elevation range (Zrange) in the 21 analysed valleys, representing the difference in the maximum elevation (Zmax) measured on the valley edge and the minimum elevation (Zmin) at the valley bottom, vary from 4.26 to 24.52 m (mean 11.28, std deviation 5.70). The most considerable elevation differences were recorded for the Koprzywno area, where the Zrange value for the valley Koprzywno 17 is 24.52 m. The lowest value of the Zrange was recorded for the valley Piaski Pomorskie 14 (4.26 m). Low values of the Zrange parameter correspond to a slight variation in the relief of the studied valleys, which may indicate that the role of overland flow in shaping the relief of the valleys was small [84].

The valley side slopes (Smax, Smin) are varied, as the maximum gradient range from 27.85° to 52.04° and the average (Smean) from 11.95° to 34.62°. The highest valley side slopes were recorded for the Koprzywno area (minimum 44.27°, maximum 52.04°) and the smallest for the Buślarki area (minimum 27.85°, maximum 45.27°). Undoubtedly, this is the result of the relief energy of the areas within which the valleys occur. The most considerable differences were observed for the Koprzywno and Rymierowo test areas, reaching over 100 m, while the smallest differences were found in the Buślarki area, amounting to just over 40 m. Similarly, the highest valley side slopes were recorded for the Koprzywno area and the lowest for the Buślarki area. These topographical conditions are conductive to gully erosion processes.

The mean aspect (Amean) parameter values of the valleys indicate the eastern direction in the Buślarki test area, the south direction in the Piaski Pomorskie test area, the southwest direction in the Koprzywno test area, and the southeast direction in the Rymierowo test site. These values are in line with the field observations and aspect maps calculated using the model.

The TRI (TRImax, TRImean, TRImin) values are low and show little variability. This means that the terrain ruggedness in the area of the valleys is low. The side of the valleys are smooth and devoid of irregularity. These results are also determined by the resolution of the DEM data.

The Topographic Wetness Index (TWImax, TWImin) values range from −0.15 to 17.19. The standard deviation of the average (TWImean) and minimum TWI (TWImin) values within all the valleys is small. The close-up values of the TWI index may indicate a similar process of overland flow formation in the analysed forms.

Hypsometric integral (HI) values are within the range of 0.36–0.68, with a mean value of 0.52. The maximum recorded value of 0.68 occurs in the Rymierewo test area and the lowest 0.36 in the Piaski Pomorskie test area. It should be noted that 5 of the 21 analysed valleys show a HI index value of above 0.6, which indicates that they are young forms [81]. The other valleys are in the mature stage of development. Three valleys may be included in the older forms (HI 0.36–0.39). The calculated HI parameter values confirm the shapes of the hypsometric curves (Figure 9a–d). There are 5 convex hypsometric curves that are typical of a youthful stage, 13 s-shaped curves that are related to a middle stage, and 3 concave curves are indicative of an old stage.

The Volume (V) parameter includes the surface of the forms and the elevation range within the forms. Based on this holistic parameter, the results are presented according to the DEM data source (

Table 7. Volumetric parameter of the erosion-denudation valleys.

Test Area	Valley	Volumetric Parameter
		V Volume (m3) Based on the Field Data	V Volume (m3) Based on the LiDAR Data	V Error (%)
Buślarki	1	5077.81	5020.97	−1.12
	2	243.10	212.30	−12.67
	3	553.97	489.15	−11.70
	4	212.52	154.13	−27.47
	5	3614.32	3728.28	3.15
	6	600.72	557.98	−7.11
	7	1222.37	1229.12	0.55
	8	2931.81	3035.10	3.52
Piaski Pomorskie	9	1584.98	1532.78	−3.29
	10	11,637.23	7257.00	−37.64
	11	1864.29	1506.53	−19.19
	12	1331.84	1187.95	−10.80
	13	1215.01	1219.49	0.37
	14	1212.70	920.31	−24.11
	15	3674.79	2376.38	−35.33
Koprzywno	16	3732.90	3769.40	0.98
	17	16,833.17	16,809.18	−0.14
	18	3244.13	3116.72	−3.93
Rymierowo	19	12,908.40	13,967.69	8.21
	20	2146.14	2392.79	11.49
	21	232.69	305.49	31.29

). The obtained values indicate the significant diversification of the volume from 212.52 to 16,833.17 m³ based on the field data and from 154.13 to 16,809.18 m³ based on the LiDAR data. Within individual research sites, valleys with both small and large volumes of denuded material were identified. The highest volume values were recorded for valley 17 in the Koprzywno test area. The volume of the remaining two valleys, 16 and 18, is close to the average value, as calculated both on the basis of field data and LiDAR data. The lowest volume values were recorded for the Buślarki test area. Among the eight mapped valleys of this area, seven have values below the average, which allows us to state that they are the smallest forms among all the studied valleys.

The differences in the results for the volumetric parameter were at an average level of 12% in comparing the two data sources (Table 7). It should be stated that the calculated volume values based on the LiDAR data are, in most cases, underestimated in relation to the measurements made directly in the field. They are on average 15% lower than those calculated based on field digital elevation models. The valleys determined based on the LiDAR data are shallower in relation to their depth measured in the field. More considerable differences should be noted for the valleys that were particularly elongated, narrow, and deep. It should also be noted that in the case of eight forms, the volume values calculated from the LiDAR data were higher in relation to the field volume. On average, the values that were overestimated differed by 7% from the field data. The most significant differences, in this case, were visible in the Rymierowo area. These valleys are located directly on the shoreline of the lake. They also have the greatest width compared with the other forms.

Comparing the results due to the interpolation method applied to the overlay DEM creation, it can be seen that the average error was 5% (

Table 8. Volumetric parameter according to the overlay DEM interpolation method.

Test Area	Valley	VNaturalNeighbor	VIDW	VIDW Error (%)	VKRIGING	VKRIGING Error (%)	VSPLINE	VSPLINE Error (%)	VTREND	VTREND Error (%)	VTIN	VTIN Error (%)
Buślarki	1	5077.81	5202.73	2.46	5228.72	2.97	5321.72	4.80	5963.83	17.45	4756.05	−6.34
	2	243.10	219.90	−9.54	242.39	−0.29	236.67	−2.64	278.84	14.70	225.38	−7.29
	3	553.97	527.22	−4.83	552.86	−0.20	546.52	−1.34	640.18	15.56	598.36	8.01
	4	212.52	183.74	−13.54	209.65	−1.35	215.67	1.48	228.26	7.41	227.58	7.09
	5	3614.32	3569.90	−1.23	3604.20	−0.28	3633.83	0.54	3886.42	7.53	3658.12	1.21
	6	600.72	604.46	0.62	610.57	1.64	663.95	10.53	653.87	8.85	616.90	2.69
	7	1222.37	1123.79	−8.06	1216.11	−0.51	1260.27	3.10	1298.92	6.26	1168.92	−4.37
	8	2931.81	3060.27	4.38	3008.40	2.61	3044.85	3.86	3021.89	3.07	2942.23	0.36
Piaski Pomorskie	9	1584.98	1555.28	−1.87	1584.70	−0.02	1613.41	1.79	1827.09	15.28	1576.29	−0.55
	10	11,637.23	9623.34	−17.31	11,295.39	−2.94	11,711.67	0.64	10,893.41	−6.39	11,921.41	2.44
	11	1864.29	1678.19	−9.98	1867.58	0.18	2083.01	11.73	1535.24	−17.65	1912.48	2.58
	12	1331.84	1269.38	−4.69	1310.65	−1.59	1341.30	0.71	1284.05	−3.59	1313.75	−1.36
	13	1215.01	1184.28	−2.53	1194.35	−1.70	1241.64	2.19	1242.05	2.23	1205.81	−0.76
	14	1212.70	1197.52	−1.25	1229.15	1.36	1426.48	17.63	1316.61	8.57	1154.62	−4.79
	15	3674.79	3703.62	0.78	3625.00	−1.36	4294.24	16.86	3549.13	−3.42	3738.13	1.72
Koprzywno	16	3732.90	3716.31	−0.44	3720.15	−0.34	4053.43	8.59	3611.89	−3.24	3739.19	0.17
	17	16,833.17	16,249.21	−3.47	16,611.01	−1.32	18,291.73	8.66	17,519.12	4.07	16,622.86	−1.25
	18	3244.13	3059.52	−5.69	3248.73	0.14	3169.06	−2.31	3431.49	5.78	3253.14	0.28
Rymierowo	19	12,908.40	11,086.76	−14.11	12,009.05	−6.97	15,025.12	16.40	15,020.46	16.36	12,728.45	−1.39
	20	2146.14	2226.96	3.77	2182.22	1.68	2228.42	3.83	2279.44	6.21	2146.22	0.00
	21	232.69	247.00	6.15	229.44	−1.40	273.31	17.46	272.62	17.16	199.29	−14.35

). The greatest differences are visible in the case of the Trend, Spline, and IDW methods, which excludes them from being used in the adopted calculation procedure. The other methods give very similar results, but the analysis of the spatial distribution of values (Figure 7) supports the use of the Natural Neighbor method.

5. Discussion

5.1. Applicability of the MAT Toolbox for the Morphometric Assessment of the Concave Landforms

The toolbox was tested on an example of 21 erosion-denudation valleys. A large variety of individual morphometric parameters within the same test area was observed, which may indicate that the analysed forms are at different stages of development. Indeed, the formation of the valleys was predisposed to the shape of the area within which the valleys occurred. The development of the landforms was not synchronous on the examined test areas due to a few factors. The most important include: the location of the valleys in geomorphologically different units (morainic uplands, outwash plains), differences in the elevation range of the ridges zones, their different longitudinal and latitudinal orientation, and as a result, differences in the aspect of their slopes. The hypsometric integral (HI) parameter values also indicate different stages of development of the valleys. It should be noticed that for genetically different landforms, such as gullies, different phases of form development were observed [2,8,87]. To test the applicability of the developed procedure on genetically different landforms, the MAT toolbox was applied to the 11 channel heads. The input data (DEM, landforms boundaries) used for calculations were obtained from the research conducted by Mazurek [88]. Channel heads determined by groundwater outflows assume the form of alcoves or niches that dissect slopes to form the headwater depressions [89]. The values of the geometric and hypsometric characteristics calculated using the MAT toolbox were consistent with the range presented by Mazurek [88], which confirms that it is universally applicable. The results are presented in the form of tables, maps, and a graph in the Supplementary Materials (Figures S1 and S2, Tables S1 and S2).

The morphometric parameters of the examined forms may be the basis for their typology and thus may be the basis for the identification of groups of landforms with similar morphometric features.

5.2. Advantages and Limitations of the MAT Toolbox

The methodology presented in the study is an approach based on the automatic extraction of the morphometry of concave landforms. The calculation is based on a minimal set of input data: the form boundaries and the digital elevation model (DEM) of the forms. The morphometric assessment using the presented toolbox can be performed for any type of concave form, regardless of its scale and shape. The calculations did not take into account the parameters related to the drainage network, which are characteristic of the analysis of river basins. This allows the use of algorithms for dry forms, where there are no watercourses, such as dry erosion-denudation valleys, sinkholes, or glacial cirques.

Calculations of the basic parameters are based on the geometry of the form boundary. The length of the form is determined by the minimum bounding box. This is a simplified method, but the form polygon can have a very complicated geometry, unlike in other tools of this type, such as for example the ACME tool [43], which is adapted to glacial cirques with simple geometry. The main watercourse is also not needed for calculations. On the other hand, the width of the form based on the surface-to-length ratio allows for the average value of the entire form to be obtained, and not only on selected cross profiles.

Calculations of the hypsometric parameters are based on the geometry of the form boundary and the digital elevation model. The Hypsometry module contains the tools, which are widely used, but they are not available in the ArcGIS Desktop system tools. These are the following functions: Terrain Ruggedness Index (TRI), Topographic Wetness Index (TWI), hypsometric integral (HI), hypsometric curves, and volume of the concave landform. The results of the mentioned parameters were successfully validated using the independent tools. The TRI parameter values were consistent with the results of the Terrain Ruggedness Index algorithm from the GDAL library available in the QGIS software. The TWI results corresponded with the calculations using the Topographic Wetness Index (TWI) function available in the SAGA GIS software. An important parameter of the presented toolbox compared with other tools of this type is the hypsometric integral (HI) and hypsometric curve graph. The results were successfully validated using the CalHypso script [42]. The volume parameter was validated using the algorithm developed by Báčová et al. [46]. The results were very similar. Minor differences resulted from the different overlay DEM interpolation methods.

The computational time of the MAT toolbox was examined. Calculations of the geometric parameters of the 21 valleys with the mean area of almost 1800 m² and the mean perimeter of 216 m were made in 4.12 s. Calculations of the hypsometric parameters of the 21 valleys using DEM of 0.5 m resolution were made in 12 min 36 s.

For comparison, the calculation of only one parameter, the volume, based on the same input data in the rillStats tool [46], took 1 h 4 min 14 s. When performing calculations using the MAT tool, it should be taken into account that the time of the hypsometric calculations is dependent on the size and spatial resolution of the input DEM data and also on the computing capabilities of the user computer. High-resolution data may take a longer time to compute.

The toolbox may be improved with new functionalities in the future. An important issue to resolve is the automatic detection of the form boundaries. The implementation of such a function would allow for the elimination of potential errors that may occur on the part of the user during the vectorisation procedure.

In conclusion, the important methodological advantages of the MAT toolbox compared with other tools can be highlighted:

(1)

Comprehensive morphometric assessment using only one tool in one software;

(2)

Minimal set of input data required;

(3)

Automatic extraction of the parameters which are widely used, but not available in the ArcGIS Desktop system tools: Terrain Ruggedness Index (TRI), Topographic Wetness Index (TWI), hypsometric integral (HI), hypsometric curves, and volume of the concave landform;

(4)

User-friendly interface (automatic names of the output files; there is no need to set additional parameters by the user);

(5)

Short computational time;

(6)

Automatic calculation of the descriptive statistics and the output raster layers of the spatial parameters.

5.3. Impact of Data Quality and Resolution

The results obtained with the use of the Morphometry Assessment Tools toolbox depend on the quality of input data. The polygon layer of the form’s boundary may be prepared with the use of various methods, such as the manual vectorisation or semi-automatic extraction of objects. If the form’s boundary is processed manually, there may be slight differences in its location depending on the accuracy and perception of the person who performs the digitisation. Báčová et al. [46] investigated the influence of the position and shape of eight different polygons of vectorised rills boundaries on the volume estimation result. The results showed that the calculated material volume differed depending on the area and shape of the polygons, but the differences were as low as ±4%, which was considered a good result.

The result of the calculations of individual hypsometric parameters may be influenced by the DEM data source, their spatial resolution, and their elevation accuracy. In this study, calculations were based on elevation data collected in the field using the total station measurements and based on the airborne laser scanning data. In the case of data obtained using LiDAR methods, errors may be related to the type of the scanner used and its calibration, the density of scanning points, the type of vegetation in the study area, and the interpolation algorithms used at the stage of data post-processing. In the conducted study the maximum difference in the hypsometric assessment was in the range of 19–38% in a few cases based on the LiDAR data compared with the reference data (Table 7). Based on the acquired LiDAR data specification [66] and the literature on the application of these data in the area of northwest Poland [90], it should be stated that the analysis of the forested small landforms based on these LiDAR data may be limited mainly due to fallen trees and piles of brushwood. These elements may be misinterpreted as the ground [90]. Both high and low levels of vegetation can be found within the studied valleys. Tree branches and windthrows in the bottom of the forms are characteristic in the valleys, where higher values of LiDAR data errors were recorded (valley nos. 4, 10, 11, 14, 15). The volume in these cases was significantly underestimated. The second reason for the differences in the results of the calculations of morphometric parameters is the relatively highly inclined valley side. The underestimation of the volume of gullies based on the LiDAR data was observed also by Perroy [22]. Measurements performed in the field using the total station allowed for the precise filtration of points in terms of their location on the ground, and their acquisition density was higher (7 per m²) than that in the LiDAR data (4–6 per m²). This contributed to a more accurate representation of the relief within the analysed forms.

To determine the influence of the spatial resolution of the used DEMs on the accuracy of the obtained results, analyses were conducted on different DEM cell sizes. The original LiDAR DEM with a resolution of 0.5 m was resampled to a resolution of 1, 2, 3, 5, and 10 m. Because the analysed valleys are very small forms, the calculation based on the coarser cell size of 10 m was possible only for the four biggest valleys (nos. 1, 10, 17, and 19 with an area of more than 3000 m² and mean width about 30 m). This threshold is related to the procedure for determining the slope, TWI, and TRI layers, where each cell is calculated based on a 3 × 3 cell neighbourhood. The results of the selected hypsometric parameters for four valleys are presented in Figure 10. The analysis was also performed for the nine medium size valleys (with an area of more than 1000 m²) using the resolution up to 5 m and for the eight small size valleys (with an area less than 1000 m²) using the resolution up to 3 m. The results are included in the Supplementary Materials (Figures S3 and S4).

It can be seen that the mean elevation (Figure 10a) increases slightly as we move to coarser resolution from 1 m until 10 m. For the other parameters, the differences are larger. Relatively, the greatest changes should be noted for the TRI (Figure 10e) and Slope (Figure 10c) parameters, which strongly depend on the DEM resolution. The same trend of parameter differences is seen for medium and small valleys (see Supplementary Materials (Figures S3 and S4). It is also worth noting that the differences in parameter values are small in the 1–3 m resolution range. The differences increase substantially when we compare the 1 m and 10 m resolution DEMs. The results indicate the loss of details by resampling the higher resolution DEM to coarser resolution. The results also suggest that if high-resolution DEM is available, it should be used in the MAT toolbox analysis instead of low-resolution DEM.

6. Conclusions

This paper presents tools that consist of sets of algorithms designed and implemented in the ArcGIS software. The functions gathered in the toolbox enable the comprehensive morphometric evaluation of any type of dry concave landform. The assessment parameters include the geometry, hypsometry, and volume of forms, calculating a total of 26 metrics. The results obtained with the toolbox were discussed with the example of 21 dry erosion-denudation valleys located in the young glacial area of the Polish Lowlands.

The calculated characteristics allow for inferences to be made about the stages of development of particular forms, their geomorphological processes, and their classification. This tool can be useful in geomorphological analyses, especially to calculate the hypsometric curves and estimate the volume of material erosion in cases where the original form surface is unknown. The recognised morphological types of erosion-denudation valleys show diversified sizes and reliefs. They represent forms that developed in the edge zones of morainic upland and outwash plains. The morphometric characteristics of the valleys record the stages of their development, thus providing information on the functioning of the processes and the transformations of the relief in the Late Glacial and Holocene periods.

Calculations were based on elevation data measured in the field and LiDAR data. The representation of the relief based on LiDAR data could be disturbed by the presence of vegetation, as well as by the shape of the forms, often narrow and deeply incised. The DEM input to the tool can also come from other remote sensing sources, such as terrestrial laser scanning (TLS) or UAV orthophotos.

The tool can be used for various types of concave terrain forms, regardless of their size and genesis. Morphometric evaluation with the use of an algorithm can be performed for natural and anthropogenic forms, such as open-cast mines or excavations.