Agricultural Expansion in Mato Grosso from 1986–2000: A Bayesian Time Series Approach to Tracking Past Land Cover Change

: Landsat 5 has produced imagery for decades that can now be viewed and manipulated in Google Earth Engine, but a general, automated way of producing a coherent time series from these images—particularly over cloudy areas in the distant past—is elusive. Here, we create a land use and land cover (LULC) time series for part of tropical Mato Grosso, Brazil, using the Bayesian Updating of Land Cover: Unsupervised (BULC-U) technique. The algorithm built backward in time from the GlobCover 2009 data set, a multi-category global LULC data set at 300 m resolution for the year 2009, combining it with Landsat time series imagery to create a land cover time series for the period 1986–2000. Despite the substantial LULC di ﬀ erences between the 1990s and 2009 in this area, much of the landscape remained the same: we asked whether we could harness those similarities and di ﬀ erences to recreate an accurate version of the earlier LULC. The GlobCover basis and the Landsat-5 images shared neither a common spatial resolution nor time frame, But BULC-U successfully combined the labels from the coarser classiﬁcation with the spatial detail of Landsat. The result was an accurate ﬁne-scale time series that quantiﬁed the expansion of deforestation in the study area, which more than doubled in size during this time. Earth Engine directly enabled the fusion of these di ﬀ erent data sets held in its catalog: its ﬂexible treatment of spatial resolution, rapid prototyping, and overall processing speed permitted the development and testing of this study. Many would-be users of remote sensing data are currently limited by the need to have highly specialized knowledge to create classiﬁcations of older data. The approach shown here presents fewer obstacles to participation and allows a wide audience to create their own time series of past decades. By leveraging both the varied data catalog and the processing speed of Earth Engine, this research can contribute to the rapid advances underway in multi-temporal image classiﬁcation techniques. Given Earth Engine’s power and deep catalog, this research further opens up remote sensing to a rapidly growing community of researchers and managers who need to understand the long-term dynamics of terrestrial systems.


Introduction
Remote sensing is a crucial tool for observing and monitoring changes in land use and land cover across large areas, allowing users to quantify changes across large regions using satellite images, LIDAR images, and aerial photographs [1,2]. The Landsat series is the most complete and longest running satellite imagery available (1972-present). Its opening to the public a decade ago has allowed access to a wealth of data for researchers [2][3][4] and opened up new avenues of research [5]. Given the enormous data volume now available, however, it has remained a substantial challenge to quickly and efficiently process Landsat images and create time series maps, particularly with legends that are specific to a given project's needs.
Time series of land use and land cover (LULC) allow researchers to identify and track changes in land use patterns and to study the effects of land use in hydrological processes, ecosystem functions, and biogeochemistry [3,6]. By providing opportunities to look reliably into the past [7][8][9][10][11], time series have increasingly been used to provide models with better baselines to make estimates of deforestation rates, land use, and effects of policy decisions [6,7,12,13]. While time series classifications, particularly of forest loss/gain, are being produced and are of high interest, customizing a time series for a given project's needs is difficult due to the time and expertise required to create classifications and the daunting complexity of making the results coherent between iterations [14].
In this study, we apply the Bayesian Updating of Land Cover: Unsupervised method (BULC-U), the recent extension of the BULC algorithm to accommodate unsupervised classifications of Landsat satellite imagery [15], to interpret imagery across the time period 1986-2000, and to create an annual time series of major LULC categories of Forest, Deforested Land, and Shrubland/Grassland at 30 m spatial resolution. With extensive field experience in northeastern Mato Grosso, Brazil (Figure 1), we chose to study that region's noteworthy changes in land cover and land use patterns during the late 1980s and 1990s.
Remote Sens. 2020, 12, x FOR PEER REVIEW 2 of 18 enormous data volume now available, however, it has remained a substantial challenge to quickly and efficiently process Landsat images and create time series maps, particularly with legends that are specific to a given project's needs. Time series of land use and land cover (LULC) allow researchers to identify and track changes in land use patterns and to study the effects of land use in hydrological processes, ecosystem functions, and biogeochemistry [3,6]. By providing opportunities to look reliably into the past [7][8][9][10][11], time series have increasingly been used to provide models with better baselines to make estimates of deforestation rates, land use, and effects of policy decisions [6,7,12,13]. While time series classifications, particularly of forest loss/gain, are being produced and are of high interest, customizing a time series for a given project's needs is difficult due to the time and expertise required to create classifications and the daunting complexity of making the results coherent between iterations [14].
In this study, we apply the Bayesian Updating of Land Cover: Unsupervised method (BULC-U), the recent extension of the BULC algorithm to accommodate unsupervised classifications of Landsat satellite imagery [15], to interpret imagery across the time period 1986-2000, and to create an annual time series of major LULC categories of Forest, Deforested Land, and Shrubland/Grassland at 30 m spatial resolution. With extensive field experience in northeastern Mato Grosso, Brazil (Figure 1), we chose to study that region's noteworthy changes in land cover and land use patterns during the late 1980s and 1990s.

