A GIS method for volumetric assessments of erosion rills from digital surface models

ABSTRACT Modern non-contact methods for data acquisition are becoming widely used for monitoring soil erosion and for assessing soil degradation after rainfall events. Photogrammetric methods are especially favored to obtain a detailed and precise digital surface model (DSM) of the surveyed area. This paper introduces the algorithm and its Python implementation as a tool for ArcGIS software, which makes efficient automatic calculations of the volume of erosion rills or gullies. The input parameters are a DSM, and the rill edge polygon. The method was tested on an artificially created rill, where the result acquired using presented method was compared to the real volume. The comparison showed that the algorithm may underestimate the volume by 10–15%. In addition, the influence of the position of the rill edge polygon was tested on two DSMs of erosion rills. The resulting volumes of the rills, calculated on the basis of eight different edge polygons, varied by 5%. The algorithm also automates interpolation of the surface prior to erosion, which simplifies its usage in firstly monitored regions. The algorithm can also be used for volumetric analyses in other research areas and it is made available as a supplement of the publication.


Introduction
In recent years, remote sensing has become widely used for monitoring soil erosion. For detailed monitoring of damage in individual catchments, or in agricultural fields, the most efficient data are supplied by very highresolution satellite images (Desprats et al., 2013), by conventional aerial photographs (Fulajtár, 1999;Maugnard et al., 2014), and nowadays also by images obtained using Unmanned Aerial Vehicles (UAVs) (d Oleire-Oltmanns, Marzolff, Peter, & Ries, 2012;Eltner, Baumgart, Maas, & Faust, 2014). Remote sensing data are suitable mainly for identifying erosion objects and for monitoring their evolution in time. A digital elevation model can be obtained from this data (especially from UAV images) by photogrammetry, which is useful for volume analysis of soil erosion effects.
Many researchers have compared contact methods (direct, destructive) and noncontact methods for rill erosion measurement: their accuracy, their time and cost requirements, and their efficiency. According to Jester and Klik (2005), major disadvantages of contact elevation measurement devices are the deformation of the original soil surface profile and limited resolution (vertical 1-2 mm, horizontal 20-25 mm). The results of their study showed that contact methods are not very precise, and that pin meter measurements were the slowest to perform of all the methods that they compared. As far as indirect methods were concerned, high resolution and high precision (sub millimeter range) were obtained with a laser scanner, while for the photogrammetric method the measurement uncertainty was approximately 1-2 mm. Photogrammetry provided the fastest data acquisition. However, DEM generation from photos was slow, needed extensive expert knowledge, and required expensive hardware and software. In addition, studies have shown that the use of laser techniques is very cumbersome, (García Moreno, Saa Requejo, Tarquis Alonso, Barrington, & Díaz, 2008;Jester & Klik, 2005), and the method is more expensive than photogrammetry and direct methods (Castillo et al., 2012). Castillo et al. (2012) concluded that photogrammetry has promising field applications in soil erosion studies. Nevertheless, authors claimed that if cost and time requirements are considered as well as accuracy, direct methods may still be an optimum approach for large-scale studies. In contrast, according to Smith and Vericat (2015), Structure from Motion (SfM) photogrammetry can be implemented easily across a particularly wide range of scales including the landscape scale. Jester and Klik (2005) concluded that each measurement technique has its own field of application.
The advantage of modern indirect methods is acquisition of high-resolution data, such as point clouds, which enable the generation of digital surface models (DSMs) and orthophotos of observed objects, such as erosion rills and gullies for further volumetric analysis. The general approach for calculating volume changes is to subtract the DSMs created for two different points in time, so-called DEMs of Difference (DoD) approach (Wheaton, Brasington, Darby, & Sear, 2010). This method was used for example by Wells et al. (2016) for quantifying rill and ephemeral gully channel erosion together with sheet erosion on agricultural fields, measurement field plots, and in a laboratory setting, using the photogrammetric technique. Time series of photogrammetric DSMs were also used by Gessesse et al. (2010) for quantifying the change in tillage-induced microtopography due to inter-rill and rill erosion on field plots. Field-scale quantification of soil loss was carried out by Martínez-Casasnovas et al. (2002) using DSMs generated from measurements by a total station before and after a heavy rainstorm.
In natural conditions it is usually not possible to get information about the surface condition before an erosion event. In order to quantify the "starting volume" of the rill or gully for the first monitoring date, it is necessary to reconstruct the (potential) original former surface prior to erosion by creating some kind of 3D cover plate (d´Oleire-Oltmanns et al., 2012). DʼOleire-Oltmanns et al. (2012) and Peter, d'Oleire-Oltmanns, Ries, Marzolff, and Hssaine (2014) used the same method for recreating the former surface and for calculating the volume of the rills. Their method consists of the following main steps: digitize the rill edge in 3D using stereoscopic glasses, create a 3D-polygon from the 3D edge line in GIS software, transform the 3D-polygon to raster format, and calculate the height difference between this raster and the rill DSM.
Our study also deals with volumetric assessments of erosion objects using very high-resolution DSMs, when no knowledge about the former surface is available. A Python tool implemented in ArcGIS software was developed for automatic volume calculations. The inputs to the tool are the DSM of a rill (or another object) and a polygon created around the edge of the rill. The method was verified on an artificially created rill, where the volume result calculated using the tool was compared with the real volume. In addition, the influence of the position and the area of the rill edge polygon were tested on two erosion rills. Eight different edge polygons were created around the rills, and a volume calculation was performed for each polygon.

