Forest Canopy Fuel Loads Mapping Using Unmanned Aerial Vehicle High-Resolution Red, Green, Blue and Multispectral Imagery

Abstract

Canopy fuels determine the characteristics of the entire complex of forest fuels due to their constant changes triggered by the environment; therefore, the development of appropriate strategies for fire management and fire risk reduction requires an accurate description of canopy forest fuels. This paper presents a method for mapping the spatial distribution of canopy fuel loads (CFLs) in alignment with their natural variability and three-dimensional spatial distribution. The approach leverages an object-based machine learning framework with UAV multispectral data and photogrammetric point clouds. The proposed method was developed in the mixed forest of the natural protected area of “Sierra de Quila”, Jalisco, Mexico. Structural variables derived from photogrammetric point clouds, along with spectral information, were used in an object-based Random Forest model to accurately estimate CFLs, yielding R² = 0.75, RMSE = 1.78 Mg, and an average ${Bias}_{rel}$ = 18.62%. Canopy volume was the most significant explanatory variable, achieving a mean decrease in impurity values greater than 80%, while the combination of texture and vegetation indices presented importance values close to 20%. Our modelling approach enables the accurate estimation of CFLs, accounting for the ecological context that governs their dynamics and spatial variability. The high precision achieved, at a relatively low cost, encourages constant updating of forest fuels maps to enable researchers and forest managers to streamline decision making on fuel and forest fire management.

1. Introduction

Forest fires play a crucial role in ecological disturbances in environments covered by vegetation, serving as natural shapers of forest ecosystems [1]. Three basic elements are required for a forest fire to occur: an ignition agent, oxygen, and forest fuels. Forest fuels encompass materials that have the potential to burn when exposed to an ignition source and consist of both dead organic matter and living vegetation [2]. According to the United States Department of Agriculture [3], forest fuel components can be grouped as ground (duff), surface (litter, herb, woody, and biomass within 2 m above ground), and canopy fuels (shrub, trees, and biomass above the surface). Forest fuel components are spatially distributed in three dimensions based on their inherent spatial variability and they are rarely correlated. Moreover, these components dynamically interact with the environment through complex biophysical processes [2]. Therefore, the necessity for accurate and updated fuel load maps becomes imperative for effective forest fuel management, demanding periodic measurements incurring high economic costs.

Wildfires are classified into three categories: ground, surface, and crown fires [4]. A crown fire is characterized as a fire typically initiated by a surface fire that ascends into the forest canopy, spreading through it. This type of fire proves challenging to control, often leading to catastrophic effects and demanding human intervention [5]. Reducing crown fire risk requires accurate description of canopy fuel loads (CFLs). Moreover, canopy fuels exert a significant influence on the characteristics of other fuel layers as they undergo constant changes, depending on the environment, weather, disturbances, and vegetation phenology. These factors condition fuels’ production, deposition, and decomposition [2]. Thus, CFLs synthesize many tree characteristics in an area with load, i.e., biomass, as its main indicator [6,7].

Recently, fire behavior models have included crown fire behavior predictions, allowing forest managers to make accurate decisions. Predictions are formulated to provide guidance for operational fire prevention and suppression responses, with a focus on prioritizing areas for hazardous fuel reduction management [8]. Several studies around the world have estimated CFLs, mainly through field measurements [9].

Thus, the estimation of CFLs is often conducted through field measurements, which are known to be highly expensive, time consuming, and labor intensive. Consequently, the literature frequently refers to the use of small sampling plots with limited re-measurements [10]. This approach is also difficult to apply in large areas and errors caused by human conditions tend to increase [11]. Prichard et al. (2019) [12] suggest plot-size units within the range of 0.04 to 0.10 ha to characterize and classify forest fuels, while Chávez-Durán et al. (2014) [13] and Ortíz-Mendoza et al. (2017) [14] have used plots of 0.05 ha. However, the small size of the plots did hamper the accurate estimation of CFLs, producing errors due to edge effects and the limited capacity to capture landscape variability [15].

