Temporally-Consistent Annual Land Cover from Landsat Time Series in the Southern Cone of South America

Abstract

The impact of land cover change across the planet continues to necessitate accurate methods to detect and monitor evolving processes from satellite imagery. In this context, regional and global land cover mapping over time has largely treated time as independent and addressed temporal map consistency as a post-classification endeavor. However, we argue that time can be better modeled as codependent during the model classification stage to produce more consistent land cover estimates over long time periods and gradual change events. To produce temporally-dependent land cover estimates—meaning land cover is predicted over time in connected sequences as opposed to predictions made for a given time period without consideration of past land cover—we use structured learning with conditional random fields (CRFs), coupled with a land cover augmentation method to produce time series training data and bi-weekly Landsat imagery over 20 years (1999–2018) across the Southern Cone region of South America. A CRF accounts for the natural dependencies of land change processes. As a result, it is able to produce land cover estimates over time that better reflect real change and stability by reducing pixel-level annual noise. Using CRF, we produced a twenty-year dataset of land cover over the region, depicting key change processes such as cropland expansion and tree cover loss at the Landsat scale. The augmentation and CRF approach introduced here provides a more temporally consistent land cover product over traditional mapping methods.

1. Introduction

Land cover across the globe has been shaped by climate over centuries and, in turn, has also influenced climate patterns [1]. In recent decades, land cover change from land use has negatively impacted ecosystems worldwide [2,3]. Commodity-driven deforestation rates, for example, remain high [4,5] as global consumption of agriculture and agroforestry products continues to grow [2,6,7]. Consequently, forests across the tropics have been cleared for agriculture [8,9,10], resulting in emissions of $2.6 {GtCO}_2 {yr}^{-1}$ [11]. While tropical regions such as the Amazon [12,13], Congo basin [14,15], and the Malay Archipelago [16,17,18] have received substantial attention, dry woody and grassland systems are increasingly being transformed and negatively impacted by land use [19,20,21,22]. Unfortunately, data on land change in extra-tropical systems of the Global South are lacking, even though these systems are particularly vulnerable to changing climate [2]. Timely estimates of land use and land cover (LULC) across these evolving landscapes are required to support a variety of efforts in the land system science community, such as deforestation monitoring, environmental impact assessments, and land conservation. In this paper, to address data shortcomings over one region in the Global South, we focus on the South American Southern Cone using 20 years of satellite-based land cover estimates. In doing so, we aim to better characterize and understand the environmental impacts and consequences of land use across the region.

The South American continent has undergone some of the planet’s fastest rates of forest loss and agricultural expansion over recent decades [23,24,25,26,27,28]. However, most LULC studies in South America have focused on the Amazon biome, leaving temperate regions on the continent largely understudied. The imbalance in regional foci is perhaps best exemplified by differences in research that has been conducted in the Amazon and Chaco forests [12,13,21,29,30,31]. The Chaco, which is the second largest forest system in South America, is understudied relative to the Amazon, even though deforestation rates are as high there as anywhere on the planet [21,24,30,32,33,34]. Further, studies that focus on soybeans in Brazil far outnumber those focused on other key soy-producing Southern Cone nations, such as Argentina, Bolivia, Paraguay, and Uruguay. Recently, however, interest and awareness regarding the impact of human activities and LULC on land areas in the Southern Cone over the last several decades has begun to increase [26,28,35,36]. Grassland and shrubland regions have received less attention, in general, thus leaving unanswered questions about LULC across large swaths of the Southern Cone [37,38,39]. As a result, key land change processes since the turn of this century in this region are largely undocumented.

As is now well known, land cover plays an important role in the delivery of ecosystem services and in the regulation of local and global biogeochemical cycles [40,41,42]. Time series of satellite imagery provide valuable information to monitor and enhance understanding of these processes. Over the Southern Cone, improved land cover estimates would also contribute toward better understanding the ecological and social impacts of forest clearance in the Chaco region [43,44]; small- and large-holder agricultural dynamics in Eastern Paraguay [45]; forest management in Misiones, Argentina [46]; shrub encroachment across Patagonia and in the Andes foothills [47]; and the impacts of grassland-to-cropland conversion across the Pampas and Espinal [48,49].

The LULC community now has Landsat imagery from 1972 to present available to monitor and map land change processes. Image classification is still a challenging endeavor, however, as illustrated by the plethora of studies and land cover datasets produced throughout the world [50,51,52,53,54,55,56,57]. To address challenges involved in image classification, satellite data acquisitions and processing methods have evolved to support temporally consistent, large-scale analyses of terrestrial events [58,59,60,61]. In particular, crop- [62,63,64,65,66,67,68,69,70,71,72,73,74] and forest-related studies [75,76] have benefited from new time series methods, analysis-ready data products, and the launches of Landsat 8 and 9 and Sentinel-2A and -2B. With increasingly longer historical records of Landsat imagery, land cover classification methods that address time series continue to evolve.

