Unveiling spatial and temporal heterogeneity of a tropical forest canopy using high-resolution NIRv, FCVI, and NIRvrad from UAS observations

Recently, remotely sensed measurements of the near-infrared reflectance (NIRv) of vegetation, the fluorescence correction vegetation index (FCVI), and radiance (NIRvrad) of vegetation have emerged as indicators of vegetation structure and function with potential to enhance or improve upon commonly used indicators, such as the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI). The applicability of these remotely sensed indices to tropical forests, key ecosystems for global carbon cycling and biodiversity, has been limited. In particular, fine-scale spatial and temporal heterogeneity of structure and physiology may contribute to variation in these indices and the properties that are presumed to be tracked by them, such as gross primary productivity (GPP) and absorbed photosynthetically active radiation (APAR). In this study, fine-scale (approx. 15 cm) tropical forest heterogeneity represented by NIRv, FCVI, and NIRvrad and by lidar-derived height is investigated and compared to NIRv and EVI using unoccupied aerial system (UAS)-based hyperspectral and lidar sensors. By exploiting near-infrared signals, NIRv, FCVI, and NIRvrad captured the greatest spatiotemporal variability, followed by the enhanced vegetation index (EVI) and then the normalized difference vegetation index (NDVI). Wavelet analyses showed the dominant spatial scale of variability of all indicators was driven by tree clusters and larger-than-tree-crown size gaps rather than individual tree crowns. NIRv, FCVI, NIRvrad, and EVI captured variability at smaller spatial scales (∼ 50 m) than NDVI (∼ 90 m) and the lidar-based surface model (∼ 70 m). We show that spatial and temporal patterns of NIRv and FCVI were virtually identical for a dense green canopy, confirming predictions in earlier studies. Furthermore, we show that NIRvrad, which does not require separate irradiance measurements, correlated more strongly with GPP and PAR than did other indicators. NIRv, FCVI, and NIRvrad, which are related to canopy structure and the radiation regime of vegetation canopies, are promising tools to improve understanding of tropical forest canopy structure and function.

1 Introduction

Important spatial and temporal heterogeneity in structurally complex and species-rich tropical forests is not well characterized. Many factors contributing to this heterogeneity, including varying microclimate, light conditions, topography, crown structure, and patterns of tree mortality and regeneration, can produce high variability in carbon fluxes, ultimately affecting coarse-scale gross primary production (GPP) measurements in forests (e.g., Xu et al., 2015; Guan et al., 2015; Morton et al., 2014; Bohlman and Pacala, 2012; Laurance et al., 2012; Clark et al., 2008; Huete et al., 2008). Improving knowledge of tropical forest dynamics at multiple scales is crucial to monitoring and predicting resilience of tropical ecosystems and productivity under climate change (Liu et al., 2021; Clark et al., 2017; Laurance et al., 2012; Malhi, 2012; Wright, 2010; Saatchi et al., 2010; Lewis et al., 2009). Remote sensing (RS) measurements have been employed to uncover vegetation patterns of structure and productivity from local to global scales, often with a focus on filling gaps in knowledge regarding variation and uncertainties in GPP estimates (e.g., Jung et al., 2011; Glenn et al., 2008; Huete et al., 2002; Ryu et al., 2018; Yang et al., 2017; Jiang et al., 2008; Zhao et al., 2010; Heinsch et al., 2006; Running et al., 2004; Turner et al., 2003). Yet, the spatial mismatch between satellite data (e.g., 30 m to 1 km pixel resolution), which provide observations across large extents at repeat intervals, and site-specific plot level data (e.g., 0.1–1 ha) is in part responsible for the uncertainties in GPP estimates (Gelybó et al., 2013; Zhang et al., 2020). A way to solve this problem is to acquire high-spatial-resolution and high-temporal-resolution data that can capture fine-grained heterogeneity of tropical forests (Clark et al., 2017; Mitchard, 2018; Saatchi et al., 2011; Lewis et al., 2009). Unoccupied aerial systems (UASs) with hyperspectral imaging sensors offer an opportunity to collect tropical forest canopy data at high spatial resolution and which could address unknowns related to the high heterogeneity of tropical forests.

Traditional reflectance-based indices (RIs) using RS data, such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), are known to capture structural changes that are coincident with changes in GPP. RIs have provided optical methods using RS to track GPP via the light use efficiency (LUE) model (Monteith, 1977; Yuan et al., 2014; Medlyn, 1998). In the most commonly used formulation of the LUE model for RS, GPP is

