Statewide Forest Canopy Cover Mapping of Florida Using Synergistic Integration of Spaceborne LiDAR, SAR, and Optical Imagery

Abstract

Southern U.S. forests are essential for carbon storage and timber production but are increasingly impacted by natural disturbances, highlighting the need to understand their dynamics and recovery. Canopy cover is a key indicator of forest health and resilience. Advances in remote sensing, such as NASA’s GEDI spaceborne LiDAR, enable more precise mapping of canopy cover. Although GEDI provides accurate data, its limited spatial coverage restricts large-scale assessments. To address this, we combined GEDI with Synthetic Aperture Radar (SAR), and optical imagery (Sentinel-1 GRD and Landsat–Sentinel Harmonized (HLS)) data to create a comprehensive canopy cover map for Florida. Using a random forest algorithm, our model achieved an R² of 0.69, RMSD of 0.17, and MD of 0.001, based on out-of-bag samples for internal validation. Geographic coordinates and the red spectral channel emerged as the most influential predictors. External validation with airborne laser scanning (ALS) data across three sites yielded an R² of 0.70, RMSD of 0.29, and MD of −0.22, confirming the model’s accuracy and robustness in unseen areas. Statewide analysis showed lower canopy cover in southern versus northern Florida, with wetland forests exhibiting higher cover than upland sites. This study demonstrates the potential of integrating multiple remote sensing datasets to produce accurate vegetation maps, supporting forest management and sustainability efforts in Florida.

1. Introduction

Southern U.S. forests are among the most productive ecosystems globally in terms of biomass production and carbon sequestration, representing a significant portion of the U.S. terrestrial carbon sink [1]. In addition to supporting over 60% of the nation’s timber production [2], these forests harbor a highly diverse range of endemic species [3]. However, these forests are frequently exposed to natural disturbances, particularly hurricanes, posing significant threats to their structure, recovery processes, and sustainability [4,5,6]. Moreover, climate change is exacerbating these challenges by increasing the frequency and severity of hurricanes [7] as well as other disturbances such as wildfires and droughts [8,9,10], further amplifying their impact on forest dynamics. These disturbances can cause immediate structural damage and initiate long-term ecological shifts, emphasizing the urgent need for comprehensive research on forest recovery, resilience, and adaptation strategies [11].

Forest canopy cover serves as a valuable metric for evaluating the distribution and resistance of forest ecosystems [12,13,14]. It not only provides a means of evaluating immediate damage after disturbances but also offers critical insights into the recovery trajectories of forests over time [15]. For instance, Senf et al. [16] demonstrated that canopy cover can serve as an indicator for predicting the duration required for ecosystems to recover to their pre-disturbance conditions. Furthermore, integrating historical canopy cover data with records of extreme weather events allows researchers to develop predictive models that forecast the potential effects of future climatic disturbances [17]. Such predictive insights are essential for informing land use planning and management strategies, thereby enhancing the resilience of these vital natural habitats against ongoing climatic challenges [18].

Nevertheless, the effectiveness of such strategic efforts is currently hindered by the limited accuracy of existing canopy cover maps [19]. Canopy cover maps derived from optical sensors lack detailed information on the vertical structure of forest canopies [20]. Specifically, sensors such as Sentinel-2 provide comprehensive wall-to-wall visual coverage but fail to capture essential metrics such as canopy height and density, which are critical for understanding forest recovery and resilience [21]. While optical sensors deliver valuable spatial coverage, their limited precision restricts their usefulness for monitoring both current and future forest dynamics, particularly in the face of disturbances like hurricanes and wildfires [22,23,24].

Airborne and terrestrial light detection and ranging (LiDAR) sensors offer a promising approach to mapping forest structure, including canopy cover, forest height, and tree density [25,26]. However, despite its superior precision, the application of LiDAR is typically restricted to smaller-scale projects due to its high acquisition costs and the scarcity of ground-based (in situ) data needed for calibration and validation [20]. These constraints hinder its widespread use for regional forest monitoring, making it less accessible for continuous, wall-to-wall assessments of forest dynamics at broader scales.

Spaceborne LiDAR sensors, particularly NASA’s Global Ecosystem Dynamics Investigation (GEDI) mission, launched to the International Space Station (ISS) in 2018 [27,28], present new opportunities to monitor forest canopy cover changes over time [27]. GEDI addresses some of the limitations of traditional sensors by providing detailed vertical structure measurements at a global scale. This system, designed to capture forest vertical structure within a ~25 m diameter footprint, offers an innovative means of estimating canopy cover across diverse forest ecosystems [27,29,30].

Although GEDI provides global forest canopy cover estimates [27], it does not generate multi-temporal wall-to-wall maps of either forest structure or aboveground biomass density (AGBD) at the fine spatial resolution required to effectively support forest management and conservation in southern U.S. forests. To overcome this limitation, GEDI can be integrated with other sensors, such as optical Landsat and Sentinel Harmonized (HLS) [31,32] and synthetic aperture radar (SAR) data [30,33], like Sentinel 1 [34,35,36] to produce high-resolution, multi-temporal forest canopy cover maps [30]. This multi-sensor fusion approach has been successfully applied to map forest height [37,38] and AGBD [37,39] in other regions, but has not yet been used for mapping canopy cover, either in Florida or elsewhere.

Despite the significant impact of hurricanes on Southern U.S. forests, detailed studies examining how such disturbances affect forest canopy cover remain scarce. This knowledge gap highlights the urgent need for targeted research to understand the implications of these disturbances on canopy cover and to inform adaptive management strategies aimed at enhancing forests’ resilience and sustainability. By focusing on Florida—a state frequently impacted by hurricanes—this research will lay the groundwork to provide a comprehensive analysis of canopy cover dynamics, establishing a framework that can be applied to other forest ecosystems globally.

The primary objective of this study was to develop a statewide forest canopy cover map for Florida for the year 2021, leveraging a multi-sensor fusion approach using GEDI, Sentinel-1 GRD, and HLS datasets. This integrated framework enabled precise mapping, offering detailed insights into the distribution of canopy cover across different regions. The outcomes of this research provide a valuable tool for shaping local adaptive management strategies. For example, counties can implement policies to preserve urban tree cover, enhance reforestation efforts, or improve stormwater management practices to mitigate the effects of canopy loss [40]. By delivering accurate, up-to-date information on canopy cover, this study enables more informed decision-making to support the management and conservation of forest resources. As the southern USA faces escalating risks from larger and more intense hurricanes and other climate-induced disturbances, this research enhances our understanding of forest resilience and recovery patterns. Ultimately, this study will contribute to broader efforts aimed at ensuring the sustainability of southern U.S. forests in the face of ongoing climatic challenges, providing a framework that can be applied to both local and global forest conservation and management initiatives.

2. Materials and Methods

2.1. Study Area

Our study area included all forested areas in the state of Florida (Figure 1). Florida hosts a diverse array of forest types which are categorized by the Florida natural areas inventory (FNAI) into wetland and upland coniferous, hardwood, and mixed forests. This broad classification helps in assessing and managing the state’s diverse forest ecosystems [41]. Among conifers, prominent species include loblolly pine (Pinus taeda L.), longleaf pine (Pinus palustris Mill.), shortleaf pine (Pinus echinata Mill.), and slash pine (Pinus elliottii Engelm.), which collectively account for approximately 45% of the forested landscape [42]. About 42% of Florida’s forests are cultivated and managed for commercial timber production [2]. On the hardwood side, Florida supports a highly diverse range of species, including red maple (Acer rubrum L.), sweetgum (Liquidambar styraciflua), live oak (Quercus virginiana Mill.), turkey oak (Quercus laevis Walter), American beech (Fagus grandifolia), and pignut hickory (Carya glabra), among others [43].

Florida is characterized by a hot and humid climate, with temperatures often exceeding 32 °C for up to six months each year. The state also experiences high levels of relative humidity, typically exceeding 50%. Annual precipitation is among the highest in the nation, with an average annual rainfall of about 1524 mm. However, temporal precipitation patterns are highly variable, leading to frequent floods in some years and droughts in others [44]. Topography is shaped primarily by the processes of deposition and erosion, influenced by fluctuations in sea level and various erosion mechanisms. Known for its notably flat landscape, the state’s highest point is Britton Hill in Walton County, which reaches an elevation of only 105 m above sea level [45]. The composition of Florida’s soils varies widely across the state, with differences in sand, clay, and organic content, as well as the depth of soil overlying the limestone bedrock. These variations are crucial in determining the types of plant communities that thrive in different regions [45].

2.2. Remote Sensing Data Collection and Processing

2.2.1. NASA GEDI Data

