Landsat sub-pixel land cover dynamics in the Brazilian Amazon

The Brazilian Amazon land cover changes rapidly due to anthropogenic and climate drivers. Deforestation and forest disturbances associated with logging and ﬁres, combined with extreme droughts, warmer air, and surface temperatures, have led to high tree mortality and harmful net carbon emissions in this region. Regional attempts to characterize land cover dynamics in this region focused on one or two anthropogenic drivers (i.e., deforestation and forest degradation). Land cover studies have also used a limited temporal scale (i.e., 10–15 years), focusing mainly on global and country-scale forest change. In this study, we propose a novel approach to characterize and measure land cover dynamics in the Amazon biome. First, we deﬁned 10 fundamental land cover classes: forest, ﬂooded forest, shrubland, natural grassland, pastureland, cropland, outcrop, bare and impervious, wetland, and water. Second, we mapped the land cover based on the compositional abundance of Landsat sub-pixel information that makes up these land cover classes: green vegetation (GV), non-photosynthetic vegetation, soil, and shade. Third, we processed all Landsat scenes with < 50% cloud cover. Then, we applied a step-wise random forest machine learning algorithm and empirical decision rules to classify intra-annual and annual land cover classes between 1985 and 2022. Finally, we estimated the yearly land cover changes in forested and non-forested ecosystems and characterized the major change drivers. In 2022, forest covered 78.6% (331.9 Mha) of the Amazon biome, with 1.4% of secondary regrowth in more than 5 years. Total herbaceous covered 15.6% of the area, with the majority of pastureland (13.5%) and the remaining natural grassland. Water was the third largest land cover class with 2.4%, followed by cropland (1.2%) and shrubland (0.4%), with 89% overall accuracy. Most of the forest changes were driven by pasture and cropland conversion, and there are signs that climate change is the primary driver of the loss of aquatic ecosystems. Existing carbon emission models disregard the types of land cover changes presented in the studies. The twenty ﬁrst century requires a more encompassing and integrated approach to monitoring anthropogenic and climate changes in the Amazon biome for better mitigation, adaptation, and conservation policies.


Introduction
The primary interface with the atmosphere is the land cover, composed of plants, topsoil, rock, water layers, and built-up infrastructure.It is also where people and a wide variety of species live.Mapping and monitoring land cover dynamics are essential to characterizing biological and geophysical attributes for hydrological and geochemical modeling of energy balance and carbon exchange with soil and atmosphere (DeFries et al., 1999;Claussen et al., 2001;Roberts et al., 2003;Jones, 2008;De Sousa-Neto et al., 2017).It is also essential for sustaining biodiversity and the growing global socio-economic demands (Krausmann et al., 2013).
While several global land cover mappings (GLCMs) exist, these maps can rarely be directly applied to regional studies of land cover dynamics, given their high level of uncertainty (Congalton et al., 2014;Pendrill and Martin Persson, 2017;Galiatsatos et al., 2020), and thematic classification infrequently satisfies regional users (Grekousis et al., 2015).The existing GLCMs do not provide land cover change estimates such as deforestation and vegetation regrowth to support policy and decision-making for conservation and planning.Anthropogenic and climatic change are the main factors influencing the land cover dynamic in the Amazon biome (Marengo et al., 2018), and the GLCMs do not provide adequate accuracy and consistency to assess these factors (Nemani and Running, 1995).
Satellite imagery is vital to mapping and characterizing land cover dynamics.The 50-year Landsat data archive (Wulder et al., 2022) and the new cloud computing capability, notably by Earth Engine (Gorelick et al., 2017), have driven the proliferation of medium spatial resolution (i.e., 30 m pixel size) land cover maps at global, continental, and regional scales.The most common approach to mapping land cover using Earth Engine is based on the pixel spectral characterization of the land attributes from the spectral bands and spectral and textural indices with machine learning classification algorithms, mainly random forest, according to a systematic literature review (Pérez-Cutillas et al., 2023).Objectbased classification has also brought attention to remote sensing map makers (Ma et al., 2017).With the novelty of cloud computing, time-series analysis (Friedl et al., 2022), and deep learning (Brown et al., 2022), image classification algorithms for handling large data sets are emerging.The MapBiomas initiative focuses on mapping at the biome scale to deal with their intrinsic biophysical and climatological characteristics.MapBiomas uses annual synthetic spectral bands and indices to obtain less cloudy (i.e., < 5.0%) Landsat mosaics to map annual land cover and land use (Souza et al., 2020), an approach replicated in several countries.
Several studies have also demonstrated that the compositional information within the Landsat pixel is valuable in characterizing land cover classes (Roberts et al., 1993;Adams et al., 1995;Asner, 2009).The pixel composition can be estimated from the multispectral bands of Landsat sensors in terms of the abundance of green vegetation (GV), non-photosynthetic vegetation (NPV), soil, and shade, which make up the land cover classes in the Amazon biome (Souza et al., 2005).For example, grasslands are composed of a mixture of GV, NPV, and soil uniquely different from dense forests (with high shade and GV and no soil) (Asner and Heidebrecht, 2002).Sub-pixel information is also essential to detect and characterize disturbances leading to ecosystem degradation (Veraverbeke and Hook, 2013;Salih et al., 2017).Therefore, pixel compositional information can potentially separate ambiguous land cover classes (e.g., pasture, bare soil, rock, and impervious surface), improving land cover mapping and characterizing ongoing climate change in land cover interactions (Song et al., 2018).
Most studies in the Amazon biome have focused on deforestation and drivers of forest degradation to characterize land cover dynamics (e.g., Souza et al., 2013;Bullock et al., 2020b;Potapov et al., 2020) or cover a few Landsat scenes and a short temporal scale (Roberts et al., 2002).Another essential demand is using land cover maps to estimate changes, such as deforestation, vegetation regrowth (Nunes et al., 2020), and water dynamics (Souza et al., 2019), directly from the land cover map series.
The objectives of this study are 4-fold: i) improve the detection and mapping of the main land cover classes of the Amazon biome using Landsat sub-pixel compositional information; (ii) characterize the annual land cover change in forest, non-forested vegetation, and aquatic ecosystems; (iii) assess how the land cover change drivers varied between 1985 and 2022; and (iv) estimate annual loss and gain among the land cover types during this period.This study is organized as follows: Section 2 presents the methods applied to process the Landsat time series to map and characterize land change.Section 3 presents the results, focusing on the annual land cover change and land change drivers.Then, we discuss the mapping uncertainty and summarize the findings of future research steps.