$$\begin{array}{}\text{(1)}& \mathrm{GPP}=\mathrm{APAR}x\mathit{\epsilon },\end{array}$$

where APAR is the absorbed photosynthetically active radiation and ε is the efficiency with which the target vegetation converts the radiation to carbon (Gamon, 2015; Yuan et al., 2014; Running et al., 2004). APAR is derived from

$$\begin{array}{}\text{(2)}& \mathrm{APAR}=\mathrm{PAR}xf\mathrm{PAR},\end{array}$$

where PAR is the incoming photosynthetically active radiation and fPAR is the fraction of absorbed PAR. RIs commonly used in the LUE model of GPP as well as direct proxies for GPP are NDVI and EVI, because of a strong relationship to fPAR (Springer et al., 2017; Morton et al., 2016; Gamon et al., 2015; Porcar-Castell et al., 2014; Glenn et al., 2008; Gao et al., 2007; Huete et al., 2002; Zarco-Tejada et al., 2013). NDVI and EVI are typically used as proxies on seasonal timescales. When used to examine changes on shorter timescales, they have been multiplied by photosynthetically active radiation (PAR) to account for changes in radiation (incoming, absorbed, and scattered) which better align with GPP changes (Springer et al., 2017; Yuan et al., 2014). However, RIs alone have often not shown enough sensitivity to capture more fine-scale or rapid changes in vegetation, such as those in tropical forests, and questions linger about the ability to track green-up with RIs in evergreen regions (Liu et al., 2021; Yang et al., 2018a; Lee et al., 2013; Xu et al., 2015; Morton et al., 2014; Samanta et al., 2010; Sims et al., 2008).

Recently, three emerging vegetation indicators have been shown to track with GPP more closely than traditional RIs. These indicators are the near-infrared reflectance of vegetation (NIRv) (Badgley et al., 2017), the fluorescence correction vegetation index (FCVI) (Yang et al., 2020), and the near-infrared radiance of vegetation (NIRvrad) (Wu et al., 2020). Because they exploit additional information from the NIR region of the spectrum, NIRv, FCVI, and NIRvrad do not saturate in dense canopies or suffer the same level of contamination from senesced vegetation and soils as traditional RIs (Baldocchi et al., 2020; Badgley et al., 2017). Additionally, these indicators require only moderate spectral resolution data and are similarly straightforward to measure and calculate as RIs, making them accessible in a broad range of studies. Therefore, NIRv, FCVI, and NIRvrad could be employed as valuable indicators of canopy structure and function (Badgley et al., 2017, 2019; Dechant et al., 2020).

NIRv is the product of NDVI and the total near-infrared scene reflectance (NIR). NIRv from moderate spectral resolution satellite imagery and field spectrometers has been shown to empirically track both measured and modeled GPP globally, although with highest uncertainties in the tropics. The NIRv ∼ GPP relationship holds at monthly to seasonal timescales presumably due to co-incident changes in canopy phenology, light capture and scattering, and GPP (Badgley et al., 2017, 2019; Dechant et al., 2020). FCVI, derived from radiative transfer theory rather than an empirical relationship, is calculated from RS data by subtracting the reflectance in the NIR from the reflectance in the visible range (Yang et al., 2020). Yang et al. (2020) demonstrated that FCVI tracked GPP and solar-induced fluorescence (SIF; a radiance-based indicator of GPP) by capturing structure and radiation information from a vegetated canopy in field experiments with crops and in numerical experiments. Yet FCVI showed differences from NIRv due to exposed soil within the vegetated study areas. In previous studies, FCVI and NIRv were similar for dense green canopies where soils have less of an impact, but this has not yet been tested in the tropics (Wang et al., 2020; Badgley et al., 2019; Dechant et al., 2020). The product of NDVI and the NIR radiance, called NIRvrad, was proposed as a proxy for GPP on half-hourly and daily timescales. In contrast, NIRv and FCVI track changes on longer timescales (Wu et al., 2020; Dechant et al., 2020; Baldocchi et al., 2020; Zeng et al., 2019). Because the radiance of NIR accounts for incoming radiation at short timescales, NIRvrad has tracked GPP and SIF on half-hourly and diurnal scales as well as seasonally in crops and, to a limited extent, natural grass and savanna ecosystems (Dechant et al., 2020; Baldocchi et al., 2020; Zeng et al., 2019; Wu et al., 2020).