The GEDI sensor provided the foundational estimates of canopy cover for this study. For this analysis, GEDI Level 2B (L2B) data [27] collected from 1 April to 30 April 2021, were used. This dataset was derived from waveforms at each 25 m diameter footprint and included the key biophysical metric of canopy cover, calculated from the directional gap probability profile of the L1B waveform [27]. To focus exclusively on forested areas, the Dynamic World dataset [46] was employed as a mask layer, isolating only tree-covered pixels, including wildland–urban interface areas. Additionally, to mitigate the influence of roads on the dataset, a buffer zone of 100 m on each side of the road was established using spatial data for all identified roads in Florida, obtained from the Florida Department of Transportation [47]. This precaution ensured the exclusion of areas affected by roads, preventing potential data distortion. Each selected footprint was filtered using a quality flag (=1) to ensure data integrity, retaining only measurements that met the highest quality standards. To further enhance the reliability of the measurements, we focused on footprints from nocturnal full-power beam transmissions, which are known to reduce atmospheric noise and interference [29].

To evaluate the spatial dependency between GEDI footprints, a semivariogram analysis was conducted to quantify the similarity of data points over increasing distances. A spatial lag distance of 1.24 km between footprints was determined as the range where spatial correlation diminished, ensuring an optimal distribution of the points [48]. This choice was based on the semivariogram curve, which showed a clear plateau at this distance, indicating minimal spatial correlation beyond this range. As a result, we selected footprints that were at least 1.24 km apart, because points closer than this distance would exhibit high similarity and thus provide less additional information. The selected canopy cover footprints were then grouped into 20 classes according to the percentage of canopy cover, ensuring a broad representation of canopy cover from low to high. Samples of footprints within each class were selected, thus capturing the full range of variability in canopy cover across the study area, rather than focusing on just a narrow subset (e.g., exclusively high or low canopy cover). After applying all filtering criteria, the final dataset consisted of approximately 9500 valid footprints, which spanned the entirety of the study area, ensuring comprehensive coverage for the analysis.

2.2.2. Airborne Laser (ALS) Scanning Data

ALS data used for external validation of the forest canopy cover maps were acquired from the United States Geological Survey (USGS) for the state of Florida [47]. This dataset was part of the comprehensive Florida Peninsula LiDAR collection and Hurricane Michael supplemental collection, with data collected between November 2018 and April 2020. The ALS dataset had a resolution of 2.5 feet (0.762 m), offering high-precision elevation data for the entire state of Florida.

Flight missions captured detailed point cloud data, with attributes such as flight dates and other metadata accessible through the Florida Geographic Information Office. These flight dates were crucial for selecting the parks used for validation (Apalachicola National Forest, Green Swamp Wilderness Preserve, and Babcock Ranch Preserve). Notably, these parks had ALS data collected from the end of 2019 until April 2020, aligning closely with the GEDI footprint collection season. This minimized the temporal discrepancies and allowed a more direct comparison between the data fusion model outputs and the ALS data, ensuring that differences in canopy metrics were not affected by changes in forest structure over time.

The canopy cover metric from ALS was computed using a height threshold of 1 m on a 30 m spatial grid, using the lidR package [49] according to Equation (1), aligning with the resolution of the canopy cover model. The ALS-derived canopy cover layer was then used to generate mosaics for each external validation site.

CC_ALS = (NFR_{HTS}/TFR) × 100

(1)

where CC_ALS is the ALS-derived canopy cover (%), NFR_{HTS} is the number of first returns above the 1 m height threshold (HTS = 1), and TFR is the total number of first returns.

2.2.3. Sentinel-1 GRD Data

Sentinel-1A, a C-band Synthetic Aperture Radar (SAR) satellite launched by the European Space Agency, is equipped with dual-polarization capabilities: vertically transmitted and vertically received (VV), and vertically transmitted and horizontally received (VH). The satellite operates at a center frequency of 5.405 GHz, providing a spatial resolution between 10 and 40 m, depending on the imaging mode. With a 12-day revisit cycle, Sentinel-1A offers consistent coverage, making it ideal for temporal analyses over large areas [50].

We utilized Sentinel-1 GRD (Ground Range Detected) data, a processed data product derived from Sentinel-1A. The GRD product is radiometrically calibrated, terrain-corrected, and projected onto a map coordinate system, providing backscatter coefficients (sigma0) that are consistent and reliable for analysis. Data from the Interferometric Wide (IW) swath mode with a spatial resolution of approximately 30 m (aligning closely with the GEDI footprint size) were sourced from the GEE platform [51]. The dataset spanned the period from 1 March to 30 April 2021, ensuring full temporal coverage during the study period.

Pre-processing of the GRD data on GEE included radiometric calibration, terrain correction, and conversion to backscatter coefficients. A seamless mosaic was produced by combining all 75 available Sentinel-1 images into a composite, capturing the entire area of interest without temporal gaps. Additionally, four vegetation indices were derived from Sentinel-1 GRD data, including the Radar Vegetation Index and Co-polarization Ratios 1, 2, and 3. These indices enhanced the analysis of vegetation structure and condition, contributing to the robustness of the study (

Table 1. Sentinel-1 GRD SAR bands and the corresponding vegetation indices, including their respective acronyms and associated wavelengths or equations.

Metric	Acronym	Equation/Wavelength	Reference
Sentinel-1 GRD
Vertical–vertical; vertical–horizontal polarization	σVV and σVH	5.405 GHz	[52]
Radar Vegetation Index	RVI	$\sqrt{{\displaystyle \frac{\mathsf{\sigma }\mathrm{VV}}{\left(\mathsf{\sigma }\mathrm{VV}+\sigma \mathrm{VH}\right)\times \frac{\mathsf{\sigma }\mathrm{VV}}{\mathsf{\sigma }\mathrm{VH}}}}}$	[53]
Co-polarization Ratio	COPOL	${\displaystyle \frac{\mathrm{V}\mathrm{V}}{\mathrm{V}\mathrm{H}}}$	[54]
Co-polarization Ratio 2	COPOL 2	$\left(\mathrm{V}\mathrm{V}-\mathrm{V}\mathrm{H}\right)-\left(\mathrm{V}\mathrm{V}+\mathrm{V}\mathrm{H}\right)$	[55]
Co-polarization Ratio 3	COPOL 3	${\displaystyle \frac{\mathrm{V}\mathrm{H}}{\mathrm{V}\mathrm{V}}}$	[56]

).

2.2.4. Harmonized Landsat and Sentinel-2 Dataset

The NASA Harmonized Landsat and Sentinel-2 (HLS) project integrates surface reflectance (SR) data from the Operational Land Imager (OLI) aboard Landsat 8 and the Multi-Spectral Instrument (MSI) on Sentinel-2 [57]. By combining data from these two sensors, HLS provides global land observations every 2 to 3 days at a spatial resolution of 30 m, enhancing temporal and spatial monitoring capabilities. The HLS observations are processed using sophisticated algorithms that ensure seamless integration between the OLI and MSI sensors. These algorithms include atmospheric correction, cloud masking, spatial co-registration, and spectral normalization. The HLS dataset encompasses over 10 million km², including data from Landsat 8 (since 2013), Landsat 9, and Sentinel-2 (since 2015) [58].

For this study, we used the HLSL30 data product with 30 m spatial resolution covering the period from 1 March to 30 April 2021. This time frame was selected to align with the GEDI footprint collection period, ensuring temporal consistency during the fusion process. To maintain data quality, we applied the Dynamic World filter to focus only on forested areas and applied additional cloud filtering to include only images with less than 30% cloud cover, using quality flags provided in the HLS data product to exclude cloud-contaminated pixels. Additionally, we addressed potential temporal misalignments by selecting the images closest in time to the GEDI acquisitions within the specified period, ensuring that the spectral data matched the conditions under which the GEDI data were collected.

We applied the Dynamic World filter to focus exclusively on forested areas, and road and water masks were applied to exclude non-forested regions and mitigate the potential influence of anthropogenic and aquatic features on our model results. These preprocessing steps ensured that the final dataset was both temporally and spatially aligned with the GEDI data and suitable for analysis.