Remote sensing offers a sound alternative to complement field measurements, especially over large and remote areas. The integration of spaceborne remote sensing with field inventory data has prompted the development of several techniques for estimating CFLs [16]. Multispectral passive satellite images have been demonstrated to be an effective and accessible tool. Passive sensors detect energy emitted or reflected from the objects and the environment. However, their primary limitation lies in their restricted sensitivity to forest structures, including fuel load, attributed to saturation at relatively low biomass levels. Active sensors provide their own source of energy to observe objects and are better suited for CFL estimation due to their higher sensitivity in capturing forest structures. In particular, Light Detection and Ranging (LiDAR) technology captures the vertical structure of forest data, allowing sharp precision and three-dimensional model generation [17]. Several studies have shown the suitability of LiDAR data to estimate canopy fuel properties from different platforms including terrestrial systems [18] and airborne [17] and spaceborne applications [19].

Several studies have proved the potential of combining active and passive remote sensors for the estimation of CFLs. Jiang et al. (2022) [20] used LiDAR and Sentinel-2 images with different Machine Learning techniques. Moran et al. (2020) [21] combined LiDAR data with Landsat through the Gradient Boosting Machine technique. García et al. (2017) [17] combined airborne LiDAR and Landsat OLI Data to extrapolate forest canopy fuel properties using Least-Squares Support Vector Machine techniques. A drawback associated with the previous studies is their coarse resolution. The spatial coverage of field plots corresponds to either a single pixel or a fraction, making it impractical to precisely estimate the variability within each plot [22]. Moreover, the overall variance tends to decrease as the pixel size increases, resulting in elevated determination coefficients that do not represent high fuel variability in the field [23].

Technological advances have paved the way for the emergence of platforms capable of transporting multiple types of sensors. LiDAR sensors have been attached to multiple airborne platforms, including conventional aircrafts and, recently, to Unmanned Aerial Vehicles (UAVs) [24]. Nevertheless, despite LiDAR’s advantages, its utilization remains expensive and the availability of UAV LiDAR sensors is still low. Instead, several algorithms have recently been developed to process conventional Red, Green, Blue (RGB) images through digital photogrammetry techniques for three-dimensional point clouds generation [25]. Furthermore, UAVs can carry sensors, allowing the acquisition of high spatial resolution images, temporal flexibility, and reliability [26].

Current advances in data collection equipment, processing techniques, and computing capability enable the development of novel methods for forest fuel management [27]. Moreover, data from diverse remote sensors as UAV high-resolution RGB and multispectral imagery, combined with field measurements and artificial intelligence techniques, are promising elements for appropriate the monitoring and management of forest fuels. Here, we present a novel approach: a cost-effective and easily replicable method which enables the monitoring of CFLs.

The aim of this paper was to devise a method for mapping the spatial distribution of CFLs in mixed forests. This was achieved by geospatial data: UAV RGB and multispectral imagery. The working hypothesis proposed that the devised method enables the consistent mapping of CFLs, considering their variability and three-dimensional spatial distribution. The results of this study offer several advantages, including optimizing time and costs involved in CFL assessment compared to traditional fieldwork. This approach enables researchers and forest managers to expedite decision-making processes for forest fire and fuel management.

2. Materials and Methods

2.1. Study Area

This research was conducted in the natural protected area “Sierra de Quila”, located in west–central Mexico, in the state of Jalisco (Figure 1); the altitude ranges from 1350 to 2540 m asl. “Sierra de Quila” hosts many flora and fauna species that maintain important ecological processes enclosed in 15,192.50 ha [28]. Vegetation is mainly composed of mixed temperate forest with the following species: Pinus douglasiana, Pinus devoniana, Quercus resinosa, and Quercus obtusata. Potential fire regime is characterized by frequent surface fires of low severity [29]. However, a devastating crown fire occurred in 1986, resulting in both forest fighter and ecosystem property losses [30]. The study area is also prone to natural and human disturbances, such as insect outbreaks, firewood production, and timber harvesting. Previous studies have classified the area according to Homogeneous Response Areas (HRAs), where field information collected in situ can be extrapolated to areas with similar characteristics [31].

2.2. Materials, Data and Methods

The estimation of the spatial distribution of CFLs involved three key stages: data collection, data processing, and CFL mapping. The data collection phase encompassed both direct tree measurements in permanent sampling plots and data from UAV remote sensors. The processing phase involved estimation of both CFLs from field data and UAV imagery processing including digital photogrammetry methods, multispectral analysis, and three-dimensional data segmentation. Finally, the mapping phase comprised the application of Machine Learning techniques, specifically Random Forests, to estimate CFLs from remotely sensed data (Figure 2).

2.2.1. Field Data