Readily available UAS-based hyperspectral sensors are capable of robust measurements of NIRv, FCVI, and NIRvrad at ultra-high spatial scales, i.e., tens of centimeters or less. In this regard, UAS-based data have the potential to improve our understanding of tropical forest structure and function over a range of scales that are poorly resolved by other RS platforms. Here, we use high-spatial-resolution UAS measurements to characterize spatial and temporal variation in a semi-deciduous tropical forest canopy during the dry season and compare commonly used spectral indices (NDVI and EVI) to newer vegetation indicators (NIRv, NIRvrad, and FCVI) by (i) examining correlations between GPP and vegetation indicators using mean values across the canopy throughout the day, (ii) evaluating the distribution of fine-spatial-resolution values (∼ 15 cm) across the canopy and examining changes in this spatial variation throughout the course of 2 d, and finally (iii) identifying the dominant spatial scale driving variation across our 10 ha study region.

2 Materials and methods

2.1 Study area

Barro Colorado Island (BCI), Panama, is a 1560 ha island (approximately 15 km²) in Gatun Lake, which was formed by the construction of the Panama Canal. The Smithsonian Tropical Research Institute manages the preserved area specifically for research. This semi-deciduous moist tropical forest receives approximately 2640 mm mean annual precipitation and has a mean temperature of 26 ^∘C with a dry season from approximately January through April (Detto et al., 2018). There is high species diversity, with approximately 500 tree species, approximately 60 species per hectare, and about 6.3 % of trees at > 30 cm diameter at breast height (dbh) (Bohlman and O'Brien, 2006; Condit et al., 2000). The UAS and ground measurements were focused on an area approximately 10 ha within the footprint of an eddy covariance tower near the center of the island (9.156440^∘, −79.848210^∘).

2.2 Data collection

The GatorEye Unmanned Flying Laboratory is a hardware and software system built for sensor fusion applications, and which includes hyperspectral, thermal, and visual cameras and a lidar sensor, coupled with a differential GNSS, internal hard drives, computing systems, and an inertial motion unit (IMU). Hardware and processing details, as well as data downloads, are available at http://www.gatoreye.org. The GatorEye Unmanned Flying Laboratory flew 13 missions on 30 and 31 January 2019 over the forest canopy within the eddy covariance tower footprint at an average height of 120 m above ground level (AGL) and at 12 m s⁻¹ (Fig. 1). In this study, we used radiometrically calibrated flight transects from the Nano VNIR 270 spectral band hyperspectral sensor (Headwall Photonics, Fitchburg, MA, USA), which covered approximately 1 ha per flight within the eddy covariance (EC) footprint in this study. The Nano sensor spectrally samples at approximately 2.2 nm and 12 bit radiometric resolution from 400 to 1050 nm. The frame rate was set to 100 fps, with an integration time of 12 ms and provided a pixel resolution of approximately 15×15 cm. The Nano sensor was calibrated to radiance by the manufacturer before the field campaign, and pixel drift was removed by dark image collection, which was corrected for during the conversion from digital number to radiance. The hyperspectral transects were equally subset for each flight in ENVI + IDL (Harris Geospatial, Boulder, CO). Each flight resulted in 1920 transects of approximately 400 m length composing three blocks discretized in 2500 data points. Simultaneous lidar was collected using a VLP-32c ultra puck (Velodyne, San Jose, CA), which was processed to a 0.5×0.5 m resolution digital surface model (DSM).

Turbulent fluxes and meteorological variables were measured from a 40 m eddy covariance (EC) flux tower (Fig. 1). The eddy covariance system includes a sonic anemometer (CSAT3, Campbell Scientific, Logan, UT) and an open-path infrared ${\mathrm{CO}}_{\mathrm{2}}/{\mathrm{H}}_{\mathrm{2}}\mathrm{O}$ gas analyzer (LI7500, LI-COR, Lincoln, NE). High-frequency (10 Hz) measurements were acquired by a data logger (CR1000, Campbell Scientific) and stored on a local PC. Other measurements made at the tower include air temperature and relative humidity (HC2S3, Rotronic, Hauppauge New York) and photosynthetically active radiation (PAR; BF5, Delta-T Devices, UK). EC data were processed with a custom program using a standard routine described in Detto et al. (2010). GPP was derived from daytime values of net ecosystem exchange (NEE) by adding the corresponding mean daily ecosystem respiration obtained as the intercept of the light response curve (Lasslop et al., 2010). Due to a power issue, EC data were not available during the 30 January flights; so only 31 January GPP data were available.

An HandHeld 2 Pro spectroradiometer (HH2; Malvern Panalytical, Boulder, CO) fitted with a diffuse cosine receptor was used on the ground in full sun at the forest edge to record incoming irradiance on 30 and 31 January 2019 (∼ 3 nm FWHM and spectral sampling at 1 nm). HH2 irradiance was resampled to match the Nano hyperspectral sensor and used to calculate reflectance. A calibrated reference tarp was placed in full sun at the forest edge, and the UAS flew over and recorded the tarp each UAS flight. Reflectance was calculated separately using the HH2 and tarp data and resulting reflectance values compared as a method to vicariously cross-calibrate reflectance from the hyperspectral data (<7.0 % difference for all data in the study). In addition, PAR was calculated with the HH2 data and compared to the tower-mounted PAR measurement (approximately 1.5 km apart) to help understand any differences in the sky conditions during flight times. PAR differences across the site for each flight time for the duration of flights (approximately 10–15 min in length each) ranged between 4.0 % and 10.3 %.

2.3 Vegetation indicators

We calculated NDVI and EVI as Tucker (1979), Huete et al. (2002), and Rouse et al. (1974):

$$\begin{array}{}\text{(3)}& \mathrm{NDVI}={\displaystyle \frac{R_{\mathrm{770}-\mathrm{800}}-R_{\mathrm{630}-\mathrm{670}}}{R_{\mathrm{770}-\mathrm{800}}+R_{\mathrm{630}-\mathrm{670}}}}\end{array}$$

and

$$\begin{array}{}\text{(4)}& \mathrm{EVI}={\displaystyle \frac{\mathrm{2.5}\left(R_{\mathrm{770}-\mathrm{800}}-R_{\mathrm{630}-\mathrm{670}}\right)}{R_{\mathrm{770}-\mathrm{800}}+\mathrm{6}\times R_{\mathrm{630}-\mathrm{670}}-\mathrm{6}\times R_{\mathrm{460}-\mathrm{475}}+\mathrm{1}}},\end{array}$$

where R is reflectance and the subscripts indicate wavelengths. Here, we used the averages of 770–800 nm for NIR, 630–670 nm for red reflectance, and 460–475 nm for blue band reflectance and normalized to reduce noise.

We further calculated the near-infrared vegetation index NIRv as

$$\begin{array}{}\text{(5)}& \mathrm{NIRv}=\mathrm{NDVI}\times R_{\mathrm{770}-\mathrm{800}},\end{array}$$

where R_770−800 is the NIR reflectance (Badgley et al., 2017). The fluorescence correction vegetation index (FCVI) was calculated from spectral data by subtracting the reflectance in the visible range (R_400−700) from the NIR reflectance (Yang et al., 2020) as follows

$$\begin{array}{}\text{(6)}& \mathrm{FCVI}=R_{\mathrm{770}-\mathrm{800}}-R_{\mathrm{400}-\mathrm{700}}.\end{array}$$

The near-infrared radiance of vegetation (NIRvrad) was calculated similarly to the NIRv, except NDVI was multiplied by the radiance, rather than reflectance, from the NIR region (R_770−800) (Wu et al., 2020) as follows:

$$\begin{array}{}\text{(7)}& \mathrm{NIRvrad}=\mathrm{NDVI}\times R_{\mathrm{770}-\mathrm{800}}.\end{array}$$

2.4 Data analysis

A workflow summarizing data analyses is provided in Fig. 1. We examined mean values across the canopy over the course of 1 d by creating a diurnal time series of scatterplots of the tower-based PAR data, tower-based GPP data, and means of all spectral vegetation indicators, on 31 January 2019, and ran comparisons using Pearson's correlation coefficients to examine correlations. Results are provided in Sect. 3.1 and Fig. 2. At fine spatial scales, i.e., pixel sizes of ∼ 15 cm, we created density plots, calculated the coefficient of variation (CV), and calculated the means of all vegetation indicators (NDVI, EVI, NIRv, FCVI, NIRvrad) for each flight to compare spatial and temporal variability. Results are provided in Sect. 3.2 and Fig. 3. To determine which spatial scales dominate the variability of each vegetation quantity, we ran power spectrum wavelet analysis using code created in the MATLAB programming language (MathWorks, Natick, Massachusetts). For each vegetation quantity and each flight, and for the lidar elevation model representing canopy height, we computed the Morlet wavelet power spectrum of individual transects (Torrence and Compo, 1998). All power spectra from the wavelet analysis were normalized to unit variance. An ensemble power spectrum for each vegetation indicator was created by averaging across all the transects of each flight and then across flights. We then compared the power spectra for each vegetation indicator and lidar data to compare the spatial scales at which the quantities captured variability as well as the spatial scale at which the lidar-based elevation model captured variability. Results are provided in Sect. 3.3 and Fig. 4. For illustration purposes, Fig. A1 is an example of two synthetic signals generated with the fractal Brownian motion algorithm and different level of noise-to-signal ratio (Signal A and B, respectively, Fig. A1) and the corresponding power spectra which decay differently at smaller spatial scales (Power Spectra, Fig. A1). Initial UAS data processing was carried out in Interactive Data Language (IDL) and Environment for Visualizing Images (ENVI) (Harris Geospatial, Boulder, CO). Other analyses, including graphical illustrations, were carried out using the R open-source environment with libraries dplyr, ggplot, and tidyverse (R Development Core Team, 2010; Wickham et al., 2018; Wickham, 2017, 2016) and MATLAB R2019a (MathWorks, Natick, Massachusetts).

Figure 1 Summary of methods. Diagram representing discrete flight times for UAS and near-continuous EC-estimated GPP (far left). Platforms and instrumentation (blue) consisted of the Malvern Panalytical HandHeld 2 Pro spectroradiometer (HH2), the GatorEye Unmanned Flying Laboratory, and the tower at Barro Colorado Island (BCI). Data collected included irradiance, hyperspectral, lidar, eddy covariance system (EC), and photosynthetically active radiation (PAR). Calculations made were PAR with the HH2 (PARHH2), normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), fluorescence correction vegetation index (FCVI), the near-infrared vegetation index (NIRv), the near-infrared radiance of vegetation (NIRvrad), the digital surface model (DSM), gross primary productivity (GPP), and PAR from the PAR sensor on the tower (PARtower). An overview of the data analysis at each scale is provided in the right of the diagram.