Methods . Study area and landsat coverage
The Amazon is the world's largest tropical rainforest, spanning the borders of eight South American nations.Its environmental functions are crucial for preserving the equilibrium of the world's climate, its vast biodiversity, and the homeplace of the original Amazonian people.Our study area, the Brazilian Amazon Biome (Figure 1), comprises more than 421 million hectares, or 50% of the Pan-Amazon extension.The Amazon biome occupies 49.5% of the nation compared to Brazilian land.
We selected all Landsat images from Collection 1 Tier 1 in the Google Earth Engine (Gorelick et al., 2017) catalog with <50% cloud cover.More than 86k Landsat scenes were used to build the Amazon annual mapping between 1985 and 2022.We assessed the number of Landsat scenes available per year and identified three distinct periods of imagery availability (Figure 2).The first was from 1985 to 1998 when the Landsat 5 TM sensor was the only operational.In that period, we found almost 15k available images, with an average of 1,070 per year.The second period was between 1999 and 2011, when Landsat 5 and 7 operated simultaneously, duplicating the total number of scenes compared with the previous period.In the second period, we selected 34k Landsat images (an average of 2,638 per year).In 2012, we observed a drop in image availability because Landsat 5 stopped acquiring data, with only Landsat 7 operational, reducing the number of scenes to 1,753 images this year.
The last period of Landsat scene availability was between 2013 and 2021.We had Landsat 7, 8, and 9 available in this period, raising the number of images per year relative to the previous periods, reaching an average of 3,358 Landsat scenes per year.For 2021, we used the Landsat 9 Collection 2 Tier 1 (n = 67 scenes) because Collection 1 data was unavailable.Including Landsat Collection 2 in the Landsat time series does not affect the annual classification results because the scenes are processed individually.In 2022, we reached almost 5k Landsat scenes (Figure 2), which may inaugurate a new period of image availability for Landsat 9.  We also assessed how the number of Landsat scenes varies across the study area (Figures 1, 2).The southern region showed 400-600 Landsat scenes over the entire period (i.e., 1985-2021; Figure 1), whereas the northern Atlantic region had fewer scenes available (<100).The central area had approximately 300 scenes, and the western and central northern areas had approximately 200. .

Land cover classes
Defining a classification scheme remains challenging for remote sensing of terrestrial land cover mapping (Yang et al., 2017).Land use is related to how people interact with land surface types (i.e., classes).In contrast, land cover relates to the properties (biophysical and chemical) of the Earth's surface (Lambin et al., 2001).We selected the most relevant land cover classes in our mapping efforts to measure change in the Brazilian Amazon, including forest, flooded forest, shrubland, natural grassland, pastureland, cropland, outcrop, bareland and impervious, wetland, and water (Table 1).We decided to include pastureland and cropland, which are associated with land use, because their extent is significant in the Brazilian Amazon, and both have distinct spectral, compositional, and phenological characteristics from natural grassland and are the main drivers of change in the region.Detailed descriptions of the land cover classes are presented in Table 1.

