AnisoVeg: anisotropy and nadir-normalized MODIS multi-angle implementation atmospheric correction (MAIAC) datasets for satellite vegetation studies in South America

The AnisoVeg product consists of monthly 1 km composites of anisotropy (ANI) and nadir-normalized (NAD) surface reflectance layers obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor over the entire South American continent. The satellite data were preprocessed using the multi-angle implementation atmospheric correction (MAIAC). The AnisoVeg product spans 22 years of observations (2000 to 2021) and includes the reflectance of MODIS bands 1 to 8 and two vegetation indices (VIs), namely the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI). While the NAD layers reduce the data variability added by bidirectional effects on the reflectance and VI time series, the unique ANI layers allow the use of this multi-angular data variability as a source of information for vegetation studies. The AnisoVeg product has been generated using daily MODIS MAIAC data from both Terra and Aqua satellites, normalized for a fixed solar zenith angle (SZA = 45^∘), modeled for three sensor view directions (nadir, forward, and backward scattering), and aggregated to monthly composites. The anisotropy was calculated by the subtraction of modeled backward and forward scattering surface reflectance. The release of the ANI data for open usage is novel, and the NAD data are at an advanced processing level. We demonstrate the use of such data for vegetation studies using three types of forests in the eastern Amazon with distinct gradients of vegetation structure and aboveground biomass (AGB). The gradient of AGB was positively associated with ANI, while NAD values were related to different canopy structural characteristics. This was further illustrated by the strong and significant relationship between EVI_ANI and forest height observations from the Global Ecosystem Dynamics Investigation (GEDI) lidar sensor considering a simple linear model (R²=0.55). Overall, the time series of the AnisoVeg product (NAD and ANI) provide distinct information for various applications aiming at understanding vegetation structure, dynamics, and disturbance patterns. All data, processing codes, and results are made publicly available to enable research and the extension of AnisoVeg products for other regions outside of South America. The code can be found at https://doi.org/10.5281/zenodo.6561351 (Dalagnol and Wagner, 2022), EVI_ANI and EVI_NAD can be found as assets in the Google Earth Engine (GEE; described in the data availability section), and the full dataset is available from the open repository https://doi.org/10.5281/zenodo.3878879 (Dalagnol et al., 2022).

1 Introduction

Anisotropy is defined as the departure from Lambertian scattering (isotropic), which is caused by the physical structure of media through which photons pass. Because most land covers are not Lambertian (isotropic), the surface reflectance measured by satellite sensors varies with the view zenith angle (VZA), view direction (backward or forward scattering), and solar zenith angle (SZA; Galvão et al., 2011). This is especially valid for images acquired over vegetated surfaces by large field-of-view (FOV) instruments such as the Moderate Resolution Imaging Spectroradiometer (MODIS; Bhandari et al., 2011). MODIS has a wide swath scanning ±55^∘ from nadir on board the Terra and Aqua satellites. For example, a reflected signal coming from the backward-scattering direction of MODIS under a large VZA and close-to-zero relative azimuth angle (RAA) between the satellite and Sun (Sun behind the platform) is generally higher than that coming from the nadir (VZA = 0^∘) or forward-scattering direction (platform facing the Sun at RAA = 180^∘). Moreover, the SZA also varies seasonally and across geographical locations, affecting the number of shadows in the surfaces observed by satellites (Galvão et al., 2013). Such view illumination effects are dependent on the land cover types, and their magnitude relates to differences in biophysical properties of the vegetation (Galvão et al., 2004; Sims et al., 2011). Therefore, the vegetation anisotropy can be seen antagonistically as sources of noise and biophysical information in the time series analysis of vegetation indices (VIs) calculated from MODIS. As a source of noise, one may consider that the reflected signal toward the large FOV satellite sensors varies with distinct view illumination geometries of data acquisition over the same surface. As a source of information, one may highlight that the anisotropy is land-cover-type dependent, showing spectral variations that may be associated, for instance, with changes in vegetation structure across different forests.

