ANALYSIS OF GULLY EROSION IN A CATCHMENT AREA IN OLIVE GROVES USING UAS PHOTOGRAMMETRY TECHNIQUES

An approach based on images captured by means of an Unmanned Aerial System (UAS) and Structure from Motion – Multi Video Stereo (SfM-MVS) photogrammetry is presented for the study of the gully erosion in a catchment area of about 16 ha. The study area is located in the province of Jaén (SW Spain) where the main land use is the olive groves. Three UAS flights have been made in April 2019, May 2019 and February 2020 with an average GSD of 2-3 cm. The image processing has been carried out using some field surveyed GNSS ground control points and the RTK positioning of the UAV with errors lower than 0.03 m and 0.04 m in XY and Z, respectively. Then DSMs and orthophotographs were obtained with a resolution of 0.10 m and 0.05 m, respectively, and DSM of differences (DoDs) were calculated, with an uncertainty of about +0.15 m. Finally, the DoDs were analysed in the GIS, in order to calculate the height differences and volumes between the three flights. The analysis has showed the areas affected by the gully erosion processes that correspond mainly to ground descents due to significant channel erosion and mass movements at the steep sidewalls. Considering the balance between depletion and deposition processes, the average height differences is about -0.03 m and the volume is about -42 m3, that is the depletion predominates. The main changes have been detected in the second period (May 2019-February 2020), while in the first one (April 2019-May 2019) they are practically insignificant and limited to small areas.


INTRODUCTION
Globally, soil erosion is one of the most disturbing phenomena of environmental degradation (Borrelli et al., 2013) that might be increased significantly in coming decades due to the acceleration processes caused by global warming (Yang et al., 2003). Thus, many studies on soil water erosion have been addressed at multiple spatial and temporal scales. In the plot scale, in which processes of laminar and rill erosion predominate, empirical measurements (Merrit et al., 2003) or different variants of the RUSLE model (Gómez et al., 2003) have been applied. Some authors (Poesen et al., 2003) demonstrate that concentrated flow processes such as gully erosion could explain between 50% and 90% of total erosion at basin scales. Therefore, gully erosion is now one of the most important and studied soil degradation processes (Castillo and Gómez, 2016), ranging from small ephemeral gullies to large permanent gullies (Poesen et al., 2003).
Images and orthoimages can be used for surface characterization, such as the estimation of lengths, widths and densities of the gully systems, and the analysis of their evolution over time (Hayas et al., 2017). In addition, image based methods permit the generation of digital elevation models (DEMs) for the estimation of gully depths and volumes as well as their variations. The availability of models at different epochs leads to multitemporal analysis of gully systems by means of the calculation of differential models (DEMs of Differences or DoDs). All of these image-based methods are often combined with point capture techniques such as LiDAR (Koci et al., 2017;Fernández et al., 2020), GNSS (Brasington et al, 2003;Rumsby et al., 2008) or conventional surveying (Castillo et al., 2012;Hayas et al., 2017).
Specifically, UAS are well suited for very high resolution and precise surveys in areas of about 0.01 to 100 km 2 . UAVs are very suitable for intermediate scales between terrestrial techniques (GNSS, close range photogrammetry and TLS) and aerial or space surveys (conventional aerial photogrammetry, LiDAR and VHR satellite imagery) but keeping low or moderate costs and allowing high temporal resolution studies. Most of the current studies using multicopter UAS are of centimeter resolution (d'Oleire et al., 2012;Stocker et al., 2015;Koci et al., 2017;Yang et al., 2019). After the image processing by means of SfM-MVS techniques, digital elevation models (both DSMs and DTMs) of high resolution and precision are obtained. Based on these DEM, morphometric measurements and volumetric estimations can also be performed (d'Oleire et al., 2012;Stocker et al., 2015;Yang et al., 2019).
A crucial aspect is the accuracy of the data collected with the different techniques. In some cases the accuracy can be estimated from errors (usually the RMS) of the image orientation at ground control and/or check points (Stocker et al., 2015;Fernández et al., 2016;Wang et al., 2016;Koci et al., 2017). Other analyses are based on the comparison of the height of profile points or DEMs, with respect to points measured with a more accurate method such as the TS/GNSS/TLS (Stocker et al., 2015;Koci et al., 2017). Finally, methods based on the comparison of repeated measurements in samples of points of the same surface can be found (Fernández et al., 2016). This paper aims to describe a methodology for the identification and quantification of gully erosion in a period near to one year based on UAS photogrammetry. The study area is an active gully in an olive grove of the province of Jaén (southern Spain). For this purpose, three UAS flights were planned and executed under similar conditions. Images were processed using the GCPs and the RTK positioning of the UAV. Then, for each flight, the corresponding DSM and orthophotographs were obtained and the DSMs of differences (DoDs) were also calculated. The DoDs were analysed in a GIS environment, in order to calculate the height differences and volumes, as well as their corresponding rates along the period considered. Finally, this methodology has been validated by means of the estimation of errors and uncertainties from the obtained models.