Study Area
The state of Mato Grosso in Brazil ( Figure 1) has undergone extensive deforestation in past decades, particularly for cattle ranching and soy cultivation [16][17][18]. To test the accuracy of our algorithm to produce the time series, we chose a study area in eastern Mato Grosso covering approximately 2 × 10 4 km 2 (166 km × 121 km) and encompassing four municipalities: Alto Boa Vista, Querência, Ribeirao Cascalheira, and Canarana, centered near 51.884°W, 12.601°S. The study area includes examples of the major tropical land covers: in the south and east, the native vegetation is predominantly savanna (a mixture of dense drought-deciduous forest and grasslands); in the north and west, the native vegetation is broadleaf evergreen forest, characteristic of the southern Amazon; the far eastern portion of the study area contains the Bananal seasonal grassland/wetland complex; and throughout are areas that are now deforested and converted to crops and pastures. To map the

Study Area
The state of Mato Grosso in Brazil ( Figure 1) has undergone extensive deforestation in past decades, particularly for cattle ranching and soy cultivation [16][17][18]. To test the accuracy of our algorithm to produce the time series, we chose a study area in eastern Mato Grosso covering approximately 2 × 10 4 km 2 (166 km × 121 km) and encompassing four municipalities: Alto Boa Vista, Querência, Ribeirao Cascalheira, and Canarana, centered near 51.884 • W, 12.601 • S. The study area includes examples of the major tropical land covers: in the south and east, the native vegetation is predominantly savanna (a mixture of dense drought-deciduous forest and grasslands); in the north and west, the native vegetation is broadleaf evergreen forest, characteristic of the southern Amazon; the far eastern portion of the study area contains the Bananal seasonal grassland/wetland complex; and throughout are Remote Sens. 2020, 12, 688 3 of 18 areas that are now deforested and converted to crops and pastures. To map the basic LULC changes of the region, we were interested in distinguishing the changes among three categories: (1) Forest; (2) Deforestation; and (3) Shrubland/Grassland mix.

BULC-U Algorithm
The BULC-U algorithm is an iterative algorithm developed in Google Earth Engine [19]. Adapted from the BULC algorithm [20,21], BULC-U has recently been applied to the GlobCover 2009 data set, in order to refine the spatial resolution of the 300 m GlobCover set to that of the Landsat imagery which it is fed [15]. BULC-U operates on unsupervised classifications from any source, and both BULC and BULC-U can run across very large regions in Earth Engine [14]. It was a relatively straightforward exercise to refine the GlobCover 2009 data set using 2009 Landsat data; because they were contemporaneous, there was substantial agreement between the data sets despite resolution differences. In this study, we asked if the same process can employ the GlobCover 2009 set to render a LULC time series both at a finer scale and two decades earlier.
To create the time series, we first ran BULC-U in reverse time order, meaning that the latest classifications (from the year 2000) were fed to BULC-U first. The remaining classifications were then processed by BULC-U by comparing each to the GlobCover data set, following Lee et al. (2018). This produced a relatively rough 1984 LULC classification that, with the exception of some unwanted noise, looked very much like a single classified Landsat image's classification from 1984; we saved that to the Earth Engine asset manager. After that spin-up, we ran BULC-U in the second, forward-running stage of the process, comparing each unsupervised LULC classification against the putative 1984 classification, creating a time series that stabilized by 1986, and ran to the end of 2000.