3 Results and discussion

3.1 Diurnal trend in spectral vegetation indicators, PAR, and GPP

The degree to which remote sensing vegetation indicators represent changes in GPP depends largely on canopy-structure-dependent light absorption and scattering processes, that is, exploiting relationships between a remote sensing vegetation quantity, PAR or APAR, and GPP. Figure 2 shows GPP, PAR, and the mean value of each vegetation quantity at each flight time over the course of 31 January, the day on which we had overlapping data between the UAS and eddy covariance system (Fig. 2a–d). Additionally, Pearson correlation coefficients among mean NIRv, FCVI, NIRvrad, EVI, and NDVI for each flight time and the GPP and PAR values at the flight times are shown in Fig. 2d. NIRv is significantly and strongly positively correlated to both FCVI (r=0.9, p<0.001) and EVI (r=0.9, p<0.01). NIRvrad is the only vegetation quantity with a significant correlation to PAR and GPP, with a strong positive relationship (0.9 and 0.81, respectively, p values <0.05; Fig. 2d). Mean NIRvrad values also have the greatest relative diurnal change among the vegetation indicators (Fig. 2c and d). These results demonstrate that a shared correlation of NIRvrad and GPP to PAR results in mean NIRvrad tracking diurnal changes in GPP to a greater degree than NIRv, FCVI, NDVI, or EVI, because NIRvrad takes incoming radiation into account, whereas the other vegetation indicators do not. The ability of NIRvrad to track APAR is notable alone. However, our evidence supports the proposed use of NIRvrad as a proxy for changes in GPP on short timescales – albeit based on only 1 d of data. NIRvrad is a more practical proxy of GPP than SIF in the sense that a separate instrument to measure PAR is not needed (Wu et al., 2020; Zeng et al., 2019). Given that the relationship between NIRvrad and GPP depends on PAR, it is unclear if the association between NIRvrad and GPP would weaken during the wet season when low light or diffuse light conditions are more common (Berry and Goldsmith, 2020).

