Mapping the spatial distribution of stand age and aboveground biomass from Landsat time series analyses of forest cover loss in tropical dry forests

Abstract

Spatial information on the timing of forest cover loss is important to identify and map stand age, which is a key factor driving the recovery of carbon pools and can also be used to estimate aboveground biomass (AGB) based on its relationship with stand age. Here, we estimated the spatial distribution of stand age and AGB of young forest (<20 years) in three types of tropical dry forest in the Yucatan peninsula using Landsat NDVI (normalized difference vegetation index) time series from 2000 to 2020. We estimated AGB based on chronosequence data and compared these results to reference field data and estimations obtained from remote‐sensing studies. The overall and user accuracy of the age map was high (95.7–99.9% and 87.35–98.5% respectively). However, lower producer accuracy values (from 31.2 to 67.2%) suggest an underestimation of the extension of young forests. We found a greater extent of young forests in the semi‐deciduous and deciduous forests compared to the semi‐evergreen ones. Mean AGB estimated from stand age (53.1 Mg ha⁻¹) was lower than that estimated from remote‐sensing studies (67.5 to 95.2 Mg ha⁻¹). These results indicate that spatial information of forest age can be accurately assessed from Landsat time series, and that the combination of stand age with chronosequence data can reduce the overestimation of AGB of recovering forests commonly found in remotely sensed data.

Introduction

Habitat loss, and particularly deforestation, is one of the most important environmental problems in the Anthropocene. Deforestation alters biogeochemical cycles, including the carbon cycle, thereby compromising the potential for carbon sequestration and climate change mitigation (Portillo-Quintero et al., 2014). It also leads to the loss of biological diversity (Malhi et al., 2008), which is crucial for the resilience of forest landscapes in the age of global change (Thompson et al., 2009). The loss of forest cover is most extensive in tropical areas due to the conversion of forest to agriculture or urban settlements (Keenan et al., 2015), especially in tropical dry forests (Allen et al., 2017; Bhaskar et al., 2018), which are responsible for approximately half of the carbon storage capacity in tropical forests worldwide (Keith et al., 2009). However, forests can regenerate after a disturbance event and, depending on the intensity and type of disturbance, attributes such as species richness (Norden et al., 2009) and aboveground biomass (AGB) can recover and approach pre-disturbance conditions within a few decades (Chazdon, 2008; Poorter et al., 2016). Other attributes such as species composition recover at much slower rates or might never reach primary forest levels (Chazdon, 2008; Rozendaal et al., 2019). Therefore, the time since the last disturbance (year of latest forest cover loss) can be related to stand age, and since AGB is strongly related to successional age (Read & Lawrence, 2003), it is possible to use this information to estimate the spatial distribution of AGB in recovering forests.

The spatial distribution of AGB has been estimated at multiple scales using a wide variety of remotely sensed products and methods (Avitabile et al., 2016; Hernández-Stefanoni et al., 2020; Rodriguez-Veiga et al., 2016; Saatchi et al., 2011). However, there is a clear tendency for remote sensing methods not only to underestimate AGB in mature forest due to a saturation of remotely sensed data, but also to overestimate the AGB in young forests (Hernández-Stefanoni et al., 2020; Rodriguez-Veiga et al., 2016). Forest age is a key factor driving the recovery of carbon pools, and thus, carbon balance in tropical forests (Coursolle et al., 2012; Gao et al., 2016). Therefore, the estimation of the spatial distribution of young forest could accurately reflect the changes in vegetation attributes such as AGB that take place during early succession. In addition, this information could help to locate recovering forests, which have a large carbon sequestration potential (Poorter et al., 2016).

The Yucatan peninsula holds the second largest extent of continuous forested land in Latin America after the Amazon. The forests here are characterized by presenting a gradient of forest cover types which transition from deciduous tropical dry forests in the drier northern part of the Peninsula to semi-deciduous and semi-evergreen forests in the more humid south east (Dupuy et al., 2015). Forests in the Yucatan Peninsula have had a long and unique history of land cover changes which varies spatially leading to a mosaic of patches of forest in different stages of recovery (stand age), and other land uses. Additionally, there is a high incidence of natural disturbances such as hurricanes and forest fires, which are also important drivers of forest dynamics (Mascorro et al., 2016). Obtaining spatially explicit information of the year of the most recent forest cover loss in these ecosystems could provide a proxy for stand age. In particular for disturbances that are patchy and generally of small spatial and/or temporal scale, such as those associated with swidden agriculture, extraction of forest products, hurricanes and fires, which characterize most of the forests of the Yucatán Peninsula. This could help to identify areas of young recovering forests, thus aid conservation efforts for, and landscape-level management of these forests. However, obtaining exact information on the spatial distribution of stand age can be a challenging task especially for tropical forests, where field-based information of stand age is often obtained from interviews to local informants, and it is often unavailable for large extents.