Our key objective in this work is to estimate annual land cover, and in doing so, develop a methodological framework for large-area estimates over long time spans. In this study, we use 20 years (1999–2019) of Landsat imagery to estimate annual maps of land cover across the South American Southern Cone region, combining multiple sources of field data, a dynamic time series smoother [77], and temporally-dependent land cover estimates. To estimate land cover transitions, we use conditional random fields CRFs [78,79] with a first-order conditional framework. The CRF modeling approach used here presents, to our knowledge, one of the most geographically expansive and longest time series applications of structured learning for land cover mapping in South America, while the land cover estimates presented here capture annual spatial and temporal dynamics at an unprecedented scale for the region.

2. Motivation for Sequential Land Cover Predictions with Structured Learning

The use of time series for land cover mapping has increased with better access to both data (e.g., Landsat) and high-performance computing. The conventional approach to estimate land cover over time is to treat each time step independently [9,24,80,81]. While recent studies have implemented post-classification probabilistic stabilization techniques [82,83,84], they still rely on existing land cover datasets as a starting point. Another novel approach for temporal land cover mapping is the continuous change detection and classification (CCDC) algorithm [85]. The CCDC method models seasonal variation in vegetation over time and uses estimated harmonic coefficients as classification predictors. This method works well in certain landscapes, such as boreal and tropical forests, where vegetation dynamics are often consistent over time and land cover change results in a pronounced satellite signal [76,86]. However, in more complex landscapes, such as dry savannas, vegetation cycles are not guaranteed to behave in a consistent, periodic manner, even when no land cover change has occurred. Moreover, many cropland systems include diverse mixtures of crop species and, therefore, do not follow predictable crop rotations.

In this study, we present a method in which temporal dependencies are incorporated into land cover predictions by directly modeling conditional dependencies across time series rather than post-inference of classified maps or by estimating land cover change events (i.e., breakpoint detection) from satellite signals. We present a framework to estimate temporally-consistent land cover over long time spans using structured learning of reconstructed satellite signals. Specifically, to estimate annual land cover we use linear-chain conditional random fields (CRFs) (see Section 5.4 for details) fit to weekly Landsat time series (described in Section 5.1 and in [77]), and apply this method to map annual land cover across 20 years in the Southern Cone of South America.

3. Study Area

The Southern Cone is a region in South America that typically refers to Argentina, Chile, and Uruguay. In this study, we also include Paraguay and the Bolivian lowlands because of their importance in regional LULC (Figure 1). In total, the included area within these five countries is approximately 4,972,837 km². Placing this in a continental perspective, the region contains the largest forest system (the Chaco) in South America outside of the Amazon, the largest contiguous area of remnant Atlantic forest in Misiones, Argentina, as well as an estimated 18% (230,849 km²) of the planet’s harvested soybean area in 2017 [23].

There is considerable regional complexity in climate and vegetation types across the Southern Cone. The region’s climate ranges from tropical in the north, desert in the northern Andes and the south, Mediterranean along the Pacific Ocean, to temperate in the eastern grasslands. The northeast, dominated by grassland and croplands, is characterized by high mean annual precipitation (>1200 mm) and high mean annual temperatures (15–20 °C) [87,88]. The south and northwest are dominated by grasslands and shrublands and are characterized by relatively low mean annual precipitation (150–500 mm) and mean annual temperatures ranging from 0–12 °C [87,88]. The Mediterranean climate in Chile supports a mixture of grasslands, shrublands, and woodlands. Lastly, this study area covers an elevation gradient from 20 m above sea level in the Pampas to more than 3000 m in the Andes.

4. Data

4.1. Satellite Data

We acquired Landsat TM, ETM+, and OLI Collection 1 (C1) surface reflectance data from the US Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center on-demand Science Processing Architecture (ESPA). We defined each 12-month period in this study as starting and ending on 1 July, corresponding to seasonality in the Southern Hemisphere. As outlined in [77], we padded the time series end years by three months to ensure data continuity. Therefore, we selected all Landsat scenes from 1 April 1999 to 1 November 2019, using the same cloud cover thresholds and pixel angle generation described in [77]. We processed and stored all Landsat bands and pixel angle layers into non-overlapping, 150 km × 150 km grid cells, gridded to the South America Albers equal-area conic projection. Data density was generally higher in the middle latitudes of the study region compared to the northern and southern locations (Figure A2).

4.2. Field Data

Using georeferenced samples from multiple sources, we compiled a large set of land cover polygons across the study region from (1) field campaigns in central and northern Argentina during the summers of 2017 and 2018 conducted by researchers at Boston University; (2) field campaigns in northern Argentina and central Bolivia during the winter of 2014; (3) field campaigns in 2016 and 2017 conducted by the Instituto Nacional de Investigación Agropecuaria (INIA) of Uruguay; (4) a field campaign in 2018 conducted by the Instituto Nacional de Tecnología Agropecuaria (INTA) of Argentina; (5) multi-year field campaigns over 2014–2018 conducted by the Buenos Aires Bolsa de Cereales (Grains Exchange); and (6) land cover data collected on-screen via image interpretation by analysts at Boston University.