GlobCover Base Image
The research community has produced several planetary-scale LULC classifications using an accumulation of images near a nominal date of interest, e.g., [22,23]. One such planetary-scale data set is GlobCover 2009, a global classification made by interpreting imagery from the MERIS sensor with 300 m pixels, with an overall accuracy of 58% for 20+ LULC classes across the globe [24]. With its 2009 nominal date, the classification distinguishes among 22 possible classes of LULC. In this study area, the GlobCover 2009 s categories were dominated by a few classes: Rainfed Croplands (6.2%); Mosaic Cropland/Vegetation (14%); Mosaic Vegetation/Cropland (13.6%), Closed to Open Broadleaved Evergreen or Semi-deciduous Forest (32.4%); Closed Broadleaved Deciduous Forest (12.1%); Closed to Open Shrubland (17%); and mosaic classes. As outlined in Lee et al. [15], we remapped pixels of most mosaic types (particularly, the Mosaic Cropland/Vegetation type) to be NoData, meaning that their type was unknown a priori. As a result, part of BULC-U's function was to reveal the non-mosaic classes contained in the~100 30 m pixels contained within each 300 m GlobCover pixel. Additionally, as in Lee et al. [15], Cropland/Mosaic Vegetation mosaic was the usual label chosen in GlobCover for what could be seen in finer-scale Landsat images as being Cropland or Pasture, and so it was remapped to the Deforested Land class for forming the "base image", the best initial estimate of LULC in the study area. BULC-U would then use the spectral characteristics of these pixels to inform its creation of the series from the much earlier period.
GlobCover's coarse spatial resolution was especially strong at identifying large areas of intact Forest (e.g., north and east of the intersection of A2 and B3 in Figure 2), large expanses of Shrubland/Grassland (e.g., D5 in Figure 2), and areas where Deforested Land (in this region, loss of Forest for either cropland or pasture) dominated the land cover across large areas (C2 in Figure 2). GlobCover was less successful in capturing fine-scale LULC, however, and it had a 69.1% overall accuracy [15] within the study area.

Imagery
We identified a time series of 58 relatively clear (less than 10% cloud cover) Landsat 5 satellite images between 1984 and 2000, from between the dry season months of May and September (Table  1), for an average of about 3.5 images annually for each study pixel. In these images, Thematic Mapper bands four, five, and seven (near-IR, near-IR and mid-IR) were the clearest and used for the segmentation and classification described below.

Imagery
We identified a time series of 58 relatively clear (less than 10% cloud cover) Landsat 5 satellite images between 1984 and 2000, from between the dry season months of May and September (Table 1), for an average of about 3.5 images annually for each study pixel. In these images, Thematic Mapper bands four, five, and seven (near-IR, near-IR and mid-IR) were the clearest and used for the segmentation and classification described below.

Segmentation & Classification
Using the segmentation and classification tools in ArcGIS 10.3, we converted the 58 multi-band Landsat 5 satellite images into 58 unsupervised classifications, which we refer to as "Events", following the BULC and BULC-U work flow [15,20] (Figure 3).

Segmentation & Classification
Using the segmentation and classification tools in ArcGIS 10.3, we converted the 58 multi-band Landsat 5 satellite images into 58 unsupervised classifications, which we refer to as "Events", following the BULC and BULC-U work flow [15,20] (Figure 3). Using the Segment Mean Shift tool [25] in ArcGIS 10.3, we identified features in each image that were spectrally distinct and spatially contiguous, with each class being spectrally uniform with little speckling for each class. After identifying the segments, each band was re-computed using the mean value within the segment as the value for all pixels within the segment. Next, we performed unsupervised classification of each segmented image using the ISODATA [26] unsupervised classification tool. Through trial and error, we identified that using 20 unsupervised classes produced spectrally and spatially distinct classifications. The number of pixels in each class was large enough that BULC-U could relate them to the coverage in the GlobCover 2009 image, but small enough that each unsupervised class tended to be dominated by just one of the LULC classes tracked through the study. Using the Segment Mean Shift tool [25] in ArcGIS 10.3, we identified features in each image that were spectrally distinct and spatially contiguous, with each class being spectrally uniform with little speckling for each class. After identifying the segments, each band was re-computed using the mean value within the segment as the value for all pixels within the segment. Next, we performed unsupervised classification of each segmented image using the ISODATA [26] unsupervised classification tool. Through trial and error, we identified that using 20 unsupervised classes produced spectrally and spatially distinct classifications. The number of pixels in each class was large enough that BULC-U could relate them to the coverage in the GlobCover 2009 image, but small enough that Remote Sens. 2020, 12, 688 6 of 18 each unsupervised class tended to be dominated by just one of the LULC classes tracked through the study.