Remote sensing offers methods for obtaining extensive historical data that have been used to overcome this challenge. The Landsat archive provides a global archive of 30-m resolution imagery from March 1984 to date, with a temporal resolution of 16 days (USGS 2020). A number of land cover change detection algorithms have been developed, taking advantage of the temporal extent of the Landsat Archive, to locate areas where a loss of forest cover has occurred by constructing time series of images and identifying changes in forest cover from time x to times x + 1, x + 2, and so on. Changes in forest cover are usually tracked by a measure of ‘greenness’, identified by vegetation indices, in ‘time x’ to sudden decreases or absences of ‘greenness’ in times x + 1, x + 2, and so on for the same location. However, in seasonally dry tropical ecosystems differentiating these changes can be challenging due to the confounding effect of seasonal fluctuations in ‘greenness’ associated with leaf abscission during the dry season.

Two common change detection algorithms are the Vegetation Change Tracker (VCT) (Huang et al., 2010), developed by NASA and the LandTrendr (Landsat-based Detection of Trends in Disturbance and Recovery) (Kennedy et al., 2010) developed by Oregon State University. These algorithms identify changes in land cover across large areas. However, these algorithms fail to capture intra-annual variability, an important characteristic in seasonal environments such as tropical dry forests, as they accept only one image per year. Therefore, even if all images are acquired at the same time of the year inter annual variation in leaf phenology can bias estimates of land cover change in seasonal forests. Moreover, global estimations of forest cover loss such as that of Hansen et al. (2013) often rely on generic definitions of forest without accounting for differences in forest types and plantations, which can lead to systematic bias due to differences in structure and seasonal patterns. The Breaks for Additive Seasonal Trend (BFAST) algorithm, is a change detection algorithm able to process dense time series constructed with diverse vegetation indices using all available imagery within a year, thus enabling the accurate characterization of seasonal changes (Verbesselt et al., 2010). The algorithm fits a harmonic model with a trend and a seasonal component for each of the pixels within the images in the time series, and then flags the pixels whose estimated values deviate from the model predictions. This method shows great potential for producing highly accurate spatial estimations of forest loss, and thus also of forest age, in highly seasonal tropical areas where small-scale swidden agricultural practices are common (Smith et al., 2019).

The objectives of this research were twofold. First, we aimed to produce accurate estimations of the year of the last forest cover loss from 2000 to 2020 as a measure of stand age in young forest stands (<20 year old) in three tropical dry forest landscapes of the Yucatan Peninsula. Second, to estimate AGB in these landscapes based on stand age and chronosequence data (relating both variables) and to compare these estimations with those of previous remote sensing studies. We estimated land cover changes in different types of tropical dry forest with slightly different land-use histories using the BFAST algorithm (Verbesselt et al., 2010). Swidden agriculture is the predominant land use in the semi-deciduous and deciduous forests, whereas cattle ranches and selective logging are more common in the semi-evergreen forest. Therefore, we expected the timing and extension of young forest to differ across sites with a greater extension of young stands in the drier forest sites, due to the cyclic nature of swidden agriculture, which involves a relatively short (< 20 years) fallow period. We also test whether information on the spatial distribution of forest age could reduce the overestimation of AGB in young forest that has previously been found in different remote-sensing studies (Cartus et al., 2014; Hernandez-Stefanoni et al., 2020; Rodríguez-Veiga et al., 2016) by capturing forest recovery time. We expected a reduction in the overestimation of AGB in young forest stands (< 20 years old) based on BFAST-derived stand age, compared to direct estimates from satellite imagery in previous studies.

Materials and Methods Study area

Three 3,600 km² square areas were selected in the Yucatan Peninsula (Fig. 1). The first area is located in ‘El Palmar’ (EP) State Reserve in the NW part of the peninsula in the State of Yucatán. This is the driest site with mean annual precipitation ranging from 500 to 700 mm per year. The predominant vegetation type is deciduous tropical dry forest (>80% of trees shed their leaves during the dry season) characterized by small stature (8−10 m) tree species (Flores & Espejel, 1994; García & Contreras, 2011). The development of henequén (Agave fourcroydes) plantations in the area and their subsequent abandonment and the practice of traditional swidden agriculture has resulted in a mosaic of forest patches of different successional age (Wyman et al., 2007). The second site ‘KK’ is located in the middle portion of the peninsula and it is named after a private reserve ‘Kaxil Kiuic Biocultural Reserve’. The dominant vegetation type is medium size (10–15 m) semi-deciduous tropical dry forest (50–75% of trees shed their leaves during the dry season). Mean annual precipitation ranges from 900 to 1100 mm. The most common land use in this area is swidden agriculture. Precipitation at these two sites is highly seasonal with a long dry season typically from November to April and a wet season from May to October when most of the precipitation occurs. The third site is ‘FCP’ located in the middle portion of the state of Quintana Roo. The dominant vegetation corresponds to semi-evergreen tropical forest, with taller trees (15–25 m) and a lower proportion (<25%) of deciduous tree species. This area is also characterized by the presence of patches of low seasonally flooded forest. Mean annual precipitation ranges from 1100 to 1300 mm per year, making this the most humid of our field sites (Fig. 1). Land use in this site includes selective extraction of commercially important woody species, pastures for cattle (Armenta-Montero et al., 2020), and small-scale swidden agriculture (Ellis et al., 2017b; Mascorro et al., 2016).