Finally, we computed thirteen vegetation indices derived from the HLS dataset to further assess vegetation structure and forest health (see

Table 2. HLS spectral bands and corresponding vegetation indices, including their respective acronyms and associated wavelengths or equations.

Metric	Acronym	Equation/Wavelength	Reference
HLS
Blue, green, red, near-infrared, short-wave infrared, short-wave infrared 2	B, G, R, NIR, SWIR1, SWIR2	0.490 µm, 0.560 µm, 0.665 µm, 0.842 µm, 1.610 µm, 2.190 µm	[57]
Normalized Difference Vegetation Index	NDVI	${\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}\mathrm{ }-\mathrm{ }\mathrm{R}}{\mathrm{N}\mathrm{I}\mathrm{R}\mathrm{ }+\mathrm{ }\mathrm{R}}}$	[59]
Enhanced Vegetation Index	EVI	$2.5\times \left({\displaystyle \frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+6\times \mathrm{R}-7.5\times \mathrm{B}+1}}\right)$	[60]
Soil-Adjusted Vegetation Index	SAVI	$1.5\times \left({\displaystyle \frac{\mathrm{NIR}-\mathrm{R}}{\mathrm{NIR}+\mathrm{R}+0.5}}\right)$	[61]
Modified -Adjusted Vegetation Index	MSAVI2	${\displaystyle \frac{\left(2\times \mathrm{NIR}+1\right)-\sqrt{{\left(2\times \mathrm{NIR}+2\right)}^2-8\times \left(\mathrm{NIR}-\mathrm{R}\right)}}{2}}$	[62]
Linear Spectral Unmixing	LSU	-	[63]
Soil	FSOIL	$0.140\times \mathrm{B}+0.160\times \mathrm{G}+0.220\times \mathrm{R}+0.390\times \mathrm{NIR}$ $+0.450\times \mathrm{SWIR}1\times 0.270\times \mathrm{SWIR}2$
Water	FWATER	$0.070\times \mathrm{B}\times 0.039\times \mathrm{G}+0.023\times \mathrm{R}+0.031\times \mathrm{N}\mathrm{I}\mathrm{R}+0.011\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+0.007\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2$
Vegetation	FVEG	$0.086\times \mathrm{B}+0.062\times \mathrm{G}+0.043\times \mathrm{R}+0.247\times \mathrm{N}\mathrm{I}\mathrm{R}+0.109\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+0.039\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2$
Simple Ratio Index	SRI	${\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{R}}}$	[64]
Normalized Difference Water Index	NDWI	${\displaystyle \frac{\mathrm{G}+\mathrm{N}\mathrm{I}\mathrm{R}}{\mathrm{G}-\mathrm{N}\mathrm{I}\mathrm{R}}}$	[65]
Green Chlorophyll Index	GCI	$\left({\displaystyle \frac{\mathrm{NIR}}{\mathrm{G}}}\right)-1$	[66]
Wide Dynamic Range Vegetation Index	WDRVI	${\displaystyle \frac{\left(0.1\times \mathrm{N}\mathrm{I}\mathrm{R}-0.2\times \mathrm{R}\right)}{0.1\times \mathrm{NIR}+0.2\times \mathrm{R}+0.5}}$	[67]
Global Vegetation Moisture Index	GVMI	${\displaystyle \frac{\left(\mathrm{N}\mathrm{I}\mathrm{R}+0.1\right)-\left(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1-0.02\right)}{\left(\mathrm{N}\mathrm{I}\mathrm{R}+0.1\right)+\left(\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+0.02\right)}}$	[68]
Chlorophyll Vegetation index	CVI	$\mathrm{NIR}\times \left({\displaystyle \frac{\mathrm{R}}{{\mathrm{G}}^2}}\right)$	[69]
Clay Minerals Ratio	CMR	${\displaystyle \frac{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1}{\mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2}}$	[70]
Kernel-Normalized Difference Vegetation Index	KNDVI	$\tanh \left({\left({\displaystyle \frac{\mathrm{N}\mathrm{I}\mathrm{R}-\mathrm{R}}{2\mathsf{\sigma }}}\right)}^2\right)$	[71]

for the complete list of indices).

2.2.5. Ancillary Imagery

Critical terrain metrics such as elevation, slope, and aspect were extracted from the NASADEM_HGTv001 dataset [72] and used in this study. NASADEM, an enhancement of earlier DEMs, utilizes refined Shuttle Radar Topography Mission (SRTM) data. The SRTM project was a collaborative effort involving NASA, the National Geospatial-Intelligence Agency (NGA), the German Aerospace Centre (DLR), and the Italian Space Agency (ASI). The reprocessing of original SRTM radar data incorporated advanced algorithms and additional elevation data from ASTER and ICESat’s Geoscience Laser Altimeter System (GLAS), significantly reducing data voids and improving accuracy [57,58]. NASADEM provides global elevation data with a resolution of 1 arc-second, covering all land areas between 60°N and 56°S latitude, representing approximately 80% of Earth’s land surface [72].

The inclusion of terrain metrics such as elevation, slope, and aspect help to account for the impact of topography on vegetation growth and spatial distribution, which is essential for understanding canopy composition and structure [73]. Although Florida’s elevation range is relatively narrow (0–105 m), even slight variations can influence vegetation patterns by affecting temperature, moisture availability, and atmospheric pressure [74]. For instance, coastal and near-coastal slopes, though they may be subtle, play a role in water drainage and retention, in turn affecting vegetation density and distribution [74,75]. Additionally, aspect influences sunlight exposure and moisture retention; north-facing slopes may retain slightly more moisture than south-facing ones even in low-elevation areas [76]. These metrics allowed the model to capture terrain-driven variations in canopy cover effectively. Geographic coordinates such as the longitude and latitude of each GEDI pixel were also used as predictor variables.

2.2.6. Spatial Transformations

For each studied variable, we calculated the maximum, minimum, mean, and standard deviation using a kernel-based spatial reclassification method within a 3 × 3 pixel moving window (

Table 3. Spatial transformations applied in this study, along with their corresponding acronyms and mathematical equations.

Spatial Transformation	Acronym	Equation	Reference
GLCM
Angular Second Moment	Metric + “_ASM”	${\sum }_{\mathrm{i}}{\sum }_{\mathrm{j}}{\left\{\mathrm{p}\left(\mathrm{i},\mathrm{j}\right)\right\}}^2$	[77]
Contrast	Metric + “_CONT”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}\mathrm{p}\left(\mathrm{i},\mathrm{j}\right){\left(\mathrm{i}-\mathrm{j}\right)}^2$
Correlation	Metric + “_COR”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}\left(\mathrm{ij}\right)\mathrm{p}(\mathrm{i},\mathrm{j})-{\mathsf{\mu }}_{\mathrm{x}}{\mathsf{\mu }}_{\mathrm{y}}/{\mathsf{\sigma }}_{\mathrm{x}}{\mathsf{\sigma }}_{\mathrm{y}}$
Variance	Metric + “_VAR”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}{\left(\mathrm{i}-\mathsf{\mu }\right)}^2\mathrm{p}\left(\mathrm{i},\mathrm{j}\right)$
Inverse Difference Moment	Metric + “_IDM”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}(1/(1+{\left(\mathrm{i}-\mathrm{j}\right)}^2\left)\right)\mathrm{p}(\mathrm{i},\mathrm{j})$
Sum Variance	Metric + “_SVAR”	$\sum _{\mathrm{i}=2}^{2{\mathrm{N}}_{\mathrm{g}}}{(\mathrm{i}-\mathrm{s}\mathrm{e}\mathrm{n}\mathrm{t})}^2{\mathrm{p}}_{\mathrm{x}+\mathrm{y}}\left(\mathrm{i}\right)$
Sum Average	Metric + “_SVG”	$\sum _{\mathrm{i}=2}^{2{\mathrm{N}}_{\mathrm{g}}}{\mathrm{i}\mathrm{p}}_{\mathrm{x}+\mathrm{y}}\left(\mathrm{i}\right)$
Sum Entropy	Metric + “_SENT”	$-\sum _{\mathrm{i}=2}^{2{\mathrm{N}}_{\mathrm{g}}}{\mathrm{p}}_{\mathrm{x}-\mathrm{y}}\left(\mathrm{i}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left\{{\mathrm{p}}_{\mathrm{x}+\mathrm{y}}\right(\mathrm{i}\left)\right\}$
Entropy	Metric + “_ENT”	$-\sum _{\mathrm{i}}\sum _{\mathrm{j}}\mathrm{p}\left(\mathrm{i},\mathrm{j}\right)\log \left(\mathrm{p}\left(\mathrm{i},\mathrm{j}\right)\right)$
Difference Variance	Metric + “_DVAR”	Variance of ${\mathrm{p}}_{\mathrm{x}-\mathrm{y}}$
Difference Entropy	Metric + “_DENT”	$\sum _{\mathrm{i}=0}^{{\mathrm{N}}_{\mathrm{g}-1}}{\mathrm{p}}_{\mathrm{x}-\mathrm{y}}\left(\mathrm{i}\right)\mathrm{l}\mathrm{o}\mathrm{g}\left\{{\mathrm{P}}_{\mathrm{x}-\mathrm{y}}\left(\mathrm{i}\right)\right\}$
Informal Measures of Correlation 1	Metric + “_IMCORR1”	$(\mathrm{HXY}-\mathrm{HXY}1)/\max \{\mathrm{XH},\mathrm{XY}\}$
Informal Measures of Correlation 2	Metric + “_IMCORR2”	${\left(1-\exp \left[2\left(\mathrm{H}\mathrm{X}\mathrm{Y}2-\mathrm{H}\mathrm{X}\mathrm{Y}\right]\right)\right)}^{0.5}$
Maximum Correlation Coefficient	Metric + “_MAXCOR”	$\sum _{\mathrm{k}}{\displaystyle \frac{\mathrm{p}\left(\mathrm{i},\mathrm{k}\right)\mathrm{p}\left(\mathrm{j},\mathrm{k}\right)}{{\mathrm{p}}_{\mathrm{x}}\left(\mathrm{i}\right){\mathrm{p}}_{\mathrm{y}}\left(\mathrm{k}\right)}}$
Dissimilarity	Metric + “_DISS”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}\left|\mathrm{i}-\mathrm{j}\right|\mathrm{p}(\mathrm{i},\mathrm{j})$
Inertia	Metric + “_INERTIA”	$\sum _{\mathrm{i}}\sum _{\mathrm{j}}{\left(\mathrm{i}-\mathrm{j}\right)}^2\mathrm{p}\left(\mathrm{i},\mathrm{j}\right)$
Shade	Metric + “_SHADE”	$\sum _{\mathrm{i}=0}^{{\mathrm{N}}_{\mathrm{g}-1}}{\left\{\mathrm{i}+\mathrm{j}-{\mathsf{\mu }}_{\mathrm{x}}-{\mathsf{\mu }}_{\mathrm{y}}\right\}}^3\mathrm{p}(\mathrm{i},\mathrm{j})$
Cluster Prominence	Metric + “_PROM”	$\sum _{\mathrm{i}=0}^{{\mathrm{N}}_{\mathrm{g}-1}}{\left\{\mathrm{i}+\mathrm{j}-{\mathsf{\mu }}_{\mathrm{x}}-{\mathsf{\mu }}_{\mathrm{y}}\right\}}^4\mathrm{p}(\mathrm{i},\mathrm{j})$
Moving window
Maximum	Metric + “_MAX”	$\underset{\mathrm{i},\mathrm{j}\in \mathrm{win}}{\max }\mathrm{F}(\mathrm{i},\mathrm{j})$	[80]
Minimum	Metric + “_MIN”	$\underset{\mathrm{i},\mathrm{j}\in \mathrm{win}}{\min }\mathrm{F}(\mathrm{i},\mathrm{j})$
Mean	Metric + “_MEAN”	$1/\mathrm{n}\sum _{\mathrm{i},\mathrm{j}\in \mathrm{win}}\mathrm{F}(\mathrm{i},\mathrm{j})$
Standard Deviation	Metric + “_STD”	$\sqrt{1/\mathrm{n}\sum _{\mathrm{i},\mathrm{j}\in \mathrm{win}}{\left(\mathrm{F}\right(\mathrm{i},\mathrm{j})-\mathsf{\mu }(\mathrm{x},\mathrm{y}))}^2}$