Validation
To understand BULC-U's ability to build from the GlobCover 2009 classification to track land cover change and stability in different decades, we assessed the classifications in the three specified categories, which could be easily observed in imagery with relatively little evaluator confusion. For estimation of area and map accuracy, the standard practices derived by Stehman [27] and Olofsson et al. [28] were used. This included creating a probability sample and calculating the inclusion probability of each point that was selected for inspection as part of the validation. We stratified a collection of reference points using the BULC-U map for 1986: 53 points from the Deforested Land stratum, 71 from the Forest stratum, and 100 from the Shrubland/Grassland stratum. The accuracy of the BULC-U time series was assessed in detail for the years at either end of the time series: 1986 and 2000. An evaluator unassociated with the project evaluated each point's reference category using a clear Landsat image from 1986 and from 2000. We used the equations of Stehman [27] to calculate unbiased estimates of the overall accuracy and each LULC category's area, user's accuracy, and producer's accuracy, along with estimates of error in these quantities.

BULC-U Classification Sequence
The 1986-2000 BULC-U land-cover time series indicates a landscape in substantial transition, from one dominated by Forest and Shrubland/Grassland to a more even mix among Deforested Land, Forest, and Shrubland/Grassland (Figures 4-6). BULC-U produced a time series of 58 land-cover classifications, dating from 1986 ( Figure 4) to 2000 ( Figure 5). At either end of the sequence, BULC-U successfully captured large changes in the landscape, such as the establishment of contiguous fields in sectors D1 and E1, A1, and A2. It was also able to capture finer-scale deforestation in sectors B2 and C2, as well as in A4 and A5. Much of the Forest category, particularly in 1986, occurred in large, spectrally steady patches that were easy to capture with the unsupervised classifications that feed BULC-U. When visually inspecting the results at either end of the sequence, it was rare to find errors in any of the categories in which, for example, an area recognizable as having been deforested had been missed by BULC-U and incorrectly labeled as Shrubland/Grassland or Forest.
The total percentage of the Deforested Land class increased substantially during the study period across the study area, for an overall increase of 3400 km 2 (13% to 32% in the unbiased estimator of Stehman 2014) and came principally from the Forest class. The Forest class experienced a nearly equivalent decline (55.0% to 38%). The unbiased estimate of the Shrubland/Grassland category was effectively steady throughout the study: from 32.1% to 30.2% of the area. Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 18

Time Series Accuracy
Following the protocol of Stehman [27] for unbiased estimation of accuracy and area, the time series starting point (1986) and ending point (2000) had an overall accuracy of 89% and 80%, respectively (Figure 7). This marked a substantial improvement over the GlobCover accuracy of below 70%, even in the same simple categories (Figure 2). Some of this improvement is likely due to

Time Series Accuracy
Following the protocol of Stehman [27] for unbiased estimation of accuracy and area, the time series starting point (1986) and ending point (2000) had an overall accuracy of 89% and 80%, respectively ( Figure 7). This marked a substantial improvement over the GlobCover accuracy of below 70%, even in the same simple categories (Figure 2). Some of this improvement is likely due to an improved ability to resolve LULC types at the finer 30 m resolution, given the inherent challenge of the relatively coarse 300 m dataset to represent fine-scale patterns. More importantly, the BULC-U process labeled about 100 times more points (as the spatial resolution was increased tenfold) more accurately in an earlier, lesser-mapped time period.
the misclassified true-deforestation points for close inspection, two factors mitigate the Deforested Land accuracy values that were lower than expected. First, the reference protocol for labelling Deforested Land was relatively vague and did not encourage the evaluator to look at multiple dates to detect deforestation. There may be, at least, some true Deforested Land pixels that were mislabeled by the evaluator as Shrubland/Grassland. Second, around half of the misclassifications in the 1986 map were, unluckily, at the edge of land clearly being actively managed, open land whose use was difficult to discern. The evaluator was not permitted to skip edge points in lieu of finding another point. As a result, though the producer's accuracy values would suggest that a large amount of deforestation was missed on both dates, we were unconvinced; instead, it appeared to us that the accuracy values imperfectly summarized that category. Coupled with the relatively high estimated error around the values for the Deforested Land and Shrubland/Grassland classes, this may be symptomatic of the known challenges of accurately estimating uncertainties in relatively small classes [29].