The field campaigns consisted of a mix of systematic and opportunistic sampling of land cover and agricultural land use (categories described in Section 5.5), but were designed to collect samples of agricultural land use. Georeferenced samples (collected in the field using GPS) were converted into polygons and manually delineated to agricultural field boundaries by interpreters. For samples collected in cultivated areas, to ensure that each field corresponded to a single crop, we interpreted boundaries from high-resolution imagery in Google Earth (GE) along with annual Landsat composites to illustrate within-field spectral homogeneity. After investigators with local knowledge reviewed field polygons for quality assurance, the field dataset included 11,464 field sample polygons of land cover and crop type.

To complement field samples, we also collected on-screen samples of land cover classes using a global sampling grid provided by the UN Food and Agriculture Organization (FAO). The grid consists of 10 km × 10 km grid cells derived from a systematic sampling design based on each latitude and longitude intersection [89]. These additional samples were necessary to provide full geographic coverage and thematic representation across the Southern Cone. Within each grid cell, we attempted to digitize three polygons with class labels defined for every available GE image acquisition over time. Samples were subjectively selected and interpreted based on time series of Landsat spectral reflectance and annual Landsat composites along with GE high-resolution imagery and the interpreters’ local knowledge. The total number of land use/cover polygons was 46,532 across all sources and years. Figure A1 illustrates the distribution of samples across the study area.

5. Methodology

5.1. Image Preprocessing and Time Series Reconstruction

We adjusted surface reflectance values for bidirectional reflectance distribution function (BRDF) effects using the c-factor method [59,77,90]. Using the BRDF-corrected surface reflectance data, we calculated the two-band enhanced vegetation index (EVI2) [91] and the woody index (WI) [92] for every Landsat image as

$$EVI2=2.5\frac{NIR-red}{NIR+2.4red+1}$$

(1)

and