Where P_(i,j) represents the probability of co-occurrence of pixel pairs with gray levels i and j within a specified window. This is derived from the normalized GLCM, where each entry (i,j) reflects the likelihood of gray levels i and j appearing together in a given spatial relationship. i,j are indices corresponding to the gray tone levels in the GLCM, and also serve as coordinates in the matrix, where each entry at position (i,j) indicates the co-occurrence probability for gray tones i and j. Ng is the number of distinct gray levels in the quantized image, which is derived from reducing the original image’s gray value range to a fixed number of levels for analysis. The mean (μ) and standard deviation (σ) of P_(i,j) summarize the average co-occurrence probability and its variability, respectively [80]. px(i) represents the probability of gray level i appearing as the first value in the co-occurring pixel pairs; py(j) represents the probability of gray level j appearing as the second value in the co-occurring pixel pairs; HX and HY are entropies of px and py; px + y is the sum distribution, representing the probability of the sum of gray levels i + j; px − y refers to the difference distribution, representing the probability of the difference between gray levels i − j; H_XY is the joint entropy of the GLCM; and t quantifies the randomness or uncertainty in the entire co-occurrence matrix. H_XY1 is the cross-entropy based on row and column marginal probabilities; H_XY2 is defined as the cross-entropy based on the joint marginal probabilities of px + y; p refers to the marginal probabilities of gray levels; and S_ent is the sum entropy, measuring the randomness of the sum of gray levels [76]. The notation “i, j ∈ win” denotes that the indices i and j are within a specified window “win” in an image. The term “win” specifies the boundary or the extent of this local neighborhood, and the operations (like max, min, or average) are conducted exclusively over the pixels that fall within this window. F_(i,j) represents a function applied to each pixel (or pixel pair) located at position (i,j) in the image.

). Additionally, we employed the Gray-Level Co-Occurrence Matrix (GLCM) to extract statistical features from the images’ grayscale values for each band and vegetation index, following methodologies established previously [77]. Incorporating the bands, vegetation indices, ancillary data, kernel moving window, and GLCM features, our dataset expanded to comprise a total of 256 predictor variables. These variables were then extracted and associated with the corresponding GEDI footprints to effectively model canopy cover. To ensure compatibility with GEDI geolocation uncertainty [78,79], all predictive metrics were calculated to a spatial resolution of 30 × 30 m. This resolution aligned with the Sentinel-1 and HLS datasets and minimized the potential impact of spatial mismatches when integrating multi-source data.

2.3. Canopy Cover Modeling and Assessment

Random forest (RF) [81] regression modeling was employed to create canopy cover models. RF modeling was chosen for its robustness, accuracy, and transferability, characteristics that are particularly beneficial for handling large and complex datasets typical of remote sensing data [82]. This approach allowed effective handling of non-linear relationships and interactions between variables. Initially, sampling was carried out, during which 1000 samples containing spatially independent GEDI-derived canopy cover estimates were selected as the response variables, and image-derived metrics were selected as predictor variables for model calibration. The dataset was then divided into training and testing subsets, with 70% of the data allocated for training and the remaining 30% used for testing. This split helped to maximize learning from the data while ensuring the model’s generalizability to unseen data (Figure 2).

An RF model was implemented using the RandomForestRegressor from the python sklearn.ensemble [83] library (version 1.6.0), configured with 250 trees to balance model complexity and computational efficiency. To ensure that only the most relevant features were included, the SelectFromModel method was applied to the feature set to rank the features based on their contribution to reducing the model’s impurity (i.e., how well a feature contributed to reducing variance or classification error). Features that met or exceeded the mean importance score were retained, reducing dimensionality and enhancing the model’s accuracy and interpretability. After selecting the most informative features, their importance scores were recorded for further analysis.

To evaluate model robustness and quantify uncertainty, bootstrapping was conducted with 100 iterations. In each iteration, the model was trained on the stratified resampled subset of the data, with replacement, allowing for variations in the sample while preserving overall model accuracy. The model’s performance was assessed using key metrics, including the coefficient of determination (R²), mean difference (MD), and root mean squared difference (RMSD), calculated for both absolute and relative terms (Equations (2)–(6)). These metrics provided a comprehensive understanding of the model’s accuracy and its capacity to predict canopy cover across the study area.

$$R^2={\displaystyle \frac{{\left[\sum _{i=1}^n\left(y_i-\overline{\mathrm{y}}\right)\left({\widehat{y}}_i-\overline{\widehat{y}}\right)\right]}^2}{\sum _{i=1}^n{{(\mathrm{y}}_i-\overline{\mathrm{y}})}^2\sum _{i=1}^n{\left({\widehat{\mathrm{y}}}_i-\overline{\widehat{y}}\right)}^2}}$$

(2)

$$RMSD=\sqrt{{\displaystyle \frac{{{\sum }_{i=1}^n\left({\widehat{y}}_i-y_i\right)}^2}{n}}}$$

(3)

$$RMSD\%={\displaystyle \frac{RMSD}{\widehat{y}}}\times 100$$

(4)