Three square sampling permanent plots 1 ha (100 × 100 m) each were established in the study area. The plot’s corners were marked with 2 cm accuracy, using a Ruide Total Station RTS-833 georeferenced with a Topcon GR-5 Global Navigation Satellite System (GNSS). Reference coordinate system was the Universal Transverse Mercator Zone 13 North (UTM 13N) (Figure 1).

In each plot, all living trees with normal diameter (ND) greater than 5 cm [32] were measured and tagged in May 2022 (Figure 3. Field measurements and tree canopies overview). For each individual tree, its genus and species was recorded; the relative spatial location of each tree within the plot was mapped using a Suunto Tandem compass^® (Suunto Oy, Vantaa, Finland) and a Truper^® tape measure (Truper, Mexico City, Mexico). ND at 1.30 m above the ground was measured using a diameter tape Forestry Suppliers^® (Forestry Suppliers Inc., Jackson, MS, USA). Crown diameter was determined with a Truper^® tape measure (Truper, Mexico City, Mexico); total height was measured with a Suunto clinometer^® (Suunto Oy, Vantaa, Finland).

CFLs for each tree were estimated using the allometric equations proposed by Rojas-García et al. (2015) [33], as described in

Table 1. Allometric equations used for the estimation of CFLs, from Rojas-García et al. (2015) [33]. Normal Diameter (ND) and Total Height (TH).

Species	Equation
Pinus devoniana	$0.182{·ND}^{1.936}$
Pinus douglasiana Pinus lumholtzii	$0.084·{ND}^{2.475}$
Pinus oocarpa	$0.058·{{(ND}^2·TH)}^{0.919}$
Quercus laeta	$0.0333·{ND}^{2.6648}$
Quercus candicans Quercus coccolobifolia Quercus obtusata Quercus resinosa	$0.0342·{ND}^{2.7590}$
Quercus rugosa	$0.0890·{ND}^{2.5226}$
Arbutus tessellata Arbutus xalapensis	$\left(0.3764·{ND}^2\right)-\left(2.3146·ND\right)-1.9106$

. In addition, the following structural variables were estimated by tree genus by plot: mean height, sum of basal area, sum of ground crown cover, and mean square diameter [34]. Previous studies have shown that differences in tree spatial distribution patterns impact the accuracy of modeling the spatial distribution of above-ground biomass [35]. Therefore, in order to assess the effect of spatial patterns on CFL estimation, Moran’s Index was used; this index ranges from 1 (strong clustering) to −1 (strong dispersion), with 0 indicating a random pattern with no spatial autocorrelation. Tree location and CFLs were used to determine spatial patterns, according to Equation (1) [36].

$$I=\frac{n}{W}\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{W}_{i,j}{z}_i{z}_j}{{\sum }_{i=1}^n{z}_i^2}$$

(1)

where $I$ is Moran’s Index; $n$ is number of observations; $W$ is the sum of weights $W_{i,j}$, $z_i$ is the difference of attribute in position $i$ with respect to its mean $(X_i-\overline{X}$), and $z_j$ is the difference of attribute in position $j$ with respect to its mean $(X_j-\overline{X}$).

2.2.2. Remote Sensing Data

UAV data were collected using a Sensefly Ebee^® drone (AgEagle Aerial Systems Inc., Wichita, KS, USA), equipped with a Parrot Sequoia+^® multispectral camera (Parrot, Paris, France) in September 2022. The sensor records spectral information in four bands: green (530–570 nm), red (640–680 nm), red edge (730–740 nm), and near-infrared (770–810 nm). Flight height was set to 212 m above the ground, resulting in a spatial resolution of 20 cm. Twenty-five images were captured per plot with 80% forward overlap and 60% side overlap. Prior to each flight, a target calibration panel was used to radiometrically correct the images [37]. Each multispectral flight was followed by a collection of RGB images, using a senseFly S.O.D.A 3D^® camera (AgEagle Aerial Systems Inc., Wichita, KS, USA) (Figure 3D–F), using the same flight height above the ground, which resulted in a spatial resolution of 7 cm. Forward and side overlap was 60% and 70%, respectively, with twenty-two images captured per plot. For all flights, a set of ground control points (GCPs) were established along with flight plans through eMotion 3.5.0^® software using grid flight missions for all flights [38]. The temporal lag between field and remote sensing data collection was based on tree species phenology [39,40].

2.2.3. Photogrammetric Point Clouds Generation