Fine-Scale LULC Time Series
Although the proportion of Deforested Land increased roughly linearly through time across the entire study area (Figure 6), a finer-scale tracking of land-cover change shows that this linear increase was composed of differing rates of conversion to Deforested Land in time and space. The area that later formed the town of Querência (Figures 5 and 6, near the edge of sectors B2 and C2) expanded considerably through the study period, and was captured well by BULC-U with respect to contemporaneous Landsat imagery (Figure 8)   Given the apparent high visual agreement between the BULC-U classifications and evidence from Landsat, it was surprising to see several low per-class user's and producer's accuracy values indicated for the 1986 and 2000 classifications (Figure 7). For the Forest class, there were few omission errors in either date, and commission errors in the class were mostly when true Shrubland/Grassland was misclassified as Forest. Shrubland/Grassland per-class accuracies were more variable, with true Shrubland/Grassland much more likely to be misclassified as Forest by BULC-U than as the Deforested Land class. For pixels mislabeled as Shrubland/Grassland on BULC-U maps, they were much more likely to be true Deforested Land than Forest on both dates. The low producer's accuracies of Deforested Land were surprising and, at first, concerning that the model failed to capture large amounts of true forest loss in the region. On both dates, true deforestation was very rarely misclassified as Forest; instead, where a Deforested Land pixel was mislabeled in BULC-U, it was nearly certain to have been mislabeled as Shrubland/Grassland. On re-inspecting the Landsat imagery from 1986 and 2000, it was difficult to find locations that we would indicate on Landsat imagery as deforested land but that were classified as Shrubland/Grassland or Forest. After isolating the misclassified true-deforestation points for close inspection, two factors mitigate the Deforested Land accuracy values that were lower than expected. First, the reference protocol for labelling Deforested Land was relatively vague and did not encourage the evaluator to look at multiple dates to detect deforestation. There may be, at least, some true Deforested Land pixels that were mislabeled by the evaluator as Shrubland/Grassland. Second, around half of the misclassifications in the 1986 map were, unluckily, at the edge of land clearly being actively managed, open land whose use was difficult to discern. The evaluator was not permitted to skip edge points in lieu of finding another point. As a result, though the producer's accuracy values would suggest that a large amount of deforestation was missed on both dates, we were unconvinced; instead, it appeared to us that the accuracy values imperfectly summarized that category. Coupled with the relatively high estimated error around the values for the Deforested Land and Shrubland/Grassland classes, this may be symptomatic of the known challenges of accurately estimating uncertainties in relatively small classes [29].

Fine-Scale LULC Time Series
Although the proportion of Deforested Land increased roughly linearly through time across the entire study area (Figure 6), a finer-scale tracking of land-cover change shows that this linear Remote Sens. 2020, 12, 688 11 of 18 increase was composed of differing rates of conversion to Deforested Land in time and space. The area that later formed the town of Querência (Figures 5 and 6, near the edge of sectors B2 and C2) expanded considerably through the study period, and was captured well by BULC-U with respect to contemporaneous Landsat imagery (Figure 8)  The area that became Querência began 1986 almost entirely covered by Forest, and the amount of Deforested Land doubled from 5.2% of the surrounding area (half the regional average at that time) to 10.5% between 1986 and 1988 ( Figure 9). Agricultural expansion continued, but slowed until 1995 (Figure 9), after which it accelerated substantially between 1996 and 1998 (Figure 10), ending the period having 32.6% Deforested Land, a value higher than the regional average and more than six times the 1986 proportion. The area that became Querência began 1986 almost entirely covered by Forest, and the amount of Deforested Land doubled from 5.2% of the surrounding area (half the regional average at that time) to 10.5% between 1986 and 1988 ( Figure 9). Agricultural expansion continued, but slowed until 1995 (Figure 9), after which it accelerated substantially between 1996 and 1998 ( Figure 10), ending the period having 32.6% Deforested Land, a value higher than the regional average and more than six times the 1986 proportion.
Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 18 Figure 9. Land cover in Querência during early part of the study period: 1986-1990.
The growth of Deforested Land in the town of Querência occurred around a center point, especially during the early boom ( Figure 9). The majority of the growth was in smaller fields which were created sporadically around the center, notably including a large field deforested in 1988. Large Forest patches still occupied a majority of the area in 1990 but the continuous Forest was starting to become isolated in patches due to ongoing Deforested Land expansion. By 1996, large patches of Forest had become even more isolated ( Figure 10) and smaller patches had become isolated and fragmented. By 1998 (Figure 10), many remaining small forest patches had disappeared. The pattern continued to 2000, as Deforested Land further encroached on Forest.