$$WI=\left\{\begin{array}{cc}0\hfill & \mathrm{if}(red+NIR)>0.5\hfill \\ 1.0-\frac{red+NIR}{0.5}\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(2)

respectively. The EVI2 was chosen because of its proven utility for the detection of agricultural phenology [93]. The WI was chosen for its non-agricultural use, particularly in dry savanna woodlands, as a large swath of the study area contains similar biophysical characteristics to the region where WI was first tested (i.e., Queensland, Australia) [92].

To generate consistent data across all years, we used a dynamic temporal smoothing algorithm (DTS) [77] to reconstruct weekly time series of EVI2 and WI. The DTS produces 60 values per 12-month cycle.

5.2. Land Cover Classification Scheme

Following land cover definitions by the USGS Multi-Resolution Land Characteristics Consortium (MRLC) National Land Cover Database (NLCD), we designed a classification scheme of 8 land cover categories for the study region (

Table 1. Land cover classes and definitions.

Class Label	Description
Cropland	Managed lands for production of annual and perennial crop species (excluding tree crops); this class represents row-crop agriculture, such as maize, soybeans, wheat, and rice.
Water	Open-water bodies, such as lakes and rives.
Developed	Urban areas, built-up structures, and roads.
Barren	Barren land, snow, or ice.
Trees	Tree cover with $>=30\%$ canopy cover (natural or managed); this class includes agroforestry such as pine and eucalyptus plantations.
Shrubland	Shrubs or cactus cover with $>=30\%$ canopy cover.
Grassland	Herbaceous grassland (natural or managed) and savannas with <30% canopy cover; this class includes pastureland used for livestock grazing.
Wetland	Seasonal wetlands, including prolonged flooding (e.g., Pampas).

). The field data collection was designed to record individual crops. Thus, for this study we aggregated all crop samples into a single cropland class and all agroforestry samples into the tree cover class (Figure 2).

5.3. Features Used for Land Cover Predictions

For land cover mapping, we used the EVI2 time series and seasonal time series metrics calculated from the weekly EVI2 and WI time series (

Table 2. Features used for vegetation indices EVI2 and WI ($VI$) to fit classification models for each grid cell. * The weekly time series were generated for EVI2 and WI and used to calculate aggregate metrics. However, only the EVI2 weekly values were used as predictive features. ** The maximum and minimum slopes were computed with a moving window of 1.5 month width.

Metric	N	Description
$EVI2$	60	* The weekly gap-filled time series
$V{I}_{q50}$	4	Median of time series each quarter
${VI}_{\mu }$	1	Mean of the time series
${VI}_{p50}$	1	Median of the time series
${VI}_{CV\mu }$	1	Coefficient of variation with ${VI}_{\mu }$
${VI}_{CVp50}$	1	Coefficient of variation with ${VI}_{p50}$
${VI}_{min\left(m\right)}$	1	** Minimum slope
${VI}_{max\left(m\right)}$	1	** Maximum slope
${VI}_{p5}$	1	5th percentile of the time series
${VI}_{p10}$	1	10th percentile of the time series
${VI}_{p25}$	1	25th percentile of the time series
${VI}_{p75}$	1	75th percentile of the time series
${VI}_{p90}$	1	90th percentile of the time series
${VI}_{p95}$	1	95th percentile of the time series
${VI}_{argmax}$	1	Date of maximum time series value

). We used the time series of both indices to compute aggregate metrics shown in Table 2. However, we excluded the WI time series from the set of predictive features to reduce the computational load on the classifier. In total, we generated 94 image layers (60 EVI2 values + 17 (x2) aggregate metrics) that we used as predictive features (X).

5.4. Land Cover Estimates with Conditional Random Fields

A conditional random field (CRF) is a discriminative model used to predict sequences of data [78]. A CRF models the conditional probability of labels and variables over a joint distribution and makes no assumptions about variable independence. In fact, CRFs are designed to estimate the most likely sequence of state transitions assuming codependency. Because of these model properties, CRFs are well suited to remote sensing problems with temporal codependencies such as land cover change. However, CRFs are more widely used for natural language processing [94] than in the field of remote sensing. Structured learning for remote sensing applications has recently increased, with several studies using CRFs to address spatial and temporal dependencies in predictive modeling of land use and land cover [95,96,97].

To fit a CRF classifier, class conditional likelihoods from a given set of labeled data (e.g., land cover) and predictive features are maximized by modeling a sequence of states (in our case, across time). A model with first-order transition dependencies, called a linear-chain CRF, is defined by the sequential dependence of the current state observation i and the previous state observation ($i-1$) as

$$P\left(y\right|x,\omega )=\frac{exp\left\{{\sum }_{i=1}^n{\sum }_{j=1}^m{\omega }_j{f}_j(x,i,y_{i-1},y_i)\right\}}{{\sum }_{y{}^{\prime }\in Y}exp\left\{{\sum }_{i=1}^n{\sum }_{j=1}^m{\omega }_j{f}_j(x,i,{y{}^{\prime }}_{i-1},{y{}^{\prime }}_i)\right\}}$$

(3)

where y is a sequence of labels and x is a sequence of variables. The feature function, $f_j(x,i,y_{i-1},y_i)$, is a set of characteristics that describes each predictive variable, based only on the current and previous state. Feature functions are designed by the user and, therefore, can be customized. In the model fitting phase, the feature weights ($\omega$), which are initialized randomly in $P\left(y\right|x,\omega )$, are derived by iteratively maximizing the class conditional likelihoods of all states (n) and features (m) by $argmaxP\left(y\right|x)$ of a regularized log-linear function. The numerator in Equation (3) is used to normalize the class conditional probabilities over all possible labels (in this study, Y is equal to the land cover classes listed in Table 1).

5.5. Land Cover Sample Design

We used a distance-based search to collect land cover samples around each 150 km × 150 km grid cell, thus generating a different set of samples and model for each grid cell. We illustrate this process for a single grid cell (Figure 3). First, we locate all land cover samples within 250 km from a grid cell center. If a grid cell target sample pool is not met, we expand the search from 250 km up to 500 km, in 100 km intervals, until sufficient samples have been collected. The radial expansion provides flexibility over regions of sparse data or with rare land cover classes such as ‘developed’. After a sample pool around a grid has been collected, they are not yet ready to use to train a CRF. The sample pool shown in Figure 3 consists of various dates, depending on the data collection source (Section 4.2). However, this sample pool does not contain repeat data, which is pivotal in order to train a sequential-based classifier such as CRF. We describe our solution to this data shortcoming in the following paragraph (Section 5.6).

5.6. Temporal Augmentation

One challenge of supervised learning with structured models is the collection of repeat training data. A CRF, for example, requires a sequence of target labels to optimize the transition weights during the fitting phase. Yet, collecting repeat labels, especially over 20 years, is difficult because the task is laborious and consistent satellite acquisitions are not guaranteed to be available. We address this impediment by temporally augmenting pseudo land cover sequences, which we then use as our model training data. More specifically, we generate training sequences using the field data from Section 4.2. Figure 4 illustrates the random selection to generate a single augmented pseudo time series sequence. The field data are distributed over the timeframe but there are no cases of repeat visits. Therefore, a temporally augmented sequence consists of samples drawn at random from a sample pool. The pseudo sequences consist of real data, but are arranged in a synthetic manner. Thus, we defined likely land cover transition scenarios to help guide the augmentation into realistic land cover transitions (

Table 3. Land cover transitions used in the temporal augmentation of pseudo time series sequences. The first column describes the land cover transitions. Single land covers represent stable sequences. The ↔ symbol indicates the sequence was generated in both directions. The second column describes the land cover change that the sequence is intended to represent. The third column shows the percentage breakdown of the land covers in the sequence. As an example, the trees-to-cropland sequence represents clearance for agriculture when sorted by trees → cropland, and cropland abandonment and regrowth when sorted by cropland → trees. In the trees → cropland sequence, the first 70% of the sequence is trees (i.e., 17 years) followed by 30% (3 years) of cropland samples, and vice versa for cropland → trees.

Transition	Land Cover Simulation	Percentage of Sequence
Cropland	Stable cover	100
Developed	Stable cover	100
Trees	Stable cover	100
Shrubland	Stable cover	100
Grassland	Stable cover	100
Water	Stable cover	100
Wetland	Stable cover	100
Barren	Stable cover	100
Cropland ↔ Grassland	Agricultural shift	70/30
Cropland ↔ Barren	Fallow	80/20
Trees ↔ Shrubland	Degradation/Regrowth/Greening	70/30
Trees ↔ Grassland	Clearance/Regrowth/greening	70/30
Shrubland ↔ Grassland	Clearance/Abandonment	70/30
Shrubland ↔ Cropland	Clearance/Abandonment	70/30
Trees ↔ Cropland	Clearance/Abandonment	70/30
Wetland ↔ Grassland ↔ Cropland	Seasonal wetlands	33/33/33

). We repeat the framework in Figure 4 using Table 3 as a guide to determine which land cover data pools to randomly draw from. The land cover transitions characterized in this study (Table 3) are representative of the most common landscape change processes in the Southern Cone—deforestation and forest degradation for ranching [21,24,98] and commodity crop agriculture [32,33,45,99,100,101], agricultural intensification [24,25,26,28], afforestation [102,103,104,105], and woody encroachment and forest regeneration [47]. In total, we generated 2500 random sequences for each grid cell, regardless of the sample count within the grid cell sampling radius.

5.7. CRF Model Fitting and Land Cover Estimates

For each 150 km × 150 km grid, we fit a CRF classifier using an exponential transformation of the time series features X (Table 2) as

$$\overrightarrow{X}=[(exp\left(X\right)-exp\left(0\right))/(exp\left(1\right)-exp\left(0\right))]5$$

(4)

which was applied to stretch feature weights and to help maximize class separability. Then, for each grid cell we optimized $P\left(y\right|\overrightarrow{X},\omega )$, where our feature functions $f_{j\in m}$ were the 94 $\overrightarrow{X}$ values. Experimentation showed that the exponential transformation from a [0,1] to [0,5] dynamic range improved CRF predictions on test samples. We used the limited-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm with L1 and L2 regularization [106,107], optimized with the sklearn-crfsuite Python interface [108] to the crfsuite software package [109]. We applied the grid cell’s fitted CRF model to the full timeframe (i.e., 1999–2019) for each pixel, with CRF class probabilities for each two-year timeframe. For each two-year timeframe, we used the land cover class with the maximum CRF probability as the class label.

5.8. Map Assessment

Guidelines for best practices in land cover change assessment stipulate that unbiased estimators of reference data be used instead of pixel counting [110,111,112,113]. Following these recommendations, we set up a stratified random sample design based on unbiased area estimates with Landsat pixels as the sampling units. We chose to assess four time periods of the map outputs (2000–2005, 2005–2010, 2010–2015, and 2015–2018) rather than the full time series, considering trade-offs between balancing uncertainty reporting and time. We consider this a reasonable compromise given the total area of the study area and the length of the time period [76]. For the purposes of map assessment, we aggregated shrubland, developed, barren, water, and wetland classes into a class called “other”. Then, we assessed stable tree, grassland, cropland, and other strata, as well as all the change strata representing the dynamics among these classes. We used a target standard error of 1.5%, resulting in 587, 596, 608, and 529 samples for each time period, respectively (Figure 5). For the last time period between 2015 and 2018, this target corresponds to 6666 km² of trees to grassland change and approximately 1164 km² of grassland to trees change.

The assessment of land cover change in these four time periods was executed using independent high-resolution imagery and medium-resolution composites that were interpreted by analysts with local knowledge and expertise in the Southern Cone region. Specifically, reference data were interpreted using the Augmented Visual Interpretation tool available in Collect Earth via the Open Foris platform [114], which enables evaluation using high-resolution imagery from Google Earth along with time series and composites of MODIS and Landsat in Google Earth Engine.

6. Results

Our study generated annual maps of land cover for the Southern Cone of South America at 30 m spatial resolution from 1999/2000 to 2018/2019 (Figure 6; annual change Figure 7, Figure 8, Figure 9 and Figure 10). In addition, we calculated area estimates of land cover with 95% confidence intervals for four approximately 5-year time periods for stable land categories and conversions

Table 4. User’s and producer’s accuracy (%) for four different time periods assessed in this study.

	2000–2005		2005–2010		2010–2015		2015–2018
Land Cover	Users	Producers	Users	Producers	Users	Producers	Users	Producers
Stable cropland	61	68	66	84	49	89	56	82
Stable other	70	92	84	91	84	92	76	89
Stable grassland	53	78	46	70	59	81	61	84
Stable trees	80	86	76	87	69	77	72	83
Cropland to other	100	6	100	3	100	8	67	3
Cropland to grassland	100	14	94	15	90	28	86	15
Cropland to trees	100	1	100	1	100	4	67	0
Other to cropland	67	5	75	7	100	13	100	6
Other to grassland	64	10	65	7	84	21	64	4
Other to trees	100	4	0	0	93	15	0	0
Grassland to cropland	92	41	76	35	81	31	100	14
Grassland to other	55	6	72	9	69	24	50	7
Grassland to trees	33	5	94	19	100	19	100	2
Trees to cropland	100	25	100	28	100	30	100	7
Trees to other	71	7	67	21	83	5	100	2
Trees to grassland	79	18	83	46	100	34	60	7

and

Table 5. Stratified unbiased area estimate and standard error of the error-adjusted area confidence interval (SE) for four different periods assessed in this study.

	2000–2005		2005–2010		2010–2015		2015–2018
Land Cover	Area (km2)	SE (km2)	Area (km2)	SE (km2)	Area (km2)	SE (km2)	Area (km2)	SE (km2)
Stable cropland	153,211	47,651	179,796	47,739	148,476	45,228	198,657	59,345
Stable other	1,340,525	119,452	1,645,780	114,346	1,620,578	110,114	1,519,434	128,984
Stable grassland	569,891	95,466	528,926	100,957	567,818	94,171	614,316	95,434
Stable trees	1,070,482	106,175	971,476	104,825	982,111	122,832	996,839	116671
Cropland to other	47,142	35,275	36,491	32,484	36,215	27,364	44,362	36,986
Cropland to grassland	39,728	29,899	62,484	35,931	67,793	36,533	65,208	40,049
Cropland to trees	58,390	38,591	36,117	30,562	10,397	14,219	48,457	40,886
Other to cropland	40,934	31,285	42,742	31,706	40,695	30,875	57,780	42,970
Other to grassland	130,780	56,644	97,214	49,063	89,557	46,749	161,125	70,127
Other to trees	81,218	51,998	82,304	53,449	130,063	61,859	93,922	57,074
Grassland to cropland	135,812	44,177	114,638	44,586	71,004	36,965	54,471	36,740
Grassland to other	143,438	60,963	91,720	49,054	57,485	35,089	44,365	35,223
Grassland to trees	56,982	41,168	89,217	44,966	134,505	56,847	72,051	48,347
Trees to cropland	29,715	22,355	25,014	19,910	13,727	12,026	24,380	26,416
Trees to other	125,544	57,870	27,574	26,855	49,169	40,565	48,448	40,364
Trees to grassland	69,920	40,419	64,891	31,069	77,341	40,286	55,509	41,197

).