Study area
The study area, with an extension of approximately 16 ha, is located in the western part of the province of Jaén (Andalusia, Spain), at a distance of about 25 km from the province capital ( Figure 1). It has an altitude of between 430 and 465 m and an average slope of 8.8°. It is located within the natural region of the Eastern Guadalquivir river basin. From the geological point of view, this basin is made up of the Guadalquivir Units (Pérez-Varela et al., 2017), a set of materials of diverse lithology, tectonically intercalated in loamy-clay sediments from the Miocene age ( Figure 1). The predominant lithologies are Triassic lutites, evaporites and carbonates, as well as Cretaceous-Paleogene marls and clays. Specifically, in the study area, Triassic lutites and sandstones with subsidiary amounts of carbonates and gypsum outcrop, as well as alluvial and colluvial deposits of Quaternary ( Figure 1).
The area corresponds to an active gully stretch in a catchment area affected by an intense erosion, both laminar and gully, in addition to other superficial processes such as landslides (Fernández et al., 2016;2020;Carpena et al., 2017). Some sections of the gully area studied, which sometimes affect rural roads and paths are shown in Figure 2.

Figures 2a and 2b
show clearly the V-shape at a tributary gully of the main gully and soil slabs falls at the steep sidewalls, respectively. Moreover, upstream erosion is affecting to a rural path, which can cause its collapse in near future as it can be seen in the photograph sequence from January 2016 to 2018 in Figure 2c.

Materials
The methodology is based on UAS photogrammetry techniques. The UAS is a DJI Phantom 4 with a RTK module integrated that provides positioning of centimeter accuracy that reduces dramatically the number of surveyed ground control points for the photogrammetric orientation. This also allows improved flight security with a flight range up to 20 minutes and capture of information for post processing kinematic (PPK). The camera is a DJI FC6310R (20 mpx and 0.0024 mm pixel size) equipped with a wide angle 8.8 mm lens.
Besides the UAS, with its own GNSS for vehicle positioning, other GNSS equipment (LEICA SYSTEM 1200+) has been used for ground control and check points (GCP/CHK) measurement.

Methodology
The methodology is based on UAV photogrammetry techniques that has been used in previous works of the research group (Fernández et al., 2016;Cardenal et al., 2019), although adapted to the erosion studies (Fernández et al., 2020). It can be summarized in the following steps: 1. Image acquisition and field work. 2. Image processing and orientation. 3. Generation of DSMs and orthophotographs. 4. Delimitation of gully areas and estimation of horizontal displacements. 5. Calculation of DSMs of differences (DoDs). 6. Estimation of height differences and volumes between models.
2.3.1 Image capture and field work. For this study, three UAS flights of very high resolution have been made on 03-April-2019, 14-May-2019 and 13-February-2020 (Table 1).
The flight planning was made once recognized the terrain, with the DJI desktop flight planning software (DJI GS RTK). Given the dimensions of the study area, approximately 850 m x 200 m, a GSD of 2-3 cm was selected in order to optimize the image resolution, the number of images and the flight time. This implied a flying height above terrain between 90-120 m, within the limits of Spanish regulations on the use of Remotely Piloted Aircraft Systems (RPAS). The photogrammetric projects consisted in a block of vertical images organized in three strips following the outline of the gully with end and side laps of 80% and 70% respectively. All flights were similarly planned. (Table  1 and 2, and Figure 3). For the first flight (April 2019) a set of 30 ground control (GCPs) and check points (CHK) was measured by means of differential GNSS with centimeter accuracy to check over the RTK image orientation proccess. The ground point network had a well-distributed pattern, all around the unstable area but also some of them inside the area following conventional distribution of GCP networks for aerial triangulation (Kraus, 2007). The GCP/CHK were artificial circular shaped targets (printed on PVC foam board) and targets sprayed with reflective paint on the ground surface and a cardboard template ( Figure 3). Some additional well-defined points were measured in unequivocally identifiable features, such as agricultural buildings and shed roof corners. After validating the block orientation on the first flight, some additional points were measured in the images to be used as second order ground control/check points for the May-2019 and April-2020 flights. Red triangles indicate the location of the ground control points. This 4-GCP network and the 3-strips flight pattern are the same for all flights. Circular targets in PVC foam board (b) and sprayed on ground with reflective green paint (c) used as GCP and CHK points in the first flight.