Figure 1. (A) Location of three study sites in the Yucatan Peninsula: EP: El Palmar; KK: Kaxil Kiuic; FCP: Felipe Carrillo Puerto. Sites are located in the three most widespread forest cover types of the Yucatan Peninsula. (B) Location of the Yucatan Peninsula.

Satellite Imagery acquisition and preprocessing

Multispectral Landsat imagery from January 1985 to Aug 2020 in Landsat Collection 1, for sensors Landsat TM, ETM+, and OLI, were obtained using the NASA Earth Explorer search engine. The search was filtered for daytime images with less than 70% cloud cover. All surface reflectance bands, NDVI, and ‘pixel qa’ (quality assessment) layers were selected. We downloaded 552 images for the semi-deciduous and deciduous sites (path 20, row 46) and 654 for our semi-evergreen site (path 20, row 45).

Landsat scenes were cropped to the extent of our study area and preprocessed in R (package ‘bfastSpatial’, Verbesselt et al., 2010, 2012). Then, images were filtered using Landsat ‘pixel qa’ (pixel quality assessment) layer to produce images with cloud cover under 70% for our study areas. Finally, the filtered images were stacked using the ‘timestack’ function, which constructs a raster of stacked, filtered NDVI scenes organized by date. This stack is the main input for the ‘bfmSpatial’, the main function which performs the time-series analysis. For detailed code and other parameters see Supplementary material S1.

NDVI was used due to its ability to differentiate photosynthetically active vegetation from other elements in the landscape by identifying the presence of chlorophyll A. Therefore, this index, which is obtained from the ratio between the sum of and the difference between red and near-infrared bands, is able to identify pixels with active vegetation. Other advantages of the NDVI are its efficiency in solving some limitations of using individual bands such as shadowing and topographic effects.

Algorithm implementation: Iterative monitoring of forest cover loss

Breaks for Additive Seasonal Trends (BFAST) algorithm was applied using the ‘bfmSpatial’ function from ‘bfastSpatial’ package in R (Verbesselt et al., 2010, 2012). The algorithm uses every pixel with cloud-free observations throughout all scenes in a time series and fits a regression model with a seasonal and a trend component (Verbesselt et al., 2010). This feature enhances the robustness of the time series and reflects seasonal systems in a more accurate way (Verbesselt et al., 2010). A first order-harmonic model was used for a better reflection of the trajectory of a pixel throughout a time series in a highly seasonal system where changes in precipitation and phenology are a normal component of the system’s functionality (Devries et al., 2015).

An ordinary least squares (OLS) residual-based moving sum (MOSUM) test, was selected to test whether the image pixels can be flagged for breakpoints. When the test indicated a significant deviation of observations from the model predictions (p < 0.05), the breakpoints were flagged, producing a layer with breakpoints, or pixels containing the date where the change occurred (Verbesselt et al., 2010). The algorithm was applied in an iterative way following Devries et al. (2015) by selecting a short monitoring period (1 year) using each year as the monitoring start year (start = 2000, 2001, and so on) and its subsequent year (end = 2001, 2002, and so on) as the end of the monitoring period. Due to restrictions in the availability of Landsat imagery in the early dates of Landsat archive this algorithm also used all historical imagery available from 1985 to the ‘start’ date to calibrate the model (history = “all”).

The ‘bfmSpatial’ function produces three to six layers of information from the timestack. For a detailed description, see Verbesselt et al. (2010). We used the first layer ‘breakpoints’ (which contains all pixels flagged where the algorithm found significant deviations from the model fit within the monitoring period) and the second layer ‘magnitude’ (the intensity of change found relative to the OLS-MOSUM model fit) to produce a new layer containing the magnitudes in all the pixels where a breakpoint occurred (Supplementary material S1). With this approach, we ensure that all data from previous years are used in order to track small-scale changes and that more than one breakpoint can be identified in any given pixel, so we were be able to detect changes that took place in the same pixel in each year throughout the length of the time series.

Magnitude threshold selection

‘Magnitudes’ represent the difference between the observed and predicted NDVI values from the harmonic model and are standard output from the bfmSpatial function. Negative magnitudes have been found to be related to a reduction in tree cover in previous analyses using BFAST (Devries et al., 2015; Verbesselt et al., 2012). Therefore, a calculation of the magnitude threshold tailored to our study area improves the sensitivity of the analysis to target forest loss. We used binomial-logistic regression between deforested sites (1) and no deforestation sites (0) for selecting a forest cover loss threshold following Devries et al. (2015). We used 440 points placed systematically within confirmed forest cover loss and no change locations to obtain the magnitude value with a 50% probability to result from forest cover loss. A dataset with the selection of confirmed forest loss pixels was obtained using the locations where magnitudes <0.06 were obtained from the BFAST algorithm, and confirmed deforestation was identified by visual interpretation in QGIS 3 and Google Earth Pro. This analysis was carried out by the ‘glm’ function (family = ‘binomial’) in R.