Discussion
The BULC-U algorithm fused Landsat data with the GlobCover 2009 land cover classification, a moderately accurate (~50%) single year land cover classification, to create a high-accuracy time series classification for a two-decade period. Using automatically created unsupervised classifications of imagery from the Landsat 5 satellite, the BULC-U algorithm tracked changes and updated the running classification as new satellite images were inputted, creating a time series of classifications that retained good accuracy at the year 2000 end date. The growth of Deforested Land in the town of Querência occurred around a center point, especially during the early boom ( Figure 9). The majority of the growth was in smaller fields which were created sporadically around the center, notably including a large field deforested in 1988. Large Forest patches still occupied a majority of the area in 1990 but the continuous Forest was starting to become isolated in patches due to ongoing Deforested Land expansion. By 1996, large patches of Forest had become even more isolated ( Figure 10) and smaller patches had become isolated and fragmented. By 1998 ( Figure 10), many remaining small forest patches had disappeared. The pattern continued to 2000, as Deforested Land further encroached on Forest.

Discussion
The BULC-U algorithm fused Landsat data with the GlobCover 2009 land cover classification, a moderately accurate (~50%) single year land cover classification, to create a high-accuracy time series classification for a two-decade period. Using automatically created unsupervised classifications of imagery from the Landsat 5 satellite, the BULC-U algorithm tracked changes and updated the running classification as new satellite images were inputted, creating a time series of classifications that retained good accuracy at the year 2000 end date.

Deforestation Expansion
The GlobCover 2009 product has provided impressive, much-needed information on land use at a snapshot in time, but what is missing from such static LULC maps is a fine-scale record of the history of change. Although Landsat data has been available for many years, constructing high resolution time-series with legends tailored to a particular use for anything but small areas has been too difficult because of the expertise required. This time series in Mato Grosso revealed a steady increase in the conversion of natural landscapes between 1986 and 2000, with roughly 13% of the study area in Deforested Land category in 1986 and nearly 32% in 2000. Almost all of that increase came at the expense of the Forest category. Since the mid-1970s, Brazil's Amazon-Cerrado agricultural frontier has been among the planet's largest and most important zones of deforestation, agricultural expansion, and agricultural intensification. Through policy, applied research, and infrastructure development, Brazil has transformed itself from a food importer to one of the world's agricultural powerhouses [18,30]. Amidst this growth in agricultural output, more than 750,000 km 2 of land was deforested in the Brazilian Amazon alone [31].
Class proportions in the BULC-U changed in a coherent manner throughout the time series and did not flip unstably or unpredictably between states during the study period ( Figure 6). The realistic trajectories of the classes and tracking with Landsat imagery suggest that the time series of BULC-U classifications was largely accurate in its progression of land cover change. The endpoints of the time series (1986 and 2000) were validated and have high fidelity with the spatial patterns apparent on contemporaneous Landsat images (Figures 4 and 5). Together, this suggests that, while it was impractical to individually validate each of the 58 classifications independently for this study, the intermediate steps plausibly have similar accuracy compared to each of the end dates. This is due, in part, to the mechanics of BULC-U: in tracking a vector of class probabilities in each pixel, the classifications are buffered by design against inconsistencies in a given Event classification, with transient errors or unusual images less likely to cause dramatic shifts in the BULC-U classification.
Although end-of-year classifications are shown here for simplicity, BULC-U produces complete classifications for each of the 58 irregularly spaced time steps, updating the estimates with the introduction of each new piece of evidence from Landsat. Because the classifications were at a sub-annual time step, the dynamics of deforestation and change occurring in land cover can be quantified at an even finer temporal granularity (about semi-annual) than that shown here. The BULC-U algorithm was able to track deforestation events in the study area soon after they occurred, with forests cleared for new fields and the expansion of existing fields (e.g., Figure 9). This type of time series information can be valuable, providing, for example, knowledge of when a municipality or state crossed a given threshold of deforestation, which may have implications for the local climate [32] or for regional compliance with Brazilian Forest Code legislation. Furthermore, knowing the time since deforestation occurred in a location or region can help in understanding the subtler ecological changes that may be occurring, such as longer-term CO 2 emissions from soils in pastures and croplands [33,34] or biodiversity in a fragmented landscape [35].