Figure 2 Diurnal time series smoothed with a LOESS filter of (a) GPP; (b) PAR; (c) NIRvrad; (d) NIRv, FCVI, NDVI, and EVI; and (e) comparisons of quantities using Pearson correlations color indicates strength of relationship: ^*: p value <0.05, ${}^{**}$: p value <0.01, and ${}^{***}$: p value <0.001.

3.2 Tropical forest canopy variation

Spatial distributions and the coefficient of variation (CV) of all pixels of NIRv, FCVI, and NIRvrad are generally similar to one another and show considerable variation spatially across the canopy and temporally over the course of a day and across days (Fig. 3a–c, Table A2). NIRv, FCVI, and NIRvrad distributions are distinct from EVI and NDVI (Fig. 3a–e, Tables A1 and A2). NIRv, FCVI, and NIRvrad have the highest CV at each flight time (between 39.78 % and 91.54 %, Table A1), followed by EVI (between 20.24 % and 37.24 %, Table A2), and NDVI varied the least at any flight time (between 9.83 % and 12.82 %, Table A2). For some indices, mean values across the canopy fail to capture extreme high (NIRv, NIRvrad, and FCVI) or low values (NDVI) during morning and afternoon hours. This pattern suggests “hot” and “cool” spots of activity related to heterogeneity in forest structure and low sun angles. In previous studies, the directional effects on NIRv have been examined on coarse spatial scales (i.e., satellites) and have been proposed as a means of improving understanding of NIRv agreement to GPP (Hao et al., 2021; Dechant et al., 2020; Baldocchi et al., 2020; Zhang et al., 2020). Our results demonstrate that NIRv, FCVI, and NIRvrad capture fine-grained heterogeneity of this tropical forest canopy, which was obscured by EVI and NDVI (Fig. 3a–e). NIRv and NIRvrad use NDVI; thus, by definition, NIR is the largest contributing factor to the heterogeneity captured (Fig. 3a, c, and e). While NIRv and NIRvrad distributions are generally similar, they diverge in the afternoons when PAR declines, which likely is why NIRvrad is better correlated with GPP. EVI variability was higher than NDVI variability but lower than that of NIRv, FCVI, and NIRvrad, indicating that EVI has a different level of sensitivity to viewing geometry and canopy components (potentially understory), light absorption, and scattering regime of the canopy than the other indices (Table A1and Table A2). We also show empirically that NIRv and FCVI are virtually the same in a dense tropical forest presumably due to both indices similarly representing the radiation regime of the tropical forest canopy, i.e., light capture and scattering, in conditions with little background soil, supporting the predictions of earlier studies (Dechant et al., 2020; Zeng et al., 2019; Yang et al., 2018b; Wu et al., 2020).