Production of forest-cover loss map and validation

For each year, a raster clumping analysis was used in order to filter single-pixel changes, which are likely to derive from noise in the time series (Devries et al., 2015; Verbesselt et al., 2012). We selected a minimum mapping unit of 0.5 ha as it is the smallest extension commonly used for swidden agriculture in the area (Garcia-Frapolli et al., 2007; Mascorro et al., 2016) for which images were clustered using ‘Rook’s case adjacency’ to join clusters of pixels with a minimum extension of 0.5 ha (‘areasieve’ in R).

All images were overlayed starting from the most recent year (2020) to the oldest year (2000) using ‘r.patch’ function in QGIS (GRASS GIS 7.2) to obtain a single raster containing all the years of possible deforestation. This method was selected to show the last forest-cover loss estimated in all pixels in our study area, which was assumed to be the year 0, when vegetation regrowth started. Forest age was calculated as the time since the last forest-cover loss event. Land use changes from forest to agriculture or urban areas were not considered in this study, since they imply more permanent changes that do not affect our estimates of stand age in young (<20 year old) forests. These areas were masked using a land cover map available for each site (Hernandez-Stefanoni et al., 2020). The total area was calculated as the sum of the area occupied by pixels flagged as forest cover loss in all sites.

A limitation on the validation of cover-change maps is the low availability of reference imagery, an issue that is most accentuated in the tropics where cloud cover persists for large parts of the year (Asner, 2001). For this reason, two points in the time series were selected according to the availability of reference imagery. Reference datasets are comprised of GeoEye historical imagery (0.5 m resolution) for years 2008 through 2015; RapidEye imagery (5 m) available for the year 2012; and Landsat Tree Cover Product (30 m) (Sexton et al., 2016) for the years 2010 and 2015. Previous work has indicated that the BFAST algorithm identified 81% of the area of forest cover loss with a maximum difference of 18 months in tropical dry forests of the Yucatan peninsula (Smith et al., 2019). Therefore, in this study we considered 1.5 years as a maximum temporal lag between observed and estimated change.

In each validation year, the proportion of area of forest cover loss and no change classes was calculated. We selected 500 samples using Cochran’s sample size formula (Cochran, 1963) with a stratified sampling design (Olofsson et al., 2013, 2014). Samples were allocated designating a minimum of 100 samples for the smallest class size (forest cover loss) and 400 in the no change class. In this process, each validation unit (sample) corresponds to a Landsat pixel. Each pixel in the validation dataset was located in the reference dataset and examined in the same year, before and after the suspected event. We considered a change event as confirmed wherever forest cover was present in the year before and absent during the estimated year of change. This event change is a result of removal of forest cover and the percentage of tree cover, derived from Landsat Tree Cover product, below <10% as a deforestation threshold, as opposed to no change when the categories remained the same or tree cover surpassed 10%. Masks for the end of the monitoring period were derived from Landsat Tree Cover product using R.

In order to assess the accuracy of the changes identified by the BFAST algorithm, we constructed a ‘confusion matrix’ (Congalton, 1991) and calculated three measures of accuracy using proportions of areas of change vs. no change rather than pixel counts (Card, 1982) as suggested in Olofsson et al., 2013. Overall accuracy, which is derived from the ratio of the correctly identified forest cover loss area by the total area within the matrix. Producer’s accuracy, the probability of a pixel being correctly classified, was estimated as the proportion of correctly classified area within a category divided by the total area in that category based on reference data. Finally, the User accuracy, the probability that a pixel within a class represents that category, was estimated as the total of correctly classified proportion in a category divided by the total area classified in that same category. The accuracies adjusted by area were calculated according to the formulas provided by Olofsson et al. (2014). Per year, variations in the number of validation sites are due to constraints in the availability of GeoEye 0.5 m resolution imagery for some regions, locations for validation points are shown in Supplementary material S2. Following Olofsson et al. (2013), adjusted areas of forest change and 95% confidence intervals were calculated for each validation year.

Aboveground biomass estimation and comparison with other studies

Chronosequence data for 276 plots in the semi-deciduous forest (Hernandez-Stefanoni et al., 2011) and 86 plots in the semi-evergreen forest (Miranda-Plaza, 2014) were used to estimate AGB of young forest in this study. These chronosequence data consist of plots of different ages located in the study area. The sampling design in both surveys was stratified to cover the whole range of forest ages over a period of 80 and 100 years in semi-deciduous and semi-evergreen forests correspondingly. In both vegetation types, the sampling units consisted of 200 m² circular plots, where all individuals were identified to species level, and height and DBH (diameter at breast height—1.3 m) were recorded. Stand age for every plot was obtained by interviewing local old landowners. More details of sampling and plot designs are in Hernandez-Stefanoni et al. (2011) and Miranda-Plaza (2014). We calculated AGB in both datasets by using the pan-tropical allometric equation of Chave et al. (2005) for trees whose DBH exceeded 10 cm. Meanwhile, for trees with DBH lesser than 10 cm we used two allometric equations, one developed in the semi-deciduous forest near the KK area and the other in the semi-evergreen forest near the FCP site by Ramírez et al. (2017) and Guyot (2011) respectively. The three equations consider wood density values for the different species. These wood density values were obtained from Hernandez-Stefanoni et al. (2020).