6.1. Accuracy Assessment

The overall accuracy of the Southern Cone map was 69% for 2000–2005, 73% for 2005–2010, 73% for 2010–2015, and 71% for 2015–2018. Class-specific user’s and producer’s accuracies for the four assessment time periods were variable (Table 4). Higher accuracies were obtained for stable land cover categories, especially for tree cover, as well as for larger map classes such as the combined “other” category, which included wetland, water, shrubland, and barren classes. Distinction between grasslands and croplands, on the other hand, was difficult, and as a result the user’s and producer’s accuracies for those two classes were lower compared to tree cover (Table 4). Confusion between these two classes was due to spectral similarities and the fact that they were often collocated, sometimes on the same field during different years. There was also map confusion with classes that represented transitional situations (e.g., shrubland encroachment on grasslands).

Producer’s accuracies were generally low and user’s accuracies were generally high for conversions between classes. Producer’s accuracies for transitions between 2010 and 2015 had higher accuracies compared to the other time periods and so did transitions among tree cover and the three other categories (i.e., cropland, grassland, and other) (Table 4). User’s accuracies were higher compared to producer’s accuracies—they ranged from 33% to 100% but were generally above 60%. In many cases, the commission error is very low (or even nonexistent) while the omission error is quite high, indicating that conversion to croplands or tree cover loss, for example, were omitted from the correct class in the final map (Table 4). In other words, while the maps do not include a lot of incorrect classifications, they sometimes failed to identify areas of tree cover loss and conversion to croplands. Given these results, our results likely represent an underestimate of how much the subcontinent has changed since 1999.