Image processing and orientation.
The images were processed and oriented by means of SfM/MVS techniques with Agisoft Metashape. The first flight of April 2019 is the reference flight. This flight was oriented with both GCP and CHK points in order to establish the camera station RTK accuracy and a second order GCP/CHK point network valid for the next flights. The RTK camera positions were corrected with data from a GNSS base station set up in the area. Since the UAS camera was a no metric camera, a self-calibration was performed during the block adjustment. In order to avoid errors propagation, mainly in the vertical component and the focal length, some GCPs were necessary. Since there were 30 fieldsurveyed points, it was possible to explore several GCP networks. Finally, a network with four GCP located at the ends of the strips was selected ( Figure 3). The remaining points were considered as CHK points. This implies a cost saving technique, since conventional aerotriangulation methods require a higher number of GCPs (Kraus, 2007).
The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Table 2 shows the orientation errors of the reference flight. These errors, estimated at the CHK points, were in the order of the image GSD (2.5 cm). Table 2 also shows the errors in camera location (RTK cam) after the block adjustment. After validating the orientation of this reference flight, some additional well-defined and stable points distributed throughout the block were measured in the images. These new points were second order points transferred to the next flights. Therefore, a common reference system (CRS) is kept for next flights (May 2019 and February 2020) and besides there is no necessity of additional field surveyed GCP/CHK for those flights.
These additional transferred points were used for orientation of the further flights, which were also RTK surveys. The GCP/CHK networks for these flights were the same as in the reference first flight, i.e. four GCP at the ends of the strips and the rest of the points well-distributed all around the area as CHK points. The number of GCPs and CHK of each flight as well as the results of the alignment process are included in the

Generation of DSMs and orthophotographs.
After image block orientation of flights, digital surface models (DSMs) were generated from the dense point clouds. DSMs resolution is 0.10 m, about three-four times the GSD of images. Although in conventional photogrammetry the resolution of models are several times the GSD of the images, the new global dense matching techniques used by current SfM-MVS software allow the reduction of this proportion, taking into account also the dimensions of the area and the characteristics of the flights. Besides, following previous works (Fernández et al., 2016;2020), the uncertainties of the DSMs were established in twice or three times the RMS errors estimated on GCP/CHK, that leads in this case to a value of about 0.10 m.
As in previous studies (Fernández et al., 2016;2020), in this paper the erosion processes was monitored with DSMs instead of DTMs because of the study area has a high density of vegetation in some sectors with grass, scrubs and bushes. Thus, using conventional tools for point clouds classification and filtering could only just remove partially the vegetation.
Besides, stereo-model edition using photogrammetric workstations was not performed because it did not ensure good results and it would have been time consuming.
Next, orthoimages for each campaign were generated with a resolution of 0.05 m. Given the RMSxy error of the orientation process and the image resolution, the uncertainty for horizontal measurements was established in 0.05 m. Finally, both products were exported as raster files to be incorporated into GIS analysis. Figure 4 shows the DSM of the first survey and the orthoimages of all surveys.