Biomass is affected by tree growth, recruitment and mortality processes during succession. However, models to predict AGB changes over time can be obtained from chronosequence data (Chazdon et al., 2006). Therefore, we fitted AGB vs stand age functions using regression analysis. We created two models, one for the semi-deciduous forests, (1) and the other for the semi-evergreen forests (2), the graphs are shown in Supplementary material S3.[Image Omitted. See PDF][Image Omitted. See PDF]

Subsequently, these two models were used to estimate the spatial distribution of AGB of young forest (1 to 20 year old) in the study area using the forest age maps obtained previously. The first model was applied to the deciduous and semi-deciduous forests, while the second model was applied in the semi-evergreen forest. The AGB map was produced by applying the AGB vs age functions obtained previously in R (Supplementary material S4).

Previous studies have indicated that BFAST algorithm is able to distinguish forest cover loss with a maximum difference of 18 months between actual and estimated time of change in tropical dry forests of the Yucatan peninsula (Smith et al., 2019). To account for errors derived from temporal lags between observed and estimated forest cover change, we produced forest AGB maps accounting for a temporal lag of 1.5 years, and calculated the differences between the original map and that with the temporal lag as a measure of error.

Finally, the mapped AGB values of young forest in this study were compared to previous maps of AGB in the same area (Cartus et al., 2014; Hernandez-Stefanoni et al., 2020; Rodriguez-Veiga et al., 2016) and reference field data (24 field plots with AGB < 50 Mg ha⁻¹ from the National Forest Inventory of the young forest in the studied area). We compared the estimated biomass values using mean AGB values and 95% confidence intervals in R (R Core Team, 2015).

Results Magnitude of threshold for estimating forest cover loss

Binomial-logistic regression between sites of confirmed forest loss (deforestation = 1) and forest permanence (no change = 0) sites showed that a pixel with a magnitude value of −0.061 had 50% of probability to belong to an actual deforestation event (Fig. 2).

Figure 2. Binomial-logistic regression between magnitude and confirmed deforestation events. Although some magnitudes will overlap, the majority of negative magnitudes beyond our threshold belong to forest cover loss area derived from confirmed deforestation events.

Spatial distribution of forest age

The spatial distribution of forest age is shown in Figure 3. Our results indicate that the greatest extents of young forests can be found in the semi-deciduous forest site, followed by the deciduous site and much smaller extents of young forests can be found in the semi-evergreen site.

Figure 3. Spatial distribution of forest age grouped by categories of 5 years in a fraction of the three study sites (A) deciduous, (B) semi-deciduous, and (C) semi-evergreen forests.

Temporal patterns of forest cover loss showed generally increasing trends over the study period (2000–2020) in deciduous and semi-deciduous forest, with the largest extent and highest variability corresponding to the semi-deciduous site (Fig. 4). The semi-evergreen site showed a much smaller area of forest cover loss.

Figure 4. Total area of forest cover loss (ha) across time (Year) of change in our study sites. An overall increase in the extent of cover loss areas is observed in the semi-deciduous and deciduous sites compared to the semi-evergreen site.

Validation of forest cover changes detected

Accuracy measures calculated differed between sites and years. Overall accuracy ranged from 95.7 to 99.9 and user’s accuracy ranged from 87.3 to 98.5 (Table 1). Standard errors were low for the overall accuracy ranging from 0.01 to 0.98 as well as user accuracy, ranging from 0.02 to 3.2. However, producer accuracy in the change class was considerably lower in deciduous forests (31.2–45.8), and semi-deciduous forests (57.4–67.2) sites compared to the semi-evergreen forests (66.8–99.9). Producer accuracy showed a wide range of standard errors, ranging from 1.07 to 72.8 in the change class. Our lowest accuracy, in all three metrics, was found in the deciduous site (Table 1). Confusion matrices obtained from validation datasets are included in Supplementary material S4.

Table 1 Overall, user and producer accuracies ± standard errors in percent. Columns represent our map estimation and rows represent reference imagery.

Site	Year	Category	Producer’s accuracy	User’s accuracy	Overall accuracy
Deciduous	2008	Cover Loss	45.9 ± 72.8	87.3 ± 0.02	98.8 ± 0.51
		No change	99.9 ± 0.52	98.9 ± 1.7
	2015	Cover Loss	31.2 ± 24.7	93.9 ± 2.1	95.7 ± 0.98
		No change	99.9 ± 1.04	95.8 ± 1.2
Semi- deciduous	2006	Cover Loss	57.4 ± 29.1	92.6 ± 3.2	99.2 ± 0.39
		No change	99.2 ± 0.40	99.3 ± 1.2
	2015	Cover Loss	67.2 ± 11.6	95.2 ± 1.8	96.9 ± 0.72
		No change	99.7 ± 0.78	97.04 ± 0.95