6.2. Comparison against Independent Sources of Data

The annual estimates of land cover allowed us to track trends over two decades in cropland area across the Southern Cone region. To complement the accuracy assessment presented in Table 4 and Table 5, we compared our results against data from the Ministerio de Agricultura, Ganadería, y Pesca (MAG) of Argentina [115] (Figure 11). We included all row crops reported by MAG, excluding banana, canary grass, jojoba, lime, mandarin, orange, pomelo, tea, tung, and yerba mate. Contrary to MAG’s estimates of planted area, our cropland estimates are of area cover over a 12-month period and, therefore, do not account for double cropping (i.e., over each pixel, we map crop cover, regardless of the number of crop cycles). As a result, our land cover area estimates of cropland will typically be lower than planted (or sown) area statistics reported by MAG. To transform our crop area estimates into planted area, we estimated the number of crop peaks, or cycles, over a 12-month period and multiplied the peak count by our cropland mask. For this comparison, we summed planted area by pixel counting and not by the unbiased area estimates reported in this study.

In general, our planted area estimates were correlated with MAG estimates in departments and provinces with more than 1000 km² of estimated planted crop area (Figure 11). We tended to underestimate planted crop area compared to MAG estimates in administrative units with less crop area and in drier regions along the southern Pampas. A prime example of this discrepancy is La Pampa province, shown by the brown dots in Figure 11. At every time period, La Pampa had the largest difference between our planted area estimates and MAG’s planted area. Similarly, southern Buenos Aires departments also had noticeably higher estimates by MAG.

