sensors Article Relief Modeling in the Restoration of Extractive Activities Using Drone Imagery Erick Russell 1 , Joan-Cristian Padró 2 , Pau Montero 3 , Cristina Domingo-Marimon 4 and Vicenç Carabassa 1,3, * 1 Universitat Autònoma de Barcelona, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain 2 Departament de Geografía, Edifici B, Universitat Autònoma de Barcelona, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain 3 CREAF, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain 4 Grumets Research Group, CREAF, Edifici C, Universitat Autònoma de Barcelona, E08193 Bellaterra (Cerdanyola del Vallès), Catalonia, Spain * Correspondence: v.carabassa@creaf.uab.cat; Tel.: +34-93-581-3355 Abstract: In the field of mine engineering, a cross-section topographic survey is usually carried out to perform volumetric calculations of earth movement in order to restore areas affected by extractive activities. Nowadays, Remote Sensing and Geographical Information System (GIS) technologies make it possible to perform the same work by using indirect methods such as images obtained by photogrammetric flights. In this context, Unmanned Aerial Systems (UAS) are considered a very convenient option to develop mapping projects in short periods of time and to provide quality geospatial information such as Digital Elevation Models (DEM) and orthophotos of centimetric spatial resolution. In the present study, this approach has been applied in a gravel extraction area to obtain data for estimating the filling volume of material required for the restoration of the relief (DEM(r)). The estimation of the DEM(r) is later used to calculate a difference of height values (DEM(r)-DEM) that will serve as a variable in the basic operation of volume calculation. The novelty of the presented method is the simulation of a relief adapted to the surrounding morphology, including the derived channel network and the visibility impact, improving what would be a simple clogging. Likewise, the Citation: Russell, E.; Padró, J.-C.; generation of 3D models allows visualizing a new morphological structure of the relief. The proposed Montero, P.; Domingo-Marimon, C.; approach, based on GIS tools, allows analyzing water flow connectivity integration of the DEM(r) Carabassa, V. Relief Modeling in the with the environment and estimating potential landscape impacts from the main focuses of a visual Restoration of Extractive Activities basin, both of which are key aspects of restoration modeling that are not always properly addressed. Using Drone Imagery. Sensors 2023, 23, 2097. https://doi.org/10.3390/ Keywords: mines; drones; photogrammetry; interpolation; modeling; DEM; water flow; visual basin s23042097 Academic Editor: Francesca Antonucci 1. Introduction Received: 19 January 2023 Mining activities are very important for the provision of raw materials and energy [1]. Revised: 1 February 2023 However, this activity generates changes in the morphology and structure of the landscape, Accepted: 3 February 2023 altering the natural conditions of the area and directly and indirectly influencing the Published: 13 February 2023 conditions of its closest environment. It has environmental impacts, such as loss and pollution of soils, vegetation cover, erosion, unstable slopes, detachments of loose material, and alteration of the flow of groundwater and surface water, among others [2]. Copyright: © 2023 by the authors. These impacts are the reason why organizations such as the United Nations recom- Licensee MDPI, Basel, Switzerland. mend improving the management of environmental impacts and, thus, complying with This article is an open access article Sustainable Development Goal 12 of Responsible Production and Consumption [3]. In this distributed under the terms and sense, Spanish legislation requires companies to implement restoration plans that apply conditions of the Creative Commons corrective measures to prevent, reduce, and rehabilitate the land affected by the extractive Attribution (CC BY) license (https:// activity. These plans aim to return the land to a satisfactory state to the greatest degree creativecommons.org/licenses/by/ possible, guaranteeing a balance of the natural space before the exploitation operations 4.0/). Sensors 2023, 23, 2097. https://doi.org/10.3390/s23042097 https://www.mdpi.com/journal/sensors Sensors 2023, 23, 2097 2 of 15 (Royal Decree 975/2009). A restoration plan should involve the recovery of environmen- tal functions and services lost during the exploitation phase, such as the regulation of the hydrological cycle, the recycling of organic matter, and the production of biomass or biodiversity so that we can talk about ecological restoration [4]. Beyond the work of creating new soils or planting forest species, the restoration of these functions depends in great measure on the correct modeling of the terrain [5]. In this context, the calculation of the filling volume necessary for correctly modeling the terrain has been addressed by topographic surveys in the field (direct methods), which are used to learn the physical characteristics of the land surface. The topographic surveys result in vector entities of type point (height levels) or lines (contour lines) that are used to generate triangulation models, such as the Triangulated Irregular Network (TIN). This model presents a simple structure based on the connection of vertices of the nearest neighbors to form triangular networks, obtaining a representation of the morphology of the terrain [6]. TIN models are commonly used in engineering work to calculate volumes from cross-sections because it is an easy method to understand and implement. However, its results are usually not accurate as it assumes that the cross-sectional area varies linearly with the length [7]. Other interpolation methods, such as Inverse Distance Weighting (IDW) [8], splines [9], or kriging [10] are used for mineral resources. TIN, IDW, splines, and kriging interpolation tools [6,11,12] are available in common GIS software (MiraMon [13], QGIS [14], ArcMap [15]) toolbox obtaining reasonable results in this and previous studies [16,17]. Nowadays, there are other mechanisms and techniques that allow studying the changes generated in the territory due to mining activity with higher precision than interpo- lating classical topographic surveys. For instance, spatial remote sensing offers numerous advantages, complementing other conventional means of observation. One of the most outstanding contributions is its ability to follow dynamic processes, making it an invaluable source for studying changes that occur on the earth’s surface [18]. On the other hand, the accelerated development of technology has allowed the refining of small platforms that integrate airborne optical sensors in order to capture remote information in short periods of time, with good spatial resolution and low cost [1]. This technology has already been applied to quarry explorations [19]. Similarly, Unmanned Aerial Systems (UAS) have become a useful source for monitoring and analyzing morphometric changes generated by extractive activities through detailed photogrammetric flights, from which orthophotos and Digital Elevation Models (DEMs) of centimeter-level resolution are obtained [20]. In this case, the variation of extracted volume can be determined through the detection of very accurate changes in altitude between the current state of the terrain (once the extractive activity has ended) and the original state (prior to the start of mining). In the absence of a DEM prior to the start of mining, a possible solution is to make an estimate based on the interpolation of height references of unaltered locations during the mining activity, as has been done in erosion studies [16,17]. The identification and digitalization of these points are performed by expert visual interpretation supported by the use of orthophotos. The corresponding altitude value is added as an attribute according to the DEM generated. However, some of those methods are more complex to use and not applicable to steep slopes like those that are dominant in the vast majority of the restorations carried out in the context of open mining. The main objective of this study is to propose an approach to estimate a relief model that simulates the filling surface (Digital Elevation Model of the filling volume of material (DEM(r))). This model allows an approximate quantification of the volume required to restore the gravel pit based on the UAS-extracted terrain model and point interpolation techniques from witnesses’ points that try to reproduce the original morphology in order to have a better approximation of the volume to be filled. The main novelty of the proposed methodology is to analyze the integration of this filling model with the morphological structure of its environment, focusing on the connection with the water network and subsequent analysis of the impact of the visual basin. A secondary objective is to provide a subsequent analysis of the impact of the visual basin. A secondary objective is to provide a GIS tool to the readership capable of reproducing the methodology described (in this case, implemented in ArcMap 10.7). Sensors 2023, 23, 2097 3 of 15 2. Study Area, Materials, and Methods 2.1. Study Area GIS tool Thetostudy the readership capable area is located inof reproducing the the methodology region of Montsià, described municipality (in this case, of Masdenverge, implemented in ArcMap 10.7). province of Tarragona—Catalonia (Northeast of Spain) (Figure 1A,B). The extractive ac- tivity called Belén EA20040002 is located at the geographic coordinates Latitude 40°43′49″ 2. Study Area, Materials, and Methods N and Longitude 0°32′54″ E, between 53 and 87 m.a.s.l. The company Construcciones y 2.1. Study Area Ayudas Ferroviarias, S.L. dedicated to the extraction of gravel is operating it [21]. The study area The official is located boundaries in Belén of the the region gravelofpit Montsià, municipality were downloaded of the from Masdenverge, Hipermapa province of Tarragona—Catalonia (Northeast of Spain) (Figure 1A,B). The extractive viewer, the catalog of layers of extractive activities of the Government of Catalonia activity [21]. called Belén EA20040002 is located at the geographic coordinates Latitude 40 ◦ 430 49” N and According to the official boundaries, the extraction area is 2.71 ha, which is replaced by a Longitude ◦ 0 54” E, between 53 and 87 m.a.s.l. The company Construcciones y Ayudas larger area0of323.24 ha after assessment by photo-interpretation (see Section 2.3.1) (Figure Ferroviarias, 1). S.L. dedicated to the extraction of gravel is operating it [21]. Figure 1. Figure 1. (A) (A)Location Locationofofthe thestudy studyarea areaatat the the autonomous autonomous precinct precinct level. level. (B)(B) Location Location at the at the prov- province ince level. (C) Location at the municipality level. (D) Orthophoto view of the Belen level. (C) Location at the municipality level. (D) Orthophoto view of the Belen waste dump. waste dump. At the The geological official level,of boundaries the thesubstrate is generally Belén gravel pit werecomposed downloaded of sedimentary material from the Hipermapa and rocks (grayish limestones, gray, black and red shales, and blackish viewer, the catalog of layers of extractive activities of the Government of Catalonia [21]. peats) [22]. Ac- cording to to According thethe visual analysis official of the time boundaries, series 2000–2021 the extraction orthophotos area is 2.71 ha, which[23], it can be is replaced byin- a ferredarea larger thatofthe exploitation 3.24 activities by ha after assessment began between 2003 and photo-interpretation (see2004 because Section 2.3.1)in(Figure previous1). yearsAtnothe outstanding geological change level, theofsubstrate earthwork is in the study generally area wasofvisualized. composed sedimentary In material contrast, and 2005rocks (grayish images show limestones, gray, black very characteristic andtones bright red shales, and of bare blackish ground andpeats) [22]. Accord- the establishment ing to the for of access visual analysis heavy of the vehicles. timethat From series 2000–2021 point orthophotos on, the extraction of [23], it canmodified material be inferredthe that the exploitation activities began between 2003 and 2004 because in natural environment, mainly the topography and the vegetation cover, initially in the cen- previous years no traloutstanding and northern change of earthwork sectors. In 2010, theinentire the study area was authorized visualized. area was already In contrast, exploited2005and images during show 2011, very characteristic the extraction bright tones operations at theofcentral bare ground and the sectors and northern establishment finished,of result- access for ingheavy in the vehicles. growth ofFromsmallthat pointand bushes on,grasslands. the extraction of material In 2021, the wholemodified area wasthe natural already environment, undergoing a mainly processthe topography and of abandonment, the vegetation pending the startcover, initially actions of restoration in the central (Figure and 2). northern sectors. In 2010, the entire authorized area was already exploited and during 2011, the extraction operations at the central and northern sectors finished, resulting in the growth of small bushes and grasslands. In 2021, the whole area was already undergoing a process of abandonment, pending the start of restoration actions (Figure 2). Sensors 2023, 23, x FOR PEER REVIEW 4 of 16 Sensors 2023,23, Sensors2023, 23,2097 x FOR PEER REVIEW 44 of of 15 16 Figure 2. Multitemporal aerial orthophotos of the period 2000–2021 [23]. Figure2.2.Multitemporal Figure Multitemporalaerial aerialorthophotos orthophotosof ofthe theperiod period2000–2021 2000–2021[23]. [23]. 2.2. Materials 2.2. 2.2.Materials Materials 2.2.1. Drone 2.2.1. Drone Data 2.2.1. DroneData Data A DJI A DJI Phantom Phantom 44 Pro Pro drone multicopter multicopter (four-rotor) with with an RGBRGB 1″ CMOS CMOS sensor A DJI Phantom 4 Prodrone drone multicopter(four-rotor) (four-rotor) withan an RGB1” 1″ CMOSsensorsensor (DJI) [24] (DJI) [24] hasbeen beenused usedforfor capturing images by planning and and executing a photogram- (DJI) [24]has has been used capturing for images capturing by planning images by planning executing a photogrammetric and executing a photogram- metric(Table flight flight (Table 1). The Rural of Agents of Catalonia carried out the photogrammetric metric flight1).(Table The Rural 1). The Agents Catalonia Rural Agents carried out of Catalonia the photogrammetric carried flight on out the photogrammetric flight 25 on October 25 October 2021, at 1 pm (CEST) approximately. The flight altitude was 100 m, flight on 252021, at 1 2021, October pm (CEST) at 1 pmapproximately. The flight (CEST) approximately. altitude The was 100 was flight altitude m, with 100 m,a with frontala frontal overlapoverlap of 80% of 80% and a and side a lapside of lap 70%, of 70%, obtaining obtaining a pixela pixel spatialspatial resolution resolution of 5 of cm. with a frontal overlap of 80% and a side lap of 70%, obtaining a pixel spatial resolution of 5 cm. photogrammetric After After photogrammetric and SfM andsoftware SfM software post-processing post-processing (see Materials (see Materials section), section), a DSMaa 5 cm. After photogrammetric and SfM software post-processing (see Materials section), DSM (Digital (Digital Surface Surface Model), Model), a DEM, a DEM, and anand an orthophoto orthophoto were obtained. were obtained. DSM (Digital Surface Model), a DEM, and an orthophoto were obtained. Table 1. 1. Characteristics of of the drone drone used(Phantom (Phantom 4 Pro)[24]. [24]. Table 1.Characteristics Table Characteristics ofthe the droneused used (Phantom44Pro) Pro) [24]. Width Width 5472 5472 px Width 5472px px High High 3648 3648 px px High 3648 px Resolution Resolution 20 mpx 20 mpx Resolution Bands 20 mpx R,G,B Bands R,G,B Bands Focal R,G,B f/2.8–11 Focal f/2.8–11 Focal f/2.8–11 2.2.2. 2.2.2. Lidar Lidar Data Data 2.2.2.Lidar Lidar Data Imaging Detection and Ranging) data were downloaded from the official Lidar (Laser (Laser Imaging Detection and Ranging) data were downloaded from the offi- mappingLidarlocal (Laser Imaging agency, Detection Institut and Ranging) Cartogràfic i Geològic data dewere downloaded Catalunya from [25]. It[25]. the offi- corresponds cial mapping local agency, Institut Cartogràfic i Geològic de Catalunya It corre- cial to mapping the second local agency, Institut Cartogràfic i Geològic de Catalunya [25]. It corre- sponds to the coverage of Catalonia second coverage carried carried of Catalonia out in 2016. out inLidar 2016.points are distributed Lidar points are distrib- in sponds toformat LAZ/LAS the second coverage by blocks of 2 kmof Catalonia × 2 km which carried out in 2016. guarantees Lidar that 95% of points arearea the block distrib- has uted in LAZ/LAS format by blocks of 2, 2 km × 2 km which guarantees that 95% of the block auted in LAZ/LAS minimum densityformat of 0.5 by blocks of points/m 2 km × 2 km resulting in awhich DEM of guarantees 2 m pixelthat size95% of the block (DEM-LIDAR). area has a minimum density of 0.5 points/m22, resulting in a DEM of 2 m pixel size (DEM- areasix The has a minimum available density quadrants ofof 0.5 points/m Lidar data were, resulting mosaicked, in and a DEM of filtering noise 2 m pixelwas sizeapplied (DEM- LIDAR). The six available quadrants of Lidar data were mosaicked, and noise filtering LIDAR). by The six eliminating available values of noisequadrants categoriesof1,Lidar 7, 11, data and 13were mosaicked, according to the and tablenoise filtering of categories was applied by eliminating values of noise categories 1, 7, 11, and 13 according to the table and wasbasic description applied of LAZvalues by eliminating files [25,26]. of noise The interpolation categories of and 1, 7, 11, Lidar13points was performed according to the table of categories and basic description of LAZ files [25,26]. The interpolation of Lidar points using the neighborhood statistics method, in which the mean value of of categories and basic description of LAZ files [25,26]. The interpolation of Lidar points all the points that was performed using the neighborhood statistics method, in which the mean value of all are waswithin a radius performed of the using 2 mneighborhood is assigned to statistics each cell.method, Likewise, in the the in which interpolation mean value tool, a of all the points that are within a radius of 2 m is assigned to each cell. Likewise, in the interpo- selection the points of that classes are2within and 8 (ground a radius points) of 2 m iswas applied, assigned toin order each toLikewise, cell. eliminate inhigher values the interpo- lation tool, a selection of classes 2 and 8 (ground points) was applied, in order to eliminate of ground lation tool,level. A2m a selection of× 2 m pixel classes 2 andsize DEM with 8 (ground a coverage points) was applied, of 12inkm × to order 12 eliminate km was higher values of ground level. A 2 m × 2 m pixel size DEM with a coverage of 12 km × 12 obtained, whichofcovered higher values groundthe surrounding level. A 2 m × 2area of the m pixel gravel size DEM pit. Lidar with point cloud a coverage data of 12 kmwere × 12 Sensors 2023, 23, x FOR PEER REVIEW 5 of 16 Sensors 2023, 23, 2097 5 of 15 km was obtained, which covered the surrounding area of the gravel pit. Lidar point cloud data wereusing processed processed using[27] LAStools LAStools [27]For software. software. the GISFor the GIS analysis as analysis as well well as for as for the the modeling, modeling, ArcMap ArcMap 10.7 10.7 andsoftware and MiraMon MiraMon software have have been used. been used. 2.3. 2.3. Methods Methods The Themethodology methodologyappliedapplied to to estimate estimate aa restored restoredfilling fillingmodel modelthat thatserves servesto toperform perform volume calculations and analyze its integration into the environment through volume calculations and analyze its integration into the environment through the analysis the analysis of water flowflow of water andand its impact on the its impact onvisual basinbasin the visual is structured in 4 stages: is structured (1) Photogrammetric in 4 stages: (1) Photogram- flight, (2) Data generation, (3) Processing and analysis, and (4) Presentation metric flight, (2) Data generation, (3) Processing and analysis, and (4) Presentation of results of re- (Figure 3). Each stage is described in detail hereafter. sults (Figure 3). Each stage is described in detail hereafter. Figure3. Figure 3. Methodology Methodology workflow, workflow,from fromStage Stage11 (photogrammetric (photogrammetric flight), flight), Stage Stage 22 (Data (Data generation), generation), Stage 3 (Processing and analysis), and Stage 4 (digital mapping results). Stage 3 (Processing and analysis), and Stage 4 (digital mapping results). 2.3.1. 2.3.1. Photogrammetric Photogrammetric Flight Flight Agisoft Agisoft PhotoScan PhotoScan [28][28] photogrammetric photogrammetric and and Structure Structure from fromMotion Motion(SfM) (SfM)software software has has been used to process the flight data, obtaining a Digital Surface Model (DSM) and been used to process the flight data, obtaining a Digital Surface Model (DSM) and aa Digital Elevation Digital Elevation Model (DEM) of 15 cm spatial resolution as well as an orthophoto ofof5 15 cm spatial resolution as well as an orthophoto 5cm cmspatial spatial resolution. resolution. ForFor the the drone-derived drone-derived product product georeferencing, georeferencing, corresponding corresponding coor- coordinates between dinates between the drone the drone images images and and the official the official mapping mapping agencyagency images images [25] were [25] were used used as control as control points; points; this allowed this allowed an accurate an accurate fittingfitting of the of the drone drone orthophoto orthophoto and theand the official official agency orthophoto, and also between the drone-derived agency orthophoto, and also between the drone-derived DEM (DEM-DRONE) and theDEM (DEM-DRONE) and the official official agency agency DEM DEM [16]. [16]. ThisThis georeferencing georeferencing method method is crucial is crucial for for further further matching matching be- between the DEM-DRONE and the surrounding DEM-LIDAR. The tween the DEM-DRONE and the surrounding DEM-LIDAR. The imagery post-processing imagery post-processing workflow workflowused usedisisdefined definedin inseveral severalprevious previousstudies studiescarried carriedout outbybythe theauthors authors[29–32]. [29–32]. 2.3.2. Data Processing and Analysis with GIS (I): Exploited Area Boundaries vs. 2.3.2. Data Processing and Analysis with GIS (I): Exploited Area Boundaries vs. Official Boundaries Official Boundaries The first analysis required for an extractive activity restoration is to evidence if the The first area real exploited analysis required is located forthe inside an area extractive activity conceded restoration to extraction. is to evidence Boundaries of theifreal the real exploited area is located inside the area conceded to extraction. Boundaries exploited area were photo-interpreted performing a manual digitization of the contour of of the real exploited the area area that were was photo-interpreted truly affected using the performing a manual digitization drone orthophoto. of the The result was contour of inconsistent the area that was truly affected using the drone orthophoto. The result was with the official limits of the gravel pit, the latter being much smaller. Therefore, a new inconsistent with the official demarcation waslimits of the defined graveltopit, in order the latter achieve being a better much smaller.ofTherefore, representation a new the relief of the demarcation gravel wasmatching pit. This defined in order to analysis achieve afor is essential better representation quantifying the realof affected the reliefarea of the to gravel pit. This matching analysis is essential restore, and for environmental administration. for quantifying the real affected area to re- store, and for environmental administration. Sensors 2023, 23, x FOR PEER REVIEW 6 of 16 Sensors 2023, 23, 2097 6 of 15 2.3.3. Data Processing and Analysis with GIS (II): Generation and Validation of a Restoration 2.3.3. Relief Model Data Processing and(DEM(r)) Analysis with GIS (II): Generation and Validation of a Restoration Often, Relief there is Model (DEM(r))cartographic information on the study area before the no available beginning Often,oftherethe exploitation is no available phase. This was the cartographic case with the information on theBelén gravel. study areaTherefore, before theit was necessary beginning of the to exploitation estimate a DEM for this phase. Thisspecific was themomentum, case with the called Belén time 0 (T0),Therefore, gravel. to model itthe projected was necessary restored terrain aDEM(r). to estimate DEM for The strategy this specific assumes momentum, that thecalled borderline time 0of(T0), the ex-to tractive basin is the only altitude reference from the past. The model the projected restored terrain DEM(r). The strategy assumes that the borderline of reference values for gener- ating the the DEM(r) extractive basin were is theextracted from reference only altitude the overlap of the from the contour past. Theofreference the gravel pit with values for DEM-DRONE generating resulting the DEM(r) from were the drone extracted fromphotogrammetric the overlap of the contour flight. Given of the that gravel thepitdrone with flight is georeferenced DEM-DRONE resulting fromusingthe control drone points extracted from photogrammetric flight.theGivenDEM-LIDAR that the drone data, the flight method is assures the georeferenced integration using of DEM(r) control points through extracted from thethe contour DEM-LIDARline interface, data, with the im- the method mediatethe assures surroundings integration of represented DEM(r) throughby the the DEM-LIDAR. contour line This way, the interface, with filling model was the immediate surroundings represented by the DEM-LIDAR. This way, the spliced with the current configuration of the adjacent relief, thus avoiding abrupt changesfilling model was spliced with the current configuration of the adjacent relief, thus avoiding abrupt changes in slope. in slope. A A continuous surface surface DEM(r) DEM(r)in inaagrid gridformat formatofof1515 cmcm resolution resolution waswas generated generated by by interpolating interpolating z-coordinates z-coordinates (altitude) (altitude) located located on theon the dividing dividing ridgeridge of aof a basin. basin. How- However, ever, as in as theinDEM-LIDAR, the DEM-LIDAR, the interpolation the interpolation of the of the points points was also wasperformed also performed using theusing the neigh- neighborhood statistics method, obtaining borhood statistics method, obtaining a 2 m pixel size grid. a 2 m pixel size grid. The The DEMDEM of of the the restored restored relief relief (DEM(r)) (DEM(r)) obtained obtained as as aa first first result result is is aa model model thatthat simulates the filling surface required. However, the authors advise simulates the filling surface required. However, the authors advise that it can present that it can present some imperfections some imperfections (Figure(Figure 4A) because only altitude 4A) because valuesvalues only altitude corresponding correspondingto the contour of the to the contour gravel pit were of the gravel pitinitially chosen, chosen, were initially withoutwithout considering higher altitude considering higher values altitudewithinvaluesthe area. within The lack of references inside the area can generate a filling model the area. The lack of references inside the area can generate a filling model presenting presenting protrusions corresponding to the highertoaltitude protrusions corresponding the highervalues. In order altitude values.to improve In order to theimprove final reliefthe and finalthe re- restoration lief and the restoration model of the morphology of the land, the method includeshave model of the morphology of the land, the method includes new points that new been points added at thebeen that have maximumadded height at the locations maximumwithin heightthe gravel pit locations (in this within the study gravelcase 35 pit (in new points) based on the higher peak values of the actual drone this study case 35 new points) based on the higher peak values of the actual drone derived derived DEM-DRONE (Figure DEM-DRONE 4). (Figure 4). Figure 4. DEM(r) generated from (A) interpolated points only at the gravel limit and (B) adding points at maximum altitude locations inside the gravel pit. Sensors 2023, 23, x FOR PEER REVIEW 7 Sensors 2023, 23, 2097 7 of 15 Figure 4. DEM(r) generated from (A) interpolated points only at the gravel limit and (B) ad points at maximum altitude locations inside the gravel pit. The integrationThe integration of the DEM(r) withof the DEM(r) withwas its environment its made environment by mergingwas this made by merging DEM(r) DEM(r) with the DEM-LIDAR. It is worth noting a key step with the DEM-LIDAR. It is worth noting a key step to success in this integration: the to success in this integrat the control points to georeference of the drone photogrammetric control points to georeference of the drone photogrammetric flight were obtained from the flight were obtai fromon information based thethe information Lidar data;based on the therefore, Lidar there data; is no therefore, step there isboth effect between no step effect betw DEMs. To complete theboth DEMs. Toofcomplete integration the integration the DEM(r) and DEM-LIDAR,of the DEM(r) and DEM-LIDAR, the rasters were adaptedthe to rasters w adapted the same grid at to the 2 m pixel same size, andgrid at 2 mwas a mosaic pixel size, and to performed a mosaic was integrate performed both to integrate b files (Figure S1). For practical purposes, this result was called the Integrated Digital Elevation Model Digital files (Figure S1). For practical purposes, this result was called the Integrated (DEM(i)). vation Model (DEM(i)). 2.3.4. Data 2.3.4. Data Processing andProcessing and Analysis Analysis with GIS (III):with TotalGIS (III):Volumetric Filling Total Filling Volumetric Calculatio Calculation To perform theTo perform the volumetric volumetric calculation of calculation of filling filling required required to restore theto restore gravel pit,the gravel pit the classic volumeclassic volume formula formula (1) was (1) at applied wastheapplied at theApixel pixel level. finallevel. sum A of final sum all the of all the va values obtained obtained will be will be performed performed to compute tototal the compute the of volume total volume of filling. filling. V = Σ(a2 × h) V = Σ(a × h) 2 (1) where V is the volume, a is the pixel resolution (15 cm or 255 cm ) and h is the Dig 2 where V is theHeight volume, a is the Model Difference pixel resolution (DHdM) [31], (15 obtained cm or 255 cm2the from ) and h is thebetween difference Digital DEM(r) Height Difference Model DEM-DRONE. (DHdM) [31], obtained from the difference between DEM(r) and DEM-DRONE. The operation was performed on a raster calculator that operates at the pixel le The operation wasinperformed resulting on a raster a new spatially detailedcalculator that operates map informing about theat height the pixel level, (h) requ attribute resulting in a new spatially detailed map informing about the height attribute (h) required to reach the fill level at each cell. This attribute is then transformed to the volume of fil to reach the fill(m level at each per 3) required cell.pixel This to attribute restore is then the transformed morphology to the of the volume terrain of filling each p by multiplying (m3 ) required per attribute by its own area (in this study case, at 15 cm × 15 cm the pixel pixel pixel to restore the morphology of the terrain by multiplying each area is 0.0255 attribute by its(Figures own area (in this × 15iscm 2) 5 and S1). study A totalcase, at 15ofcm volume filling the pixel obtained byarea is 0.0255 adding up the mvalues of e (Figure 5 and Figure pixel. S1). A total volume of filling is obtained by adding up the values of each pixel. Figure 5. (A) Calculation of the height difference (DHdM) value between the restored DEM (DEM(r) and the drone-acquired DEM-DRONE. (B) Calculation of the volume needed to obtain the DEM(r), obtained by multiplying the pixel area (in this study case 0.255 m2 ) and the DHdM pixel value. Sensors 2023, 23, 2097 8 of 15 2.3.5. Data Processing and Analysis with GIS (IV): Water Flow Modeling This application is specifically designed for the creation of hydrological modeling, taking into account the maximum and minimum levels of the DEM(i) input file, generating a model that simulates a drainage structure without losing the surface continuity of the relief. To analyze the integration of the estimated filling model with its environment, it is considered appropriate to simulate the spatial distribution of the surface runoff of the area to be restored, understanding that the runoff is conditioned by the morphology of the terrain. The direction of the decline of each cell is calculated through the altitude values of DEM(i), taking as a reference the method of D8 [32], which considers that a cell has the option of the direction of its 8 closest neighbors, choosing the one with the lowest value. According to the number of times a cell converges towards another, a weight value is assigned, and the cells with the highest number of flow accumulation (higher weight) are classified as channels of the water flow network. At the scale of this study case, cells with an accumulated flow weight value > 50 were considered water stream channels; cells with values < 50 were assigned null values; cells of weight 0 identified the ridges (see Figure S2 in Supplementary Materials). 2.3.6. Data Processing and Analysis with GIS (V): Visual Basin Modeling In order to analyze the landscape impact generated by the restored relief of an extrac- tive activity, it is important to evaluate the visual scope that is perceived from different observation points. The calculation of the visual basin was based on DEM(i). Points of maximum height levels within the limit of the gravel pit were classified as observation points, taking a buffer of 4 km and considering elements such as roads and urban centers. ArcMap 10.7 visibility tool was used to calculate the visual basin [33], which determines the raster surface locations visible to a set of observer features or identifies which observer points are visible from each raster surface location. The result is a binary map assigning values 1 to visible cells and values 0 to non-visible cells. 3. Results The first result derived from the drone imagery is the identification of discrepancies between the real exploited area boundaries and the official boundaries. This is due to the extraordinary detail of the drone imagery lower than 10 cm, which allows for a very accurate photo-interpretation. The second important result is obtaining a restoration relief modeling (DEM(r)) based on the DEM-DRONE heights of the exploitation contour boundaries and peak point heights inside the exploitation basin. This result allows the calculation of a volumetric estimation for the filling restoration. A third result is obtaining the integrated DEM(i) based on the merging of the DEM(r) with the official mapping agency DEM (DEM-LIDAR). This result allows the calculation of the resulting hydrological channel network and the visual impact. Let us look at each result in detail. 3.1. Volumetric Calculation The quantitative analysis of the inner extractive basin is important for both the propri- etary benefits and the environmental restoration. The volumetric calculation, computed as the difference (DHdM) between DEM(r) and DEM-DRONE, gives a 15 cm × 15 cm sampling grid for the entire pit area. Each pixel area, multiplied by the height difference value, offers extremely detailed volumetric data in spatial terms. A total of 1,428,467 pixels of 15 cm × 15 cm were classified as filling pixels and as a function of the filling volume in order to identify the severity of the restoration needs (see Figure S3). The total volume required to restore the relief of the Belen gravel pit corresponds to 124,941 m3 . The area and volume required to restore the gravel pit are distributed as shown (Table 2 and Figure 6). volume in order to identify the severity of the restoration needs (see Figure S3). The tota volume required to restore the relief of the Belen gravel pit corresponds to 124,941 m3. Th area and volume required to restore the gravel pit are distributed as shown (Table 2 and Sensors 2023, 23, 2097 Figure 6). 9 of 15 Table 2. Volumetric calculation by the intensity of volume class. Table 2. Volumetric calculation by the intensity of volume class. Classes m3 Area m2 Volume m3 0 Classes 209.97 m3 Area 209.97 m2 0.00 Volume m 3 0.00 0–0.05 0 10,401.93 209.97 209.97 0.00 12,546.50 0.00 0–0.05 10,401.93 12,546.50 0.05–0.10 10,391.96 0.05–0.10 10,391.96 34,645.90 34,645.90 0.10–0.15 5862.53 Total landfill area 31,384.08 Total backfill Total backfill volume 0.10–0.15 5862.53 Total landfill 31,384.08 volume 0.15–0.20 0.15–0.20 3857.92 3857.92 32,140.49 area 32,140.49 29,642.78 29,642.78 124,941.22 124,941.22 0.20–0.25 1226.48 11,887.03 0.20–0.25 1226.48 0.25–0.30 399.67 11,887.03 4834.93 0.25–0.30 399.67 4834.93 Figure6.6.Spatial Figure Spatial distribution distribution of volume of the the volume required required to restore to restore thegravel the Belén Belénpit. gravel pit. ItItisisimportant to present important the GIS to present thetool GISthat toolautomatizes the generation that automatizes of the reliefofre-the relie the generation stored DEM(r) and compares it with the current relief DEM, obtaining restored DEM(r) and compares it with the current relief DEM, obtaining a raster of differences a raster of differ (DHdM) and a statistics report (Figure 7 and Supplementary Materials). ences (DHdM) and a statistics report (Figure 7 and Supplementary Materials). Sensors 2023, 23, 2097 10 of 15 Sensors 2023, 23, x FOR PEER REVIEW 10 of 16 Figure 7. Volume calculation workflow diagram in ArcMap 10.7. The reader can find the tool in the Figure 7. Volume calculation workflow diagram in ArcMap 10.7. The reader can find the tool in the Supplementary Materials (Volumetric_Tool.zip). Supplementary Materials (Volumetric_Tool.zip). 3.2. Water Flow Modeling Results 3.2. Water Flow Modeling Results A very dense water network has been obtained, not only within the restored model A very dense water network has been obtained, not only within the restored model but also integrated into the surrounding relief. but also integrated into the The surrounding method relief. that the estimated relief restored model presents good allows observing The method results allowssince observing that the estimated the water flow is relief restored not interrupted modelthe and/or because presents good by the flows generated Sensors 2023, 23, x FOR PEER REVIEW results since the water flow is not interrupted and/or because the flows 11 of 16 by the generated morphology follow a course and are integrated into the main network (Figure 8). morphology follow a course and are integrated into the main network (Figure 8). Figure 8. Figure 8. Drainage Drainagepattern patternintegrated intointo integrated the environment. the environment. 3.3. Visual Basin Modeling Results Results of the analysis of the visual basin of the gravel pit indicate that it has a visual impact of 17% on areas located within the 4 km buffer, so 83% of the watershed has no visual contact with the gravel pit. Two urban centers (Campredo and Masdenverge) have Sensors 2023, 23, 2097 11 of 15 3.3. Visual Basin Modeling Results Results of the analysis of the visual basin of the gravel pit indicate that it has a visual impact of 17% on areas located within the 4 km buffer, so 83% of the watershed has no visual contact with the gravel pit. Two urban centers (Campredo and Masdenverge) have Sensors 2023, 23, x FOR PEER REVIEW 12 of 16 also been identified within the impacted area, as well as two secondary roads and a railroad (Figure 9 and Table 3). Figure9.9.Visual Figure Visualbasin basinwith withaaradius radiusof of44km kmgenerated generatedfrom fromtwo twomaximum maximumpoints pointsidentified identifiedatatthe the northern (76 m) and southern (88 m) ends. northern (76 m) and southern (88 m) ends. Table 3. Results obtained from the visual basin. Table 3. Results obtained from the visual basin. Type Area (m2) Percentage (%) Type Area (m2 ) Percentage (%) Visible 8,642,592 m2 17 Visible 8,642,592 m2 2 17 Not visible 41,621,176 m 83 Not visible 41,621,176 m2 83 Total Total 50,263,768mm 50,263,768 2 2 100 100 4. Discussion 4. Discussion Drones offer very good quality products, such as centimeter-level resolution DEMs, Drones offer very good quality products, such as centimeter-level resolution DEMs, which are very useful for terrain morphology studies. At the same time, it is important to which are very useful for terrain morphology studies. At the same time, it is important have time series data, since from this information, it is possible to analyze dynamic pro- to have time series data, since from this information, it is possible to analyze dynamic cesses such processes suchasaschanges changesand andland landuse useininthe theterritory territory [34,35]. [34,35]. Specifically, Specifically, in in applications applications oriented toward the restoration of extraction areas, it could be considered oriented toward the restoration of extraction areas, it could be considered a disadvantage a disadvantage nottotohave not havedata dataprior priortotothe thestart startofof operations operations inin thethe gravel gravel pit,pit, sincesince information information regard- regarding the physical conditions of the territory is lost. To overcome this lack of information, itit ing the physical conditions of the territory is lost. To overcome this lack of information, isispossible possibleto toestimate estimatethe therelief reliefmorphology morphologyusing usingterrain terrainreferences referencesand andinterpolation interpolation methods [16,17]. methods [16,17]. Themethodology The methodologyproposedproposedin inthis thisstudy studyconsiders considersgenerating generatingaasimplesimplerelief reliefmodel model estimatedfrom estimated fromelevations elevationsobtained obtainedfromfromthethecontour contourof ofthe theextraction extractionarea areaas asaabasic basicproxy proxy for the initial estimation of filling volume in geomorphological restoration. for the initial estimation of filling volume in geomorphological restoration. However, the However, the method overpasses the problem of inner elevations in the gavel method overpasses the problem of inner elevations in the gavel basin. When superposing basin. When superposing the DEM(r) with the DEM-DRONE (current state of the relief), some elevations higher than the estimated model were identified, due to the existence of altitude values higher than the filling model. The analysis of this group of values revealed that they correspond to slightly altered areas that maintain topographic conditions of the original relief, so it is Sensors 2023, 23, 2097 12 of 15 the DEM(r) with the DEM-DRONE (current state of the relief), some elevations higher than the estimated model were identified, due to the existence of altitude values higher than the filling model. The analysis of this group of values revealed that they correspond to slightly altered areas that maintain topographic conditions of the original relief, so it is essential to identify these maximum elevations and add them to the interpolation method in order to estimate a filling model that completely covers the exploited area. The application of these fills can ensure the morphological integration of the gravel pit to its immediate surroundings, reducing external dynamic effects such as erosion, landslide, and rock fall, among others [36]. The volumetric calculation by indirect methods could have better results than those obtained through conventional methods, since the proposed methodology works with information at a pixel level of 15 cm spatial resolution, compared with the volume calculation method by cross sections that averages contiguous areas [7]. The integration of the modeled relief restored DEM(r) with the official mapping agency surrounding DEM (DEM-LIDAR) is crucial for the correct analysis. The integrated DEM (DEM(i)) was assured by georeferencing the drone flight using control point coordinates obtained from the DEM-LIDAR, and thus, the merging of both DEMs was good enough to avoid errors in the hydrological analysis. If this integration is wrong, it may present errors at the extremes, over the contour area, showing step pattern effects in the channel network. As an alternative to photo-interpreted control points from the official maps, in situ topographic work can improve the delimitation of the study area and capture altimetry data of time-invariant points to refine the interpolation model and obtain a finer integration of the filling model. The detection of water flow is an element of analysis that corroborates the integration of the model in the landscape; it also confirms that the estimated filling model is morpho- metrically close to the desired relief for restoring the gravel pit [31]. The water runoff path identification and the Strahler order [37] of each channel point out the main water erosion threat locations and allow tracking of the lixiviated particles. Finally, the analysis of the visual basin has shown that the integration of the UAS- extracted terrain model with the Lidar data can quantify, with high accuracy, the potential areas with views over the open pit mine. This is therefore a useful, fast, and accurate tool for decision making regarding the visual impact generated by extractive activities. 5. Conclusions Geomorphological restoration modeling of extractive activities has multiple solutions. This paper presents an approach based on the calculation of the volume of filling material needed for the restoration, based on information from drones and GIS processing. Nev- ertheless, the main contribution is the expansion of the scope out of the exploited area by analyzing the integration of the relief restoration model with its surrounding relief. The clue of the method resides in georeferencing the drone-derived Digital Elevation Model using photo-interpreted control points obtained from the official agency Digital Elevation Model, thus assuring a fine merging between both. Results show that the method improves the relief modeling of the restored space and, especially, contributes to the combined study of the modeled relief and the current relief context. Likewise, the modeled relief integration with the environment is corroborated by analyzing the water flow behavior. Although this method is still under users’ evaluation (users can test the provided tool (see Supplementary Materials)), this approximation to the final morphology is valid for estimating in an easy and very precise way the amount of filling material necessary to restore the morphology of the extraction pit and for modeling of its consequences to the environment. The method contributes to decision-making criteria for the treatment and morphological restoration of these types of spaces since the information has a very good spatial resolution, and the mathematical operations are carried out at a very high-resolution pixel level. Further research on the analysis using more variables to guarantee its effectiveness is needed, which will add more certainties to the models. Sensors 2023, 23, 2097 13 of 15 Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/s23042097/s1, Figure S1: Integration of DEM-LIDAR and DEM (r); Figure S2: Reclassified raster with flow accumulation values >50 and null values; Figure S3: Classi- fication of raster in cut and filling categories. The following file was also added: Volumetric_Tool.zip with the ArcMap 10.7 tool to calculate the restored DEM with sample data. Author Contributions: Conceptualization, J.-C.P., P.M., E.R. and V.C.; methodology, E.R., J.-C.P. and P.M.; software, E.R.; validation, J.-C.P., P.M. and V.C.; formal analysis, E.R.; investigation, P.M.; resources, V.C.; data curation, V.C. and P.M.; writing—original draft preparation, E.R.; writing— review and editing, J.-C.P., C.D.-M., P.M. and V.C.; visualization, C.D.-M.; supervision, J.-C.P., C.D.-M. and V.C.; project administration, V.C.; funding acquisition, V.C. All authors have read and agreed to the published version of the manuscript. Funding: This research was funded by the Government of Catalonia through the projects “Re- search and innovation in the process and restoration of extractive activities” (CREAF 01012015) and “GRUMETS” (SGR2021-00554) and by the European Union LIFE Preparatory project “New Life for Drylands”, grant number LIFE20 PRE/IT/000007. Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: All drone and airborne orthophotos, DEMs, shapefiles and code will be made available on request to the correspondent author’s email with appropriate justification. Acknowledgments: We are very grateful to the reviewers for their valuable comments that helped to improve the paper and to the RPAS Unit of the Rural Agents Corp (Government of Catalonia) for operating the drone flight. Conflicts of Interest: The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results. Abbreviations The following abbreviations are used in this manuscript: DEM Digital Elevation Model DEM(i) Integrated Digital Elevation Model DEM(r) Digital Elevation Model of the filling volume of material (restoration relief) DEM-DRONE Digital Elevation Model derived from drone DEM-LIDAR Digital Elevation Model derived from lidar DHdM Digital Height Difference Model [DEM(r)-DEM-DRONE] DSM Digital Surface Model GIS Geographical Information Systems IDW Inverse Weighted Distance SfM Structure from Motion software TIN Triangulated Irregular Network UAS Unmanned Aerial Systems References 1. Ren, H.; Zhao, Y.; Xiao, W.; Hu, Z. A review of UAV monitoring in mining areas: Current status and future perspectives. Int. J. Coal Sci. Technol. 2019, 6, 320–333. [CrossRef] 2. Ramli, M.; Thamrin, M.; Asrafil, M. Analysis of Soil Erosion on Mine Area. IOP Conf. Ser. Mater. Sci. Eng. 2020, 875, 012052. [CrossRef] 3. United Nations Environment Programme. Sustainability Reporting in the Mining Sector. Current Status and Future Trends. Available online: https://www.unep.org/resources/report/sustainability-reporting-mining-sector (accessed on 1 January 2023). 4. Young, R.E.; Gann, G.D.; Walder, B.; Liu, J.; Cui, W.; Newton, V.; Nelson, C.R.; Tashe, N.; Jasper, D.; Silveira, F.A.O.; et al. International principles and standards for the ecological restoration and recovery of mine sites. Restor. Ecol. 2022, 30, e13771. [CrossRef] Sensors 2023, 23, 2097 14 of 15 5. Carabassa, V.; Ortiz, O.; Alcañiz, J.M. RESTOCAT Instruccions d’aplicació. Protocol per a l’avaluació de la Restauració d’activitats Extractives. Available online: https://mediambient.gencat.cat/web/.content/home/ambits_dactuacio/empresa_i_produccio_ sostenible/restauracio_dactivitats_extractives/restocat/documents/restauracio_restocat.pdf (accessed on 1 January 2023). 6. Mitas, L.; Mitasova, H. Spatial Interpolation. In Geographical Information Systems: Principles, Techniques, Management and Applications; Longley, P., Goodchild, M., Maguire, D., Rhind, D., Eds.; John Wiley & Sons: New York, NY, USA, 1999; pp. 481–492. 7. Casanova, L. Topografía Plana. Available online: http://www.serbi.ula.ve/serbiula/libros-electronicos/Libros/topografia_ plana/pdf/topografia.pdf (accessed on 1 January 2023). 8. Achilleos, G.A. The Inverse Distance Weighted interpolation method and error propagation mechanism—Creating a DEM from an analogue topographical map. J. Spat. Sci. 2011, 56, 286–304. [CrossRef] 9. Longley, P.A.; Goodchild, M.; Maguire, D.J.; Rhind, D.W. Geographic Information Systems and Science, 3rd ed.; John Wiley & Sons: New York, NY, USA, 2010; p. 560. 10. Du, D.; Yan, S.; Yang, F.; Zhu, Z.; Song, Q.; Yang, G. Kriging Interpolation for Evaluating the Mineral Resources of Cobalt-Rich Crusts on Magellan Seamounts. Minerals 2018, 8, 374. [CrossRef] 11. Hutchinson, M.F.; Stein, J.A.; Stein, J.L.; Xu, T. Locally Adaptive Gridding of Noisy High Resolution Topographic Data. Available online: https://www.mssanz.org.au/modsim09/F13/hutchinson.pdf (accessed on 1 January 2023). 12. Hutchinson, M.F.; Xu, T.; Stein, J. Recent Progress in the ANUDEM Elevation Gridding Procedure. In Proceedings of the Geomorphometry 2011: Five days of Digital Terrain Analysis, Redlands, CA, USA, 7–11 September 2011. 13. Pons, X. MiraMon. Geographical Information System and Remote Sensing Software, Version 8.01b; Centre for Ecological Research and Forestry Applications (CREAF). Bellaterra. 2004. Available online: http://www.creaf.uab.cat/miramon/Index_usa.html (accessed on 1 January 2023). 14. GIS Development Team. QGIS Geographic Information System. Open Source Geospatial Foundation Project. Available online: http://qgis.osgeo.org (accessed on 1 January 2023). 15. ESRI. ArcGIS Desktop: Release 10.7; Environmental Systems Research Institute: Redlands, CA, USA, 2011. 16. Padró, J.-C.; Cardozo, J.; Montero, P.; Ruiz-Carulla, R.; Alcañiz, J.M.; Serra, D.; Carabassa, V. Drone-Based Identification of Erosive Processes in Open-Pit Mining Restored Areas. Land 2022, 11, 212. [CrossRef] 17. Carabassa, V.; Montero, P.; Alcañiz, J.M.; Padró, J.-C. Soil Erosion Monitoring in Quarry Restoration Using Drones. Minerals 2021, 11, 949. [CrossRef] 18. Chuvieco, E. Teledetección Ambiental: La Observación de la Tierra Desde el Espacio, 3rd ed.; Ariel: Madrid, Spain, 2010; p. 592. 19. Shahriari, H.; Honarmand, M.; Mirzaei, S.; Padró, J.-C. Application of UAV photogrammetry products for the exploration of dimension stone deposits: A case study of the Majestic Rose quarry, Kerman province, Iran. Arab. J. Geosci. 2022, 15, 1650–1666. [CrossRef] 20. Padró, J.-C.; Muñoz, F.-J.; Planas, J.; Pons, X. Comparison of four UAV georeferencing methods for environmental monitoring purposes focusing on the combined use with airborne and satellite remote sensing platforms. IJAEO 2019, 75, 130–140. [CrossRef] 21. Gencat. Activitats extractives. Expedients amb les àrees d’afecció de les Activitats Extractives a Catalunya a febrer de. 2014. Available online: https://sig.gencat.cat/visors/hipermapa.html (accessed on 1 January 2023). 22. Institut Cartogràfic i Geològic de Catalunya Mapa geològic comarcal 1:50.000. Available online: https://www.icgc. cat/Administracio-i-empresa/Descarregues/Cartografia-geologica-i-geotematica/Cartografia-geologica/Mapa-geologic- comarcal-1-50.000 (accessed on 1 January 2023). 23. Open ICGC. Connector QGIS Open ICGC. Available online: https://www.icgc.cat/en/Downloads/Open-ICGC-plugin (accessed on 1 January 2023). 24. DJI. Phantom 4 Pro. Available online: https://www.dji.com/es/phantom-4-pro (accessed on 1 January 2023). 25. Institut Cartogràfic i Geològic de Catalunya. Dades LIDAR. 2016. Available online: https://www.icgc.cat/Descarregues/ Elevacions/Dades-lidar (accessed on 1 January 2023). 26. Graham, L. The LAS 1.4 Specification. PERS 2012, 78, 93–102. Available online: https://www.asprs.org/wp-content/uploads/20 10/12/LAS_Specification.pdf (accessed on 1 January 2023). 27. LAStools. “Efficient LiDAR Processing Software” (Version 141017, Academic). Available online: http://rapidlasso.com/LAStools (accessed on 1 January 2023). 28. Agisoft LLC. Agisoft PhotoScan User Manual Professional Edition, Version 1.4. Available online: https://www.agisoft.com/pdf/ photoscan-pro_1_4_en.pdf (accessed on 1 January 2023). 29. Padró, J.-C.; Carabassa, V.; Balagué, J.; Brotons, L.; Alcañiz, J.M.; Pons, X. Monitoring opencast mine restorations using Unmanned Aerial System (UAS) imagery. STOTEN 2019, 657, 1602–1614. [CrossRef] [PubMed] 30. Padró, J.-C.; Muñoz, F.-J.; Ávila, L.; Pesquer, L.; Pons, X. Radiometric Correction of Landsat-8 and Sentinel-2A Scenes Using Drone Imagery in Synergy with Field Spectroradiometry. Remote Sens. 2018, 10, 1687. [CrossRef] 31. Carabassa, V.; Montero, P.; Crespo, M.; Padró, J.C.; Pons, X.; Balagué, J.; Brotons, L.; Alcañiz, J.M. Unmanned aerial system protocol for quarry restoration and mineral extraction monitoring. JEMA 2020, 270, 110717. [CrossRef] [PubMed] 32. Jaramillo Baltra, R.; Padró, J.-C. Generación de cartografía a partir de imágenes captadas con dron de ala fija, asociada a proyectos hidráulicos fluviales. Int. Rev. Geogr. Inf. Sci. Technol. 2020, 26, 93–117. [CrossRef] 33. Jenson, S.K.; Domingue, J.O. Extracting topographic structure from digital elevation data for geographic information-system analysis. PE&RS 1988, 54, 1593–1600, USGS Publications Warehouse. Sensors 2023, 23, 2097 15 of 15 34. Rotnicka, J.; Dłużewski, M.; Dabski, ˛ M.; Rodzewicz, M.; Włodarski, W.; Zmarz, A. Accuracy of the UAV-Based DEM of Beach– Foredune Topography in Relation to Selected Morphometric Variables, Land Cover, and Multitemporal Sediment Budget. Estuaries Coasts 2020, 43, 1939–1955. [CrossRef] 35. Lindner, G.; Schraml, K.; Mansberger, R.; Hübl, J. UAV monitoring and documentation of a large landslide. Appl. Geomat. 2016, 8, 1–11. [CrossRef] 36. Carabassa, V.; Ortiz, O.; Alcañiz, J.M. RESTOQUARRY: Indicators for self-evaluation of ecological restoration in open-pit mines. Ecol. Indic. 2019, 102, 437–445. [CrossRef] 37. Strahler, A.N. Quantitative analysis of watershed geomorphology. Eos Trans. AGU 1957, 38, 913–920. [CrossRef] Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.