. Image processing
Our approach to mapping and characterizing land cover dynamics in the Brazilian Amazon Biome with Landsat time series follows five steps: (i) data preprocessing; (ii) spectral mixture analysis (SMA); (iii) land cover classification; (iv) postclassification; and (v) accuracy assessment (Figure 3).We present below the detailed methodology for each of these steps.

i) Data preprocessing
We used Landsat Collection 1 Tier preprocessed and ready for analysis (Dwyer et al., 2018), which fulfills standard geometric and radiometric quality criteria, supporting consistent time series applications (Qiu et al., 2018).Detailed information about Landsat scene selection and availability is presented in Section 2.1.Once the scenes were selected, the next step was to apply cloud and shadow masks.We used the Temporal Dark Outlier Mask (TDOM) algorithm and the Band Quality Assessment (BQA) band available in the Landsat Collection for that purpose, presented in detail in the code script (Supplemental material 1).
We also used annual Landsat cloud-free mosaics, obtained by applying a median reducer to all scenes selected every year from 1985 to 2022, as inputs to classify outcrops, bareland, and impervious classes.The Landsat quality assessment (QA) band was used to remove the shadow and cloud from each Landsat scene used to build the annual mosaics.Furthermore, we calculated the median value of the stacked bands for each year, removing the remaining cloud and shade not filtered by the QA band.We also used ancillary data to assist in mapping the outcrop class.The Global 30 meters height above the nearest drainage (HAND) (Rennó et al., 2008) dataset was used for that, combined with the sub-pixel compositional information generated with spectral mixture analysis (see Section 2.5.1 for more detail).HAND offers a worldwide depiction of elevation relative to the closest drainage point at a resolution of 30 meters, which helps detect regions with steep slopes and ridges associated with outcrops.The SRTM data are also an input to map the outcrop class.ii) Spectral mixture analysis We estimated the pure spectral composition's subpixel proportion (or fraction), named endmember of green vegetation (GV), soil, non-photosynthetic vegetation (NPV), and cloud in each pixel using a generic SMA model (Souza et al., 2005).The SMA assumes that the reflectance of each pixel in each band of the Landsat image can be modeled as a linear combination of the fraction image Fi of the endmember components and its additional residual: where ρ b is a reflectance of the n th mixed, in band b, i is the reflectance of the pure component i, in band b, F i the fraction (or proportion) of the endmember component i, and ε b the residual error of each band.The fractional values of the endmember in the pixel i sum up to one (i.e., 100%) so that, The SMA is solved to estimate the unknown fraction values F i of each endmember component within the Landsat image pixel, thus obtaining sub-pixel proportions of vegetation, NPV, soil, and cloud.The shade fraction is obtained as a complement of the sum of all fractions for each pixel.In this case, the shade endmember is set to zero (i.e., photometric shade) in all Landsat bands.The endmember values used in this study have been published elsewhere (Souza et al., 2005).Fractional information of each was also used to calculate additional SMA indices, including the normalized difference fraction index (NDFI) (Souza et al., 2005) and the canopy shade factor index (CSFI), a novel fractional index presented in this study.These sub-pixel compositional fractional and derived indices had been used in time series analysis (Bullock et al., 2020a), and in this study, we used them to classify the land cover classes targeted in this study.Detailed information to calculate the SMA fractions and indices using Google Earth Engine is presented in Supplementary material 1. iii) Land cover classification The land cover classification steps included: (i) defining the reference dataset to train, calibrate, and test the map results; (ii) applying data augmentation to train and calibrate RF models for each Landsat scene; (iii) building the RF classification model and applying it to each Landsat scene to generate annual land cover maps; (iv) applying postclassification algorithms to remove classification errors and artifacts; and (v) assessing the map accuracy.
We used the point reference data generated by the MapBiomas Project (Parente et al., 2021;Sales et al., 2021), .
/ gc. .which contains tens of thousands of sample pixels.The reference dataset was built annually from 1985 to 2021 with colocated land cover classifications over Brazil.Sample pixels were randomly selected using topography and land change proportion strata.Several trained interpreters carefully assigned land cover and land use classes to Landsat pixels through visual interpretation, spectral information from Landsat and Modis, and higher spatial resolution imagery (when available).The reference dataset classification was harmonized with the classification legend of this study (Table 1).We used all 35,000 pixels that fell inside the borders of the Amazon Biome and split them into training and calibration (10k samples; 29%) and test datasets (25k samples; 71%).Because our land cover classification approach is applied to all Landsat scenes, we had to implement a data augmentation technique to get enough samples per class to train and calibrate the RF algorithm.First, we segment the Landsat scene using the simple non-iterative clustering (SNIC) algorithm, using six Landsat bands (red, green, blue, NIR, SWIR-1, and SWIR-2).Then, we overlay the point reference data sample with the segmented Landsat scene matching the year.The image segments with more than one land cover class were excluded, leaving only segments with a single class.We assessed the level of class purity within the image segment using existing land cover maps (i.e., MapBiomas Collection 7).Finally, we generate random samples within the selected image segment (i.e., polygon region) and assign the same thematic classification to the reference data sample used to select the image segment.This process guarantees enough samples for the RF training and calibration for the majority of the thematic classes (Table 1).This approach was applied in the +86k Landsat scenes selected for the Amazon biome through the time series.
We used the random forest (RF) classifier available in Google Earth Engine (Gorelick et al., 2017) using as input data the SMA sub-pixel fractional information, which included: green vegetation (GV) soil, non-photosynthetic vegetation (NPV), cloud, shade, NDFI, and CSFI.This selection was based on the feature importance analysis, which pointed out that the fraction bands were in the top 10 list of the most important features to classify the Landsat scenes together with the NIR and SWIR bands.The number of input features selected was the default value for this parameter (i.e., mtree, which is given by the square root of the number of features as defined for each biome).The RF algorithm was set to run 50 to 100 iterations, which allows for the lowest calibration error.
The land cover classes wetland, outcrop, bareland, and impervious did not have enough reference sample data (i.e., classified pixels) for training and calibrating the RF classifier.To overcome that, we classified these classes independently and integrated them into the post-classification step (Figure 3).Wetlands can be found in three land cover types: forest, shrubland, and natural grassland (Table 1).These classes can be permanently or sporadically flooded.We randomly selected stratified samples using the MapBiomas Surface Water dataset to train and calibrate the RF model to classify wetland (Supplementary material 2).Then, we used SRTM, HAND, and SMA fractions as inputs to the RF to classify wetland.We combined the wetland binary map with the forest class to obtain the flooded forest (i.e., overlapping with wetland).The areas where shrubland and natural grassland overlapped with wetland were reclassified to the latter class.
Outcrop has a unique spectro-temporal behavior, with a stable compositional trend of soil, NPV, and GV.We empirically defined the outcrop compositional trend and used empirical rules to detect and map this class.Additionally, we included information about altitude, steep slopes, and geomorphological features such as escarpments and the top of hills using SRTM and HAND.Then, we selected stratified samples (with outcrop and no outcrop) to train and calibrate an RF model to map outcrop.