We also compared our CRF land cover estimates with an independent product (MapBiomas Chaco [116] and Pampas [117]) by generating an error matrix with overall accuracy. We used all validation points (from Section 5.8) that intersected the MapBiomas Chaco and Pampas map extents. The accuracy results from this simple comparison are shown in

Table A1. Overall accuracy (%) of two land cover products for four different time periods. These accuracy numbers were generated using the same set of validation samples used in the unbiased assessment. To compare the two maps at a comparable legend, we recoded MapBiomas Chaco maps to match our land cover legend as close as possible. The MapBiomas Chaco product contains twenty land cover categories, including multiple tree cover (all woodland classes were recoded to “trees” except “flooded natural woodlands”) and grassland (all grassland categories were recoded to “grassland” except “flooded grassland”) categories.

	2000–2005	2005–2010	2010–2015	2015–2018
CRF	90	91	93	79
MapBiomas Chaco	73	75	74	79

. Note that the numbers reported here for CRF (this study) are higher than the numbers reported for the overall accuracies across the entire study region, while the accuracies for the MapBiomas regions are lower than reported by the MapBiomas projects (consistent range of around 80%). The set of validation points used are biased toward the CRF maps, as they were generated using a stratified random sample of change estimated from this study. Additionally, Figure A5 illustrates the temporal stability of CRF for a single pixel compared to the MODIS MCD12Q1 global land cover product and MapBiomas Chaco product.

6.3. Spatial and Temporal Trends of Change in the Southern Cone of South America

Our study estimated land cover annually over twenty years from 1999/2000 to 2018/2019 and our results revealed widespread and systematic changes in the distribution of land cover in the Southern Cone of South America in the 21st century (Figure 6). Changes took place across all ecoregions of the subcontinent but they were most pronounced in the savanna ecoregions in the northeast (e.g., Humid Pampas and Uruguayan Savanna) and in the forested ecoregions in the north (e.g., Gran Chaco and Atlantic Forest). At the regional scale, noticeable net changes included loss of tree cover (Figure A4) and replacement by natural or managed grassland in western Paraguay, as well as an increase in cropland area at the expense of grasslands across the Argentine Pampas. While the regional net changes revealed striking land cover transformations, one of the key contributions of this paper was the annual estimates in land cover. For example, the expansion of cropland in the Uruguayan savannas (Figure 7), loss of tree cover in western Paraguay (Figure 8) and eastern Bolivia (Figure 9), as well as gradual cropland expansion near Santa Cruz, Bolivia (Figure 10). The CRF model produced annual land cover estimates with minimal “pixel flipping” noise from year to year, as illustrated in Figure 7, Figure 8, Figure 9 and Figure 10. Additionally, Figure A3 illustrates land cover change at a larger scale, while Figure A4B highlights the consistency of this method for accurately detecting long-term gradual change events.

In the beginning of the 21st century (from 2000 to 2005) we estimated forested area across the Southern Cone of South America at 1,070,482 km² (±106,175 km²), the area of croplands at 153,211 km² (±47,651 km²), and grasslands at 569,891 km² (±95,466 km²) (Table 5). At the end of the study period (2015–2018), forest area decreased to 996,839 km² (±116,671 km²) while croplands expanded to an area of 198,657 km² (±59,345 km²) and grasslands to 614,316 km² (±95,434 km²) (Table 5). The magnitude of deforestation was the highest at the beginning of the time period and was associated with an increase in the area of “other”, a category that included shrublands and barren ground (Figure 12). Conversely, an increase in conversion from tree cover to grasslands and croplands in the last ten years of the study period corresponded to a decrease in the area converted from forests to the “other” category (Figure 12). Throughout the 20-year study period, conversion from tree cover to grasslands was consistently higher than conversion from tree cover to croplands. Starting in 2005, conversion directly to grasslands became more widespread compared to conversion to “other” land cover (Figure 12).

7. Discussion

The Southern Cone of South America experienced large-scale land cover change over the course of the last twenty years since 2000 (Figure 6). Stable tree cover experienced the most substantial losses—it decreased to 996,839 km² (±116,671 km²) by the end of 2019 (Table 5). Forested landscapes are mostly converted to grasslands and “other” land cover types (e.g., shrublands) and only later converted to croplands (Figure 12), indicating that there is a time lag from deforestation to land utilization and that livestock production is the proximate driver of deforestation while crop production is the underlying driver [21,24]. Our results reveal noticeable deforestation in the Dry Chaco of western Paraguay (Figure 8), northern Argentina (Figure 6), and in the Chiquitania of eastern Bolivia (Figure 9 and Figure 10). Deforestation in the Southern Cone region is largely driven by agricultural production–soybean cultivation in the Atlantic Forest of eastern Paraguay and the Dry Chaco of northern Argentina [32,33,45,99,100,101] and cattle ranching in the Dry Chaco of western Paraguay and the Chiquitania of eastern Bolivia [21,24,98]. These findings are consistent with recent publications documenting large-scale deforestation in the dry woodlands of South America [24,25,28,36], a significant source of carbon emissions since the 1980s [20].