Midday distributions of NIRv, FCVI, and NIRvrad on 30 January at 12:00 and 13:30 and 31 January at 12:30 are less skewed than at other times of the day, whereas morning and afternoon distributions are skewed toward lower values, except for 31 January at 15:30 (Fig. 3a–c). On both days, when mean values peak at midday, the variation for all vegetation indicators is lowest (30 January at 12:00 CV between 47.6 and 49.2 and 31 January at 12:30 CV between 45.6 and 47.2) (Fig. 3, Table A1). The highest variability occurred in the afternoon on both days (30 January at 16:30 CV between 91.3 % and 91.5 and 31 January at 14:30 CV between 83.3 % and 83.8 % for all quantities) (Fig. 3, Table A2). At midday, NIRv, FCVI, and NIRvrad variability was low and means were high, indicating that viewing and sun geometry drive the higher and lower values during morning and afternoon. This effect is greater in the afternoon than the morning (Fig. 3, Table A2). However, a different pattern is apparent on 31 January during the 15:30 flight time when mean values increased from the 14:30 flight time means and the CV values were the lowest of any flight observations in the study, and this influence appears to be greatest on EVI. It is possible that this was due to another type of effect on illumination geometry, such as wind influencing the UAS, diffuse radiation effects, or hotspot effects.

Figure 3 NIRv (a), FCVI (b), and NIRvrad (c) density plots for each flight time on 30 and 31 January 2019. Colors of distributions indicate day.

3.3 Power spectrum analysis