Semi-evergreen	2010	Cover Loss	66.8 ± 24.4	97.9 ± 1.2	99.5 ± 0.34
		No change	99.9 ± 0.35	99.5 ± 0.71
	2015	Cover Loss	99.9 ± 1.07	98.5 ± 1.1	99.9 ± 0.01
		No change	99.9 ± 0.01	99.9 ± 0.78

The mapped area of forest cover loss, the adjusted area and 95% confidence intervals are reported in Table 2. Our mapped area of forest cover loss was within the confidence intervals of the adjusted area of forest change in all three sites. Greater uncertainty was found in deciduous and semi-deciduous sites (Table 2).

Table 2 Estimated and adjusted area of forest cover loss and its uncertainty in hectares (ha) (with 95% confidence intervals) for each site and validation year.

Site	Year	Mapped area of forest loss (ha)	Adjusted area of forest loss (±95% CI) (ha)
Deciduous	2008	4501.7	3932.5 (±3623.6)
	2015	7397.8	6946.04 (±6938.9)
Semi-deciduous	2006	4004.4	3709.9 (±2804.7)
	2015	20650.3	19653.4 (±5066.6)
Semi-evergreen	2010	3128.9	3062.8 (±2421.6)
	2015	3499.7	3447.1 (72.7)

Aboveground biomass estimated in young forest

AGB estimations for <20-year-old forests showed a mean of 66.1 (±22.8) Mg ha⁻¹ for the deciduous site, 44.0 (±23.2) Mg ha⁻¹ for the semi-deciduous site and 76.1 (±31.3) Mg ha⁻¹ for the semi-evergreen site. A greater extent of young forest cover was found in the semi-deciduous site, followed by the deciduous site, and a smaller extent of young forest was found in the semi-evergreen (Fig. 5).

Figure 5. Spatial distribution of estimated aboveground biomass of young forest in fraction of the 3600 km2 for the three study sites (A) deciduous, (B) semi-deciduous, and (C) semi-evergreen tropical forests.

Considering the possible error derived from temporal lags differences between observed and estimated forest cover change we calculated error maps and frequency histograms that are shown in Supplementary material S6. Errors ranged from zero to 12 Mg ha⁻¹ in our study sites. However, more frequent error are between 6 and 8 Mg ha⁻¹.

Our estimations of AGB of young forest were lower than those obtained by Cartus et al. (2014), Rodriguez-Veiga et al. (2016), and Hernandez-Stefanoni et al. (2020). They also were more similar to field estimates from the National Forest Inventory, although the estimated values of AGB from all studies (including this one) were significantly higher than the field reference data (Fig. 6, Table 3).

Figure 6. Mean aboveground biomass values and 95% confidence intervals of young forests estimated from Hernandez-Stefanoni et al. (2020), Rodriguez-Veiga et al. (2016), Cartus et al. (2014), this study and National Forest Inventory Field data.

Table 3 Descriptive statistics (mean, standard deviation (±SD), 95% confidence intervals, and range) of estimated aboveground biomass in Mg ha ⁻¹ estimated in this study, Cartus et al., 2014, Rodriguez-Veiga et al., 2016, Hernandez–Stefanoni et al., 2020 and field reference data.

Estimation	Mean AGB (±SD)	CI 95%	Range
Hernández-Stefanoni et al. (2020)	80.7 (±33.2)	80.7–80.8	5.5–226.5
Rodriguez-Veiga et al. (2016)	95.2 (±52.7)	95.2–95.3	2.08–197.9
Cartus et al. (2014)	68.1 (±32.8)	68.08–68.1	0.93–136.4
This study	53.06 (±28.5)	53.04–53.08	11.2–98.4
Field reference data	37.5 (±13.8)	26.4–48.6	12.9–48.9

Discussion

The main objective of this research was to obtain spatially explicit estimations of forest age using the year of latest forest cover loss obtained from Landsat NDVI time series. Then, we used this information to estimate AGB of young forest stands in three tropical dry forest sites in the Yucatan peninsula. Overall accuracy for the three stand age (latest forest-cover loss) maps ranged from 95 to 99%. These high values are comparable to those reported in other studies using BFAST to track loss of forest cover or ‘deforestation’ in tropical forests in Guanacaste, Costa Rica, and Yucatan, Mexico (Smith et al., 2019); Bolivia (Dutrieux et al., 2015) and other forest types such as montane forests in Ethiopia (Devries et al., 2015) and in southeastern Australia (Verbesselt et al., 2010). This indicates that a large proportion of the cover change area in the map is located as land change in the field. However, the semi-deciduous and deciduous sites showed low producer accuracies, indicating high omission errors. Previous work has highlighted that higher omission errors in the same dataset are an indicator of underestimation of changes in the class being evaluated (Smith et al., 2019). This error (calculated as 100 – producer’s accuracy) is higher in the deciduous site, with the highest values in 2015 (68.7%) and 2008 (54.1%) in the cover loss class. This suggests an underestimation of the extent of young forests in the study area in general, and in the deciduous forest site in particular.