Finer Scale Land Cover Changes
The BULC-U time series reveals locations where, at finer spatial scales, LULC change occurred at varying rates through the study period, even though the full study area's rate of conversion to agriculture was roughly linear through time (Figures 6 and 9). In the town of Querência, the history of agricultural expansion contributed to-but can be seen here to be distinct from-the regional expansion pattern. Conversion rates were sporadic, probably in response to large and local scale social and economic drivers, with rapid conversion occurring between 1987 and 1988 and between 1996 and 1998, and slower conversion rates between 1991 and 1995. Knowing how local and regional deforestation and conversion rates changed may help researchers to better understand the social and economic drivers of deforestation-e.g., is the price of commodities or legislation in a particular year or location reflected in subsequent deforestation rates with a time lag [30,36,37]? The high spatial and temporal resolution of the BULC-U time series make such analyses more feasible.

Forest Fragmentation
Forest fragmentation decreases biodiversity, disrupts ecosystem function, and increases forest degradation and carbon emissions [38][39][40]. The BULC-U time series shows the short-term, fine-scale patterns of deforestation and fragmentation. In areas like Querência, forest remained the dominant land cover by the end of the study in 2000, but the forest had become increasingly fragmented, especially during the latter half of the period. Larger patches were more fragmented and smaller patches in the area disappeared with the continuing expansion of agriculture ( Figure 10). This ability to identify and date patterns of forest fragmentation and agricultural expansion can be a valuable tool for conservation management, biodiversity research, land use management, and policy decisions. In possible future work, techniques such as network theory and circuit theory could be applied to the time series classification created by BULC-U to quantify decreases in connectivity and fragmentation at each time step, thereby creating greater conservation value.

Earth Engine/Generality
This study focused on the categories of Deforested Land, Forest, and Shrubland/Grassland through time, but the BULC-U algorithm's application to create time series is not limited to these specific LULC classes. Rather, users can create time series of any phenomena that can be reliably mapped from imagery or, more generally, mapped by any method with periodic updates. This means, for example, that if the changes in a specific grassland were of interest and reliably (though not always perfectly) detectable, a user of this approach could trace change and stability of that class, along with other classes of interest, through time. This generality opens possibilities for users to create LULC classifications targeted to their particular needs and could be relevant for those interested in changes beyond the more familiar forest/non-forest distinctions.
Earth Engine facilitated this research in two specific ways. First, as the pre-eminent cloud computing platform for satellite imagery, it permitted the analyses to be run with ease from any computer. Recent advances in Earth Engine's web presence mean that BULC-U could potentially be offered for exploration to a wider audience of non-expert stakeholders. Second, Earth Engine's speed permitted very rapid prototyping, in which one could quickly view the outcomes of different modeling scenarios (for example, to explore the outcome of adding cloudier images to the mix to increase the temporal resolution at the potential cost of decreased accuracy). Additionally, the mid-2018 addition of segmentation algorithms, like those used here, mean that this entire workflow could conceivably be done at the same time in a single browsing session [41,42].

Conclusions
The BULC-U technique applied here represents a novel way to accurately create sub-annual time-series land cover classifications and quantify land cover change in past time periods. Because the algorithm operates on unsupervised classifications, it permits users to map land cover change using only a base map and the statistical characteristics of unclassified satellite imagery, with little manual user involvement to classify individual events. The high-resolution time series provides more information than would otherwise be possible in coarser products with discrete large time jumps. Our application of the technique to eastern Mato Grosso illustrates the potential value of time series in land cover mapping. Eastern Mato Grosso, on average, had a large and steady decline in forest cover from 1986. But the higher spatial resolution time series produced by BULC-U clearly showed the regional and temporal nuances-rapid change took place at fine scales in individual years, springing up around growing towns, with multi-year slowdowns in conversion before another rapid growth phase. This fine-scale depiction of change, when combined with socio-economic data, can be valuable in the future for understanding how macro-scale drivers (e.g., fluctuating commodities prices, environmental legislation, or market forces) influence local land use decisions.