Power spectrum analysis was used to identify the dominant spatial scales driving variability across the canopy (Fig. 4). In Fig. 4, the area beneath the curve is proportional to the variance because it is the spectrum divided by the corresponding scale and then plotted as a function of the log of the scale (example signals and power spectra provided Fig. A1). Similar to their spatial distributions (Fig. 3), NIRvrad and FCVI are indistinguishable in their dominant scales of spatial variability (Fig. 3) (Dechant et al., 2020; Zeng et al., 2019). Power spectrum analysis shows a distinct peak around 50 m spatial scale for NIRv, NIRvrad, FCVI, and EVI, whereas NDVI peaks at approximately 90 m. The largest tree crown sizes on BCI are on the order of 20–30 m in diameter, and the most common crown sizes are between 4–10 m (Fig. A2). Thus, the spatial variability of the vegetation indicators is strongly influenced by larger forest structures, such as forest gaps and tree clusters, rather than individual tree crowns.

This larger scale of variability is also confirmed by the power spectrum of the lidar-derived canopy surface model, which displays a peak at 70 m scale, indicating that larger than tree crown scales produce the most variability in canopy height. In other words, UAS-based lidar data also show that canopy heights within a 70 m spatial scale create strong spatial features on the landscape. Vegetation indicators and the lidar canopy surface model appear less effective at capturing smaller scale differences within a canopy (leaves or leaf clumps) or among the most frequent tree crown sizes on BCI (4–10 m sunlit tree crown sizes determined by stereophotos; Fig. A2). However, the peaks in the vegetation indicators are broader than the peak in the lidar data, showing that smaller features of the canopy are still contributing to the total spatial signal in the power spectra. These results suggest that satellite data with a spatial resolution greater than ∼ 50 m may miss important variation in diverse tropical forest canopies. NDVI displays a different shape with a slower decay at small scales, indicating less distinguishable spatial structures from the canopy, and a peak shifted to the larger scales (Fig. 4); i.e., NDVI does not distinguish smaller spatial structures. At much larger scales (>100–200 m), the vegetation indicators decline smoothly, while NDVI and especially lidar show an increase in variance probably associated with topographic heterogeneity.