The total estimated area of forest loss and extent of young forest stands were larger in the semi-deciduous and deciduous forest compared with the semi-evergreen forest (Fig. 3, Table 2). This result may reflect the combined effect of multiple factors both intrinsic to the analysis, such as the accuracy of the algorithm in the different sites and the vegetation index used, and extrinsic factors, such as climatic conditions of the site, proximity to urban settlements, land use and land tenure, which differ across sites. Variability in change detection has a direct impact on the estimation of the area of forest cover loss. The results of this research show that a greater number of change pixels were found in semi-deciduous and deciduous sites compared with the semi-evergreen site. Drier forests such as deciduous and semi-deciduous, tend to have a higher probability of cloud-free observations within the time series, and therefore, a greater probability to find significant change pixels within the time series (Smith et al., 2019).

Adjusted estimates of forest cover loss area also indicate greater uncertainties in the deciduous forest site (Table 2). The amplitude of the confidence interval suggests an underestimation of the forest cover loss area in the deciduous site, although all estimates fell within the estimated 95% confidence interval. Underestimation of the area of forest cover loss in deciduous forests may be partly attributable to the properties of the vegetation index used. The sensitivity of greenness-based indices, such as NDVI, may be affected when rapid regrowth occurs after a clearing as red and near-infrared values saturate quickly (Grogan et al., 2016). Thus, if the timing of the imagery matches the regrowth period after a clearing event, a positive magnitude can be found in an area that was recently deforested. This area would be masked out of the analysis as we considered only negative change detection points as indicators of forest cover loss. This may be more accentuated in deciduous and semi-deciduous forests where swidden agriculture is more widespread than in semi-evergreen forests. Both medium scale and low impact swidden agriculture are practiced within the deciduous and semi-deciduous site, whereas the semi-evergreen site is dominated by selective logging in communally owned areas and pastures for cattle although swidden agriculture is also practiced in the area. Land tenure has played a central role in shaping the patterns of loss of forest cover in central Yucatan where private and federal property showed significantly higher forest loss rates compared to communal property (Ellis et al., 2017). The low impact practices carried out in the semi-evergreen forests have been estimated to reduce disturbance and carbon emissions in the area, outperforming protection schemes carried out in protected areas (Ellis et al., 2019) and thus preserving mature forests within the area. This could help explain the higher producer accuracy estimates and smaller area of forest loss and young stands found in the semi-evergreen forest compared to the other two forest types.

It is important to consider that drivers of forest cover loss potentially differ in extent, frequency, and magnitude among the three forest types assessed in this study. Repeated disturbances over the same area could be more extensive in semi-deciduous and deciduous forests where swidden and subsistence agricultural practices are common. Moreover, during the period of 2005 to 2010, twelve hurricanes impacted the Yucatan Peninsula (Mascorro et al., 2016). Although hurricanes strongly affect forest structure their effect on stand age is not only less clear and difficult to assess, due to the patchy nature of their impact but also, there is a confounding effect of prior land-use history. Many plant species have the capacity of resprout and regenerate (Bonilla-Moheno, 2010; Heartsill et al., 2010; Mcgroddy et al., 2013; Uriarte et al., 2004) thereby affecting the probability of detecting negative changes in the NDVI time series. In addition, forest degradation, which is not explicitly considered in this study, is likely to have an impact on forest age and AGB estimation. In severely degraded forests increased mortality could signify forest cover loss below the 10% cover threshold and could also be considered as deforestation. Therefore, degradation could have an effect on the extension of younger forests and their differences in extent among the forest sites. On the other hand, forest degradation may also impact the AGB contained in the forest and not the forest cover, therefore estimations from passive sensors may overestimate AGB in degraded forest. Accurately characterizing forest degradation could improve mapping of forest age, however, there is much room for improvement in the assessment of the effect of degradation in the loss of forest cover. In particular, active sensors such as L-band ALOS-PALSAR and LiDAR show great potential for the monitoring of tropical forest change as they provide cloud free day observations of terrestrial ecosystems and have been used for the estimation of structural attributes of forests previously (Hernández-Stefanoni et al. 2020; Mitchard et al., 2009). However, ALOS-PALSAR data are unavailable prior to 2006 and the L-band saturates at >150 Mg ha⁻¹ (Mitchard et al., 2009); and LiDAR data are unavailable for large extensions therefore could not be used exclusively to track forest change. Data fusion from active and passive sensors could improve the characterization of forest degradation. Considering the recent efforts in measuring forest degradation with active sensors have proven effective in tropical forests (Joshi et al., 2015), and in particular, with sufficient data availability, the fusion of LiDAR and SAR data could overcome limitations from previous estimations of forest loss in different vegetation types (Collins & Mitchard, 2015).