Data and methods
Close range stereo-photogrammetryrill R1 The first testing rill R1 was artificially created in fallow field conditions. The field is situated on the experimental catchment of the Býkovický stream in the Central Bohemian Region near the town of Benešov (Báčová & Krása, 2016). The rill was manually incised on the surface and shaped by a water inflow in order to approach similar geometrical parameters to a nature-formed rill.
A metal reference frame was designed for calculating the exterior orientation. The frame was provided by 40 equally distributed ground control points (GCPs). The frame size was set to 1.6 × 1.1 m. The reference frame was kept in horizontal position by four screws throughout the imaging. All photographs were taken with the same camera (Sony NEX-5N), the same settings (aperture f-8, exposure time 1/40 and less, focal length 18 mm), and the same lighting throughout the imaging. The camera was oriented vertically, with a height above ground of approximately 1.3 m and an image parallax of 0.3 m. Each image covered the whole test area and all reference points.
PhotoModeler Scanner (PhotoModeler Technologies, Canada, Vancouver) software was used for generating the DSM of the rill. First, multi-sheet camera calibration was executed, and then a dense surface point mesh was created, with the following steps: • select and add an available pair of stereo photographs, and match it automatically with an appropriate camera from the calibration library, • use the auto-marking tool to detect GCPs, process the project, and verify that a good result has been obtained by checking the largest residual and camera positions, • define the scale and the rotation definition of the coordinate system of the project, • define the target area used to generate a point mesh, • set the dense surface creating parameterse.g. the example sampling interval (1 mm in our project), the depth range, the matching region radius, and other advanced parameters, • execute the automatic dense surface point mesh, • export the point mesh to ASCII format for ArcGIS processing.

Structure from motionrill R2
The SfM method (Szeliski, 2010;Ullman, 1979) enables simultaneous computation of relative projection geometry and a set of 3D points, using only corresponding image features occurring in a series of overlapping photographs captured by a camera moving around the scene (Verhoeven, Doneus, Briese, & Vermeulen, 2012). This method was used for generating a DSM of the second testing rill R2, in this case a natural rill caused by a heavy rainstorm. The study area is located close to the first one, in the central part of the Czech Republic, near the town of Benešov, in the village of Postupice. The location was impacted by an intensive rainfall event in the spring of 2013 with a total precipitation amount of 17.3 mm and a maximal rainfall intensity of 37 mm h −1 . An extreme number of ephemeral rills and gullies originated during this event, and soil from the field was transported to the neighboring pond. The rill selected for the test reached 60 m in length and the average width was 0.3 m. Above the selected testing part of the rill, the metal frame with GCPs was placed for referencing of the model. Photos were again taken from ground, using a Sony NEX-5N camera with the same settings (aperture f-8, exposure time 1/40 and less, focal length 18 mm).
The DSM was created in Agisoft PhotoScan Professional (Agisoft LLC, Russia, St. Petersburg) software in the following steps: • align the photos in the highest accuracy level, with a key point limit of 40,000 and a tie point limit of 4000, • reference a sparse cloud by placing markers on GCPs positions and filling in the coordinates of the markers, • build a dense cloud with ultra-high quality and disabled depth filtering, • generate DSM.