A Structure-from-Motion (SfM) algorithm was implemented to obtain three-dimensional data from S.O.D.A 3D^® camera (AgEagle Aerial Systems Inc., Wichita, KS, USA), which enabled the reconstruction of the three-dimensional structure of a scene from a set of two-dimensional images [41]. This algorithm was executed using the Agisoft Metashape 1.6.4^® software, which comprises three stages: (1) pre-processing, (2) densification, and (3) output generation.

The internal and relative orientation of individual images was applied based on their metadata spatial reference. Densification involved the assignment of descriptors to key points, which were matched using nearest-neighbor approximations and random sample consensus algorithms. Triangulation methods were used to estimate three-dimensional point positions and scene geometry reconstruction [41]. Point cloud, the three-dimensional collection of points in space representing the structure of the scene, was refined through the clustering stereoscopic multi-views algorithms [42]. Subsequently, relative coordinates were transformed into absolute coordinates by GCPs, achieving horizontal and vertical average RMSE values less than 0.50 m and 0.80 m, respectively. Point clouds were processed by linear triangulation resampling [43] and results were exported as point clouds and orthomosaic at the original spatial resolution.

2.2.4. Point Clouds Segmentation

Once the photogrammetric point clouds were generated, we attempted the identification of trees using a watershed-based approach. A point height normalization was first applied using the CloudCompare 2.12.4^® software. The Cloth Simulation Filter (CSF) tool was used to filter point clouds into ground and non-ground data [44]. Ground returns were subsequently used to create a digital elevation model (DEM) and non-ground returns were used to generate digital surface models (DSM). The restricted ability of passive sensors to penetrate the canopy led to a scarcity of ground points, impeding the generation of a precise DEM. Therefore, it was necessary to use the National Institute of Statistic and Geography of Mexico (INEGI) DEM, with 15 m spatial resolution [45], to carry out point height normalization. The average ground slope was relatively constant, with values below 10% in each plot; with these values, it was feasible to use the DEM from INEGI, applying the nearest-neighbor resampling method [46] and using 1 m spatial resolution.

Subsequently, a Canopy Height Model (CHM) was generated, portraying the height of the crown trees in relation to the ground at a spatial resolution of 1 m. CHM segmentation (automatic delineation of trees in the scene) was performed using TREETOP, a Shiny-based application which was originally developed for the extraction of ecological information from airborne LiDAR data [47]. TREETOP detects individual trees and delineates crowns from CHM using the local maximum algorithm; this algorithm assumes that CHM local maxima represent the highest point of each tree [48]. Crown edges were segmented using the Voronoi Tessellation-based algorithm [49]. The segmented outputs did not align with individual trees, but they represented clusters of intermingled trees. After segmenting the tree crowns, the Canopy Volume of each segment (CV) was calculated based on the area and the height, derived from the CHM of the polygons [25], using the Python programming language [50].

2.2.5. Multispectral Analysis

Twelve vegetation indices (

Table 2. Equations used to estimate vegetation index. Reflectance of near infrared ($NIR)$; Reflectance of red ($R)$; Reflectance of green (G); Reflectance of red edge ($Rededge)$; Soil adjustment factor ($L$) = 0.5; $\alpha$ = 0.96916.

Vegetation Index	Equation	Reference
Normalized Difference Vegetation Index (NDVI)	$\frac{NIR - R}{NIR + R}$	Rouse et al. (1974) [56]
Soil Adjusted Vegetation Index (SAVI)	$\left[\frac{\left(NIR - R\right)}{\left(NIR + R + L\right)}\right]·\left(1+L\right)$	Huete (1988) [57]
Modified Soil Adjusted Vegetation Index (MSAVI)	$\frac{2·\mathrm{N}\mathrm{I}\mathrm{R} + 1 - \sqrt{{\left(2·\mathrm{N}\mathrm{I}\mathrm{R} + 1\right)}^2 - 8\left(\mathrm{N}\mathrm{I}\mathrm{R} - \mathrm{R}\right)}}{2}$	Qi et al. (1994) [58]
2-band Enhanced Vegetation Index (EVI2)	$2.5·\left(\frac{NIR - R}{NIR + 2.4·R + 1}\right)$	Jiang et al. (2008) [59]
Difference Vegetation Index (DVI)	$NIR - \alpha R$	Richardson & Everitt (1992) [60]
Green Normalized Vegetation Index (GNDVI)	$\frac{NIR - G}{NIR + G}$	Gitelson et al. (1996) [61]
Green Ratio Vegetation Index (GRVI)	$\frac{NIR}{G}$	Sripada et al. (2006) [62]
Green Difference Index (GDI)	$NIR - R + G$	Gianelle and Vescovo (2007) [63]
Green Red Difference Index (GRDI)	$\frac{G - R}{G + R}$	Gianelle and Vescovo (2007) [63]
Red edge normalized difference vegetation index (NDVIre)	$\frac{NIR - Rededge}{NIR + Rededge}$	Gitelson and Merzlyak (1994) [64]
Red edge simple ratio (SRre)	$\frac{NIR}{Rededge}$	Gitelson and Merzlyak (1994) [64]
Datt4	$\frac{R}{G·Rededge}$	Datt (1998) [65]

), together with homogeneity texture for each spectral band, were estimated; the moving window size considered in the analyses was 5 × 5 pixels. Spectral indices were used to estimate canopy fuel properties, including canopy fuel load [16]. Vegetation indices and textures were generated by multispectral orthomosaics [51]. One of the most efficient and widely used texture methods is the Gray Level Co-occurrence Matrix (GLCM) [52], whose features are based on neighboring pixel pair occurrence. This method produced different indicators such as contrast, correlation, energy, entropy, or homogeneity, among others. Homogeneity texture is widely used for its great importance in vegetation data due to its high level of detail [53]. Vegetation indices provide spectral insights linked to vegetation photosynthetic activity, vigor, status, and cover [54], and all are considered as valuable predictors of ecosystem features. To reduce the number of variables and to avoid collinearity of vegetation indices and texture variables, a principal components analysis (PCA), along with its factorial matrix, was carried out [55].