An important goal of this study was to estimate the spatial distribution of AGB of young forests using the time since last forest cover loss as a proxy of stand age and chronosequence data to relate stand age to AGB. We also compared our AGB estimates to those from previous remote sensing studies such as Hernandez-Stefanoni et al. (2020), Rodriguez-Veiga et al. (2016), Cartus et al. (2014), and reference field data from the National Forest Inventory (CONAFOR). Even though our estimations were significantly higher than field estimated AGB, we obtained lower estimations compared to Hernandez-Stefanoni et al. (2020), Cartus et al., 2014 and Rodriguez-Veiga et al., 2016 in young forest stands. One factor which may have contributed to our overestimation is that we do not have chronosequence data from deciduous forest, instead we used the chronosequence from semideciduous forest to estimate AGB in both deciduous and semi-deciduous forests. Semi-deciduous forest has a relatively taller canopy (10–15 m) than deciduous forest (8–15 m) and also higher AGB. Therefore, lower AGB values from deciduous forest may not be captured by chronosequence data from the semi-deciduous forest. In addition, the chronosequence datasets show a larger spread in the semi-deciduouous site, this may be due to the assignation of forest age which is done through interviews to local informants (Hernandez-Stefanoni et al., 2011; Miranda-Plaza, 2014).

Previous studies aimed at estimating AGB with remote sensing products have found that algorithms such as random forest tend to overestimate the AGB for younger stands (Hernández-Stefanoni et al., 2020; Rodriguez-Veiga et al., 2016; Zhao et al., 2016). A number of factors contribute to this overestimation such as the sensors’ sensitivity to changes in seasonality (Cartus et al., 2014), the algorithm used (Rodriguez-Veiga et al., 2016) and the sensor’s sensitivity to the presence of non-forest areas, canopy gaps, and clearings (Hernández-Stefanoni et al., 2020). The comparison of AGB estimates from this study and previous ones, suggests that including information on the loss of forest cover as a proxy for stand age through time series of reflectance indices such as NDVI as used in this study might decrease this overestimation.

Current work has highlighted the importance of the sequestration potential of tropical dry forests (Chazdon et al., 2016), and particularly young forests (Cook-Patton et al., 2020). However, the comparison with previous AGB mapping efforts from Hernández-Stefanoni et al. (2020), Rodriguez-Veiga et al. (2016) and Cartus et al. (2014) shows that the carbon density in these areas might be currently overestimated. Here we estimate a lower forest biomass (~53.1 Mg ha⁻¹) in comparison to Hernández-Stefanoni et al., (~78.7 Mg ha⁻¹), Rodriguez-Veiga et al., (95.2 Mg ha⁻¹) and Cartus et al., (67.6 Mg ha⁻¹).

The semi-deciduous and deciduous forests of the Yucatan peninsula would contain a lower recovery potential (~44 Mg ha⁻¹) and (~66.9 Mg ha⁻¹) respectively, than the semi-evergreen forest (~76.1 Mg ha⁻¹). In addition, semi-deciduous and deciduous forests have larger areas in a stage of recovery compared to semi-evergreen forests. This would represent the lower carbon sequestration potential of these two vegetation types for the region. Overall, the forests in the Yucatan Peninsula show lower AGB values than other dry forests in Mexico (Corona-Nuñez et al., 2018), South America such as the Guanacaste region (Smith et al., 2019) and central Africa (Xu et al., 2017), however having high deforestation rates, the Yucatan Peninsula is part of the REDD+ early actions program. Therefore, information on the carbon density of biomass of these forests, along with accurate spatial estimations of forest cover loss could be useful for implementing payments of ecosystem services, which could significantly contribute to the decrease in deforestation in the Yucatan, and for obtaining estimations of greenhouse gas emissions from deforestation for this area.

Conclusions

Estimations of the timing and extent of forest cover loss are important as they provide information on the spatial distribution of stand age. This is particularly challenging in seasonal forests, such as the tropical dry forests of the Yucatan, where the detection of forest loss from remote sensors can be confounded by seasonal and inter-annual variation in leaf phenology. BFAST provided estimations of forest loss from significant changes in NDVI time series in tropical forests of the Yucatan Peninsula allowing to map forest cover loss with high overall and user’s accuracy. Although producer’s accuracy was considerably lower in semi-deciduous and deciduous forests indicating an underestimation of the forest loss area in these forest types, adjusted area of forest loss indicates that the area of change found in all three sites lies within the confidence intervals.

Information on the timing and extent of forest cover loss can be used in conjunction with chronosequence data to estimate AGB and identify areas of young recovering forests where potential carbon sequestration is greater in tropical forests, and thus aid conservation efforts to protect the ecosystem services of young forests in tropical areas. Moreover, spatially explicit estimations of forest age, or the time since the latest forest cover loss event obtained from time series of satellite imagery show potential in reducing uncertainties associated with the overestimation of the AGB in young tropical forests from remotely sensed data.

Acknowledgments

The authors thank Erik Lindquist, and the USFS Silvacarbon for their guidance in the implementation of the BFAST algorithm. We also thank Fernando Tun-Dzul and Filogonio May Pat for their support during field work. The first author would like to thank CONACyT for the Ph D. scholarship awarded.

Funding Information

This research was financed by Ecometrica LTD and the United Kingdom Space Agency as part of the project Forests 2020.