Rill volume calculation method
A Python tool was developed for rill (gully) volume calculation in ArcGIS (Esri, California, Redlands) software. The script utilizes the "arcpy" Python framework and enables automated calculation of the rill volume. The essential input data are the DSM of the rill in a raster format and the rill edge polygons. As more polygons can be processed within a single execution an attribute field holding the polygon identifier must be chosen.
The rill edge polygon can be created manually over underlying orthophoto, or it can be visible from DSM or contours, if the rill is sharply cut and the sides of the rill are steep. Automatic rill edge detection is also possible, for example by using the Canny edge detector (Canny, 1986), which was applied by Eltner et al. (2014) in her approach for the automatic rill extraction. The polygon edge is then used to estimate the surface prior to erosion event. The polygon needs to be placed carefully not to reach inside the assessed rill but it does not need to be placed exactly on the edge of the rill. A mild buffer outside the rill is optimal.
The calculation of the volume includes following steps ( Figure 1): First, the (planar) input polygons are assigned an elevation from the DSM used. Optional densification of the polygon vertices prior to processing can be chosen in the tool settings as the elevation is assigned to the vertices of the polygon (see Figure 2) where the value represents the maximum distance between two vertices on an edge. The vertices with elevation values are then used to generate a triangulated irregular network (TIN) surface model which replaces the unknown surface definition before the storm eventthe "original surface". The TIN surface is then converted to a raster model in order to have an identical format with the rill DSM. The two raster models (the original surface and the DSM of the rill) are then subtracted and the differential raster is generated.
The DoD is limited by a threshold value (pre-set to zero) to exclude areas that are higher than the original surface. Those negative or, more precisely, opposite values can occur near the edge of the raster (Figure 3). These opposite values originate due to errors during the creation of the edge polygon, mainly when the vertices of the polygon are placed on local depressions ( Figure  4), or if there is an insufficient number of vertices, so that the surface around the rill is not sufficiently representative. Another possible cause of this phenomenon is that there is some object that is placed higher than the bottom of the rill, such as stick or stone situated at the edge. In these places, the original surface is situated lower than the rill DSM, and the subtraction result has a value with an opposite sign to the numbers in the rest of the rill. These border opposite values are not desired for the calculation of the final volume and are therefore removed by the threshold operation.
The difference threshold can be set to non-zero values to account for homogenous compaction of the soil, sheet soil erosion, or uncertainty in the DSM data. The threshold limited DoD serves as the rill extent definition (mask) for further processing. The mask is used as a zone definition for calculation of maximum and average depth within the rill. The volume of the rill is calculated by multiplying the mean depth by the extent area. Optionally, the rill extent mask can be saved as polygon with some additional simplification settings provided (see Figure 2). Intermediate results of the whole procedure can be maintained for later error check by setting so in the dialog window prior to the tool's execution.

Validation of the method
The rill volume calculation method itself was validated in two parts. First, the method was tested on an artificial rill, and the results were compared with the real results obtained with the use of knowledge about the original surface from the period before the rill was formed.
The second part involved testing the influence of the position and the shape of the boundary polygons, which are important inputs into the tool.
Calculating the volume of the artificial rill To verify the rill volume calculation method, an artificial rill, referred to as R1, was formed manually in an agricultural field. The surface before the formation of the rill was photographed with the reference frame to obtain the original surface ( Figure 5). In the next step, the rill was dug out and the reference frame was situated in the same place, using screws which were fixed into the ground. Then, the rill was photographed again for use in creating the rill model ( Figure 5).
A DSM of both surfaces (before and after the formation of the rill) was generated using PhotoModeler Scanner software, as described in section " Test rill R1". The volume of the rill was calculated by subtracting the DSMs. For a comparison with the results calculated using the tool, subtraction was performed only for a specific area, i.e. for the polygons presented in the following section.
Influence of the position and the shape of the rill edge polygons The influence of the position and the shape of the rill edge polygons were tested on both rills. Eight polygons with different shapes and positions were manually created for each rill. Their positions should represent the most frequent options for locating the polygons, e.g. accurate location right behind the edge of the rill (P1 in rill R1), less precise location including fewer vertices (P2-P7 in rill R1), and also a distant position (several centimeters from the edge of the rill) of the polygon (P8 in rill R1).
The volume of the rills was calculated using the Python tool, and then cross sections passing through the original surfaces, formed above the rills, were performed.