2.2.6. CFL Spatial Distribution

Segments extracted using TREETOP were aligned with field-acquired CFL data to derive reference values for CFLs. Descriptive statistics were estimated for all variables (i.e., CV, vegetation index, and texture metrics) for each segment. Differences among plots were tested by the non-parametric Kruskal–Wallis test as homoscedasticity was not fulfilled [66,67]; according to the Levene test, there was no homogeneity of variances among plots (p-value < 0.05).

The spatial distribution of CFLs was estimated by the Random Forest (RF) algorithm [68] using Python programming language [50]. RF is one of the most efficient Machine Learning algorithms [69], which implements a regression technique through bagging and random subspace methods [70]. The CFL was considered the response variable and CV, vegetation index, and textures as the explanatory variables. Segment metrics of the three plots were pooled and split into training data (70% of the sample) and validation data (30% of the sample). To train the model, 400 decision trees, bootstrap resampling, and square error were used [71]. The importance of each explanatory variable was evaluated through the mean decrease in impurity index (MDI) [72]. The algorithm was calibrated and applied using the following libraries: numpy [73], pandas [74], matplotlib [75], scikit [76], and GDAL OSGeo [77] for spatial data management.

The performance of the model was evaluated through the determination coefficient (R²) and Root Mean Square Error (RMSE) using validation data [68]. The RF-trained model was also applied over the total number of segments by plot; the R² coefficient and absolute and relative CFL Bias were used to validate the accuracy across plots. The absolute Bias was calculated to evaluate the behavior of the RF-trained model on each plot and the relative Bias was used to compare plots [78]. Both estimations were performed according to Equations (2) and (3).

$${Bias}_{abs}={\sum }_i^N\widehat{{CF}_i}-{\sum }_i^N{CF}_i$$

(2)

$${Bias}_{rel}=\frac{{\sum }_i^N\widehat{{CF}_i}-{\sum }_i^N{CF}_i}{{\sum }_i^N{CF}_i}·100$$

(3)

where ${Bias}_{abs}$ and ${Bias}_{rel}$ correspond to absolute and relative Bias, respectively; $\widehat{{CF}_i}$ is the estimated CFL for each segment and ${CF}_i$ is the CFLs estimated from the field data.

Statistical analyses were carried out using the following R libraries “readr, northest, car, FactoMineR, and ggplot2” [79,80,81,82,83] available in R-project 4.3.0 [84].

3. Results

3.1. Plots’ Forest Structure