To reduce the bidirectional effects as a source of noise, a nadir-normalized dataset can be created. We can normalize the surface reflectance of the MODIS bands to a specific set of VZA and SZA, using the bidirectional reflectance distribution function (BRDF) represented by a model such as the Ross–Thick Li–Sparse (RTLS; Wanner et al., 1995). To ensure confidence in the data analysis, we can also use the multi-angle implementation atmospheric correction (MAIAC) for atmospheric correction. MAIAC is a new generation of cloud screening and atmospheric correction algorithm that uses an adaptive time series analysis and processing of groups of pixels to derive atmospheric aerosol concentration, cloud mask, and surface reflectance without typical empirical assumptions (Lyapustin et al., 2011, 2012). By mitigating atmospheric interference and advancing the accuracy of surface reflectance over tropical vegetation by a factor of 3 to 10, MAIAC offers substantial improvement over conventional products such as the MOD09 (Hilker et al., 2012). Because of the better data quality retrieval, MAIAC is also an alternative to the MCD43A4 16 d Nadir Bidirectional Reflectance Distribution Function (BRDF)-Adjusted Reflectance (NBAR) product due to the less variable seasonal signal (3 to 10 times) over evergreen forests resultant from reduced effects of sun view geometry. While the MCD43A4 NBAR product offers view illumination correction, by using the MAIAC products, one can also correct for solar illumination effects at the same time. Due to the improvements in cloud detection, aerosol retrieval, and atmospheric correction, the MAIAC algorithm provides from 4 % to 25 % more high-quality retrievals than the traditional MOD09 product, with the largest estimate being observed for tropical regions (Lyapustin et al., 2021). Studies have used MODIS MAIAC observations with nadir-normalized geometry to assess the Amazon forest structure, functioning, and impacts of environmental and climate change (Hilker et al., 2014; Wagner et al., 2017; Anderson et al., 2018; Dalagnol et al., 2018; Fonseca et al., 2019; Bontempo et al., 2020; Gonçalves et al., 2020; Zhang et al., 2021). For instance, such products have provided a reliable time series of surface reflectance data that allowed the identification of large-scale communities of bamboo species and their dynamics in the southwestern Amazon (Dalagnol et al., 2018). Last, by improving the cloud screening and minimizing BRDF artifacts in comparison to uncorrected data, the MAIAC greatly contributed to the understanding of the long-standing debate in the Amazon over the possible existence of the green-up phenomenon observed during the dry season of each year or with severe droughts (Morton et al., 2014; Bi et al., 2015; Saleska et al., 2016; Wu et al., 2018). The existence of this phenomenon has implications on the comprehension of the resilience of tropical forests to climate change.

To use the bidirectional effects as a source of information, we generate an anisotropy dataset that is dependent on land cover types and captures the variations in sunlit and shaded canopy components viewed by the sensors (Chen et al., 2003; Gao et al., 2003). The use of multi-angular information to obtain the metrics of anisotropy and extract information on forest structure was suggested 2 decades ago (Gobron et al., 2002; Diner et al., 2005). One of the early experiments exploring the use of anisotropy to extract information about vegetation structure was conducted by calculating the ratio between backward- and forward-scattering data and generating the anisotropy index (ANIX) for studying short-stature grass-type vegetation (Sandmeier et al., 1998). Other indices have been developed and validated afterwards (Schaaf et al., 2002; Lacaze et al., 2002; Chen et al., 2005; Pocewicz et al., 2007; de Moura et al., 2015; Sharma, 2021). However, this remains an understudied topic, with limited results reported in the literature, especially in tropical regions. For instance, observations from the Multi-angle Imaging SpectroRadiometer (MISR)/Terra in the backward- and forward-scattering directions facilitated the discrimination of savanna physiognomies in Brazil (Liesenberg et al., 2007). MODIS MAIAC data from both directions were also used to calculate an anisotropic VI that explained part of the large-scale photosynthetic activity in the Amazon, where higher photosynthetic activity was associated with higher anisotropy values (de Sousa et al., 2017). De Moura et al. (2015) employed a more sophisticated approach based on scattering at backward and forward view directions using multi-temporal and multi-angular observations of MAIAC MODIS and BRDF modeling. The resultant metrics of anisotropy were further validated against field and airborne light detection and ranging (lidar) observations, showing strong linear relationship with leaf area index (LAI; R²=0.70–0.88), canopy heterogeneity (R²=0.54), and photosynthetic activity (R²=0.73–0.98; de Moura et al., 2015, 2016; Hilker et al., 2017). Although showing great potential in vegetation studies, the aforementioned anisotropy metrics were never computed over larger areas of the world as is proposed in this study for South America.