FIGURE
Image processing steps and algorithms applied to map land cover dynamics in the Brazilian Amazon biome at the Landsat sub-pixel level.
Similarly, we identified the compositional trend of bareland and impervious surfaces using SMA fraction images and derived indices.We used the dynamic world (Brown et al., 2022) dataset available in the Google Earth Engine to randomly select stratified samples for the classes impervious/bareland and others.Then, RF was used to generate annual binary maps of impervious/bareland to update the annual land cover maps.Detailed information on how these classes were mapped can be found in Supplementary material 3. iv) Post-classification In this step, we generated the final annual land cover and land cover change maps using the Landsat scenes classified individually.For each year, we calculated intra-annual map metrics to measure classification stability, class changes, and classification errors, as proposed by Pontius et al. (2017).
Then, we trained and calibrated RF with the samples used to classify each Landsat scene to generate the annual land cover maps (Figure 4).More details about the annual classification are provided in Supplementary material 1.
Finally, we applied spatial and temporal filters to generate the final annual land cover map collection, which we used to measure land cover change (Figure 3).The spatial filter uses a majority-neighborhood logic to replace single pixels or a group of pixels.The mode of the neighborhood classes replaces areas smaller than 1 hectare.This approach helps to eliminate undesirable transitions deriving from classification.The temporal filter is a set of rules for non-allowed transitions applied to each image classified in a given year.For example, it was possible to remove clouds and correct non-allowed transitions.A total of 50 temporal rules, distributed in 3 groups, were used: (a) rules for cases not observed in the first We evaluated the mapping accuracy and area estimation using the good practices protocol (Olofsson et al., 2014) based on 25,000 independent test samples obtained with reference data splitting (Figure 3).Then, we estimated each class's overall accuracy and the user's and producer's accuracy with 95% confidence.We then assessed each year's map accuracy using the reference samples.Because not all samples could be utilized annually due to data cloud coverage or the inability to classify a sample pixel, the sample size varied annually.