Delimitation of gully areas and estimation of horizontal displacements.
After visual interpretation of orthoimages, the gully banks were delimited and digitised, using the corresponding tools of the GIS software (QGIS 3.0). Then, the resulting vector line files (in shp format) were converted in point files with a spacing of about 1.6 m (5391 points). These points were used to calculate displacements between the gully bank lines corresponding to the different surveys. The average displacements were also calculated. 2.3.5. Calculation of DSMs of differences (DoDs). The models have been calculated from the DSMs, which has allowed the detection of the areas that undergo vertical displacements of the ground surface between successive dates. Displacements may be negative or positive, depending on whether each model compared with a reference model lies below or above it. This allows the identification of areas of ground descent (mass depletion or erosion) or ascent (mass deposition or accumulation), respectively. The vertical uncertainties of the DoDs are estimated as follows (Brasington et al., 2000): Unc. DoD YEAR1-YEAR2 = (Unc. DSM YEAR1 2 + Unc. DSM YEAR2 2 ) 0.5 The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII-B2-2020, 2020 XXIV ISPRS Congress (2020 edition) Taking into account that the uncertainties of the DSMs are estimated in 0.10 m, the uncertainties for DoDs can be established in about 0.15 m, in the same order of previous works (Fernández et al., 2016). being ND, the digital number in the corresponding RGB bands.

Estimation of height differences and volumes between models. As mentioned, the GIS analysis of
The GLI is obtained for the orthoimages corresponding to the different dates. A threshold of 0.05 is applied to separate the areas covered and not covered by vegetation (index value higher of lower than the threshold, respectively). GLI images are combined to obtain a filter image of unconsidered areas.
After filtering the orthoimages to discard uncertainty and vegetation areas, we can calculate the regions affected by gully erosion. The calculation of the average values from the DoD (height differences) in the gully areas allows us to analyze the balance between depletion and deposition processes, estimating the general losses or gains of soil material. If the average balance is negative, the depletion processes predominate, and if it is positive, the deposition processes do. The average values of negative (depletion) and positive (deposition) height differences were also calculated, considering in each case the corresponding areas where height differences were higher than 0.15 m in absolute terms. The annual rates of depletion, deposition and balance height differences can be estimated by dividing the corresponding values by the time interval between models.
Finally, the volumetric calculations are carried out, estimating the volume balance of soil material (losses if the balance is negative or gains if the balance is positive) in a given area. Same as height differences, the depletion, deposition and balance volumes were estimated in the areas where height differences were higher than 0.15 m in absolute terms. And in the same way, the volume rates were also calculated by dividing those values by the time interval. Finally, the volume rates were transformed to mass rates by surface and time (t/ha*year), dividing by the area in ha and considering an average soil density of 1.5 t/m 3 .

Areas affected and horizontal displacements
The areas affected by the erosion processes are shown in Table  3. First, the gully area (about 1.5 ha) undergoes an increase of 0.2% in the first period (apr19-may19) and near 3% in the second (may19-feb20) and the whole period (apr19-feb20). The perimeter (about 3.4 km) increases 1.1% in the first period, 2.8% in the second one and almost 4% in the whole period.
Within the gully area, the area covered by vegetation (unconsidered area) represents 88.6%. Meanwhile the uncertainty area reach 9.7% for the first period, 4.9% for the second period and 4.8% for the whole period. Thus, the area with significant changes in the ground surface only represents 1.7% in the first period, 6.5% in the second one and 6.6% in the whole period. Within of this change area, the depletion area represents 69% in the first period, 52% in the second period and 57% in the whole period. Meanwhile the deposition area represents 31%, 48% and 43%, respectively.  Table 3. Areas affected by the processes in the gully area. Units are in m (perimeter) and m 2 (area). Table 4 shows the horizontal displacements of the gully bank lines. The average is higher in the first period (0.024 m) than in the second (0.117 m) and the whole period (0.136 m). In the same way, there is a great number of points with larger displacements in the second period than in the first one. Thus, the number of points with displacements larger than 1 m and 0.1 m are 113 and 1278, respectively, in the second period. Meanwhile, these are 11 and 330 in the first period.  Table 4. Horizontal displacements of gully bank lines.

Height differences and volumes
The height differences are shown in Table 5 Table 6. Volumes in the gully area.