The objective of this work is to present the AnisoVeg product, and how it can be used for vegetation studies. We use MODIS Collection 6 (C6) MAIAC (Lyapustin et al., 2018) monthly data (2000–2021) generated at 1 km spatial resolution for the entire South American continent, with two different types of layers, i.e., (1) nadir-normalized (NAD) data, for the surface reflectance of MODIS bands 1 to 8 and two VIs (normalized difference vegetation index, NDVI, and enhanced vegetation index, EVI), and (2) anisotropy data (ANI), calculated from the difference between backward- and forwarding-scattering estimates of bands 1 to 8 and VIs (de Moura et al., 2015). The motivations for generating this product extend from developing applications of multi-angle observations for vegetation studies to producing analysis-ready and openly available datasets of anisotropy and nadir metrics for a larger community of users. The paper is organized in several sections to present the processing steps for generating the AnisoVeg products, a brief evaluation of data products over experimental areas, and, finally, an example of its potential application in vegetation studies.

2 Methodology to compute the AnisoVeg product

2.1 Daily MODIS MAIAC surface reflectance data over South America

Daily surface reflectance data were obtained from the MODIS product MCD19A1 v006 (collection 6) for the tiles covering South America (Fig. 1). According to the MODIS traditional tiling system, these tiles ranged from 9–14 (horizontal) and 7–14 (vertical). The input data consisted of cross-calibrated surface reflectance from Terra and Aqua satellites on eight spectral bands (Table 1), with 1 km spatial resolution from 2000 to 2021 (Lyapustin and Wang, 2018; https://doi.org/10.5067/MODIS/MCD19A1.006). This product provides surface reflectance data corrected for atmospheric effects by the MAIAC algorithm and is controlled for cloud-free and clear-to-moderately-turbid conditions with aerosol optical depth (AOD) at 0.47 µm below 1.5 (Lyapustin et al., 2018). The MAIAC algorithm uses a time series approach for improved cloud filtering, amongst other filters such as surface reflectance change, in order to provide the most accurate surface reflectance estimates. The raw data were obtained from NASA's Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC; available at https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MCD19A1/, last access: 16 January 2023).

Figure 1 AnisoVeg product concept and the area of coverage. (a) Schematic representation showing the observational geometry and the processing steps for producing NAD and ANI data from MODIS and to provide information on vegetation heterogeneity and structure. (b) The visualization of the anisotropy EVI (EVI_ANI) for South America from August 2021 at 1 km spatial resolution, showing the coverage of the product in South America and the location of three sites used to demonstrate potential applications. The sites are (1) Tapajós National Forest, (2) São Félix do Xingu, and (3) Xingu Indigenous Park. Red lines indicate the country boundaries.

Table 1 MODIS spectral bands. NIR is for near-infrared, and SWIR is for shortwave infrared.

2.2 The AnisoVeg product

The AnisoVeg product consists of two main types of data, spanning from 2000 to 2021, in monthly composites at 1 km spatial resolution, with (a) the nadir-normalized (NAD) data and (b) the anisotropy (ANI) data. Each data type has 10 layers, corresponding to the MODIS bands 1 to 8, and two VIs (NDVI and EVI). Additionally, the product provides auxiliary layers of backward scattering and forward scattering, including part of the bands (description in Sect. 5).

2.2.1 The nadir-normalized (NAD) data

In order to minimize the differences in sun sensor geometry between the MODIS scenes and generate the NAD dataset, the daily surface reflectance data were normalized to a fixed 45^∘ SZA and to the nadir observation (VZA = 0^∘) using the BRDF and the RTLS model (Wanner et al., 1995). The parameters of the RTLS BRDF model are part of the MAIAC product suite (MCD19A3 product) reported every 8 d. The MAIAC algorithm detects significant land cover changes (e.g., fire and deforestation) within the 8 d period and does not use those observations for the BRDF inversion (Lyapustin et al., 2018). A minimum of three observations in the 8 d window was required to accurately model the signal. The closest RTLS parameters in time were used to normalize the daily data. The normalized bidirectional reflectance factor (BRFn) for the NAD surface reflectance (SZA = 45^∘, VZA = 0^∘, and RAA = 0^∘) was calculated using Eq. (1) as follows (Lyapustin et al., 2018):

$$\begin{array}{}\text{(1)}& \text{BRFn}=\text{BRF}\times {\displaystyle \frac{k^{\mathrm{L}}+F_{\mathrm{0}\mathrm{V}}\times k^{\mathrm{V}}+F_{\mathrm{0}\mathrm{G}}\times k^{\mathrm{G}}}{k^{\mathrm{L}}+F_{\mathrm{V}}\times k^{\mathrm{V}}+F_{\mathrm{G}}\times k^G}},\end{array}$$

where k^L, k^V, and k^G are the BRDF isotropic, volumetric, and geometric optical kernel weights, respectively. F_0V and F_0G are the BRDF kernel values for the given geometry listed in Table 2. F_V and F_G are the kernel values of the RTLS model for the specific MODIS observation, respectively (Lyapustin et al., 2018). F_V and F_G values are available at 5 km cells and were resampled to 1 km using the nearest-neighbors method to match the spatial resolution of the spectral bands. This resampling step does not create spatial artifacts in the data because the geometry changes slowly over time (Lyapustin et al., 2018).

Table 2 Viewing angle normalizations and corresponding BRDF kernel values.

We aggregated the normalized daily data into monthly composites by keeping the median values for each pixel. During the temporal aggregation, we also calculated the per-pixel number of samples (or observations) for each monthly composite, which can be used as auxiliary data to filter pixels with a low number of observations (less reliable estimates of surface reflectance). The tiles were mosaicked for the entire South American continent and then re-projected from the original sinusoidal projection to the geographic coordinates system (datum WGS 84; EPSG: 4326). The output spatial resolution corresponded to 0.0091^∘, which is approximately equivalent to 1 km in projected coordinates.

We also calculated two traditional vegetation indices, namely NDVI (Rouse et al., 1974; Eq. 2) and EVI (Huete et al., 2002; Eq. 3).

$$\begin{array}{}\text{(2)}& {\displaystyle }\text{NDVI}={\displaystyle \frac{\mathit{\rho }\text{NIR}-\mathit{\rho }\text{Red}}{\mathit{\rho }\text{NIR}+\mathit{\rho }\text{Red}}}\text{(3)}& {\displaystyle }\text{EVI}=\mathrm{2.5}\times {\displaystyle \frac{\mathit{\rho }\text{NIR}-\mathit{\rho }\text{Red}}{\mathit{\rho }\text{NIR}+\left(\mathrm{6}\times \mathit{\rho }\text{Red}-\mathrm{7.5}\times \mathit{\rho }\text{Blue}\right)+\mathrm{1}}},\end{array}$$

where ρ is the surface reflectance of a MODIS band, ρNIR is the NIR reflectance (band 2), ρRed is the red reflectance (band 1), and ρBlue is the blue reflectance (band 3). The constants in Eq. (3) (bands 6, 7.5, 1, and 2.5) represent the aerosol coefficient adjustment of the atmosphere for the red and blue bands, the adjustment factor for the soil, and the gain factor, respectively (Huete et al., 2002).

2.2.2 The anisotropy (ANI) data

For the ANI data, the daily surface reflectance data were first normalized to two viewing angles at the backward (SZA = 45^∘, VZA = 35^∘, RAA = 180^∘) and forward (SZA = 45^∘, VZA = 35^∘, RAA = 0^∘) scattering, using Eq. (1) and values from Table 2. The VZA was set to the near-hot spot (VZA = 35^∘) instead of the actual hot spot (VZA = 45^∘) to keep VZA closer to the actual range of MODIS observations across the South American continent and minimize errors coming from the extrapolation of the BRDF (de Moura et al., 2015). The standard deviation for this modeling was thoroughly investigated in a previous study and determined to be 10 % of the observed variation in anisotropy (de Moura et al., 2015). Furthermore, we aggregated the backward- and forward-scattering data temporally into monthly composites following the same procedures as before for the NAD data. We then calculated the NDVI and EVI for each viewing angle normalization. Finally, we obtained the difference between backward- and forward-scattering estimates for each of the eight MODIS bands, in addition to the NDVI and EVI, effectively generating the ANI layers as follows (Eq. 4; de Moura et al., 2015):

$$\begin{array}{}\text{(4)}& {\text{ANI}}_i={\text{Backward}}_i-{\text{Forward}}_i,\end{array}$$

where i is the spectral band or VI selected in the calculation.

2.3 Algorithm and computation

All data processing was done in R v4.0.2 (R Core Team, 2016), and the code is available from GitHub (https://github.com/ricds/maiac_processing, last access: 16 January 2023; Dalagnol and Wagner, 2022). Besides processing the AnisoVeg product from the daily MAIAC MODIS data, the code can also generate 16 or 8 d temporal composites, mosaics, and VIs. Although we focused on South America when developing AnisoVeg, the code can readily be adapted to process data for other parts of the world and generate corresponding NAD and ANI layers. Below, we provide the computer specification for anyone who wishes to process the data independently.

For the presented dataset, the computation was performed using a HP Z840 workstation with Intel Xeon CPU E5-2640 V3 (2.60 GHz; 32 cores) and 64 GB (gigabytes) RAM. The daily MODIS data for the whole South American continent from 2000 to 2021 accounted for 6.69 TB (terabytes). Processing monthly composites is computationally intensive due to loading all of the daily data for each month at once for a given tile. Thus, the main bottlenecks are for RAM and hard drive writing speed. For the workstation with 64 GB memory, the usage of 10 cores running in parallel processing was the optimal choice. The average processing time of each monthly composite for one tile was 6 min. Therefore, it took 26.2 h to process the 262 composites (March 2000 to December 2021) for each tile. Since we had 31 tiles covering South America, the total amount of time to process one viewing normalization was approximately a month (33.8 d). Consequently, the total time spent in computation was 101.5 d for processing the three viewing normalizations (nadir, backward, and forward scattering) and generating the NAD and ANI layers. Processing can also be done with less potent computers with a minimum of 16 GB RAM memory and four processing cores.

2.4 Time series availability and uncertainty

The monthly compositing process returned a time series dataset over all of the South American continent, with an average of 242±35 out of a maximum of 262 composites (period between March 2000 and December 2021) for each pixel, with some data missing due to a lack of high-quality observations (Fig. 2). Only 34.3 % of the available pixels have the full time series (262 composites). The Amazon region shows a lower mean number of samples in the time series, with an average of 231±29 composites, which can be seen in Fig. 2. This lower number of samples is due to the innate high cloud cover (Durieux et al., 2003). It is important to note that the AnisoVeg product was strictly created to analyze the land surface and does not cover waterbodies. Moreover, the period between March 2000 and June 2002 has higher amounts of missing data because it preceded the launch of the Aqua satellite. When data from both satellites (Terra and Aqua) were combined to create the product after 2002, we had a much better pixel-level data availability to produce dense time series. Although we have a dense time series across the Amazon Rainforest (Fig. 2a), the mean number of daily observations within a month for this region is relatively lower than that observed in more dry and seasonal regions of South America (Fig. 2b). Thus, we suggest using a number of sample layers as a proxy for the uncertainty in the retrieval of monthly composites to filter out pixels with a low number of samples (e.g., less than three observations per composite). The lower the number of samples 1 pixel has, the higher the uncertainty in the data analysis. Although we use the median values to aggregate observations within months and mitigate potential land cover changes, stand-replacing changes may cause inaccurate anisotropy estimates for the given monthly estimates. Hence, we advise filtering data for land use and land cover changes before using them to obtain the most accurate anisotropy estimates.

Figure 2 AnisoVeg time series availability and uncertainty over South America. (a) The number of composites in the time series representing pixel availability. The maximum number of composites in the time series is 262 for the period between March 2000 and December 2021. (b) Mean number of daily observations within a month used to create the monthly composites as a proxy for uncertainty. The maximum daily observations in a composite is 60 (twice a day every day for a month).

3 Spatial and temporal distribution of NAD and ANI data across the Amazon forest

We selected three experimental areas in the Brazilian Amazon Rainforest to show the spatial and temporal distribution of NAD and ANI data (rectangles in Fig. 1). These areas show old-growth rainforests with a distinct canopy structure and aboveground biomass (AGB) stocks. The AGB increases from semi-deciduous forests at Xingu Indigenous Park (190±19 Mg ha⁻¹) and open ombrophilous forests with lianas at the São Félix do Xingu (241±31 Mg ha⁻¹) to dense ombrophilous forests at the Tapajós National Forest (288±38 Mg ha⁻¹), as estimated by the ESA/CCI (Climate Change Initiative) AGB map from 2017 (Santoro and Cartus, 2021). These are large-scale AGB estimates and may underestimate the true AGB at higher values such as in the Tapajós site. These three sites are also expected to show different phenological dynamics because their selected pixels cover distinct phenoregions in the study reported by Xu et al. (2015).

When compared to the nadir-normalized EVI (EVI_NAD) images (Fig. 3a–c), the anisotropy EVI (EVI_ANI) data showed different spatial patterns across sites (Fig. 3d–f). While the forests over the three sites showed approximately similar EVI_NAD values (EVI_NAD≈0.50; Fig. 3a–c), they showed more variability in EVI_ANI between the Xingu Indigenous Park (EVI_ANI>0.20), São Félix do Xingu (EVI_ANI>0.24), and Tapajós (EVI_ANI>0.27) sites (Fig. 3d–f). This increase in EVI_ANI between sites goes into the same direction of the AGB gradient observed from the Xingu Indigenous Park to the Tapajós National Forest. This result may indicate different forest canopy structures that were not captured in the EVI_NAD observations but were captured by the EVI_ANI. Overall, the EVI_ANI is high over forests (0.20 to 0.30) and low over pastures and crops (less than 0.10). This means large anisotropy between the reflected energy in backward- and forward-scattering MODIS directions due to the structural complexity of forest canopies. The association between anisotropy and forest canopy structure has been previously shown for the same region in a previous work (de Moura et al., 2016).

Figure 3 The spatial distribution in August 2020 (dry season) of the nadir-normalized enhanced vegetation index (EVI_NAD) is shown in panels (a), (b), and (c) for Tapajós National Forest, São Félix do Xingu, and Xingu Indigenous Park, respectively. The corresponding results for the anisotropy EVI (EVI_ANI) are shown in panels (d), (e), and (f), respectively. The triangles plotted over panels (a–c) indicate the sites used to obtain the profiles in Fig. 4.

From the comparison of different sites (triangles in Fig. 3a), we observed that the mean EVI_NAD signal over the time period did not vary much between the selected forests, while the EVI_ANI varied greatly (Fig. 4) for Tapajós (mean EVI_NAD=0.49; mean EVI_ANI=0.27), São Félix do Xingu (mean EVI_NAD=0.51; mean EVI_ANI=0.24), and Xingu Indigenous Park (mean EVI_NAD=0.51; mean EVI_ANI=0.22). Moreover, EVI_NAD and EVI_ANI values were moderately positively correlated at Tapajós ($r=+\mathrm{0.37}$), weakly correlated at São Félix do Xingu ($r=+\mathrm{0.06}$), and moderately negatively correlated at the Xingu Indigenous Park ($r=-\mathrm{0.28}$). The EVI_NAD and EVI_ANI seasonal variability and phase correlation changes from site to site, suggesting that different canopy dynamics processes are likely being captured by the two metrics at the three sites. Understanding exactly what those effects mean for these forests is beyond the scope of this paper. However, it indicates open avenues for studying forest functioning using these products. For example, previous studies have shown that EVI_NAD metrics captured different compositions of leaf ages in the canopies of the central Amazon (Gonçalves et al., 2020).

Figure 4 Time series of AnisoVeg's MODIS enhanced vegetation index (EVI) from 2000 to 2021 for old-growth forests of (a) Tapajós National Forest, (b) São Félix do Xingu, and (c) Xingu Indigenous Park. The black line indicates the nadir-normalized signal (NAD layer), while the red line represents the EVI anisotropy (ANI layer). The profiles are the mean value of 3×3 pixels, whose locations are indicated by triangles in Fig. 3.

To demonstrate the potential of AnisoVeg for a large-scale forest structure inference, we compared the NAD and ANI data against forest height measurements from the Global Ecosystem Dynamics Investigation (GEDI) lidar sensor. We found that EVI_ANI was able to explain up to 55 % of the height variability in the Amazon forest, according to a simple linear relationship (R²=0.55; p<0.01; Fig. 5). This is a very strong predicting power for a single variable, considering a simple linear model, especially for satellite passive optical data which are often underrated for forest structure estimates in comparison to synthetic aperture radar (SAR) data. EVI_NAD was significantly but weakly associated with height variability (R²=0.16; p<0.01), reinforcing the increase in the explanation power owed to the anisotropy metrics built from multi-angle observations. The height data were derived from the GEDI lidar sensor aboard the International Space Station. They were obtained more specifically from the product GEDI L2A elevation and height metrics data version 2 (footprint size 25 m) acquired from April 2019 to October 2020 (available dates at the time of download). GEDI data were downloaded from Earthdata Cloud service system (https://earthdata.nasa.gov, last access: 16 January 2023). We selected the relative height metric at the 98th percentile (RH98), which represents the top canopy height. The selected RH98 metric was averaged over each 1 km grid cell and filtered using a threshold of greater than or equal to 50 shots per kilometer squared to have a high confidence of reliable height estimation representing the 1 km mean. The AnisoVeg data used for this comparison were based on the same time period as GEDI and filtered for EVI_NAD larger than 0.35 to exclude non-forested areas. While we only showed the plot for the strongest EVI_ANI : GEDI relationship in June 2019 (Fig. 5), the other months also showed significant (p<0.01) and strong relationships with R², ranging from 0.36 to 0.55 (mean R²=0.46). Future studies should explore relationships using ANI from different months and other indices, either alone or in combination with each other, to further understand their significance for explaining forest structure. This is important to determine how the anisotropy data can contribute to aboveground biomass and carbon estimates in conjunction with other sources of data such as those from SAR sensors.

Figure 5 Relationship between forest height (GEDI mean RH98) and two AnisoVeg layers obtained in June 2019 over the Amazon. (a) EVI_NAD and (b) EVI_ANI. The RH98 metric is the relative height at the 98th percentile, which represents the top of canopy height. A total of 7000 random matching pixels were used in this analysis (1 % of 700 000 total matching pixels available), resulting from the filtering of both GEDI and AnisoVeg data. The red line indicates the fitted line by a simple linear model.

Terrain illumination is a factor of spectral variability which can affect EVI_NAD determination and its relationship with biophysical attributes of vegetation, as shown by previous literature (Huang et al., 2010; Chen and Cao, 2012). Even at 1 km spatial resolution, EVI_ANI results of Figs. 3, 4, and 5 can be affected to some extent by terrain illumination effects observed locally at some sites. For instance, topographic effects on EVI_ANI probably occurred at the São Félix do Xingu site, where topographic roughness, observed in SRTM (Shuttle Radar Topography Mission) data (results not shown), was coincident with increased EVI_ANI values in Fig. 3e. Furthermore, even in relatively flat terrains, variations in topographic aspect (surface orientation to the Sun) can affect the EVI variability in MODIS data because of the different amounts of energy reflected in the NIR towards the sensor by inclined surfaces in the forward- and backward-scattering view directions. Such effects have been observed in southern Brazil with MODIS at 250 m spatial resolution and increased in magnitude at higher spatial resolution data obtained by other sensors (Galvão et al., 2016). Therefore, it may prove useful to include topographic variables in modeling exercises to offset these effects.

In a prospective analysis, we also explored the behavior of the two EVI AnisoVeg metrics over the Amazonian phenoregions mapped by Xu et al. (2015). The EVI_NAD and EVI_ANI monthly means over different phenoregions highlighted the strong heterogeneity of the Amazon forest (Fig. 6). For instance, the profiles showed strong differences between both metrics from January to September in a phenoregion with well-defined dry and wet seasons (phenoregion 1 in Fig. 6a at Xingu Indigenous Park). Large differences between EVI_NAD and EVI_ANI were also observed in some phenoregions without a very long dry season in the northwestern Amazon (phenoregion 5 in Fig. 6e). On the other hand, EVI_NAD and EVI_ANI showed temporal decoupling in phenoregion 3 located in the central–eastern Amazon (Fig. 6c). Overall, while the seasonality of EVI_NAD has been investigated by many studies in the past, the seasonality of EVI_ANI is something to be further explored with the support of auxiliary data (e.g., airborne lidar and field campaigns). This is important to better understand the differences in seasonal patterns between both AnisoVeg metrics.

Figure 6 Monthly means of EVI_NAD (black) and EVI_ANI (red) for nine phenoregions mapped by Xu et al. (2015) in the Amazon. The phenoregions are shown in increasing order from 1 to 9 in corresponding panels (a) to (i). They represent forests with similar seasonality and landscape structure. Solid lines and shaded areas represent the mean and 95 % confidence interval around the mean. The values were extracted from 20 years of data (from 2001 to 2021) for 100 random coordinates within each region and extracted from 3×3 windows of pixels.

Table 3 Examples of other multi-angular anisotropy indices that can be further calculated using the layers of the AnisoVeg product. Lambda represents the selected spectral band or vegetation index. N, B, and F represent nadir view normalization, backward-scattering, and forward-scattering estimates, respectively.

^* ANI is included in the AnisoVeg product. Source adapted from Sharma (2021).

4 Code and data availability

All code is available from the Zenodo repository (https://doi.org/10.5281/zenodo.6561350; Dalagnol and Wagner, 2022). The full dataset can be found at the official AnisoVeg Zenodo repository (https://doi.org/10.5281/zenodo.3878879; Dalagnol et al., 2022). The dataset was organized in compressed files (*.zip format), sub-divided by years (currently 2000–2021) and layers (bands 1–8, NDVI, and EVI) for both nadir normalization (code = NAD) and anisotropy (code = ANI). The number of sample layers (code = NO_SAMPLES) is also provided. Inside each compressed file are 12 image files (*.tif format), with one per month, except for the year 2000, which starts in March. The storage size for the whole dataset is 162.6 GB. The data have a scale factor of 10 000 to reduce the file storage size. Thus, to obtain surface reflectance values of bands or correct the range of values for indices, you should divide the layers by 10 000. The exception is the number of samples, which already shows the correct range of values (from 0 to 60 observations). It is planned to update the dataset on a yearly basis. Auxiliary data that allow the calculation of other anisotropy metrics (listed in Table 3) are included in two separate Zenodo repositories for backward (https://doi.org/10.5281/zenodo.6040300; Dalagnol, 2022a) and forward scattering (https://doi.org/10.5281/10.5281/zenodo.6048785; Dalagnol, 2022b), including the selected layers (red, NIR, NDVI, and EVI). The EVI_ANI and EVI_NAD layers were also uploaded to the GEE platform using the geeup tool v0.5.3 (Roy, 2022). They can be accessed through the GEE ImageCollection assets (projects/anisoveg/assets/evi_anisotropy and “projects/anisoveg/assets/evi_nadir), which can be found at https://code.earthengine.google.com/?asset=projects/anisoveg/assets/evi_anisotropy (last access: 16 January 2023, Dalagnol, 2022d) and https://code.earthengine.google.com/?asset=projects/anisoveg/assets/evi_nadir (last access: 16 January 2023, Dalagnol, 2022c).

5 Prospective use of the dataset

The NAD layers from the AnisoVeg product have been used in previous studies to explore the climate drivers of the Amazon forest greening (Wagner et al., 2017), the large-scale Amazon forest sensitivity to drought (Anderson et al., 2018), the structure and dominance of bamboo species in the southwestern Amazon (Dalagnol et al., 2018), the productivity in a flooded forest in the eastern Amazon (Fonseca et al., 2019), the productivity and relationship with sun-induced fluorescence over the Brazilian Caatinga biome (Bontempo et al., 2020), the relationships in the leaf age demography in central Amazon (Gonçalves et al., 2020), and the relationships between fire disturbance and SAR-based vegetation optical depth in southern Amazon (Zhang et al., 2021).

The ANI layers from the AnisoVeg product have been mainly used to characterize the Amazon forest structure properties (de Moura et al., 2015, 2016). These layers now open new avenues of investigation for vegetation, including (but not limited to) the characterization of biophysical attributes of forests, including their seasonality and trends, the assessment of changes in vegetation structure due to natural disturbances or degradation (logging, fire, and edge effects), and the evaluation of forest health and productivity (greenness and browning). We expect that this dataset contributes to upscaling studies over large areas of key forest properties, such as the AGB and canopy roughness (Foody and Curran, 1994; Saatchi et al., 2008). This information is required for dynamic vegetation models to accurately represent the carbon cycle. This dataset is not limited to the study of Amazonian forests and can be used to explore other biomes of South America, such as the Atlantic Forest, savannas (Cerrado), Caatinga, Chaco, Pantanal, and Pampas. Such studies could improve our understanding of large-scale vegetation functioning, carbon storage, and cycling. Ultimately, they can contribute to refining global ecosystem models and obtaining accurate estimates of the carbon cycle in response to climate and environmental change. Furthermore, auxiliary backward- and forward-scattering data are also available with the dataset. Beyond the use of the provided ANI layers, this effectively allows the computation of several other multi-angular anisotropy indices from the literature (Table 3). The advantage (or disadvantage) of one specific anisotropy index rather than others is not established in the literature, given the range of vegetation applications and the lack of available datasets to date. We calculated and provided only ANI, due to its demonstrated relationships with Amazon forest structure and functioning (de Moura et al., 2015, 2016; Hilker et al., 2017). However, we expect other indices, including ratios and normalized differences between the backward and forward scattering components and offer additional possibilities for tropical vegetation studies which should be explored in future studies.