$$MD={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left({\widehat{y}}_i-y_i\right)$$

(5)

$$MD\%={\displaystyle \frac{MD}{\widehat{y}}}$$

(6)

where n is the number of observations; y_i and ${\widehat{y}}_i$ are the observed and predicted values for canopy cover for the i-th observation; $\overline{y}$ and $\overline{\widehat{y}}$ are the observed and predicted mean values across all observations.

After completing all the bootstrap iterations, variable importance analysis was conducted to evaluate each predictor’s contribution to the model. During each run, variables were assessed based on their ability to enhance model accuracy, quantified by the mean decrease in impurity, which quantified how effectively a variable split the data into more homogeneous groups, reducing uncertainty or disorder in the predictions [81]. The cumulative mean decrease in impurity was then calculated for each variable across all bootstrap iterations. Final importance scores were determined by evaluating how frequently each variable was selected during the iterations, weighted by its contribution to the model’s performance. These scores were converted to percentages, ranging from 0% to 100%, to allow easier comparison and ranking. This process highlighted the most significant predictors, ensuring that the final model emphasized the variables with the greatest predictive power, thus improving both the model’s generalizability and its accuracy in predicting canopy cover.

2.4. Wall-to-Wall Canopy Cover Mapping

The previously trained RF models were initially trained as regressors to predict continuous canopy cover values. For implementation on GEE [84], these models were adapted into a classifier format, enabling canopy cover predictions at a 30 × 30 m pixel resolution across Florida. Using 100 models generated through bootstrap iterations, this classifier was applied statewide on a 30 × 30 m grid, providing detailed spatial predictions of canopy presence and facilitating precise geospatial analysis within GEE. The decision to use 100 bootstrap iterations was based on findings from previous studies, such as Rosenfeld et al. [85], which demonstrated that 100 iterations strike an optimal balance between computational efficiency and model stability. Additionally, preliminary testing in this study confirmed that increasing the number of iterations beyond 100 did not lead to significant improvements in the stability or accuracy of the model’s outputs (e.g., R², RMSD, and MD), further supporting the use of 100 iterations as a reliable and efficient choice.

The output maps from each of the 100 models were subsequently stacked to compute the mean and standard deviation across all models. This process resulted in the creation of a final upscaled map, representing the average predicted values across all 100 models, and a standard deviation map, capturing the variability or uncertainty in these predictions. These maps, as delineated by Equations (7) and (8), represent the spatial variability and consistency of the modeled data across Florida. The mean map illustrates the central tendencies of the modeled variables, while the standard deviation map identifies areas with greater uncertainty, indicating where the model’s predictions were less consistent.

$${\widehat{x}}_{i,j}={\displaystyle \frac{1}{M}} \sum _{k=1}^M{x}_{i, j, k}$$

(7)

where x_i,j,k is the value at grid cell (i,j) in the k-th map; i, j are the indices representing the position of the cell in the grid; k is the index of the map (k = 1,…,M), one for each bootstrapped model. M is the total number of bootstrap maps.

$${\sigma }_{i,j}{\sqrt{{\displaystyle \frac{1}{99}} \sum _{k=1}^M(x_{i,j,k}-{\widehat{x}}_{i,j})}}^2$$

(8)

where ${\sigma }_{i,j}$ is the average standard deviation value for cell (i,j); x_i,j,k is the value at grid cell (i,j) in the k-th map; ${\widehat{x}}_{i,j}$ is the mean value for the grid cell (i,j) across all bootstraps. The factor 99 in the denominator corresponds to the degree of freedom.

2.5. External Validation

To evaluate the performance of our model on an independent dataset, we selected three state parks in Florida for external validation: Apalachicola National Forest (ANF), Green Swamp Wilderness Preserve (GSWP), and Babcock Ranch Preserve (BRP) (Figure 1). These sites were chosen based on the availability of ALS data collected within a similar phenological season, ensuring consistency and avoiding strong discrepancies during validation. Additionally, the sites were selected for their inherent range of canopy cover, providing a robust test of the model’s ability to accurately predict areas of both lesser and greater canopy coverage. The sites are located in Tallahassee, Lake, and Punta Gorda counties, respectively. The ALS data (refer to Section 2.2.2) from these sites were employed for map validation (see Figure 2).

To ensure comprehensive validation, we first assessed the impact of including each park individually in the validation process and then evaluated the model’s performance across all three parks combined. For the validation process, which involved comparing our model-derived maps with the ALS ground-truth data, 2000 pixels measuring 30 × 30 m each were randomly chosen from each site and assessed on a one-to-one basis. This number of pixels was selected as a balance between ensuring robust statistical representation and addressing computational constraints, as processing larger pixel counts would have exceeded the platform limitations. To quantify the model’s precision, accuracy, and agreement with the observed data, we calculated key metrics, including the R² and absolute and relative (%) RMSD and MD.

In addition to validating our model with ALS data, we calculated the mean canopy cover for each individual park (refer to Equations (9)–(11)). Specifically, we estimated the mean value $\widehat{x}$ of the predicted canopy cover, which served as the central tendency of the predictions for each park. Accompanying this mean, the 95% confidence interval provided the range within which the true mean canopy cover was expected to be found. To assess the variability in these predictions, we calculated the total variance $V\left[\widehat{x}\right]$ for each park. From this variance, the standard error (SE) was computed, providing a measure of the average deviation of the estimated means from the true population mean. Additionally, the relative standard error (rSE) was calculated, relating the standard error to the mean canopy cover, offering a relative measure of uncertainty. This was particularly useful for comparing the precision of the model across different parks. These error measurements provide a comprehensive evaluation of the reliability and consistency of the predictions in different geographic contexts, enabling a detailed assessment of the model’s accuracy and precision to ensure that it can be applied effectively throughout the state’s diverse environments.

$$\widehat{x}={\displaystyle \frac{1}{N}}\sum _{g=1}^N\left({\displaystyle \frac{1}{M}}\sum _{k=1}^M{x}_{k,g}\right)$$

(9)

where $\widehat{x}$ represents the mean value for each park, N is the total number of observations in each park, M is the total number of bootstrap maps, g refers to the position of the grid cell, $x_{k,g}$ is the value at position g in the k-th bootstrap map for the respective park.

$$SE=\sqrt{V\left[\widehat{x}\right]}$$

(10)

where $SE$ is the standard error for each park and $V\left[\widehat{x}\right]$ represents the variance.

$$rSE\left(\%\right)={\displaystyle \frac{SE}{\widehat{x}}}\times 100$$

(11)

2.6. Characterization of Forest Canopy over Across the State of Florida

To calculate the canopy cover percentage for each forest type, we integrated the final average canopy cover map generated by our model with the FNAI vegetation classifications [41]. This analysis utilized the vegetation type data from the Florida Department of Environmental Protection (FDEP), obtained from their statewide land use and land cover inventory [86], allowing us to accurately assign vegetation type to the entire extension of the study area. This integration enabled the precise assignment of vegetation types across the entire study area, resulting in a detailed map illustrating canopy cover percentages by forest type. Additionally, we calculated both the mean and standard deviation of canopy cover for each forest type to further analyze the distribution and variability within these categories.

A Kruskal–Wallis test [87] was performed at the 95% confidence level to assess whether there were significant differences in canopy cover among the various forest types studied. The null hypothesis (H₀) assumed that the median canopy cover values across the forest types are equal, while the alternative hypothesis (H₁) suggested that at least one group’s median would be different, indicating a significant variation in canopy cover among the forest types. If significant differences were detected, Dunn’s post-hoc pairwise comparisons [88] with Bonferroni correction were applied to determine which specific forest types significantly differed from one another.

For our county characterization, we clipped the mean canopy cover map to the geographic boundaries of each county to calculate the mean canopy cover per county. To estimate the standard deviation for each county, we used the average standard deviation map (refer to Equation (8)), and by extracting the standard deviation values within each county’s boundaries, we computed the variance for each county using Equation (12). The standard deviation for each county was then obtained by taking the square root of the variance (Equation (13)). This process provides a measure of the uncertainty in the canopy cover predictions at the county level.

$$V_{county}={\displaystyle \frac{1}{N-1}}{\sum _{i=1}^n({\sigma }_i-\widehat{\sigma })}^2$$

(12)

where $V_{county}$ represents the variance in standard deviation per county, N is the number of extracted values of standard deviation within the county, ${\sigma }_i$ is the standard deviation for each pixel in the county, and $\widehat{\sigma }$ is the mean standard deviation across the county.

$${\sigma }_{county}=\sqrt{V_{county}}$$

(13)

3. Results

3.1. Model Performance Assessment

The mean R² of the 100 aggregated models was 0.69, indicating that, on average, 69% of the variance in canopy cover was explained by the model. The RMSD values, reflecting the average prediction errors, were 0.17 (absolute) and 44% (relative) (Figure 3). The MD values were notably low at 0.001 (absolute) and 0.34% (relative), demonstrating satisfactory bias in the predictions.

The range of R² across the models varied from 0.62 to 0.77, indicating moderate variations in model accuracy between iterations. RMSD values oscillated between 0.15 to 0.19 (absolute) and 40.2% to 49.8% (relative), highlighting the variability in error size relative to the actual canopy cover values. Additionally, the MD fluctuated from −0.025 to 0.029 (absolute) and ranged from −6.27% to 7.79% (relative), indicating changes in the direction and magnitude of the prediction bias across different model iterations.

Longitude emerged as the most influential predictor, contributing over 95% to the model importance in the RF models (Figure 4), and was used by the RF in 83 of the 100 models. Following closely, the red and red mean predictors showed significant influence, contributing to 55% of the model variance and being selected in 53 models. Latitude also had considerable impact with a 53% contribution, while red mean had 49% importance.

Lesser but still notable contributions came from green mean and swir2 min, each accounting for 31% of the model’s variability, and green min at 30%. Other predictors varied in their influence, with importance scores ranging from 29% to 10% (e.g., green, swir2 mean, evi min, red savg, and others). For clarity and ease of visualization, Figure 4 includes only those metrics that demonstrated more than 10% in the importance scoring. It is important to note that our approach to modeling canopy cover considered a total of 260 metrics (Table 1, Table 2 and Table 3), but through feature selection, only 50 variables were chosen for the final model.

3.2. Wall-to-Wall Canopy Cover Maps for Validation Sites

The final wall-to-wall canopy cover maps revealed distinct variations in canopy cover across the three validation parks. BRP exhibited a mean canopy cover of 20.76%, indicating relatively sparse vegetation in some areas. In contrast, GSWP showed a higher mean canopy cover of 43.07%, reflecting denser foliage. ANF also displayed a high mean canopy cover of 42.40% with a lower standard deviation of 9.10%, suggesting more uniform coverage. These maps enable a comprehensive evaluation of specific regions within these state parks, offering a detailed perspective beyond the limitations of relying solely on GEDI footprints (Figure 5).

When evaluating the mean canopy cover for each individual validation area, we found that BRP had the lowest mean canopy cover at 20.8%, while ANF and GSWP had higher estimates of canopy cover at 42.4% and 43.9%, respectively (

Table 4. Summary of canopy cover (%) and associated errors for each validation site. The reported canopy cover includes ±2 times the standard error to indicate the 95% confidence interval. GSWP = Green Swamp Wilderness Preserve; BRP = Babcock Ranch Preserve; and ANF = Apalachicola National Forest.

Site	Estimated Canopy Cover Mean (%) ± Standard Error	Standard Error (%)	Relative Standard Error to the Mean (%)
BRP	20.8 ± 36.6	18.3	87.9
ANF	42.4 ± 41.0	20.5	48.3
GSWP	43.9 ± 40.9	20.5	46.7

). The validation revealed that the canopy cover estimates for BRP had the lowest precision, with the highest rSE value (87.9%). In contrast, the canopy maps for ANF and GSWP exhibited greater precision, with rSE of 48.3% and 46.7%, respectively. Overall, these results underscore the variability and precision of our canopy cover estimates, with GSWP showing the most reliable results among the evaluated sites.

As shown in

Table 5. Results of external validation using ALS data. GSWP = Green Swamp Wilderness Preserve; BRP = Babcock Ranch Preserve; ANF = Apalachicola National Forest; all sites = GSWP + BRP + ANF.

Site	R2	RMSD		MD
		Abs. %	Rel. %	Abs. %	Rel. %
GSWP	0.64	0.30	0.30	−0.24	−0.24
BRP	0.61	43.05	42.12	−34.21	−34.71
ANF	0.57	0.26	0.29	−0.17	−0.22
All sites	0.70	64.15	47.31	−44.04	−36.15

, we found that using all three parks for validation yielded a higher R² value (0.70), but it also resulted in a slight increase in errors compared with BRP, with an absolute RMSD of 0.29% and an MD of −0.22%. These results reveal that incorporating a broader range of canopy cover conditions—by including all three parks—allowed the model to explain more variability in the dataset, as indicated by the higher R². However, the increased RMSD and MD indicate that this approach introduced greater complexity into the predictions, possibly due to the model having to account for more extreme or diverse canopy conditions. When the individual parks were examined, GSWP provided the best validation results, with an R² of 0.64, an RMSD of 0.30%, and an MD of −0.24%. This park was the most homogeneous and had the highest canopy cover compared with the other two parks, suggesting that our model performed better in areas with high canopy cover than in areas with lower or mixed canopy cover. Conversely, ANF produced the poorest results among all three parks, with the lowest R² value (0.57), and higher RMSD and MD at 0.30% and −0.24%, respectively. Overall, the validation results indicate underestimation in canopy cover, as evidenced by the negative MD values for all individual sites and the combined all-sites validation. When considering all sites for validation, the model underestimated canopy cover by an average of −0.22, or −36.21% relative to the actual values.

3.3. Spatial Characterization of Canopy Cover Across the State of Florida

Figure 6 presents the prediction and uncertainty of canopy cover across the state of Florida. Observations indicate a higher concentration of canopy cover in the northern regions of the state, with a noticeable decrease towards the south. Notably, areas with more canopy cover also exhibit higher uncertainty, as shown by the red and purple colors on the map. In contrast, the blue areas in the south represent regions with less canopy cover but higher uncertainty, while the yellow areas indicate places with lower uncertainty associated with their canopy cover predictions.

The canopy cover in the state ranges from as low as 2% to as high as 77%. Most of the state exhibits canopy cover between 20% and 30%, followed by a substantial portion with canopy cover ranging from 40% to 54%. Regarding the standard deviation across the 100 model iterations applied to Florida, values ranged from 0.35% to 37%. Higher standard deviations were associated with regions of dense canopy cover, indicating that the model estimated a wide range of canopy cover values for these areas. Figure 6b,c depict how longitude and latitude influence canopy cover for the state. Canopy cover increases from southern to northern latitudes. Forest canopy cover also generally increases longitudinally from east to west, except in areas with less land (from 88°W to 86°W) and regions with cloud occlusion.

Across Floridan forest types, wetland areas generally exhibit higher canopy cover compared with upland forest types. Specifically, wetland hardwood forests include the most canopy cover at approximately 46%, followed closely by mixed wetlands at 45% and wetland coniferous forests at 44%. In contrast, upland forest types display lower rates of canopy cover, with mixed uplands at 42%, upland hardwoods at 41%, and upland coniferous forests at 37% (Figure 7).

The Kruskal–Wallis test revealed a significant difference in canopy cover across all forest types (χ² = 6212.3, df = 5, p < 2.2 × 10⁻¹⁶). For instance, results for the upland hardwood forest (UHF) and upland mixed forest (UMF) did not differ significantly from each other (p = 1.000, both labeled “a”), indicating similar canopy cover. Upland coniferous forest (UCF) included significantly less canopy cover compared with other upland forests (p < 0.001, labeled “b”). Among wetland types, wetland coniferous (WCF) and wetland mixed forest (WMF) showed no significant difference from each other (p = 0.126, both labeled “c”), but were significantly different from all upland types. Wetland hardwood forest (WHF) had the highest rate of canopy cover, and this result was significantly different from all other forest types (p < 0.001, labeled “d”).

Furthermore, the standard deviation of canopy cover indicated considerable variability within these forest types. Wetland hardwood forest, which had the highest canopy cover, also exhibited the highest standard deviation 10%. Other forest types, such as wetland coniferous, upland hardwood, and upland coniferous, show standard deviation values ranging from 9% to 9.7%, indicating relatively high variability. Conversely, mixed wetlands and mixed uplands had the lowest standard deviation values at 6.9% and 8.3%, respectively, suggesting more uniform canopy cover within these types.

Our analysis revealed substantial variation in canopy cover across counties. The counties of Baker, Holmes, and Union lead with the greatest area of average canopy cover, each at 43%; following closely are Bradford, Columbia, Dixie, and Lafayette counties, each with canopy cover of 42%; at the lower end, Broward, Palm Beach, Miami-Dade, and Martin counties have ≤12% canopy cover. Notably, these counties with lower canopy cover are predominantly located in the southern part of the state. The remaining counties exhibit canopy cover ranging from 39% to 14% (Figure 8c).

An observable trend indicates that higher canopy cover is concentrated in the northern, northeastern, and western parts of the state, while central regions show moderate coverage and southern regions have significantly lower rates of tree cover. Additionally, counties such as Monroe, Jefferson, and Miami-Dade exhibit high variability in canopy cover, suggesting diverse landscapes within these areas. In contrast, counties like Gilchrist and Suwannee show minimal variation, indicating more uniform canopy distribution.

4. Discussion

This paper focuses on developing and testing a data fusion framework that combines GEDI-derived canopy cover samples with continuous data from HLS optical imagery, Sentinel-1 GRD radar, and other ancillary data to create a high-resolution, wall-to-wall canopy cover map for the entire state of Florida. While previous research has used similar fusion approaches using GEDI to estimate biomass [29,30,33,89,90] and tree height [37,91,92], this is the first study to adapt such methods specifically for canopy cover, marking a significant advancement in the field. The novelty of this study lies not only in the application of this data fusion technique but also in the reliable, high-accuracy canopy cover maps it produces. Canopy cover is a critical metric for the early detection of forest changes, as it is often the first visible indicator of forest alteration that can be observed remotely. This is particularly relevant in Florida, where forests are regularly impacted by hurricanes, making it essential to monitor canopy cover dynamics over time.

4.1. Canopy Cover Model Performance and Prediction Assessment

The RF model developed and utilized in this study demonstrated robust performance in predicting canopy cover, with an average R² of 0.69 and RMSD values ranging from 0.15 to 0.19. These results are consistent with those reported by [79] and represent an improvement over earlier studies using simpler regression models (e.g., linear regression) for similar approaches [30,93]. To date, linear regression models have been the most widely used methods for AGBD estimation; however, this approach does not effectively capture the complex nonlinear relationships between forest structure/AGB and remotely sensed data [94].

Numerous studies have indicated that algorithms based on decision trees, such as RF, can better estimate forest metrics such as AGB, height, and canopy cover [79,95,96,97,98,99]. The strength of the RF method lies in its ability to accommodate large datasets with numerous predictor variables and its robustness against overfitting, making it particularly well suited for remote sensing applications where complex interactions between variables are common [81]. In this study, the RF model effectively integrated data from diverse sources, including GEDI, HLS, Sentinel-1 GRD, and ancillary datasets, demonstrating its ability to capture the majority of the variability in canopy cover across Florida.

While the model exhibited strong performance overall, moderate variability was observed in key metrics across the bootstrap iterations. For instance, R² values ranged from 0.62 to 0.77, reflecting differences in the model’s ability to explain variance in canopy cover across subsets of data. Similarly, RMSD values ranged from 0.15 to 0.19 (absolute) and 40.2% to 49.8% (relative), highlighting variability in error magnitude relative to actual canopy cover. These fluctuations are likely to have stemmed from the diverse environmental conditions across Florida, including abrupt changes in canopy cover in wetlands or coastal areas, for example, which present additional modeling challenges in the remote sensing realm.

Despite the challenges posed by abrupt changes in vegetation and environmental variability, the narrow ranges observed in the performance metrics indicate the model’s robustness and reliability in predicting canopy cover under varying conditions. For instance, the RF model demonstrated strong performance across Florida’s diverse environmental gradients, underscoring its ability to adapt to the state’s unique and heterogeneous forest ecosystems. While this suggests potential applicability to other regions with varied ecological characteristics, further analysis is needed to confirm its ability to generalize across different landscapes. Testing the model’s spatial transferability will be essential to determine its effectiveness in large-scale environmental monitoring and management beyond the study area.

The high R² and low RMSD values suggest that our model captured the majority of the variability in canopy cover across Florida, as highlighted earlier. This has significant implications for forest management, conservation planning, and ecological research. For example, accurate canopy cover maps can help in monitoring deforestation, assessing habitat quality, and evaluating the effectiveness of conservation interventions. These results underscore the importance of using advanced methodologies that integrate diverse data sources to improve predictive accuracy.

Optical and radar data provide complementary information about forest structure and canopy cover [100,101,102], which can enhance model robustness. Comparing our results with findings from other studies that relied solely on single data sources [103,104], those studies reported lower accuracy and higher error margins, emphasizing the potential for integrating multiple data sources, such as GEDI, Sentinel-1, and HLS, to improve model performance. However, it is important to recognize that these studies focused on different forest types and target metrics, which may account for some of the observed differences in performance.

Our model identified several key predictors of canopy cover, including longitude, red mean, red, latitude, and red min. The importance of these predictors aligns with findings from other studies that have emphasized the relevance of geographic and spectral variables in canopy cover estimation [105]. Specifically, longitude explained more than 95% of model variation, being considered important in 83% of the models. This was higher than the importance scores for spectral indices alone, in agreement with previous reports [105].

The relationship between latitude and vegetation distribution has long been recognized, dating back to Von Humboldt’s observations between 1799 and 1804, and it is strongly related to climatic zones [106,107]. The spatial and temporal coordinates of data points are often neglected in ecological studies, but it has been argued [107] that these should be included as predictors when modeling the spatial distribution of organisms. Similarly, it was found [108] that incorporating coordinate metrics such as longitude and latitude enhanced model performance in predicting forest structure.

Longitude and latitude are associated with the spatial distribution of species [106,107] and are closely linked to climatic zones and soil types [109,110]. These geographic metrics help capture environmental gradients that influence canopy cover. For instance, in our study, the trend of canopy cover distribution from south to north revealed that lower rates of canopy cover were prevalent in southern counties, while more canopy cover was observed in northern counties. This spatial pattern probably explains why latitude was one of the most important metrics, accounting for more than 53% of the model’s variation.

While the longitudinal patterns are not immediately obvious when looking at Figure 6a, there are reasons why this metric appeared in more than 85 out of 100 models. Florida is a peninsula, meaning that its vegetation distribution is highly impacted by the distance from the ocean. Coastal areas are subjected to harsher conditions such as salt spray, higher winds, and frequent disturbances like hurricanes, which can limit tree growth and reduce overall canopy density [111,112]. Additionally, coastal regions are often more developed due to urbanization and tourism [113], further decreasing canopy cover. Conversely, inland areas generally experience more stable conditions that are conducive to denser and more extensive forest cover. Thus, longitude is likely to serve as a proxy for these varying environmental and anthropogenic influences across the state.

Unexpectedly, vegetation indices (e.g., NDVI, NDWI) were not the most important metrics; instead, coordinate metrics and spectral bands were. The red, red mean, and red min metrics emerged as significant predictors of canopy cover. Spectral indices, particularly those in the red band, are critical in remote sensing for assessing vegetation health and density. In particular, the red band is sensitive to chlorophyll content and can indicate the presence and condition of vegetation [114]. Studies have shown that red band reflectance is inversely related to vegetation health; healthier vegetation absorbs more red light due to a higher chlorophyll concentration [99,114]. This relationship makes the red band a valuable indicator for modeling canopy cover. The red mean metric captures the average reflectance in the red band, providing a measure of overall vegetation health across an area. Meanwhile, red min indicates the minimum reflectance value, highlighting areas with the densest and healthiest vegetation, which is crucial for accurate estimation of canopy cover. Consequently, these spectral bands might also directly reflect canopy characteristics more effectively than composite indices, leading to their higher importance in the model.

Our findings are consistent with those of other studies that have highlighted the significance of geographic and spectral variables in canopy cover modeling. For example, studies [108,109] have emphasized the role of geographic metrics such as longitude and latitude in predicting vegetation patterns. Jongman et al. [106] similarly reported that spectral bands, particularly in the red region, were essential for accurate estimation of canopy cover.

The variability in canopy cover estimation performance across validation sites can be linked to site-specific factors, including forest types and structural characteristics. For example, GSWP is dominated by wetland forests with high canopy cover and a relatively homogenous structure [115], while Babcock Ranch Preserve features a mix of pine flatwoods and cypress swamps, representing moderate canopy cover [116]. In contrast, ANF is ecologically diverse, with landscapes ranging from longleaf pine savannas to hardwood swamps, making it the most structurally heterogeneous among the sites [117,118].

Our model performed best at GSWP, where the homogenous vegetation and higher rate of canopy cover probably contributed to lower errors and higher R² values. Conversely, ANF showed the poorest performance, with lower R² values and higher errors. This result aligns with the ecological complexity of ANF, where diverse vegetation types and canopy structures pose greater challenges for accurate predictions. Additionally, ANF has been heavily impacted by disturbances such as hurricanes (e.g., Hurricane Michael in 2018), which caused significant canopy loss and altered forest composition, further complicating the model’s performance [118].

These insights not only confirm the validity of our approach but also emphasize the importance of site-specific considerations in canopy cover modeling, especially in regions with diverse forest types. The incorporation of external ALS data for validation further underscores the innovative nature of our study, as validating a statewide canopy cover map across such ecologically diverse areas is a novel and critical contribution to the field. This approach enhances confidence in the model’s generalizability while providing a roadmap for similar studies in other regions with complex forest landscapes.

4.2. Spatial Characterization of Canopy Cover

The spatial distribution of canopy cover across the state shows a clear trend from south to north. Southern counties, such as Miami-Dade and Monroe, generally have lower mean canopy cover percentages (11–17%), while northern counties exhibit higher canopy cover (35–43%). This gradient can be attributed to climatic and soil differences, with northern Florida having cooler, wetter conditions and more fertile soils conducive to dense forest growth. According to the Köppen climate classification [119], southern Florida’s climate is classified as Aw (equatorial with dry winter) and Am (equatorial with monsoon), while northern Florida’s climate is classified as Cfa (warm temperate, fully humid with hot summer), contributing to the differences in vegetation patterns.

In particular, Miami-Dade County is characterized by flat topography, sandy soils, a shallow groundwater table, and well-developed canal systems [120]. These factors contribute to the lower canopy cover observed in the region. The flat topography and sandy soils are less supportive of dense forest growth, while the shallow groundwater table and canal systems alter the natural hydrology, impacting vegetation patterns.

The variability of canopy cover within counties is particularly notable in Monroe, Jefferson, Miami-Dade, Madison, Palm Beach, and Martin, as indicated by high standard deviation values. This variability can be attributed to the heterogeneous landscapes in these counties. For instance, southern counties exhibit contrasting environments such as beach areas with little vegetation alongside sparse forests. In northern counties like Jefferson and Madison, high cloud cover in HLS images may have introduced noise, leading to greater variability in canopy cover estimates.

Land use practices and history of disturbance significantly influence canopy cover distribution [121,122,123]. Northern counties tend to have more stable forest ecosystems, partly due to experiencing fewer hurricanes compared with the southern regions, which are more frequently affected by these [124], impacting forest cover dynamics. Additionally, conservation practices vary across the state [125], with some regions implementing more effective land management strategies that promote canopy cover.

Urban development patterns also play a crucial role, as more densely populated areas typically exhibit lower rates of canopy cover due to land conversion for residential and commercial use [122]. In Florida, population distribution varies significantly across the state. According to data from the U.S. Census Bureau [126], southern counties such as Miami-Dade, Broward, and Palm Beach have significantly higher population densities compared with the more rural northern counties. Central Florida, including Orange County (home to Orlando) and surrounding counties like Seminole, Osceola, and Lake, has moderate population densities but is experiencing rapid growth due to urbanization.

Moreover, our model also indicates that wetland forests generally exhibit more canopy cover compared with upland forests. This higher rate of canopy cover in wetlands can be attributed to their unique ecological conditions, including abundant water supply and fertile soils [127], supporting dense vegetation growth. These areas are also protected by law and are often less disturbed by human activities [128], allowing for more mature forest development. On the other hand, upland forests may experience more significant water stress, particularly during dry seasons [129], which is a primary reason for their having less canopy cover. Additionally, their well-drained and less nutrient-dense soils further contribute to reduced vegetation density [130]. Upland forest types, such as longleaf pine savannas, are also adapted to regular fire regimes, which can periodically reduce canopy cover [110].

The spatial characterization of canopy cover across Florida underscores the value of having detailed, county-level canopy cover information. This granularity enables a more nuanced understanding of regional forest dynamics and the factors that drive canopy cover variability, such as climate, soil type, land use, and disturbance history. Accurate county-level data enables targeted conservation and management efforts, as stakeholders can tailor strategies to the specific ecological and socio-economic conditions of each region [131]. For example, identifying areas with less canopy cover in densely populated southern counties could inform urban greening initiatives, while understanding the higher canopy in wetland forests highlights areas of high conservation value. By providing insights into both broad trends and localized variations, county-level canopy cover information becomes a crucial tool for effective forest management, biodiversity preservation, and urban planning throughout the state. This approach aligns with the needs of forest managers, conservationists, and policymakers who require precise data to make informed decisions at both local and statewide levels.

4.3. Limitations and Future Work

The comparison of our model’s performance and the identified key predictors with the results from existing studies highlights the effectiveness of integrating GEDI, Sentinel-1 GRD, and HLS data. However, despite these strengths, our model tended to underestimate canopy cover, as indicated by the negative MD. In our case, the MD of −36.15% suggests that the model tended to predict lower rates of canopy cover than were observed in the ALS dataset. This underestimation is a common issue in remote sensing models, particularly those using RF [132].

This underestimation may have stemmed from a combination of factors. First, optical and radar sensors have certain limitations, one of which is saturation [97,133,134]. Specifically, spectral saturation occurs when spectral reflectance values are no longer sensitive to increases in biomass or leaf area index, resulting in low estimation accuracy [135]. This phenomenon has been demonstrated in several studies [136,137,138,139,140] and is reflected in our results. Despite considerable research efforts to address saturation in high-density vegetation, it remains a significant challenge in the remote sensing community [64,139,141]. However, studies have shown that the fusion of different types of sensors can significantly mitigate the saturation problem [139,142,143,144]. In our study, the integration of LiDAR, SAR, and optical data substantially improved the accuracy of our canopy cover estimates, suggesting that without this data fusion, the issue of saturation would have been much more pronounced. LiDAR, in particular, stands out because unlike optical sensors, it is not prone to saturation, enabling it to capture fine-scale details across a wide range of canopy densities [139]. This characteristic allowed LiDAR to effectively penetrate dense canopies, providing valuable structural information that enhanced the precision of our canopy cover assessments.

In addition to the challenges of underestimation and sensor limitations, several other potential sources of error may have influenced our canopy cover predictions. One important factor is geographical variability. Environmental differences between northern and southern Florida can affect model accuracy, particularly in areas with heterogeneous landscapes, such as southern counties, where abrupt, small-scale changes in canopy cover present added challenges. Another significant source of error is geolocation uncertainty. The geolocation accuracy of GEDI footprints is estimated at around 7 m at the 1-sigma level [92]. While this uncertainty was not directly modeled in our analysis, we attempted to reduce its impact by using a spatial resolution of 30 m. However, differences in resolution between ALS and GEDI data, along with our own canopy cover maps, also contributed to discrepancies during validation. Despite the use of a 30 m spatial resolution, variations in data collection and canopy cover computation between ALS and GEDI introduced potential inaccuracies, especially since the ALS data were processed with minimal modeling, making such differences more noticeable when compared with our maps.

Seasonal variations and phenological changes can further exacerbate these errors, as collecting data within a narrow time frame does not entirely eliminate the discrepancies between predicted and actual canopy cover. For example, even small differences in vegetation growth stages between the ALS and GEDI datasets could have influenced canopy measurements, as foliage density, leaf area, and canopy structure fluctuate seasonally. These phenological changes, particularly in rapidly growing or deciduous species, can cause significant variations in canopy cover between the datasets [145,146], making it difficult to align the model’s predictions with ground-truth observations. Additionally, cloud cover in some HLS images, particularly in northern counties, may have reduced the quality of the data and introduced noise, leading to greater variability in the model’s predictions. These combined factors underscore the complexity of accurately estimating canopy cover in dynamic landscapes. Differences in seasonal conditions, phenology, and data quality due to cloud interference can introduce additional sources of error, limiting the model’s ability to predict canopy cover with full precision.

While our model effectively integrates multiple data sources and performs well in explaining variance in canopy cover, it is crucial to address the identified sources of error to enhance accuracy. Future work should focus on refining data integration methods, particularly to better capture high-density canopy areas, and on mitigating the impacts of geolocation uncertainty, resolution differences, seasonal variations, and cloud cover. Despite our model’s competitive performance metrics, the observed underestimation bias underscores the need for these refinements in order to achieve greater accuracy and robustness in canopy cover prediction. Furthermore, exploring advanced machine learning techniques such as deep learning for improving prediction accuracy has shown promise in recent studies [147]. Incorporating temporal data could also provide insights into seasonal variations [148] and long-term trends in canopy cover, further enhancing the model’s utility for dynamic forest monitoring. By addressing these aspects, future models can achieve higher precision and provide more reliable data for forest management and conservation strategies.

Overall, our study contributes significantly to the field of ecological modeling by providing a robust framework for canopy cover estimation and highlighting the critical role of data fusion in enhancing the accuracy and utility of environmental models. By integrating GEDI-derived canopy cover samples with optical imagery, radar, and ancillary data, our framework not only enhances canopy cover mapping at finer resolutions but also offers a realistic pathway toward better forest management and decision-making. By applying our model, we can more accurately calculate the impacts of yearly hurricane events on Florida’s forests, assessing changes in vegetation and carbon stock over time.

Another key strength of this study lies in the rigorous external validation using ALS data, ensuring the robustness of the model and increasing its reliability and applicability to unseen areas. Additionally, by incorporating spatial predictors such as geographic coordinates alongside spectral and structural data, we have demonstrated how combining diverse remote sensing variables can significantly improve vegetation mapping. This methodology also highlights the importance of accounting for environmental gradients, which are particularly pronounced in Florida due to its unique geographic and ecological characteristics.

Florida’s forests are characterized by a remarkable diversity of ecosystems ranging from upland to wetland forests, presenting unique challenges and opportunities for forest monitoring. Our study provides valuable insights into the spatial distribution of canopy cover across these forest types, offering the potential to identify areas that are more or less resilient to disturbances such as hurricanes. These findings have practical implications for forest management, including prioritizing reforestation efforts and mitigating the impacts of canopy loss. Furthermore, the ability to map canopy cover annually establishes a baseline for monitoring hurricane-induced damage over time, contributing to adaptive management strategies and long-term forest sustainability.

5. Conclusions

In this research, we investigated the efficacy of a data fusion framework using LiDAR, SAR, and optical data to create a comprehensive, wall-to-wall map of canopy cover for the state of Florida. By incorporating GEDI, HLS, and Sentinel-1 GRD data, we successfully addressed the spatial resolution limitations of GEDI alone, achieving a detailed and extensive canopy cover map for Florida, which had previously been unavailable. This integration proved to be highly effective, marking a significant advancement in canopy cover mapping. Our findings demonstrate the effectiveness of using geographic coordinate metrics and spectral bands as key predictors for capturing environmental gradients and vegetation health, which are essential for accurate estimation of canopy cover. Notably, we observed distinct spatial patterns, with higher canopy cover in northern Florida and lower canopy cover in the southern regions. This variability reflects the diverse ecological conditions across the state and offers valuable insights for targeted conservation and forest management efforts. Identifying and addressing potential sources of error, such as geolocation uncertainty and resolution differences between datasets, will be critical for improving future models and enhancing their predictive accuracy. In summary, our model approach holds significant potential for the conservation of Florida’s forests. The framework can be adapted to predict canopy cover for different years, including post-hurricane scenarios or in the wake of other environmental disturbances like droughts. This capability makes it an effective tool for assessing forest area loss and planning conservation strategies in response to environmental changes.