DISCUSSION
The RMS errors obtained in the alignment process of the UAS flights are lower than 0.04 m in XY and 0.05 m in Z for the CHK points (Table 2). Taking the values of CHK as reference the uncertainty can be established in 0.05 m in XY and 0.10 m in Z (DSMs). These values are very similar as those estimated in previous works (Fernández et al., 2016(Fernández et al., , 2020Cardenal et al., 2019). The XY uncertainty is equal to the spatial resolution of orthoimages, so the measurements made on them can be considered as reliable. Regarding DSMs uncertainty, this allows to establish the DoDs uncertainty in +0.15 m. Thus, all vertical changes of the ground surface or height differences lower than this threshold, in absolute terms, are considered as not significant while those higher than the threshold are considered in the calculations In addition to the uncertainty area, we have had to take into account the area cover by vegetation in each survey, which have not been considered for the calculations within the gully either. This area represents almost 89% of the gully area (Table  3), so the calculations have been made in the remaining surface not covered by vegetation. The disregarding of these areas can produce a subestimation of the volumes involved.
From the calculations, we firstly can observe that the area and perimeter of the gully stretch increased about a 3-4% in the whole period analysed, although most of it occurred in the second period (Table 3). It leads to consider a moderate progress of the gully in the whole period, although greater in the second period regarding the first one. The displacements in points of the gully bank line (perimeter) also confirms this observation with a higher average displacements (0.12 m in the second period and 0.02 in the first) and a higher proportion of the point with larger displacements (in some points even greater than 1 m). In general, the increase of gully areas and perimeter and the displacement of the points of the gully banks informs about the main mechanism of erosion or depletion that is the gully walls retraction (Figure 5a). This retraction is larger in the north-western bank than in the south-eastern bank, exposed to the NW and with a greater vegetation cover.
All this agrees with that the second period presents a larger area with significant changes and a smaller uncertainty area, Moreover, most of the change area corresponds to depletion or erosion processes regarding the deposition, especially in the first period where the proportion becomes almost 70/30%. This proportion decreases to 52/48% for the second period and 57/43% for the whole period. However, it should be taken into account that change areas represents only 1.7% in the first period and about 6.5% in the second and the whole period. In addition, although the maximum and minimum values of height differences reach up to 5 m in some specific sector, the average values of depletion and accumulation heights are about 0.3 m in the first period and about 0.4 m in the second and the whole period. These moderate values of depletion and deposition as well as the limited extension of areas with significant changes lead to very low balance height differences (about -0.02 m). The conclusion is that the erosion processes are not so intense in the period studied, especially in the first one, although this is a very short interval (one month). If we compare these data with the height differences observed in the study area in previous studies (Fernández et al., 2020), the values are of the same order than those periods less active, such as 2005-2009 and 2013-2016. They are very far from the periods with higher height differences such 2009-2011 and 2011-2013 (between -0.6 and -1 m for the average of balance height differences). Besides, these average values are not very significant since they refer to areas that represent only 1.7-6.5% of the total.
The same may be observed in the analysis of volumes, where the values of depletion volume are in general higher in absolute terms than the deposition volumes that produces a negative (waste) volume. The volumes estimated are about 50 m 3 for depletion, 25 m 3 for deposition and -28 m 3 for the balance in the first period. This values increases up to 200 m 3 , 177 m 3 and -20 m 3 in the second period, respectively. In general, the volumes implied in this study are of one-two order of magnitude lower than the volumes estimated in previous years (Fernández et al., 2020) or even more if we consider the more active years 2009-2013) when volumes reach values up to 18000 m 3 for depletion and 15000 m 3 for the balance. The predominance of depletion is higher in the first period than in the second, in which the situation is more balanced. Thus, the erosion process occurs in a moderate way affecting especially to gully walls in the last years (Figure 5b), as was revealed in previous studies (Fernández et al., 2020). After 2011, the erosion mechanism is the lateral growth, unlike previously in which it was the vertical incision. However, the energy is not enough to transport the material over large distances and it is deposited in the foot of the gully walls or some meters downstream (Figure 5b and c). Taking into account the rates, for a better comparison with previous studies, the height differences range between 0.5-2.5 m/year for depletion and accumulation, being close to 0 for the balance. These values are comparable with the rates obtained for the previous years, even the most active (2009-2011 and 2011-2013), but given the limited proportion of the change areas and also the very shortness of the periods (especially the first one) this rates are not very significant. More valuable is the analysis of volume rates, in which the values range between 250-465 m 3 /year for depletion, 200-235 m 3 /year for deposition and 25-250 m 3 (negative) for the balance. These rates are in the lower range of the rates obtained for the less active periods such as 1990-1986, 2001-2005, 2005-2009 and 2013-2016; and far from the more active periods such as 2009-2011 and 2011-2013 when the rates reach values between 5000-10000 m 3 /year for depletion and balance (Fernández et al., 2020).
In fact, translating these data to mass rates, the values range between 8.5-15 t/ha*year for depletion, about 7 t/ha*year for deposition and 1-8 t/ha*year for the balance. These rates are also in the lower range of the rates found for this gully in the previous years and very much smaller than the values of the most active periods, when the rates reach values between 300-500 t/ha*year. Moreover, these rates are also lower than those obtained in other areas: 40 t/ha*year (Hayas et al., 2017); 160-430 t/ha*year (Martínez-Casasnovas et al., 2004); or 700 t/ha*year (Lane (2003).
In summary, the gully area studied has a low activity in the 1year period analysed, following the pattern of the last years (2013-2016) and some previous periods (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009). In these periods, same as in the analysed 2019-2020 period of this paper, the rainfalls were scarce and do not reach the values found in 1996-2001, 2009-2011 and 2011-2013. In these wet periods, the weekly rainfalls showed values higher than 100 mm and the monthly rainfalls around 250-300 mm. Specifically, in the period analysed, significant rainfalls about 60-80 mm were only reached in some months such as April, September, November and December 2019 (the first one related to the first period and the other three to the second period). This moderate to low activity was focused in the gully banks, with retraction through slab falls in the walls and material accumulation in their feet. In some cases, excavation of the lower part of the walls occurs (Figure 5d), but practically there is not excavation of the gully bottom. The deposition of material is located not only in the walls feet but in the gully bottom at a given distance from the origin, in relation to plane areas of presence of obstacles (vegetation, blocks, etc.).
All these details, such as the changes in the gully bank delimitation and in the microtopography can been observed due to the high resolution and precision of the products obtained by means of the UAS photogrammetry. In this sense, the suitability of the methodology applied has been proved.

CONCLUSIONS
In this paper, we have developed an approach to study the erosion processes in an active gully stretch of about 16 ha using UAS photogrammetry techniques. Data capture has been made in three UAS flights with an average GSD of 2-3 cm. The image processing has been carried out by means of SfM-MVS techniques using some field surveyed GNSS-based GCPs and the RTK positioning of the UAS. RMS errors obtained are lower than 0.04 m in XY and 0.05 m in Z for the CHK points.
Digital Surface Modes (DSMs) and orthophotographs have been obtained with a very high resolution of 0.1 m and 0.05 m, respectively. Considering the values of CHK, the uncertainty can be established in 0.05 m in XY and 0.10 m in Z (DSMs). DSM of differences (DoDs) were calculated, with an uncertainty of about ±0.15 m, and the bank gully lines were digitised on the orthophotographs. The GIS analysis of these highly accurate products has allowed the observation of some details that could not be done with other techniques.
Thus, first we can observe an increase of 3-4% in the gully area and perimeter, and displacements of bank line points with an average of 0.14 m. Moreover, a moderate depletion of some sectors of the ground surface corresponding to the gully banks is also detected, being the average of height differences of around 0.3-0.4 m in the change area that represents no more than the 6% of the gully area. The depleted volume is about 220 m 3 . Meanwhile, deposition is also observed in the lower part of the gully banks and in the gully bottom, being the average of heights differences about 0.3-0.4 m and the volumes of 175 m 3 . The balance of height differences (-0.02 m) and volumes (-40 m 3 ) is quite limited. All of these observations lead to established that the main mechanisms of the gully evolution is the retraction of the gully bank through the slab fall of the walls and the deposition of material in the walls feet and in the gully bottom after a given distance of transport. The energy is not enough to transport the materials through long distances and evacuate them out of the study area. The moderate rates of height differences and volumes (-0.04 m/year and -50 m 3 /year for the balance) are related to the scarce rainfalls of this period (April 2019-February 2020), regarding the higher activity registered in previous wet periods such as 2009-2013.
Future works will be lead to overcome the limitations of this study, especially in the improvement of the accuracy to reduce the uncertainty areas and the vegetation areas. In this sense, the use of other sensors such as laser scanners and multispectral cameras can be help in the filter and classification of these vegetated areas. Meanwhile, the better exploitation of GNSS-RTK positioning of cameras for direct orientation will allow to reduce the field work and increase the extension of study areas to make estimation of erosion in larger catchments.