Test rill R1
The DSM of the artificial rill R1 is shown in Figure 6. Eight polygons are presented, with their different shapes and positions, especially on the detailed part of the edge.
The so-called original surfaces, created within the tool based on single polygons, are shown in Figure 7. The cross sections are drawn with the same type of lines as the related polygons in previous figure. The position of the cross section line is marked by the blue dashed line in Figure 6.
The cross section includes the original natural undisturbed surface as it was before the rill was excavated. The shape of the original natural surface appears to be more curved than the estimated surfaces. However, the estimated original surfaces will always have a planar or slightly curved shape. This is because they are created as a TIN between the vertices of the polygons, which include the surface elevations on the edge of the rill, but not above the rill. It is not possible to estimate the elevations and the roughness of the natural original surface without information about the situation before the erosion event. Figure 7 contains a detail of the left border side of the cross section of R1. It can be deduced from the picture that the positioning of the polygons, or of the vertices of the polygons, influences the design of the calculated original surfaces. For example, the surfaces created from polygons P1, P2, P3, P4, and P6 are placed lower than the edge of the rill. This is because the vertices of the polygons were located in local depressions in this place. The surface created from the elevations in these depressions is consequently also lower than the edge of the rill. However, this effect, which can of course be observed in several places, is only local. It is illustrated on the right side of the rill (Figure 7), where all the surfaces are located correctly on the same level as the edge of the rill.
Another inaccuracy visible in Figure 7 is that the surface P3 is not connected with the edge of the rill. The reason is that the polygon outline segments is a straight line (shortest connection of two adjacent  vertices) and does not follow the curvature of the rill's edge. That is the reason why the polygon edge runs above the rill even if the two points at both ends of the outline segment are placed correctly on the surface. This is shown in detail in Figure 6, where the line P3 is situated above the deeper part of the DSM. In this case, one more vertex would be needed for the correction of this inaccuracy.
Unlike P3, most of surfaces in Figure 7 are drawn across the edge. This means that the polygons were  placed too far from the rill. However, this error does not influence the results, because after the DSMs have been subtracted, minor results with opposite signs are eliminated from the result for the volume (see the detailed explanation in section "Rill volume calculation method").
The results for the rill volumes, calculated via the eight different polygons, are presented in the graph (Figure 8). In numerical form, the percentage difference between the results is not more than 5%, which is good agreement. However, fact that the area of the polygon influences the volume should not be overlooked. In this case, the dependence between the volume of the rill and the polygon area is almost linear (Figure 8). This effect is also evident in the cross section of the rill (Figure 7), where the original surfaces created from the largest polygons P7 and P8 form the best cover of the rill. The larger the polygon is, the better the cover of the rill, which leads to the bigger volume of the rill. However, Figure 6 shows that polygon P8 is located far from the edge of the rill. If, for example, there was another rill or a local depression, the result would be overestimated.
The second part of the verification on the artificial rill was a comparison between the volume results and the real volume, which was calculated from the real natural original surface of the rill, obtained before the rill was manually formed.
The real rill volume was calculated by subtracting the DSM of the rill from the DSM of the surface before the rill was excavated. However, because the area of the polygon influences the results, the volume had to be calculated again through the area of eight tested polygons P1-P8 for comparability of the results.
The comparison in Table 1 shows that the volume of the rill, calculated using the tool and without knowledge of the natural original surface, is 10-15% smaller than the real volume, depending on the input polygon. This difference is caused mainly by the fact that the real natural surface was more curved, as is shown in Figure 7. However, for the calculation using the tool, no information was available about the surface before erosion. Bearing this in mind, the results are in good agreement.
In this case, the biggest volumes seem to be the most precise results. However, if the natural original surface were lower, smaller volumes would be closer to the real result. It is therefore appropriate to consider the results to be only a quite accurate estimate.