Land cover change between and
The Landsat annual map time series generated in this study allowed us to assess land cover change and the dynamics of the Brazilian Amazon biome.We present first the changes that happened between the initial date (1985) of the Landsat time series analyzed and the end date (2022) (Figure 5).In 1985, the Amazon Biome comprised 97% of natural land cover classes (forest, shrubland, natural grassland, and water).We expect some level of anthropic disturbance in the natural land cover during this period.For example, forests are affected by logging, fire, and fragmentation; fires in shrubland and natural grassland and damming water bodies affected these land cover classes by 1985.Pasture covered only 2.9% (i.e., 12.10 Mha) of the biome with a tiny relative area of cropland (400 thousand ha) (Table 2).
Astonishingly, over the 37 years of our Landsat time series analysis, we found extensive land cover change, affecting primarily the forest class.The most extensive land cover change was the conversion of forest to pastureland.We estimated that 46.6 Mha (11% of the forest in 1985) of forests were lost over the 37 years analyzed, which represents an average loss of 1.254 Mha/year (Table 2).Pastureland increased by 44.7 Mha (10.6% increase relative to 1985), which means that forest conversion to pastureland dominates in the Amazon (95.7%).The comparison of the initial and end maps of the land cover time series also shows that forest loss no longer concentrates in the so-called arc of deforestation region, which is in the southern fringe between the Amazonia and Cerrado biomes, going up north along the eastern border.In 2022, we found three new hot spots of forest loss: in the midlands of Pará, along the BR-163 road connecting Santarém-PA to Cuiabá-MT, and in the southern Amazonas (Figure 5).