One reason why vegetation indicators and lidar captured variability at spatial scales larger than the most common tree crown sizes on BCI is that canopy heights tend to be more uniform on BCI compared to other tropical forests, possibly due to wind (Bohlman and O'Brien, 2006). For example, Dipterocarpus-dominated Southeast Asian forests have emergent trees, unlike BCI, which can reach up to 60 m in height. Additionally, tree crowns on BCI tend to be more flat-topped than conical or rounded, and trees can be found clumped in similar heights, which could explain why the most often detected unit is larger than the mean of a single crown. On the other end of the spectrum, forest gaps can be larger than a single crown because treefall often affects neighboring trees.

Vegetation indicators and the lidar-derived surface model represent the spectral and structural properties most broadly of the upper canopy, and thus it is conceivable that they display similar spatial variability. However, NIRv, FCVI, NIRvrad, and EVI discriminated details at a different spatial scale from NDVI and lidar. These results parallel the variability detected in their distributions (Fig. 3 and Table A1), where NDVI patterns were distinct from the other vegetation indicators. Taken together, these results show that NIRv, FCVI, and NIRvrad have a smoother spatial pattern and peak at finer scales than NDVI, which is known to saturate at high green biomass (Zhu and Liu, 2015; Huete et al., 2002), whereas NIRv, FCVI, and NIRvrad should better correlate with aspects of photosynthetic capacity. Thus, these emerging indicators should measure finer-resolution spatial heterogeneity and should be more adept at monitoring changes in structure and function of the canopy than NDVI. Additionally, the emerging indicators can potentially disaggregate the physiological and structural component of SIF when SIF measurements are available since changes in structure of the forest coincide with changes in GPP (Wang et al., 2020; Wu et al., 2020; Yang et al., 2020; Dechant et al., 2020). Emerging indicators' heightened ability to differentiate the fine-scale spatial variability in the canopy is likely due to the influence of high upwelling of NIR from the canopy and understory, particularly in the dry season, which tends to blur the signal of the upper canopy for NDVI. Notably, EVI and NDVI, two common indicators of vegetation greenness, show differences in their power spectrum, in particular the slope of the curve for scales less than 20 m. EVI was designed to better capture vegetation changes by exploiting variability in the reflectance in the blue range, especially effective in dense green canopies. This may help explain the scale of variability in this canopy where variation in the blue may be expected to manifest, especially because deciduous crowns, which have high reflectance in blue wavelengths compared to fully leaved crowns, are present on BCI (Bohlman, 2008).

Figure 4 Ensemble wavelet power spectra for all the quantities used in this study and a lidar-derived digital surface model (DSM). Note that FCVI and NIRv are similar; thus, the NIRv curve is obscured by the FCVI. Ensembles were created by averaging the spectrum of individual transects and then averaging across flights. Note that in this representation, the spectrum divided by the corresponding scale as a function of the log of the scale, the area beneath the curve is proportional to the variance.

4 Conclusions

We examined NIRv, FCVI, and NIRvrad, emerging vegetation indicators related to fPAR of a semi-deciduous tropical forest canopy using UAS-based hyperspectral data. Our findings demonstrate that NIRvrad has greater potential to track GPP over the course of a day than the non-radiance-based indices as evidenced by a shared correlation among NIRvrad, PAR, and GPP. Thus, NIRvrad is a potential proxy for tracking GPP on short timescales without the need for separate measurements of incoming irradiance. Also, NIRv, FCVI, and NIRvrad at high spatial resolution (∼ 15 cm) unveil greater spatial and diurnal variability of BCI's tropical forest canopy versus EVI or NDVI, which may pave the way to improve our understanding of the relationship between GPP and remote sensing observations. For instance, by benchmarking changes of vegetation function and structure that underlie a GPP measurement representing the whole EC footprint, fine-scale NIRv, FCVI, or NIRvrad measurements may reveal highly differential behaviors of tropical species diurnally to seasonally. The dominant scale driving spatial variability of spectral measurements and lidar data is larger forest structures occurring on BCI, such as groups of similar trees or forest gaps. Yet, smaller, broader peaks in the power spectra of NIRv, FCVI, NIRvrad, and EVI indicate these four indices incorporate smaller scale information compared to NDVI. Taken together, the demonstrated potential to track GPP, measure spatial heterogeneity and variability, and capture forest structural characteristics of BCI open greater possibilities to examine structure and function within and across this tropical forest.

Because remote sensing advancements are making it possible to capture physiological responses of vegetation, the importance of improved techniques to examine the radiation regime, for instance estimating fPAR or APAR, can be overlooked. However, recent studies have highlighted the importance and difficulties of measuring fPAR and APAR, the strong dependence of measurements on illumination and viewing geometry, and the need for increased understanding of structure-related radiation regime information more generally (e.g., Hao et al., 2021; Dechant et al., 2020; Baldocchi et al., 2020; Rocha et al., 2021; Zhang et al., 2020). For NIRv, FCVI, and NIRvrad, inclusion of the NIR spectral region makes the emerging indices more sensitive to incoming, absorbed, and scattered radiation, which can be influenced by illumination and viewing geometry, changes in canopy leaf angles, or associated structure changes. In the case of NIRvrad, which was most strongly associated with GPP, changes in light regime and associated photosynthetic capacity can even be captured diurnally. Further use of NIRv, FCVI, and NIRvrad measurements, especially at high spatial and temporal resolution, can inform our understanding of information each captures from a canopy, as well as improving interpretation of traditional reflectance-based indices, and other remote sensing measurements, such as SIF. This study highlights the importance of understanding the incoming solar radiation, absorbed and scattered radiation, and illumination and viewing geometry of any remote sensing data, but it also encourages exploiting RS observations to improve our ability to measure structure-related light capture and scattering patterns. It is in this role that we show these measurements should be further investigated as valuable tools to improve our understanding of complex tropical forest canopies and potentially as an improved estimate of fPAR, APAR, or GPP. While this study focuses on BCI, these techniques could be applied more broadly for the purposes of defining the dominant scale of spatial variability, tracking structural changes, monitoring coincident changes in GPP or light regime, or as inputs to vegetation models of tropical forest structure and function.

Appendix A

Figure A1 Sample signals with relatively higher noise (signal A) and lower noise (signal B) and their corresponding power spectra ensemble plotted as normalized on log scale. Note the representation of the variance by area under the curve is preserved by multiplying the power (S(f)) by the frequency (f). In this way the area beneath the curve is still proportional to the variance.

Figure A2 Distribution of tree crown sizes on BCI in a sample ∼ 10 ha plot taken from digitized high-spatial-resolution stereophotos that were linked to stems in the field (Bohlman and Pacala 2012). This ∼ 10 ha plot does not coincide with the ∼ 10 ha area sampled by the UAS near the eddy covariance tower in this study.

Table A1 Mean, standard deviation (SD), and coefficient of variation (CV) of NIRv, NIRvrad, and FCVI measurements for the study.

Table A2 Mean, standard deviation (SD), and coefficient of variation (CV) of NDVI and EVI measurements for the study.