Three genus and twelve species were recorded within the plots (Table 1). The largest number of trees was recorded in P1, with 541 individuals, followed by P3, with 332 and P2 with 236. Pinus, followed by Quercus, was the dominant genus in all plots; Pinus exhibited the largest Quadratic Mean Diameter, Basal Area, Total Height, and Surface Crown Covered (

Table 3. Tree genus structural variables: Quadratic Mean Diameter (Dg), Basal Area (BA), Ground Crown Cover (GCC), Canopy Fuel Loads (CFLs), and mean value for Total Height (TH) (Standard error). To BA, GCC, and CFL, the values in each column correspond to the sum of the respective variable in each plot.

Plot	Genus	Number of Trees	Dg cm	BA m2	TH m	GCC m2	CFL Mg
1	Pinus	202	35.88	20.42	19.74 (1.06)	8078.60	134.37
1	Quercus	339	15.81	6.65	12.92 (0.51)	3913.88	29.10
2	Pinus	129	46.66	22.05	25.37 (1.52)	6451.49	120.57
2	Quercus	107	25.28	5.37	16.02 (1.17)	3537.99	30.50
3	Pinus	156	40.73	20.32	25.38 (1.21)	7097.84	75.52
3	Quercus	166	21.18	5.85	11.21 (0.80)	3123.20	25.07
3	Arbutus	9	17.92	0.23	8.79 (3.07)	131.160	0.72

). Pinus douglasiana and Quercus resinosa were the most abundant species in both genus. The highest CFL value was recorded in P1, with 163.47 Mg, followed by P2 and P3, with 151.07 Mg and 101.31 Mg, respectively. Pinus contributed the greatest amount of CFL in all plots, reaching values greater than 74% in each plot.

Distinct spatial distribution patterns were observed among the three plots (

Table 4. Spatial distribution patterns of tree individuals by genus. Moran’s Index (MI) where +1 indicates strong clustering, −1 indicates strong dispersion and 0 indicates random pattern. (*) pattern tendency with significance of 0.05.

Plot	Genus	MI	p-Value
P1	Pinus	+0.177	<0.05 *
P1	Quercus	+0.077	<0.05 *
P2	Pinus	+0.209	<0.05 *
P2	Quercus	+0.033	0.51
P3	Pinus	+0.032	0.45
P3	Quercus	+0.012	0.61
P3	Arbutus	+0.107	0.42

). In P1, both Pinus and Quercus exhibited clustering, whereas in P2, only Pinus formed clusters, with Quercus individuals displaying random distribution. In contrast, all trees in P3 were randomly distributed.

3.2. Remote Sensing Estimates

A set of multispectral orthomosaics and CHMs were obtained as a product of digital photogrammetry and point clouds processing. The mean CHM values were 17.84 m, 21.95 m, and 20.73 m in P1, P2, and P3, respectively. The segmentation analysis identified 31, 28, and 28 segments, covering 7883 m², 8104 m², and 8782 m² in P1, P2, and P3, respectively (Figure 4).

The number of trees identified using TREETOP was 511, 229, and 322 for P1, P2, and P3, respectively. These values represented a detection rate of 94.5%, 97%, and 97.3% for P1, P2, and P3. The absence of certain trees within each plot led to a minor underestimation of CFLs at the plot level when focusing solely on the fuel load of the trees correctly detected (161.02 Mg, 147.60 Mg, and 98.61 Mg in P1, P2, and P3, respectively), as indicated in Table 3. P1 had the lowest omission in terms of CFL, despite the lowest detection rate, with 1.50%, followed by P2 with 2.29% and P3 with 2.67%.

According to the factorial matrix from the PCA, four variables were uncorrelated, exhibiting a robust relationship with the first three principal components: Normalized Difference Vegetation Index (NDVI), 2-band Enhanced Vegetation Index (EVI2), Red band homogeneity texture and Red edge band homogeneity texture.

The median CFL values per segment were 4.54 Mg, 5.09 Mg, and 3.35 Mg in P1, P2, and P3, respectively (Figure 5A), with no significant differences according to the Kruskal–Wallis test: H(2) = 3.64, p = 0.16. There were statistically significant differences for CV, vegetation index, and textures among plots. For CV, H(2) = 11.96, p < 0.05, with median values of 4000.66, 6740.58, and 6590.29 m³ for P1, P2, and P3, respectively (Figure 5B). Finally, for the NDVI, EVI2, Red band texture and Red edge band texture, the results were H(2) = 56.73 with p < 0.05, H(2) = 33.00, p < 0.05, H(2) = 36.96, p < 0.05, and H(2) = 12.48, p < 0.05, respectively (Figure 5C–F).

3.3. Spatial Distribution of CFLs

