Vegetation Cover Estimation in Semi-Arid Shrublands after Prescribed Burning: Field-Ground and Drone Image Comparison

: The use of drones for vegetation monitoring allows the acquisition of large amounts of high spatial resolution data in a simple and fast way. In this study, we evaluated the accuracy of vegetation cover estimation by drones in Mediterranean semi-arid shrublands (Sierra de Filabres; Almer í a; southern Spain) after prescribed burns (2 years). We compared drone-based vegetation cover estimates with those based on traditional vegetation sampling in ninety-six 1 m 2 plots. We explored how this accuracy varies in different types of coverage (low-, moderate-and high-cover shrublands, and high-cover alfa grass steppe); as well as with diversity, plant richness, and topographic slope. The coverage estimated using a drone was strongly correlated with that obtained by vegetation sampling (R 2 = 0.81). This estimate varied between cover classes, with the error rate being higher in low-cover shrublands, and lower in high-cover alfa grass steppe (normalized RMSE 33% vs. 9%). Diversity and slope did not affect the accuracy of the cover estimates, while errors were larger in plots with greater richness. These results suggest that in semi-arid environments, the drone might underestimate vegetation cover in low-cover shrublands.


Introduction
Vegetation cover, which can be defined as the percentage of ground area covered by green vegetation, is a widely used parameter to describe the amount of vegetation in a given area [1].The estimation of this parameter is essential in the monitoring and management of arid and semi-arid regions [2,3], due, among other reasons, to its importance in key biophysical processes such as evaporation, transpiration, and photosynthesis [4].Likewise, it is one of the parameters studied to evaluate vegetation recovery after disturbances, such as fires.
Estimates of vegetation cover have traditionally been obtained by classical vegetation sampling (e.g., plot surveys [5]) or classification of digital photographs [6][7][8].Although these approaches are still used for monitoring vegetation, they have certain limitations.In particular, they are relatively subjective, time-consuming and labor-intensive, and thus seem unrealistic for carrying out accurate large-scale, medium-to long-term dynamic monitoring [9].
Rapid technological advances in remote sensing have enabled an increase in its application in numerous scientific disciplines related to natural resource assessment.In recent years, a wide variety of remote sensing satellites and sensors with different spatial, spectral, and temporal resolutions have been developed, providing high-quality spatial and temporal data at a range of scales [10].Compared to traditional methods, remote sensing allows the estimation of vegetation condition over larger areas with greater accuracy and efficiency [11], optimizing the assessment of vegetation response to practices such as prescribed burning [12,13].Given these advantages, numerous studies have used spectral indices derived from remote sensing to estimate vegetation cover [2,[14][15][16][17].Most remote sensing studies have included field validations, although some appear to be biased due to the differences in the scale of the field sampling vs. the pixel size of the information obtained from remote sensing [18].
Likewise, there has been a boom in the use of remotely piloted aircraft systems (RPAS) or drones in ecological, forestry, ecosystem restoration, and environmental resource assessment research (among other fields).This is mainly due to the fact that users can acquire large amounts of high spatial resolution data that allow detailed mapping of structural and functional attributes of vegetation [9,10,[19][20][21].As with remote sensing studies, the use of drones for vegetation cover estimation requires field validation.
This study assessed vegetation cover after prescribed burning, defined as the use of fire in a planned manner to reduce fuel load under appropriate environmental conditions [22,23], in semi-arid Mediterranean shrublands [24,25].The aims of this study were: (i) to evaluate the accuracy of vegetation cover estimation by drone in semi-arid shrublands under different prescribed burning treatments; and (ii) to analyze the effect of type of vegetation, treatment type (prescribed burning and pyric herbivory), diversity, and slope on the estimation of vegetation cover.

Study Site
The study site is located on the north side of Sierra de los Filabres (Almería, Spain) (37 • 18 35.30"N 2 • 36 9.67" W) (Figure 1a), at 5.5 km and 3.5 km from the Special Conservation Areas of Calares de Sierra de los Filabres and Sierra de Baza, respectively, which are included in the Natura 2000 Network.The climate is semi-arid Mediterranean, with arid summers and cold dry winters.It is characterized by a mean annual temperature of 13.9 • C and mean annual precipitation of 444 mm.Plurennial droughts are relatively frequent and intense.The study area has steep topography (slopes ranging from 12% to 82%) with altitudes between 1170 and 1320 m, though maximum elevations in Sierra de los Filabres exceed 2100 m.The soils in the study area are classified as eutric cambisols, eutric regosols, and chromic luvisols with lithosols, mainly developed on schists and quartzites [26].

Experimental Design
Two management treatments were carried out in our study area (11.87 ha, see Table 1): prescribed burning (PB), and prescribed burning followed by grazing (i.e., pyric herbivory [23,25], PH) (Figure 1).The prescribed burns were carried out at two different times of year: autumn and spring (Table 1).After these burns, in the PH treatment, traditional grazing was carried out with a 300-head flock of sheep, whereas the PB zone was fenced to prevent livestock grazing.One of the goals of these experiments was to evaluate the role of livestock in wildfire prevention [29,30], since livestock may decrease fuel amounts, continuity, and height, and increase fuel moisture content [22,30,31].For each treatment, a number of randomly distributed 500-m 2 circular plots were selected: four for PB and eight for PH.Within each circular plot, eight 1-m 2 subplots were established, seeking to capture the variability in vegetation cover in this shrublands (before the prescribed burns): low-, moderate-and high-cover shrublands, and high-cover alfa grass steppe.Vegetation sampling was carried out to estimate plant diversity, richness, and vegetation cover.After PB, vegetation monitoring was carried out every autumn and spring until 2021.Drone flights were also carried out over the 12 circular plots to estimate plant cover (Figure 1b).

Vegetation Monitoring
In each of the 96 subplots, the total plant and bare soil cover were recorded, as was the cover of each plant species.These data were used to determine species cover, species richness, and plant diversity, estimating the latter with the Shannon-Wiener index [32] as where p i is the relative frequency of species i in each plot.We used guides to the flora of Spain and Portugal [33] and flora of eastern Andalusia [34] as references for taxonomical identification.Data from the vegetation surveys conducted in the spring of 2021 were selected for comparison with drone data.

Drone Flights and Image Processing
A total of six drone flights (three flights for RGB and three for multispectral sensors) were made on May 2021, with a multicopter DJI Matrice 210 RTK V2 (Da-Jiang Innovations, Shenzhen, China).Four of the 12 circular plots were covered in each flight.The flights followed a simple grid.To address the potential issues of scale heterogeneity due to lowaltitude flights and steep terrain [35], and to obtain the most accurate information possible, the flight routes were planned using UgCS PRO software [36].The flights were performed at an altitude of 30 m in the central hours of the day, when the sun's rays strike nearly perpendicular to the terrain, minimizing shadows that distort the spectral information.For the planning and execution of the flights, a digital terrain model updated in 2020 (previously obtained from a low-resolution prospective drone photogrammetric flight) was incorporated into the flight planning software, enabling the same height of image acquisition with respect to the terrain and the same spatial resolution to be maintained across the study area.A total of 8900 visible and multispectral images were acquired.The longitudinal and transversal overlaps were 80% and 75%, respectively, with a 120-min total flight time at 4 m/s.At the time of flights, the sun elevation ranged from 57.6 • to 72.9 • .
Prior to the acquisition of flight images, 16 ground control points (GCPs) were installed in the study area for georeferencing.The imagery co-registration and georeferencing was undertaken based on the WGS84 (UTM zone 30 N) coordinates of sixteen GCPs scattered across the study area following the recommendations of Martínez-Carricondo et al. [37].The three-dimensional (3D) coordinates of the GCPs were measured with very accurate equipment composed of a Trimble R6 receptor working in the post-processed kinematic (PPK) mode and a wayside base station.
For the acquisition of multispectral images, a Micasense RedEdge camera was used [38] which consists of a 5-band multispectral sensor (1.2 MP resolution): blue (475 nm), green (560 nm), red (668 nm), red edge (717 nm) and near infrared (842 nm); with a field of view of 47.2 • .The spatial resolution for each spectral band was 2 cm.A Zenmuse X7 camera with a resolution of 24 MP and a spatial resolution of 0.7 cm was used to obtain RGB images [39].
Once the flights had been completed, the images were preprocessed and georeferenced, and radiometric corrections were applied.The images acquired in the flights were processed with PIX4D photogrammetry software [40] that allowed the creation of the orthomosaic and the digital surface model of the complex under study (RGB sensor) and the generation of the different reflectance bands (multispectral sensor).The quality of the resulting imagery was evaluated using the root mean square error (RMSE): Supervised spectral classification was used to distinguish the different elements present in the study area.We used a support vector machine (SVM) as a classifier since it yielded a better solution (lower misclassification rate, higher model accuracy and faster training time) than other methods tested, such as decision trees or isoclustering (data not shown).The SVM classifier can be used on segmented and raster layers, which makes it an efficient tool as it reduces the number of processing steps [41].In addition, it allows load configurations that make it less prone to generating error due to noise, which is key in the detection of elements.This allowed us to discriminate between plant material and inert material, and to subsequently classify the different vegetation cover types.To evaluate the SVM, we used the Sorensen-Dice coefficient [42] calculated as: where A and B represent the two datasets (binary images), and | • | the cardinality.This similarity coefficient represents the degree of matching between two binary images.
The |A ∩ B| represents the number of pixels common to the two sets (imagery) obtained by automatic segmentation (the SVM algorithm) and manual segmentation (ground truth).The ground truth was obtained by manually segmenting RGB images of a smaller zone within the study area.Then, a multi-class Dice score, i.e., a mean value of the Dice coefficient was obtained for each class.On the other hand, an algorithm was generated to model the terrain from the clouds of millions of three-dimensional points generated by the photogrammetry software.In this way, raster files were generated with the information on the height, surface and phytovolume of each pixel of the study area previously classified as vegetation.From this point, the vectorization of the masks allowed extraction of the information on vegetation cover for each plot.

Statistical Analysis
First, we explored the correlation between the vegetation cover estimated from drone and that obtained from field measurements.The fit between the field measurement and the drone data was evaluated using the RMSE.We performed a regression analysis between the drone-based and field measurement-based cover estimates.Then, the predictive error sum of squares (PRESS) statistic was computed as: where ŷi,−i is the fitted values of a model with the ith observation removed.The lower the PRESS statistic the better the predictive power of the model.We also computed the predictive R 2 as: where SST is the total sum of squares: where y is the mean of the response variable.
To explore the potential differences in the correlation between drone-based and field measurement-based estimates as a function of the vegetation cover in our study area before the prescribed burns, we used two approaches.On the one hand, we explored how the correlation between the two cover estimates varied for each of the cover types identified in our study area before the burns: low-, moderate-and high-cover shrublands, and highcover alfa grass steppe.To make RMSE comparable between different cover types, the normalized root mean square error (NRMSE) was calculated as: where Max and Min were the maximum and minimum value respectively of the field measurement vegetation cover for each cover type considered.On the other hand, we used a classification and regression tree (CART) algorithm [43,44] to explore whether the correlation between the two approaches for the estimation of the vegetation cover (drone and field measurement) was homogeneous across all levels of cover, or there were groups which depended on the level of vegetation cover.With this approach, the data are partitioned along the predictor axes (in our case, the vegetation cover estimated by field surveys) into subsets with homogeneous values of the dependent variable (vegetation cover estimated by drone).To avoid overfitting, the final tree was pruned using a complexity parameter of 0.075 [45,46].Additionally, to explore whether the correlation between drone-and field measurementbased estimates of the plant cover were affected by the treatments, we compared the NRMSE of the correlations for each group of subplots, i.e., subplots within PB areas and those within the PH area (see Figure 1).
Finally, we assessed how the correlation between the vegetation cover estimates based on drone and field measurements was influenced by other variables related to plant diversity (richness and Shannon-Wiener index), and topographic features (e.g., terrain slope).For this purpose, we explored the relationship between the residuals of the dronefield measurement correlation and each of the variables of interest (richness, diversity, and slope).
All analyses were performed with R software v. 4.0.2[47].We used the rpart package [45] for the CART analysis.A full list of all the packages used (generated with the grateful package [48]) is included in the Supplementary Materials.

Results
Following the geometric correction of the images using the GCPs, the RGB sensorbased images showed RMSE values of 0.018 m, 0.016 m and 0.022 in the x-, y-and zdirections, respectively.For the multispectral sensor, RMSE values obtained were 0.0045 m, 0.0062 m, and 0.041 m in the x-, y-and z-directions, respectively.
The average ground sampling distance (GSD) of the RGB images were 0.60 cm•pixel −1 and 1.97 cm•pixel −1 for the multispectral images.The classification of the training set yielded a multi-class Dice score of 0.86.The vegetation cover estimated by drones was strongly correlated with that obtained by field surveys (R 2 adjusted = 0.821; R 2 predictive = 0.817; RMSE = 9.89%) (Figure 2).The linear relationship between the two metrics was significant (p-value < 0.0001).
tions, respectively.For the multispectral sensor, RMSE values obtained we 0.0062 m, and 0.041 m in the x-, y-and z-directions, respectively.
The average ground sampling distance (GSD) of the RGB images were 0 and 1.97 cm•pixel −1 for the multispectral images.The classification of the yielded a multi-class Dice score of 0.86.
The vegetation cover estimated by drones was strongly correlated with by field surveys (R 2 adjusted = 0.821; R 2 predictive = 0.817; RMSE = 9.89%) (Figure 2 relationship between the two metrics was significant (p-value < 0.0001).This relationship varied between the different vegetation cover classes, values ranging from 9.78 to 33.18% (Figure 3a).The high-cover alfa grass st the best fit between drone-based and field measurement-based estimates cover (NRMSE = 9.78%; R 2 predictive = 0.938), while areas with lower cover (low showed the highest NRMSE values (33.18%) and the lowest R 2 predictive (0.0145 cant regression was also obtained for this vegetation cover type (p = 0.069; F This relationship varied between the different vegetation cover classes, with NRMSE values ranging from 9.78 to 33.18% (Figure 3a).The high-cover alfa grass steppe showed the best fit between drone-based and field measurement-based estimates of vegetation cover (NRMSE = 9.78%; R 2 predictive = 0.938), while areas with lower cover (low shrublands) showed the highest NRMSE values (33.18%) and the lowest R 2 predictive (0.0145).No significant regression was also obtained for this vegetation cover type (p = 0.069; Figure 3a).
The results of the variance partitioning by vegetation cover in the field showed a split into two groups of data (Figure 4), with a cut-off threshold at 36% of vegetation cover.The NRMSE value for the lower vegetation cover group (i.e., subplots with less than 36% cover, see Figure 2) was higher (NRMSE = 30.68%)than that for the group with higher cover (more than 36%; NRMSE = 23.07%).The results of the variance partitioning by vegetation cover in the field showed a split into two groups of data (Figure 4), with a cut-off threshold at 36% of vegetation cover.The NRMSE value for the lower vegetation cover group (i.e., subplots with less than 36% cover, see Figure 2) was higher (NRMSE = 30.68%)than that for the group with higher cover (more than 36%; NRMSE = 23.07%).The comparison of the estimates for the PB treatment subplots with those for PH treatment subplots yielded very similar NRMSE results (Figure 3b).
No significant pattern was observed in the residuals of the relationship between the drone-based and field-measurement-based estimates of vegetation cover with either terrain slope of the subplots or plant diversity index (Shannon-Wiener) (Figure 5).In contrast, we obtained a significant but very weak pattern for these residuals as a function of richness (p < 0.01) (Figure 5), with larger residuals for higher values of richness.The comparison of the estimates for the PB treatment subplots with those for PH treatment subplots yielded very similar NRMSE results (Figure 3b).
No significant pattern was observed in the residuals of the relationship between the drone-based and field-measurement-based estimates of vegetation cover with either terrain slope of the subplots or plant diversity index (Shannon-Wiener) (Figure 5).In contrast, we obtained a significant but very weak pattern for these residuals as a function of richness (p < 0.01) (Figure 5), with larger residuals for higher values of richness.The comparison of the estimates for the PB treatment subplots with those for PH treatment subplots yielded very similar NRMSE results (Figure 3b).
No significant pattern was observed in the residuals of the relationship between the drone-based and field-measurement-based estimates of vegetation cover with either terrain slope of the subplots or plant diversity index (Shannon-Wiener) (Figure 5).In contrast, we obtained a significant but very weak pattern for these residuals as a function of richness (p < 0.01) (Figure 5), with larger residuals for higher values of richness.

Discussion
Overall, the results of our study show a strong correlation (R 2 predictive = 0.817; RMSE = 9.89%, Figure 2) between vegetation cover values estimated from field measurements and those estimated from drone data.This finding is in line with those from studies in shrublands in other ecosystems [3,[49][50][51], reinforcing the use of the drone method as a complementary, and less costly method of estimating the attributes of the vegetation, since it provides reasonably accurate measurement of vegetation structural attributes [21].Nonetheless, caution should be exercised when using vegetation cover estimates from drone-derived data, since we observed how the accuracy of the estimates varied depending on the vegetation coverage (Figures 3b and 4).The results of our two approaches, i.e., exploring the precision of the estimates as a function of vegetation cover types, and using a variance partitioning algorithm, suggest that the accuracy of the estimates is poorest for low vegetation cover (highest NRMSE values) (Figures 3a and 4).Specifically, for low-cover shrubland, it seems that drones underestimate vegetation cover.This is probably due to very small and scattered plants being detected by visual assessment but not by drone sensors.Our results are in agreement with recent studies that indicate an overestimation of vegetation cover when using traditional methods versus drone methods, mainly explained by the presence of species which either have highly branched, diffuse canopies with highly reduced leaves or have high overlaps between individuals, thus disproportionately increasing ground-based area totals relative to UAV-based area measures [52].
The positive relationship of the estimate residuals with richness seems to indicate that estimates vary more and are less accurate in subplots with higher species richness (Figure 5).Therefore, and according to our results, in these semi-arid environments, the estimation of vegetation cover with drones underestimates the real level in low-cover shrublands, which coincide with the areas containing the greatest number of species.This result is in line with recent findings indicating the unfavorable use of drones in the detection of ultra-fine scale demographics changes and species-level detail [52].However, in semiarid shrublands, the areas with very low cover (<25%) are key because they have low fuel loads and high plant richness, which are of great importance due to their roles in nutrient cycling, water retention, and as seed sources after a fire [53,54].
Regardless of our general good correlation between vegetation cover estimated by drones and that obtained by field measurements, caution is required when sampling larger areas.In our study, we flew at an altitude of 30 m above ground level, and the survey was carried out on 12 plots of 500 m 2 .However, larger areas would require the use of drones covering larger areas, more flights, or flights at higher altitudes, with the associated image resolution loss.Breckenridge et al. [49], in a comparison study of different RPAS vehicles to estimate vegetation cover, found that fixed-wing drone allowed the assessment of larger areas, but when high-detail data were required, multicopter vehicles performed better than fixed-wing drone.The flight altitude determined the flight distance, the flight duration, the size of images generated, and the computer processing time [55].Some studies have evaluated how drone flight height affects the estimation of vegetation parameters such us forest biomass [56] or forage mass [57].Others assessed how the estimation of vegetation cover varied using several flight altitudes in steppe ecosystems using a fixed-wing drone [49].DiMaggio et al. [57] did not find large differences in forage mass estimates when using different flight altitudes, but they did not include wood plants in their assessment.Therefore, an interesting question for semi-arid mountain shrublands is to determine how the estimation of vegetation cover varies as a function of flight altitude, and then establishing the optimum flight altitude for collecting data using RPAS multicopter.Other studies have proposed a multiple-scale approach combining field data with drone data, and then scaling up with satellite information to estimate the fractional cover of vegetation over large areas [50,58,59].Those approximations have been successfully applied in the estimation of the cover of some invasive species [50,58], or in the estimation of vegetation cover in tundra areas [59]; however, testing the approaches in low-cover but high-richness species areas, such as those in our study ecosystems, would be of interest.
Drone imagery has been promoted as an alternative to replace arduous and costly field work in remote areas and in adverse weather conditions; however, for this to be possible, there must be a high concordance between results from field data and those obtained with drones [3].In light of our results, it would be interesting to adapt the cover estimation method (drone-based vs. field measurement-based) to management objectives, and the types of shrubland in our study area.For instance, the potential use of drones to detect vegetation changes in early postfire shrubland communities has recently been highlighted [52].Thus, if the aim is to analyze cover in high-cover shrublands, drones may be useful.However, if the aims are focused on monitoring changes in cover in areas with high species richness, drone methods should be complemented by field estimation methods.Our results also suggest that, in addition to the potential problems related to variations in lighting conditions between flights or shadow related issues [60], it is important to consider cover type and the scale to be monitored, particularly in semi-arid environments where high spatial heterogeneity exists.In this sense, some studies have shown good agreement between field and drone measurements, both for the estimation of vegetation cover and for the estimation of bare ground [49], the biological importance of which is crucial in semi-arid environments.Likewise, to evaluate the effect of a prescribed burn, it may be interesting to see how reliable drones are in estimating bare ground, since after a disturbance such as prescribed burns, an increase in the percentage of bare ground is expected.Provided that we are aware of their possible limitations, multispectral images obtained with drones can play a key and complementary role in vegetation monitoring, as some studies have indicated [3,20], particularly in areas with high heterogeneity, such as semi-arid shrublands.

Figure 1 .
Figure 1.Location of the study area and experimental design (a).Two treatment zones were lished: prescribed burning (PB), and pyric herbivory (PH; pyric herbivory: prescribed burn grazing).Within each area, 500-m 2 circular plots were selected (b), and within these, eight subplots were established and evaluated by field surveys and drone (RPAS).

Figure 1 .
Figure 1.Location of the study area and experimental design (a).Two treatment zones were established: prescribed burning (PB), and pyric herbivory (PH; pyric herbivory: prescribed burning + grazing).Within each area, 500-m 2 circular plots were selected (b), and within these, eight 1 m 2 -subplots were established and evaluated by field surveys and drone (RPAS).

Figure 3 .
Figure 3.Comparison of RPAS (drone)-vs.field measurement-based estimates of vegetation cover by vegetation type (a); and by treatment (b).The adjusted R 2 , normalized root mean square error (RMSE), and predictive R 2 are shown.

Figure 3 . 14 Figure 4 .
Figure 3.Comparison of RPAS (drone)-vs.field measurement-based estimates of vegetation cover by vegetation type (a); and by treatment (b).The adjusted R 2 , normalized root mean square error (RMSE), and predictive R 2 are shown.

Figure 5 .
Figure 5. Relationship between model residuals (drone-based vs. field measurement-based estimates of vegetation cover) and plant community diversity (Shannon-Wiener index) (left), plant species richness (number of species) (center), and terrain slope (in percent) (right).

Figure 4 .
Figure 4. Results of the variance partitioning.The CART algorithm identified two separate groups with a threshold of 36% vegetation cover.

Figure 4 .
Figure 4. Results of the variance partitioning.The CART algorithm identified two separate groups with a threshold of 36% vegetation cover.

Figure 5 .
Figure 5. Relationship between model residuals (drone-based vs. field measurement-based estimates of vegetation cover) and plant community diversity (Shannon-Wiener index) (left), plant species richness (number of species) (center), and terrain slope (in percent) (right).

Figure 5 .
Figure 5. Relationship between model residuals (drone-based vs. field measurement-based estimates of vegetation cover) and plant community diversity (Shannon-Wiener index) (left), plant species richness (number of species) (center), and terrain slope (in percent) (right).

Table 1 .
Characteristics of prescribed burn treatments.