While croplands sometimes expanded onto forests in the 21st century, they mostly and consistently, throughout the 20-year study period, expanded onto grasslands (Figure 12, Table 5). In total, cropland area increased to 198,657 km² (±59,345 km²) by the 2015–2018 time period. In only 5 years—from 2000 to 2005—135,812 km² (±44,177 km²) of grasslands were converted to croplands. In addition to pastureland conversion to croplands along forest frontiers, our study demonstrates that there is cropland expansion at the expense of grasslands (natural, pasturelands, and rangelands) in the Argentine Pampas and the Uruguayan Savanna (Figure 6 and Figure 7). Expansion of croplands across Argentina is especially prevalent in the provinces of Buenos Aires, Córdoba, and Entre Ríos (Figure 6 and Figure 11) and is largely driven by expansion of soybean cultivation [26,48]. Taken together, conversion of natural grasslands, pasturelands, and tree landscapes to croplands represents large-scale agricultural intensification across the subcontinent that started as early as the 1960s [24,25,26,28].

A relatively less studied land cover change phenomenon across the Southern Cone is the increase in tree cover due to reforestation and afforestation. Over the last several decades, there has been a proliferation of plantations in Uruguay [102], eastern Argentina, and Chile [103,104,105] (Figure 6). According to our results, increase in tree cover in the Southern Cone predominantly happens on land that was previously either grasslands or “other” (Figure 12). Most conversion to tree cover occurred between 2010 and 2015, with 134,505 km² (±56,847 km²) from grasslands and 130,063 km² (±61,859 km²) from “other” (Table 5). They are most likely the result of afforestation secondary growth, reforestation, and woody encroachment, respectively. Note, however, that it is common across Chile for plantations to replace natural forests, a transition that is not captured by our study because we did not differentiate between different kinds of tree cover [103,104,105].

This study represents—to our knowledge—the first application of conditional random fields to land cover change classification at a large spatial scale using two decades of Landsat observations. The CRF classifier is well suited for studies of land cover change because it learns temporal transitions in the presence of noisy time series and, therefore, helps to address spurious annual land cover estimates. Our results demonstrate that the CRF approach to land cover classification indicates promise but future studies can improve the map accuracies achieved in this analysis by addressing data collection and advancing the CRF framework. For example, data collection at “revisit sites” would strengthen temporal transition model learning (i.e., sequence-based models such as the CRF used in this study) and provide an important dataset to use for model benchmarking and validation. The need to scale and enrich the sequential modeling framework that was used here is evident, given the underestimation of some key land cover transitions. The sample design could also be tested on a sequentially-based deep learning framework. As shown in [118], large-scale land cover mapping using sequence-based deep learning architectures is becoming a reality. Lastly, the methods in this study could easily be extended and enhanced by the inclusion of Sentinel-2 data. Given that the Sentinel-2 sensors share similar radiometric properties to the Landsat series of sensors, harmonization between the two is relatively straightforward [59]. A complete replication of the methods used in this study would require a limited set of optical wavelengths, and there are no sensor restrictions to the CRF framework. Therefore, Sentinel-2 and additional bands or indices could be tested with minimal changes to the modeling framework.

Despite methodological advancements in this study, there are still challenges and limitations that we face in mapping land cover and land cover change across the large latitudinal, altitudinal, and climate gradient in the Southern Cone region of South America. For example, our map accuracies are lower than targeted for less common classes and transitions such as “other” to grassland and grassland to “other” (Table 4). These transitions are consistent with Chaco forest clearance that can be characterized by a land cover change sequence (tree to grassland to shrubland to cropland) that is hard to capture due to similar spectral response of different classes and debris that are often left on the field. The spectral similarity and spatial proximity of grasslands and croplands makes them difficult to separate, especially in the Pampas and Uruguayan Savanna. Wetlands, which are also part of the “other” class in the accuracy assessment, are challenging to discriminate since they are usually seasonally flooded and occur in grasslands and treed landscapes [119]. We also highlighted some underestimation of planted crop area (Figure 11). Many of La Pampa and Buenos Aires departments with higher MAG planted area predictions (Figure 11) are located in the Southern Pampas, where the historical precipitation gradient declines north-to-south, livestock grazing becomes more commonplace than row-crop agriculture, and the sea of maize and soybeans that defines the Argentine agricultural heartland gives way to cereals and forage crops, such as sunflowers, wheat, barley, and sorghum. A large proportion of underestimated departments were along the southern fringes of this Pampas region. Based on our analysis and collective expertise on regional agriculture in Argentina, we concluded that this was the main region where we incorrectly classified cropland with herbaceous cover (grassland/pastureland). In addition, reference data collection carries its own uncertainty due to difficulties in interpretation of land cover, especially early in the time series when high-resolution data are less available and sometimes the available imagery has lower resolution. These difficulties have been previously described in the literature and are not unique to this study [120]. Lastly, this study was fortunate to have funding for large-scale data processing on a high-performance compute environment (Boston University Shared Computing Cluster). The methods used in this study, particularly the time series reconstruction, cannot be efficiently applied over large regions without a compute cluster. Although we did not record processing compute times, replication of this work using a long time series of nearly every available Landsat image would likely require access to a high-performance compute environment.