According to validation data, our RF model achieved a strong performance with R² = 0.75 and RMSE = 1.78 Mg. At the plot level, P1 reached R² = 0.98, ${Bias}_{abs}$ = −17.12 Mg, and ${Bias}_{rel}$ = 10.63%; P2 yielded R² = 0.85, ${Bias}_{abs}$ = 28.06 Mg, and ${Bias}_{rel}$ = 19.01%; and P3 offered a performance of R² = 0.72, ${Bias}_{abs}$ = 25.85 Mg, and ${Bias}_{rel}$ = 26.21%. (Figure 6). Therefore, the modeled CFL per plot was 143.90 Mg, 175.66 Mg, and 124.46 Mg for P1, P2, and P3, respectively (Figure 7D–F).

The most significant variable was CV, at 80.67%, followed by Red band texture at 7.00%, Red edge band texture at 5.99%, NDVI at 4.15%, and EVI2 ay 2.19%.

4. Discussion

The spatial distribution of the CFL, estimated through a combination of digital photogrammetry, segmentation, multispectral analysis, and Machine Learning, yielded robust results, according to R², RMSE, and Bias. The study of CFLs in mixed forests with structural complexity and different spatial distribution patterns was conducted through the integration of various techniques. These techniques allowed for a more comprehensive analysis and understanding of the CFL in the complex nature of mixed forests of the “Sierra de Quila”.

4.1. Plots Forest Structure

The Pinus genus exhibited the highest CFL values, contributing to over 70% in each plot. This aligns with the findings of a previous study by Villavicencio-García et al. (2005) [85], which investigated the entire altitudinal gradient of “Sierra de Quila” and identified species from the Pinus genus as dominant. Typically, Pinus species occupy the upper canopy, while Quercus species tend to dominate the middle canopy, with tree crowns in close proximity or even overlapping [86]. In many regions, including areas like “Sierra de Quila”, most Pinus species in Mexico rely on fire, with these regions characterized by frequent surface fires of low severity [29,87].

In the present study, the Pinus genus presented the highest values in structural variables such as quadratic mean diameter, basal area, tree height, and vegetation cover. This result is highly relevant as values of the CFL depend on tree structural variables [33]. Previous studies by Figueroa-Navarro et al. (2023) in “Sierra de Quila” [88] obtained similar results to those achieved in the present work, with average densities of 556 trees per ha and a mean basal area of 20 m²/ha, with Pinus douglasiana being the most dominant species, regarding Pinus genus and Quercus resinosa as the most dominant species of the Quercus genus.

Tree spatial distribution patterns differed among the three plots according to Moran´s Index. Moreover, differences among the plots also occurred in terms of CFLs. Thus, the highest degree of clustering was present in P1, showing the greatest values of CFLs (median value of 16.98 kg/m²). Conversely, the greatest random pattern was observed in P3, presenting the lowest CFL values (median value of 10.71 kg/m²). These differences limit the use of linear and allometric model regression because tree spatial distribution patterns prevent generalizable models from being obtained [89,90]. Other studies have identified that mixed forests tend to present random distributions following a disturbance occurrence and they are reverted to clusters over time [91]; therefore, these may be the cause of differences in the spatial distribution patterns of trees in our study area.

4.2. Remote Sensing

Information derived from the photogrammetric point clouds provided valuable results by capturing significant three-dimensional information on trees, despite the fact that LiDAR data lack penetration ability [92]. Mixed forests in Mexico make it difficult to accurately estimate CFLs, mainly due to the high tree density that results in crowns close to each other or even overlapping. In such forest types, estimates focusing on crown segments, rather than on individual trees, provide a viable alternative [93]. CV has been a good alternative for estimates focusing on crown segments. This agrees with Riggi et al. (2021) [25], who used CV to estimate the above-ground biomass of olive plantations by linear regression models, reaching up to R² = 0.68 and RMSE = 39.19 kg. In our work, mixed forests presented greater structural complexity, with irregular spatial distribution of trees. Nevertheless, texture and vegetation indices included in the RF algorithm proved to be a valuable complement to CV, allowing the inclusion of differences in trees’ spatial distribution patterns between plots [51].