. Land cover trend analysis
This section explores the yearly land cover change over the entire time series between 1985 and 2022.We investigated the trend (Figure 6) and breaking points (Figure 7) in the time series that reveal major shifts in change among the land cover classes.We applied the Loess function to the annual land cover maps to estimate the changing trend (Figure 6).
The forest class steadily decreased between 1985 and 2022, from 362.0 Mha to 315.6 Mha, which reduced the forest cover in the Brazilian Amazon biome approximately to 75% (Table 2).When inspecting the annual forest cover change, we observed four  temporal breaking points associated with forest loss (Figure 7).The first breaking point period, between 1985 and 2004, showed the fasted rate of deforestation (1.62 Mha/year).The second period entered the deforestation control phase between 2004 and 2012, when forest loss reached a low annual rate in 2012, with an average annual rate of 0.94 Mha/year.The year 2013 marked the start of a new period of increasing deforestation that will extend until 2022, the end of our time series analysis.We also found that some years showed forest gain due to secondary regrowth [not captured in Brazil's official forest monitoring system (Valeriano et al., 2012)].From 2013 to 2022, the average forest loss was 1.04 Mha/year (Figure 6).Natural grassland and shrubland land cover classes showed a monotonous temporal trend, with some years of decrease in extent, possibly associated with droughts and fires (i.e., in 2000 and 2010).We identified two possible decadal change cycles, marked in 2004, in these land cover classes, which deserve further in-depth analysis with climate and other variables (not the subject of this study).The land cover classes associated with water land cover classes (i.e., water, wetlands, and flooded forests) showed a decreasing trend in annual extent (Figure 6).Surface water started reducing in 1998, with several years below the average annual change, notoriously 2004-2005, 2010, 2016, and 2020, possibly associated with extreme droughts (Marengo et al., 2023) (Figures 6, 7).The annual change in the flooded forest class showed a similar temporal trend as in water.However, wetlands exhibited a steeper loss after 2010, a hallmark year of consecutive dry years in the Amazon region (Parrens et al., 2023).As expected, the outcrop class showed no significant change over the 37-year-time series due to its longterm geologic stability.Our outcrop mapping approach minimized the vegetation phenological changes expected in the land cover class because we used dry season imagery mosaics to map it.
The annual time series trend analysis also revealed temporal breaking points associated with pastureland, cropland, and bareland and impervious land cover classes (Figures 6, 7).A detailed inspection of the annual area change behavior showed three temporal patterns of pastureland expansion (Figure 7).First, between 1985 and 2000, the pastureland area expanded to 22.66 Mha at an annual average rate of 1.51 Mha/year; then decelerated to 0.99 Mha/year between 2000 and 2013, and slightly increased until 2022 at 1.0 Mha/year.Cropland expansion was marginal between 1985 and 1996 at an annual average of 2,550 ha/year (Figure 6).Cropland annual area changes exhibited three distinct temporal rate patterns (Figure 7).The trend Loess function revealed that 1996 was the temporal breaking point for cropland faster expansion in the Amazon Biome, with two periods: 1996-2008 and 2008-2022.In the first period, the cropland annual average expansion rate was 158,500/year, and in the second one, 233,000 ha/year (Figures 6, 7).
Bareland and impervious are a mixed land cover class, including urban areas and infrastructure (roads and large mining infrastructure) for the impervious surfaces, and abandoned pastures with unvegetated topsoil (i.e., bareland), also found in most rural towns.We also found three breaking point periods for this land cover class: 1985-2000, with an expansion of 21,682 ha at a rate of 1,457 ha/year; 2000-2013, expanding almost 95,000 ha (at 7,271 ha/year); and 2013-2022 gaining 143,929 ha (at 15,992 ha/year) (Figures 6, 7).

. Land cover change drivers
We used the annual land cover trend analysis presented above to estimate the annual area change in each mapped class and to define temporal breaking points that indicated how fast or slow the changes happened.In the section, we used the breaking points to understand the land cover class transitions in forest and aquatic ecosystems (Figure 8).
The annual land cover change analysis identified three key breaking points along the 1985-2022 time series associated with forest loss : 1998, 2004, and 2013.These breaking points are temporal markers of deforestation, mainly due to pastureland expansion (Figure 8A).We also detected, to a smaller extent, pastureland-to-cropland conversion and pasture abandonment (followed by potentially secondary growth).We also estimated that 2 million hectares of pastures were abandoned and potentially entered into secondary vegetation growth, and that 27 thousand hectares were shifted to cropland (Figure 8A).Forest directly converted to cropland was more than four times larger (i.e., 127 thousand hectares) relative to pastureland conversion to this class.A tiny area (600 hectares) of cropland was abandoned and followed secondary growth.Interestingly, 28.7 thousand hectares of cropland returned to pastureland in this period (6.2% more than the opposite).Only 10.1 thousand hectares remained as cropland in the following period, i.e., 35% (Figure 8A).Between 1998 and 2004, it was also a period with very high annual deforestation rates, with 2004 considered the second largest annual rate of forest loss in the region, according to the official deforestation statistics in Brazil (Valeriano et al., 2012).The annual rate of pastureland abandonment also increased by 15% (2.3 Mha), and pastureland converted to cropland increased almost 20-fold (532 thousand hectares) relative to the previous one, which marks the acceleration of crop expansion in the Brazilian Amazon.We detected approximately five times as much direct forest loss to cropland (i.e., 536 thousand hectares) relative to the previous period.Cropland shift to secondary forest regrowth was 1.18 thousand hectares (Figure 8A).
The 2004-2012 period was recognized as the 'golden' period of deforestation control in the Brazilian Amazon (West and Fearnside, 2021).However, our analysis showed that the annual deforestation rate in the first 3 years of this period was not different from 1985 to 1998.Lowering deforestation started after 2008 and lasted until 2012.Cropland expansion over pastureland reached 1.0 Mha, and 153 thousand hectares were the opposite.This type of land use shift occurs in farms as a soil management practice in crop-livestock systems (Monteiro et al., 2024).The Pastureland and Cropland convertions to secondary growth forest were 2.4 Mha and 590 thousand hectares, respectively (Figure 8A).
The 2013-2022 period is notable for a lack of deforestation control in the Brazilian Amazon biome (Gatti et al., 2023).Annual forest loss rates to pastureland and cropland were 1.44 Mha/year and 44.2 thousand hectares/year (Figure 8A).The total direct forest to cropland conversion area summed 385 thousand hectares, and pastureland shifted to 3.3 Mha, establishing extensible crop fields in southern Mato Grosso (Figures 5, 8A).
The aquatic ecosystems showed two prominent temporal breaking points: 1998 and 2010 (Figure 7).Surface water was reduced by 10 thousand ha between 1985 and 2022, with distinct dynamic patterns in the breaking point periods.In the 1985-1998 period, water changed among the water-influenced land cover classes (i.e., water, flooded forest, and wetland).We found that 3.36 thousand ha of water class changed to wetland, indicating that surface water was reduced to the point that soil and vegetation mixed with water (Figure 8B).The change from water to natural grassland was 1.2 thousand ha and 0.74 thousand ha to flooded forest, indicating surface water loss.Conversely, flooded forest and wetland were covered by water in this period, with a total change of 4,62 thousand ha (96% to wetland).This may indicate that new areas were permanently flooded.Overall, ∼4.22 thousand ha of surface water (i.e., water class) mapped in 1985 was displaced,

FIGURE
Land cover class transitions in the Amazon Biome at the temporal breaking points associated with forest conversion (A) and in the water ecosystem (B).
and new surface water was gained, extending approximately 100 thousand ha.From 1998 to 2010, the change in water to wetland almost doubled to 6.42 thousand ha, relative to the previous period (Figure 8B).A similar trend was observed from water to natural grassland, with a change of 2.70 thousand ha.
We also mapped the spatial distribution of the land cover change along the temporal breaking points (Figure 9).The forestto-pastureland conversion concentrates on the arc of deforestation between 1985 and 2004, along the fringe of the Amazon biome with the Cerrado and Caatinga biomes in the southern and eastern borders, respectively.After that, pastureland replaced forest, mostly toward the core of the Amazon (Figure 9A).Cropland expansion was concentrated in the southern Mato Grosso states, with less extent in the eastern Amazon.Secondary growth vegetation occurred along the regions where pastureland is concentrated, with more concentration in the northeast of the Amazon biome.The water loss occurred along the Amazon River margins in its estuary on the Marajó-PA Island (Figure 9D).The varzea region on the northern Amazon border with Roraima also showed large areas of water loss, as did the southern border of the Amazon biome with Pantanal and in the Eastern Amazon.Water gains are mostly associated with damming for hydroelectric power, with less extent flooding in Marajó-PA Island in the estuary of the Amazon River.
The quantification, characterization, and mapping of land cover change indicated that pastureland expansion is the major change driver in the Amazon biome, followed by cropland expansion, and that water losses, especially after the 2010 severe drought followed by more frequent similar events (Marengo et al., 2018), suggest that climate change is a newer driver affecting water ecosystems in the region.

Land cover mapping uncertainty
We estimated annually the mapping accuracy between 1985 and 2021 (there was no reference data for 2022) using the independent reference test sampled (25k).The mean overall accuracy during these years was 89% ± 0.0023.In this study, we present the users' accuracy per class, which is more meaningful because it measures how each land cover class mapped agrees with the reference data (Stehman and Czaplewski, 1998) (Supplementary material 5 presents all accuracy metrics per year).
The reference sample we chose for this study was designed to train a different land cover classification scheme.We harmonized the reference samples to the Amazonia biome, which is further integrated into the final MapBiomas mapping product.The reference test sample estimated the overall accuracy with a margin of error of 5%.It is crucial to highlight that the sample size for reference data in several classes is small, as described above, requiring another sampling design, even though we decided to present the area error estimate using good practices (Olofsson et al., 2014) (Figure 10).
Comparing the land cover area estimates with pixel counting and areas adjusted using the accuracy error matrix (Figure 10) reveals that our mapping efforts to reconstruct the 37-time series presented in this study can identify the major land change drivers in the Amazon biome.The forest class showed similar annual area estimates using these two approaches, with pixel counting accounting for 1.5% less area on average.Our mapping approach underestimates the area of Pastureland by 14% (SD = 8%), with more difference between pixel counting and area adjusted in 1985-1994 (26%), reducing to 10.5% in 1995-2018, and converging (2.7%) in the last 3 years of the time series (2019)(2020)(2021)(2022).Pixel counting overestimated the cropland area along the time series relative to the area adjusted, except for 1985 and 1986, but showed a similar area expansion until 2012 (27% larger on average in 1987-2012).After 2013, these area estimates showed different trends, with the area-adjusted estimate below 400 thousand ha (Figure 10).
Our mapping approach also underestimated the area of flooded forest by 28% on average (SD = 4%) and the wetland class by 42% (SD = 3%).Natural grassland was under-mapped on average by 5% until 2011, increasing to 17% after 2012.Shrubland was overmapped with pixel counting along the entire time series at 13% (SD = 6%) but within the confidence interval (CI = 95%) in most of the years.Even though we estimated that the area adjusted to bareland, impervious, and outcrop, the results must be disregarded due to the limited number of samples needed to estimate the mapping accuracy (Figure 10).
The low mapping accuracy of the classes above did not affect the land cover change analysis because most of the changes occurred within the most predominant classes, and these classes showed higher accuracy (i.e., forest, pastureland, cropland, and water).Additionally, the area estimates of these main land cover classes showed similar results using pixel counting, and the area was adjusted with the accuracy metrics according to "Good Practices" (Olofsson et al., 2014), allowing us to infer the land cover change drivers in the Amazon Biome (presented in Section 3.4).The remaining cloud and shade after applying the TDOM algorithm did not affect the land cover mapping results because the SMA model handled them.However, we could not assess possible remaining atmospheric contamination (e.g., haze) in the Collection 1 Tier preprocessed Landsat data and its impact on the mapping uncertainty because we processed more than 86k Landsat scenes, and assessing individual scenes was unfeasible.

Deforestation estimation
Given the high accuracy of mapping the forest class (93-94%), we designed a methodology to estimate annual deforestation directly from the land cover time series.We also applied a methodology to map the annual extent of secondary growth based on the methodology proposed by Nunes et al. (2020).We combined the annual forest losses and gains to estimate the annual forest balance and the annual deforestation of old-growth forests and secondary forests ≥ 5 years old (Supplementary material 6).
The total deforestation estimate in the Amazon biome was 52.7 Mha between 1985 and 2022 (Figure 11), and the total forest loss in this period was 46.4 Mha, less than the total deforestation because the forest class includes secondary forests.Subtracting the secondary growth forest (> 5 years old, because 1-5 years old is considered fallow land) gain of 5.7 Mha in 2022, we get 47 Mha of forest loss in 2022, meaning that the Amazon Biome lost 52.7 Mha over the 37 years analyzed, and 5.7 Mha is under secondary vegetation growth.The annual deforestation method used masked out forest loss annually, and to unmask it, we had to apply another method (Nunes et al., 2020) (Supplementary material 6).Therefore, we present the annual deforestation of old-growth forests and secondary vegetation (> 5 years old).
We compared the annual deforestation rate for old-growth forests with other sources (Supplementary material 6, Figure 2).Our total deforestation estimate was 9.3% greater than the Prodes deforestation estimate (48.2 Mha) for the Amazon region from 1988 to 2022 (period overlapped with Prodes; Figure 2, Supplementary Material 6).We observed more annual differences in our estimate with the Prodes one between 1991 and 2000 and less difference between 2000 and 2012.After 2012, our deforestation estimate was 18.9% larger than Prodes.We attribute the overall differences to the spatial resolution of our maps (i.e., 30m sub-pixel) compared to Prodes, which used 6.5 ha.Our mapping approach used all useful Landsat (>50% of cloud cover) scenes, whereas Prodes uses only one scene per year.
Another innovation brought by our land cover change analysis was to estimate the deforestation of secondary growth forests (Figure 11b).The tree loss cover product (Hansen et al., 2013), which considers both old-growth and secondary forests in the deforestation estimation analysis, does not separate these two deforestation classes.The annual deforestation from the tree loss cover product showed similar results to our estimates until 2013, increasing over most of the year until 2021.The deforestation in secondary forests > 5 years old adds up to

Conclusion
In this study, we presented the land cover dynamics of the Brazilian Amazon biome between 1985 and 2022 and the associated mapping uncertainties.The novelty of our mapping approach was to process more than 86 thousand Landsat scenes to extract subpixel compositional information as input to machine learning and empirical rules to map the 10 most significant land cover classes in the biome.The annual land cover maps revealed that forest loss is mainly associated with pastureland for cattle ranching, followed by cropland expansion.It was also possible to derive annual deforestation and secondary growth vegetation maps and statistics, with estimates compatible with Brazil's official deforestation and global deforestation products.We also characterized non-forested ecosystems, including shrubland and natural grassland classes, and assessed their anthropogenic change.It was also possible to identify water loss in aquatic ecosystems, which may be associated with extreme and frequent droughts, extended summer seasons, and higher temperatures.Future steps for our research will focus on the measure of forest disturbances related to forest degradation drivers, which include selective logging, fires, and forest fragmentation, based on this study's land cover mapping efforts.The 37-year improved spatial-temporal characterization of land cover change and drivers in the Amazon, revealed in this study, offers crucial details about the threats to the Amazon biome and how human intervention and, possibly, climate change have changed the region.

FIGURE
FIGUREBrazilian Amazon and Landsat scenes of Collection during and with < % of cloud cover.

FIGURE
FIGURENumber of Landsat scenes in Collection Tier per year processed during and .

FIGURE
FIGURE Example of the intra-annual land cover classification for a subset ( × pixels) of the Landsat / , in the Santarém-PA region; (A) the Landsat scenes processed for (n = ); (B) the classification metrics used in the RF post-classification; and (C) the annual land cover classification.The cropland extent was larger in the more recent images; the area increased in the mode metric; however, crop occurrence is concentrated in small areas along the year.The reference training and calibration samples showed a higher concentration of pastureland in this area, making the RF model classify most of the areas in this land cover class.

FIGURE
FIGUREAmazon biome land cover in and .The main change was pastureland replacing forest and cropland expansion along the southern Mato Grosso border with the Cerrado biome.

FIGURE
FIGURELand cover area yearly trend during and detected by the Landsat time series at the sub-pixel level.The red line is the Loess function time series trend.

FIGURE
FIGUREAnnual area change and the Loess trend line, with uncertainty interval, during and for the land cover classes of the Amazon Biome.

FIGURE
FIGURELand cover change mapping in the temporal breaking points.

FIGURE
FIGUREPixel counting and area adjustment for the land cover classes mapped.The area adjustment for was not presented because the reference data were not available.

FIGURE
FIGURE Spatial distribution of annual deforestation between and (A) and annual deforestation rate of old-growth forest and secondary vegetation > years old (B).

7
Mha between 1992 and 2022.It is important to highlight that secondary forest deforestation happens in areas previously deforested.When considering the total deforestation (i.e., in oldgrowth and secondary forests > 5 years old), the area of forest loss increases to 63.7 Mha from 1985 to 2022.
TABLE Land cover classes mapped and used in the change detection analysis.
TABLE Summary statistics for land cover change in the Amazon biome between and .
a Forest loss includes old-growth forest and secondary growth.b The area change is associated with water loss.