Test rill R2
Erosion rill R2 is illustrated in Figure 1. Eight different polygons were again created around the edge of the rill. Their areas are shown in Table 2, together with the volume results.
The volume results calculated using the tool for eight different polygons again deviate from each other only slightly. The variation between the results is ±4 %, which is even less than for the rill R1.
The dependence between the area of the polygon and the rill volume is also apparent in this case. This may again be because some polygons are not large enough to cover the entire rill, so the resulting volume is underestimated. Conversely, some polygons may be too large and may include depressions behind the rill edge, and the volume is therefore overestimated.
For this rill, the real volume could not be calculated. The rill was not created artificially, but was formed during a natural erosion event. No previous data were availablewhich is indeed the usual situation.

Discussion
The presented methodology is a promising approach to make fast, non-invasive assessments of erosion damage at a field scale. Non-contact methods, such as photogrammetry and laser scanning, are currently the most efficient methods for obtaining precise detailed DSMs. The Python tool takes advantage of DSMs created using these methods for calculating the volume of erosion rills and other erosion objects. The calculations within the tool are fully automatic, apart from the creation of one of the inputsthe rill edge polygon. This step is closely linked to the accuracy of rill edge detection. Peter et al. (2014) concluded that the digitized rill edge is based on the individual human perception of the person who digitizes it, and in some cases this human perception can influence the accuracy of the overall results. Some authors consider the determination of the geometry of the rill, especially its width, to be a crucial issue (Casalí et al., 2006;Evans & Lindsay, 2010;Peter et al., 2014). In this study, the influence of the position and the shape of the polygons have been tested on two rills. The results showed, as expected, that the polygon, and especially its area, has a relatively big influence on the volume result: the greater the area of the polygon, the greater the rill volume. However, the results calculated on the basis of eight different polygons varied only within a range of 5%, which is a reasonably good result.
The resolution and the accuracy of the input data also influence the resulting volume. In photogrammetry, these parameters are influenced by many factors (Eltner et al., 2016): scale/distance, type of camera and its calibration, image network geometry, image-matching performance, surface texture and lighting conditions, and GCP characteristics. The source of errors is variable; however, it appears that the distance between sensor and object is a particularly important factor (Smith & Vericat, 2015). The comparison of several studies carried out by Eltner et al. (2016) indicates that with an increase in distance, the errors in models increases. However, according to Eltner et al. (2016), there is no linear trend detectable. In contrast, the quality of spatial resolution decreases with increasing sensor to surface distance almost linearly (Báčová, 2018;Mesas-Carrascosa, Notario García, Meroño de Larriva, & García-Ferrer, 2016). Similarly, as it was found in author's previous study (Báčová, 2018), the volume of assessed object is influenced nearly directly proportional to the distancethe increasing distance causes decreasing of the resulting volume.
The raster format of DSM is used in the tool, because it is favorable for both computational and time requirements. However, application of so-called 2.5D approach may result in underestimation of the rill volume as in the model every xy-coordinate correspond to only one z-coordinate, which means that undercut walls are not represented in the model. However, several authors have compared cross sections and volumes of ephemeral rills or gullies calculated from 3D and 2.5D surface models, and considering the volume, both approaches showed only small differences in results of 1.4% in maximum (Di Stefano, Ferro, Palmeri, Pampalone, & Agnello, 2016;Frankl et al., 2015). Errors though depend on individual conditions, more significant overhangs and undercuts are typical for bank gullies in sedimentary deposits and generally for gully headcuts. Moreover, the morphology of the steep sidewalls cannot be mapped detailed from an airborne perspective (e.g. using aerial and UAV photographs) (Giménez et al., 2009, Stöcker et al., 2015, which can also lead to a decision to use reduced 2.5D approach.
Besides the single rills, the tool can be used to estimate the volume of eroded material from larger areas, for example from agricultural fields or catchment areas, as well as extensive gullies. However, bigger the extent of the task is, bigger the uncertainties of the results are. The presumption that the original terrain can be estimated by "stretching a rubber sheet" within the outline polygon may be no longer valid. However, processing of greater areas may be split into a number of smaller tasks where the polygon outlines are carefully placed on undisturbed locations. Soil loss from larger areas may be then extrapolated from the results of these tasks, on the assumption that the tasks where located to the representative parts of the field.
For the field-scale areas, terrestrial methods for monitoring an erosion objects are not efficient. In recent times, UAV has become the most widely used and the most suitable instrument for field-scale monitoring. It can be used effectively and quite promptly for photographing erosion damage of fields from a relatively low height level, so it is still possible to detect rills and other erosion objects. As it was mentioned above, the resolution and the accuracy of DSMs generated from UAV photos are lower than from the terrestrial observations. The maximum flight altitude 15-20 m can be recommended for volume calculations, where, according to author's tests (Báčová, 2018), the resulting volume differ from terrestrial photogrammetry in about 5-6%. Of course, it depends also on other error sources. However, if there is no requirement for a high-precision result, such a data are sufficient for estimating the volume of eroded material from rills over larger areas. In addition, terrestrial measurements can be applied for checking, verification, or to fill up data gaps.
It is necessary to have in mind that the method described here provides a volumetric assessment of erosion objects such as rills, gullies, or caverns, but that sheet erosion or deposition of erosion material cannot be quantified by this method. For these erosion effects, long-term monitoring and knowledge of the prior ground surface are needed. The method is useful for cases when prior information is not available, which is a common situation, and it estimates only the volume of in-depth erosion objects.

Conclusions
This paper has presented a volumetric assessment of soil erosion objects using a method which was developed and implemented as a tool for ArcGIS. The assessment does not require information about the original natural surface prior to the erosion. The calculation method makes use of high-resolution DSMs as the input data. The calculation within the tool is fully automatic, except for the creation of a second inputthe rill edge polygon.
The influence of the accuracy of the digitized rill edge polygon on the resulting volume of the rill was tested on two rills. Eight polygons with different shapes were created for each rill. As expected, the influence of polygon area was proved and also the influence of the position of the vertex of each individual polygon. All inaccuracies of the polygons and their influence on the results have been presented in the results section. Generally, it can be concluded that the vertices should be placed carefully behind the edge of the rill in order to cover the entire rill. However, it is not appropriate to put the polygon too far from the edge, as this can cause additional close surface features to be included in the volume of the rill. However, the volume, calculated via the different polygons, varied within a range of 5% in both rills. This result can be considered as very good, especially on this detailed scale.
An important question for future research is to find how to detect rill edge polygons automatically. However, this task remains challenging due to difficulties in the positioning of the polygons, as we have discussed. It is necessary to find some compromise between accuracy of the results and time requirements.
The method was validated using an artificial rill. Before the rill was manually formed, the surface was photographed to create a DSM of the original natural surface. Due to these steps, it was possible to calculate the real volume of the rill, and to compare it with the results of the method we developed. The comparison showed that result from the tool was underestimated by 10-15%. This underestimation can be explained by the fact that the real natural surface was curved, whereas the tool is able to create only flat or slightly curved cover from the vertices of the polygon. However, this result is very sufficient for a calculation, or rather, for an estimate of the volume of the rill.
The main advantage of the method is that a fast and effective calculation of the volume of an erosion object can be made without knowledge about the surface before the erosion event. Naturally, it is always better to have available time-series data for long-term monitoring of soil erosion. Unfortunately, however, there is commonly no data until an erosion event comes along. The method is useful for an estimate of soil loss, even if sheet erosion is not included.
The tool implementing the method can be used for making a detailed estimate of single rills, where it is appropriate to obtain the input DSM using close range photogrammetry or TLS. A volumetric assessment of larger areas can also be made by this method, typically for agricultural fields, catchments, or extensive gullies. A DSM for these calculations can be obtained by UAV photogrammetry or by aerial laser scanning. Soil erosion is not the only field where the method can be used. For example, volume calculations of mines, excavations, potholes, and other objects can also be made. The tool is made available as an appendix of the manuscript.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This study was supported by the Ministry of Agriculture of the Czech Republic (project No. QJ1330118 Using Remote Sensing for Monitoring of Soil Degradation by Erosion and Erosion Effects and No. QK1720289 Development of automated tools for optimizing monitoring erosion of agricultural land using remote sensing methods), and by the Czech Technical University in Prague (grant No. SGS17/ 173/OHK1/3T/11 Experimental research of erosion and transport processes in agricultural landscapes).