Despite the more limited information extracted from photogrammetric point clouds as compared to LiDAR point clouds, our results showed strong agreement with field estimates of CFLs. Thus, CFL differences within the segments were less than 3% compared to the total CFLs measured in the field. This remarkable agreement could be attributed to the size of the plots used, which were notably larger than those used in other studies. This factor contributed to minimizing the border effect. Tree trunks typically deviate for a vertical orientation, resulting in a discrepancy between their geolocation and the boundaries of segments. This border effect is mostly evident in small plots, while in large plots this error tends to be compensated [15]. According to Knapp et al. (2018) [94], to attain precisions with errors below <10%, sample plots must be of at least 100 m per side. Our plots measured 100 m on each side, enabling the segmentation process and field measurements to capture the inherent three-dimensional spatial distribution of the trees in accordance with their spatial variability.

4.3. Spatial Distribution of CFLs

The RF regression model enabled an accurate description of the spatial distribution of CFLs within our plots. The trained model achieved errors of less than 30%, even when trees had different patterns of spatial distribution. The lowest errors were observed in P1, where trees showed the highest density. These contrast with the findings of Guerra-Hernández et al. (2016) [95], who suggested that higher-density plots may yield inaccurate results due to the limited penetration of photogrammetric point clouds into the canopy. In our case, we addressed this situation by incorporating multispectral data and textures, along with the segmentation approach based on surfaces rather than on individual trees. In addition, employing an external DEM with proven accuracy contributed to enhancing our estimates of the structure metrics derived from UAV data. The results in this study were similar to those reported by Maesano et al. (2022) [24], who achieved ${Bias}_{rel}$ close to 10% by combining data from LiDAR, RGB images, and RF algorithms. This confirms that our approach enables the development of robust and accurate models, even in the absence of LiDAR data.

According to the MDI, CV was the most important variable of the model, accounting for over 80%. In contrast, the texture and vegetation indices, considered together, showed importance values close to 20%. The influence of CV on accuracy was substantial, particularly in instances where trees exhibited clustered spatial distribution patterns, leading to improved model performance. Moreover, texture and spectral indices were useful in overcoming the limitations of point clouds and addressing the structural complexity inherent in mixed forests [51,54].

The developed method offers high potential to accurately estimate CFLs within the context of native mixed forests characterized by a complex spatial distribution of trees, even in the absence of LiDAR data. The proposed method for mapping CFLs’ spatial distribution not only allows for the optimization of costs and time associated with intensive field measurements [10] but also serves as a cost-effective alternative to LiDAR technology [96]. In this research, sampling plots showed relatively low ground and slope values, staying below 10%. The accuracy may be impacted in areas with higher slopes or predominantly rugged terrain; therefore, a DEM should be applied instead of a DSM. However, this situation can be solved using topographic measurements.

The field data utilized in this research belong to the same HRA of the study area, where they share a common ecological background. This context enables in situ data collection to extrapolate to areas with similar characteristics [31]. However, the robustness of the method ensures its replicability across various latitudes. Moreover, it is feasible to extend its applicability from a local geographic scope to wider areas thought upscaling techniques [17,22]. The present method uses data from field measurements, digital photogrammetry, multispectral analysis, and Machine Learning techniques, mostly through open source software; the method is also straightforward to implement. Although the method has undergone local testing, its transferability to other areas is seamless, facilitated by the easy and cost-effective collection of data across diverse geographical regions.

After collecting field measurements, the proposed method requires only a few hours for data processing. Moreover, trained Machine Learning models can be used in areas with similar ecological context, allowing forest managers to make informed decisions based on current and reliable information.

5. Conclusions

The present research introduces a methodology for mapping the spatial distribution of CFLs by employing image segmentation, multispectral analysis, and Machine Learning modeling. This approach enables the extraction of structural information, yielding commendable accuracy in CFL estimation. Mapping CFLs was consistent even when the study area presented a high level of complexity pertaining to the mixture of genera and species in the forest and the variation in patterns of the trees’ spatial distribution. Machine Learning modeling, in particular RF based on CV complemented with multispectral variables derived at the object level, enabled a high accuracy at a relatively low cost.

This method can be easily applied to drone-derived data, which is particularly important in areas with frequent disturbances, to obtain reliable and updated information for decision making. The present method offers the opportunity for CFL monitoring, which can be used to generate models of temporal dynamics and their interactions with the environment.

The investigation of the complex of forest fuels holds a pivotal position in fire research, accenting the overriding significance of studying the spatial distribution of CFLs. This exploration is essential for depicting the intricate nature of the fuel complex, comprehending its dynamics, and unraveling spatial variability within the broader ecological framework. Such insights are indispensable for researchers, academics, and forest managers alike, providing a central foundation for making informed decisions regarding the management of forest fuels and the mitigation of forest fires.