A Suite of Tools for Continuous Land Change Monitoring in Google Earth Engine

Land cover has been designated by the Global Climate Observing System (GCOS) as an Essential Climate Variable due to its integral role in many climate and environmental processes. Land cover and change affect regional precipitation patterns, surface energy balance, the carbon cycle and biodiversity. Accurate information on land cover and change is essential for climate change mitigation programs such as UN-REDD+. Still, uncertainties related to land change are large, in part due to the use of traditional land cover and change mapping techniques that use one or a few remotely sensed images, preventing a comprehensive analysis of ecosystem change processes. The opening of the Landsat archive and the initiation of the Copernicus Program have enabled analyses based on time series data, allowing the scientific community to explore global land cover dynamics in ways that were previously limited by data availability. One such method is the Continuous Change Detection and Classification algorithm (CCDC), which uses all available Landsat data to model temporal-spectral features that include seasonality, trends, and spectral variability. Until recently, the CCDC algorithm was restricted to academic environments due to computational requirements and complexity, preventing its use by local practitioners. The situation has changed with the recent implementation of CCDC in the Google Earth Engine, which enables analyses at global scales. What is still missing are tools that allow users to explore, analyze and process CCDC outputs in a simplified way. In this paper, we present a suite of free tools that facilitate interaction with CCDC outputs, including: (1) time series viewers of CCDC-generated time segments; (2) a spatial data viewer to explore CCDC model coefficients and derivatives, and visualize change information; (3) tools to create land cover and land cover change maps from CCDC outputs; (4) a tool for unbiased area estimation of key climate-related variables like deforestation extent; and (5) an API for accessing the functionality underlying these tools. We illustrate the usage of these tools at different locations with examples that explore Landsat time series and CCDC coefficients, and a land cover change mapping example in the Southeastern USA that includes area and accuracy estimates.


INTRODUCTION
Remote sensing of the land surface can contribute to our understanding of the climate system in a variety of ways. The land surface represents the boundary conditions for atmospheric circulation models, providing essential information on exchanges of mass, energy and momentum between the surface and the atmosphere (Sellers et al., 1997). The Global Climate Observing System (GCOS) has identified Essential Climate Variables (ECVs), which include a suite of land surface products related to the cryosphere, biosphere, and hydrosphere that are "needed to understand and predict the evolution of climate, to guide mitigation and adaptation measures, to assess risks and enable attribution of climate events to underlying causes, and to underpin climate services" (WMO, 2020). In this paper, we focus on one of these Climate ECVs, land cover, with particular attention to change in land cover as it has been a primary driver of increases in atmospheric CO 2 . Historically, more than half of the CO 2 added to the atmosphere by humans can be attributed to land use and land cover change, which releases large amounts of carbon, often from burning of biomass (Houghton and Nassikas, 2017). Land use and land cover change remain as important sources of greenhouse gases, and for this reason the REDD+ program was created. By the time the REDD+ program started, ∼12% of the annual carbon emissions were due to land use practices, including deforestation and forest degradation (Houghton et al., 2012;Goetz et al., 2015). The REDD+ Program, an abbreviation of Reducing emissions from deforestation and forest degradation and the role of conservation, sustainable management of forests and enhancement of forest carbon stocks in developing countries, was designed to reduce greenhouse gas emissions from the land sector by providing economic incentives. The program stipulates that developing countries will receive economic compensation for measuring, verifying, and reporting on reductions of greenhouse gas emissions from one or more of the five REDD+ activities: (1) deforestation, (2) forest degradation; (3) conservation, (4) sustainable management, and (5) enhancement of forest carbon stocks (GFOI, 2016). Thus, land cover change monitoring and reporting by participating countries is at the heart of the REDD+ Program.
Monitoring land surface change has been one of the primary uses of remote sensing over the past several decades. The most common approach has been to compare images from before and after the changes of interest to map the location, extent and kinds of change. This approach can be thought of as using satellite images as "endpoints" to monitor change. There are many examples in the literature, most of which have used data from the Landsat Program (Coppin and Bauer, 1996;Hansen and Loveland, 2012;Zhu, 2017). The Landsat Program has been collecting global data since 1972, and since 1982 a series of four satellites have collected data at 30 m spatial resolution, providing an archive of imagery valuable for monitoring land change over decades. Many other sensors have also been used for monitoring land change. The MODIS sensors on Terra and Aqua have been frequently used for monitoring change, particularly fires (Giglio et al., 2018). However, the relatively coarse spatial resolution of MODIS (250 m for the NIR and red bands and 500 m for the other "land bands") has limited its utility for monitoring land use change, which often occurs over much smaller areas. More recently, the Copernicus Program has launched the Sentinel-1 (2014 and 2016) and Sentinel-2 (2016) series of satellites, which have begun collecting data suitable for monitoring land change. Sentinel-1 is a C-Band Synthetic Aperture Radar that has four operating modes, which collect data at spatial resolutions that range from 5 to 25 m. The Sentinel-2 satellites collect visible and near infrared data at 10 m spatial resolution and shortwave infrared and several specialized bands (e.g., red edge) at 20 m. There are also spectral bands at 60 m spatial resolution, that are primarily used for atmospheric characterization.
Since the advent of free data from the Landsat Program in 2008, new methodologies for mapping land cover change based on analysis of many observations (or time series) have emerged. Rather than rely on a pair of observations before and after to detect a change, algorithms can now be used to monitor land change in a more continuous fashion . The first time-series approaches used a single observation for each year, enabling the detection of changes on an annual basis (Huang et al., 2010;Kennedy et al., 2010;Hughes et al., 2017). These algorithms were usually based on a single spectral index and were designed to find specific types of change, often change in forests. Other methods that used MODIS data (usually NDVI, the Normalized Difference Vegetation Index) composited to 16 or 32 day time periods to look for land change were developed simultaneously (Verbesselt et al., 2010). Over time, more emphasis has been placed at Landsat resolutions on using more than annual observations to better understand the timing and dynamics of land change. One approach that attempts to use all good quality observations is CCDC [Continuous Change Detection and Classification (Zhu and Woodcock, 2014)], which is intended to be a general-purpose algorithm for finding many kinds of changes in a continuous manner and mapping the land cover before and after. Parallel efforts have attempted to study and characterize land surface phenology using the information contained in time series of remote sensing observations directly, rather than deriving land cover data from them. Some examples include the study of trends in spring onset and fall phenology (Schwartz et al., 2006;Dragoni and Rahman, 2012), crop yield forecasting (Bolton and Friedl, 2013), global phenology products and intercomparisons (Ganguly et al., 2010;Moon et al., 2019), among other applications (Morisette et al., 2009;Henebry and Beurs, 2013). Land surface phenology is an important but extensive research topic that is beyond the scope of this paper. For this reason, we will focus exclusively on land cover and land cover change mapping.
The relationship between REDD+ and Landsat is essential because Landsat is frequently the most viable data option for REDD+. For example, countries are required to establish a Forest Reference Level (FREL) against which more recent emissions can be compared, often requiring multiple decades of images to contribute to the establishment of FRELs and monitoring thereafter. For many of those countries, using Landsat data is the most plausible option due to its relatively high resolution, availability since the 1980's for many regions in the world, and consistent pre-processing and archival in collections that allow the delivery of high-quality surface reflectance products (Olander et al., 2008;USGS, 2018). Additionally, the Landsat Program has improved the processing and consistency of their datasets by including atmospheric correction, cloud and cloud shadow masking and improved multitemporal registration, while providing the data at no cost to the user (Zhu and Woodcock, 2014;Vermote et al., 2016). Furthermore, the development of algorithms to harmonize Landsat and Sentinel-2 surface reflectance products is expected to augment these collections with data of similar characteristics (Claverie et al., 2018). However, providing the necessary satellite data is just the first step -subsequent steps include determining the best use for the data to monitor REDD+ activities in ways compliant with the Intergovernmental Panel on Climate Change criteria and also useful for decision-and policy-making.
To achieve compliance, effective and comprehensive monitoring in countries where it is needed the most, capacity building efforts have been sponsored and organized by various organizations including SilvaCarbon, USAID, UN-FAO, NASA-SERVIR, and GOFC-GOLD. As direct and active participants in these efforts, we have witnessed a dramatic reduction in the technical steps required to obtain end products of land cover change of adequate quality. Initial training and capacity building workflows in support of REDD+ that involved teaching participants how to register images, remove clouds and cloud shadows from datasets, and correct images for atmospheric effects have been replaced by streamlined applications or preprocessing done by the Landsat Program. These developments, combined with the growing adoption of the Google Earth Engine (GEE; Gorelick et al., 2017), have alleviated much of the preprocessing and allowed the community to instead focus on increasing precision, accuracy and comprehensiveness of remote sensing studies.
A key advantage of GEE is the close connection between data and the algorithms, both of which can be accessed through an application programming interface (API). With the exception of BFAST (Verbesselt et al., 2010), all of the temporal segmentation algorithms mentioned previously have been natively implemented in GEE (But a recent user implementation of BFAST monitor in GEE has been made available by Appel, 2020). The CCDC algorithm, originally developed in MATLAB and later implemented as a Python package by USGS-EROS as part of the LCMAP initiative (Brown et al., 2020), has been recently implemented as a temporal segmentation algorithm in GEE by Noel Gorelick, Zhiqiang Yang, and the Earth Engine Team and will be formally described by their authors in an upcoming publication (Google, 2020). Therefore, the CCDC description included here is intended only as a quick overview of the algorithm to aid in the understanding of the concepts mentioned in this paper but does not constitute the official description of CCDC in GEE.
CCDC is an "online" time series algorithm, meaning it starts at the beginning of the time series and is updated as new observations become available. Thus, at the beginning, once there are sufficient observations a time series model is fit. If new observations added to the time series do not diverge widely from the model, the new observations are added to the time series and an updated model is fit. When the new observations no longer follow the expected values from the existing model, a "break" is inserted. The time between breaks (or the beginning of the time series and the first break) define "time segments." The breaks indicate that the spectral/temporal behavior of a pixel is different than in the past, presumably due to a change in land cover or a change in condition. A break is detected when the residuals, normalized by the model root-mean-square error (RMSE), exceed a threshold for a number of consecutive observations. Parameters for defining the critical region and the number of observations to confirm a change can be tuned by the user to control the sensitivity of the change detection. The RMSE for a model is an indication of the typical variability of observations relative to that model. Details of the inner workings of CCDC are beyond the scope of this paper and can be found in Zhu and Woodcock (2014).
The CCDC algorithm is well-suited for the needs of REDD+, as it allows for continuous monitoring of change, enabling analysis of trends in land cover and even land use through time, such as changes in rates of deforestation . The time segments fit by the algorithm can be classified into land cover classes that are consistent over time, making them particularly useful for understanding land cover change. Additionally, the time segment coefficients can be used on their own to study land surface phenomena that exhibit variable trends and seasonality (Pasquarella et al., 2018). While the availability of the algorithm in GEE has removed the significant computational requirements that restricted the use of CCDC to academic environments, improving the ability of users to interact with the outputs of the algorithm is essential for a broader use of the algorithm, particularly by countries hoping to participate in REDD+.
The goal of this paper is to facilitate the use of CCDC even more by providing tools that will help users (countries in the case of REDD+) to get the benefits of processing the entire archive of Landsat data for studying land cover and change. We present a series of tools and educational materials designed to support all steps from training data to maps and area/accuracy estimates, without the need to preprocess data, write complex code or construct statistical estimators. The tools presented in this paper are designed to interact with the outputs of the Continuous Change Detection and Classification algorithm (CCDC) as implemented in GEE. We believe that the tools we have developed in the GEE platform will greatly facilitate the use of remote sensing data in support of climate change mitigation by enabling a comprehensive analysis of change on the surface of Earth while following the best practices in the field.

METHODS
The tools and code repository presented here (Arévalo and Bullock, 2020) allow the users to interact with time series of Landsat observations and the outputs of the GEE implementation of the CCDC algorithm, as well as to create products from them. The tools were implemented as Earth Engine Apps when possible and are accompanied by an API built atop the GEE JavaScript API. The tools to display and interact with time series and CCDC results for a pixel were inspired by the legacy of the TSTools plugin for QGIS2 (Holden, 2017). Documentation for using the tools and API with access to a sample dataset, as well as up to date links to access them and the source code can be found at https://gee-ccdc-tools.readthedocs.io. A curated copy of the source code suitable for suggestions, comments and contributions is hosted at https://github.com/parevalo/gee-ccdctools. Exploratory tools that make use of the Earth Engine Apps platform can be used directly in an internet browser without requiring a GEE account. In order to use other tools, and access the API and its source code, users need to have an active GEE account, which can be created here: https://earthengine.google. com/signup/.

CCDC Algorithm
The current description of CCDC is intended to provide a general understanding of the algorithm, and primarily serves the purpose of helping users understand the tools presented in this paper. In this regard, the concepts of time segments and breaks are fundamental, as well as the coefficients for the various elements in the harmonic regression models: an intercept (similar to the mean); a slope (indicative of a trend); coefficients for the sines and cosines of the 12, 6, and 3-months harmonics. Figure 1 illustrates the use of multiple harmonics to fit a model to a time series of observations with no breaks. In turn, Figure 2 shows an example of CCDC results for a single pixel where multiple breaks have been detected. All the examples presented in this paper have been obtained from running the CCDC algorithm using these Landsat bands as inputs: blue, green, red, Near Infrared (NIR), Shortwave Infrared 1 and 2 (SWIR1, SWIR2), and thermal.
Another construct made possible through the use of CCDC results is a "synthetic image." Once the harmonic regression models (i.e., time segments) have been fit by the CCDC algorithm, we can use them to estimate a value for a pixel at any time within the input data time range. This capability can be used for the creation of a "synthetic" image, which is simply the predicted value for all pixels for a given date (Zhu et al., 2015). For example, in Figure 2, an observation for a single pixel from a real Landsat image is represented by any of the blue dots, whereas the synthetic value for the same pixel would correspond to any value along the harmonic models, or their predicted value for the gaps between models. The use of synthetic images is discussed below in Section spatial data viewer.
The GEE implementation of the CCDC algorithm is data agnostic, meaning that in principle it can be applied to other multitemporal datasets other than Landsat, such as Sentinel-1 and -2. However, the algorithm was originally created to use Landsat observations, and most of the published research using CCDC has been applied to Landsat observations. For this reason, the applications shown here are demonstrated using Landsatderived results and allow the user to interact with Collection 1 Landsat imagery and the products derived from the coefficients for the time segments fit by CCDC.

Graphical Interface Tools
The graphical interface tools we present here can be categorized in three groups: 1) Visualization and exploration tools. 2) Classification and mapping tools. 3) Area estimation and accuracy tools.
The tools were designed as standalone applications that can be used independently from each other, except for the Land cover tool, which requires the output of the Classification tool. Still, certain processes can be streamlined if all the apps are used together. For example, a user may be interested in generating a land cover map for a geographical area they are not familiar with. In this scenario the visualization tools could be used to explore the spatio-temporal land dynamics in that area. Once those are better understood, the user may collect training data and use it in the classification and mapping tool to produce maps of land cover change. Finally, the user may need to quantify the accuracy of such maps, or estimate areas of change with quantifiable uncertainty, for which they may use the last set of tools. However, considering that users may have different purposes for the tools that are not exclusively related to land cover mapping, only the core usage of each tool is described below, and independent case examples are presented in the results section.

Time Series Viewers
The time series viewers are simple interactive Earth Engine Apps that allow the user to explore series of Landsat observations for any pixel on the globe where such data are available. Cloudy or shadowed observations are filtered from the time series using the quality data flags contained in the Landsat collection images. The user can select the spectral band and time range to visualize and click on any point in the time series to load the corresponding image into the map viewer. The tools also show the time segments detected on-the-fly by the CCDC algorithm. Two versions of the tool are available. The first is the simple viewer (Figure 3), which can be found in this link: https://parevalo_bu.users.earthengine. app/view/quick-tstools. In this version, the options that the user can modify are minimal (i.e., time range and spectral band to visualize), as it was designed for use in mobile devices to facilitate the analysis of the landscape history directly in the field. The tool works on any mobile or desktop browser that is supported by Earth Engine Apps and runs CCDC using the default parameters. The second version of the tool is the "advanced" time series viewer (Figure 4), created for users who want to have more control on the parameters used for running CCDC and visualization. In addition to the capabilities described for the simple viewer, it is possible to modify the parameters to run the CCDC algorithm. The image visualization parameters or the range of spectral indices may also be modified. The tool also allows adding ancillary data stored as public assets to the map. The tool can be found here: https://parevalo_bu.users. earthengine.app/view/advanced-tstools.

Spatial Data Viewer
The spatial data viewer tool ( Figure 5) expands the capabilities of the time series viewers by allowing the user to spatially explore Frontiers in Climate | www.frontiersin.org FIGURE 1 | Landsat time series (dots) for the SWIR1 band, and the corresponding time segments detected by the CCDC algorithm, using different numbers of harmonics. (A) Time segment using only one harmonic (annual harmonic), (B) two harmonics (annual and six-month harmonics), and (C) three harmonics (annual, 6-, and 3-months harmonics). The pixel corresponds to a stable forest in the Blue Hills Reservation, in the commonwealth of Massachusetts, USA. Only a single spectral band is shown here, but CCDC fits harmonic models to all bands included in the analysis, which are specified by the user.
the coefficients of the time segments and the detected changes, as well as creating various Landsat derivatives. To use the core functionality of the tool, the user needs to specify the path to existing CCDC coefficients saved as GEE assets, either as the Earth Engine image that is created as the algorithm output, or as an Image Collection that contains one or more of those images. Once loaded, the tool enables the exploration of the following: 1) Model coefficients: as explained in Section CCDC algorithm, each time segment is defined by a set of coefficients: the intercept; the slope; and the harmonic terms of the regression model. These coefficients, as well as the RMSE, and the amplitude and phase of the harmonics, can be visualized using this tool. The user can add individual coefficients or color composites of coefficients to the map. Model coefficients, and quantities derived from them like the phase and amplitude of the harmonic models, can be used in a variety of ways. The most common use is in the classification step of CCDC where time segments are used as an input to the classification process. The classification process and the tool for its implementation are discussed below in Sections classification and spatial data viewer. However, coefficients and their derivatives have intrinsic value. For example, the amplitude of the harmonic regression is indicative of the magnitude of the changes in spectral bands seasonally, which could be related to deciduousness. 2) Synthetic images: the time segments can be used to predict the surface reflectance for any point in time within the time range of the input data. Imagery that is composed of predicted surface reflectance is referred to as "synthetic" imagery (Zhu et al., 2015), and can be generated using a color combination specified by the end user. Synthetic imagery can be generated for any day within the time range of the input data even in the absence of Landsat observations. These data represent a new kind of "derived" data product, as they are based on past Landsat observations, but are not themselves observations. They can be valuable and easy-to-use images as, by definition, they are free of the effects of clouds, cloud shadows, and snow. Additionally, they can capture the seasonality of locations when produced for multiple times of year, which can improve applications like image classification (Pasquarella et al., 2018).
3) Change information: for each pixel there is a variety of change information available. If changes ("breaks") have been detected by the CCDC algorithm, there is information about the timing of the break, the magnitude of the break and the number of breaks. The magnitude of the break is calculated internally by the CCDC algorithm as the difference between the mean of the last five observations before the break and the mean of the first five observations after the break in any spectral band or index. The user can display the number of changes for a selected time period, the change with the largest magnitude (e.g., deforestation) and the time of the most recent or largest change.

Classification
For classification, we provide two tools that can be used together for land cover and/or land use classification using the CCDC model parameters as predictor variables. The tools allow the creation of categorical maps, and the testing of input parameters without writing any code. The first tool, referred to here as the Classification Tool, allows users to classify the time segments created using the CCDC function in GEE. The user provides the path to the CCDC results (stored as Earth Engine Assets), training data (stored as an Earth Engine Feature Collection), and  selects the machine learning classifier, predictor variables and study area to use. The user can select a country boundary to use as the study area or can draw the study area on the map. A list of the available predictor variables can be found in Table 1. The result is an N-band image corresponding to N time series segments. Since CCDC operates on the time series of every individual pixel, the start and end dates for the segments vary by pixel accordingly. Therefore, to extract a categorical map for a specific date, or change between multiple dates, a separate tool is required. The second tool, referred to here as the Land Cover Tool, allows users to create land cover or land use maps from the classified segments produced by the Classification Tool at any date in the study period. In addition, the Land Cover Tool allows users to specify a change layer, corresponding to conversion between one or multiple classes at a specific date to one or multiple other classes at a later date. The result is an N-band image corresponding to N-number of dates specified by the user, which can be visualized "on-the-fly" or exported as an image to an Earth Engine Asset, Google Cloud, or Google Drive (Figure 6). Unlike the tools described previously, the Classification and Land cover tools allow exporting the results as GEE assets. However, the Earth Engine App platform does not support exporting to GEE assets directly and therefore these tools cannot be distributed as "standalone" Earth Engine Apps but must be distributed as GEE scripts. For this reason, static links to the tools cannot be generated and need to be regenerated every time the code is updated. Up-to-date links to the tools can be found at https:// gee-ccdc-tools.readthedocs.io.

AREA2: Area Estimation and Map Accuracy Assessment
Land cover and land change are important climate variables, and analyses of remote sensing data are often the only feasible way to determine the extent of these variables. Traditionally, maps constructed by classification of remote sensing data were used to quantify such variables by summing map units assigned to map classes (i.e., pixel counting). However, as emphasized in the literature, areas obtained in this way may be subject to bias due to classification errors (McRoberts, 2011;Olofsson et al., 2014;GFOI, 2016). Instead, a sample-based approach in which an unbiased estimator is applied to sample data for estimation of areas of interest is recommended. In this context, "unbiased" refers to a property of an estimator: the bias of an estimatorμ of a population parameter µ is the difference between µ and the expected value ofμ over all possible samples. Expressed mathematically,μ is unbiased if Bias μ = E μ − µ = 0 (Casella and Berger, 2002, p. 330).
Designing surveys, determining reference conditions on the land surface, and selecting, constructing and using appropriate estimators are all important but complex steps and decisions. To assist users navigating such choices, we have developed a toolkit titled AREA2 (ARea Estimation and Accuracy Assessment;  (Hijmans et al., 2005) "area two"). AREA2 (Bullock et al., 2020a) provides support for the three major components of an estimation protocol: sampling design, response design, and analysis (Stehman and Czaplewski, 1998). The first component, sampling design, refers to the protocol for selecting a subset of population units (i.e., a sample) in the study area. AREA2 supports simple random (SRS) and systematic (SYS), stratified random (STR), and two-stage sampling. An important part of designing a survey is determining the sample size. With AREA2 we have implemented support for sample size determination using the variance estimators for SRS/SYS and STR designs [Equations 4.2 and 5.25 in Cochran (1977)] by specifying a target precision and then calculate the sample size required to achieve that precision.
Following the sampling design, the reference land surface condition needs recorded at each sample unit location by examining reference data. Good reference data have traditionally been commercial high-resolution data, aerial photos and field inventories but such data are costly and less than ideal for observing changes on the land surface through time. Reference data based on time series of observations have proven a powerful dataset for observing reference change conditions, especially when combined with high resolution imagery [e.g., Olofsson et al. (2016), ]. Traditionally, compiling time series of Landsat data for sample locations was demanding as whole images in a series needed downloading and stacking before a time series could be extracted. Attempts were made to mitigate these issues, the best example being TimeSync which extracted annual "best image" chips of Landsat data at sample locations . The viewer in AREA2 extracts time series of Landsat surface reflectance and various spectral indices from Google Earth Engine during a user specified time period (1982 at earliest). Clicking an observation in the series displays the associated Landsat image. The last step corresponds to the analysis, which includes protocols for estimating area, map accuracy and uncertainty (Olofsson et al., 2014). Central to the analysis is the estimator (preferably an unbiased estimator), which is a formula by which an estimate of a population parameter is calculated from the sample data (Cochran, 1977). Even though the estimator must correspond to the sampling design, there are usually more than one unbiased estimator to choose from and selecting the most efficient estimator in a given situation is often complex. Further, the implementation of the estimator (and the associated variance estimator) is not always straightforward, especially for complex designs such as two-stage sampling. AREA2 provides tutorials guiding users through the decisions, and supports the following estimators (the references provide the underlying mathematical frameworks and proof of unbiasedness): • expansion estimators for simple random and systematic sampling (Cochran, 1977). • stratified estimator for stratified random sampling (Cochran, 1977). • post-stratified estimator for simple random and systematic sampling (Cochran, 1977). • model-assisted regression estimator for stratified random sampling (Särndal et al., 2003). • ratio estimator for stratified random sampling when the map classes that defined strata are different from the map classes in the final error matrix (Stehman, 2014). • ratio estimator for use with two-stage sampling (Särndal et al., 2003).
AREA2 can be used with any categorical map, even if it was not created with the CCDC algorithm. For this reason it was created as a separate project, and the detailed documentation and tutorials are available at https://area2.readthedocs.io/.

Case Study for Land Cover Mapping and Area and Accuracy Estimation
We mapped land cover and land cover change due to expansion of built land in southeastern United States as a case study for using the classification tools and AREA2. We made use of previously developed training data from the NASA GLANCE project (http://sites.bu.edu/measures/) and reference data from LCMAP (Brown et al., 2020) to perform classification and assessment for the case study using the Classification tool. The study area consisted of the states of Florida, South Carolina, North Carolina, Louisiana, Mississippi, Georgia, Arkansas, Kentucky, and Alabama. The CCDC results used for this exercise correspond to preliminary CCDC outputs made available to us that are currently unavailable to the general public. However, a similar workflow can be carried out with any CCDC results run by a user, or by other publicly available or shared CCDC datasets. Following the steps of the Classification tool, we specified the path to the CCDC results and the training data, selected the FIGURE 6 | Screenshots of the Classification Tool used to classify CCDC segments for the country of Ethiopia (A), the Land Cover Tool which allows users to extract categorical maps, including maps of change, for any desired date or dates (B), and a subset of (B) that shows pixels classified as cropland gain due to displacement of forest (C).
Frontiers in Climate | www.frontiersin.org features to be used in the classification (All of the inputs specified in Table 1), the classifier to use (a Random Forest classifier) and exported the resulting multi-band image containing a land cover class per time segment per pixel. These results were used as inputs to the Land Cover tool to map the expansion of the built class, referred to from now on as the "New Built" class. The change period selected for the creation of the change layer was 1999-2016, but any period between the original date range of the classification results can be studied. After the change map finished loading "on-the-fly" in the tool map panel, we exported it as an asset to use in the AREA2 tool. Accuracy of the classification results were obtained using the reference sample, which was created independent on the training dataset. Since the reference sample was selected under a simple random sample design, it can be used to estimate map accuracy with AREA2 using the "direct estimator" available in the subfolder "3. Estimation" of the AREA2 repository.

Application Programming Interface
An API was built atop the GEE JavaScript API to facilitate custom operations with CCDC outputs and facilitate building the graphical tools presented here. The API currently contains the following modules: 1) Inputs: module to create Landsat, Sentinel 1 and 2 stacks for time series plotting, as well as the computation of indices derived from the optical sensors. 2) CCDC: module to manipulate the results from existing model runs done using Landsat observations, enabling the extraction of regression coefficients, the generation of synthetic images, and the extraction of time, magnitude and frequency of changes. 3) Classification: module to pre-process training data, classify CCDC coefficients and obtain land cover maps from them. 4) Dates: module to interact with dates in the formats required by the CCDC algorithm in Google Earth Engine. 5) Change: module to post-process the change results obtained from CCDC. Currently it contains functions to evaluate the presence of errors of omission and commission of change (Bullock et al., 2019).
Up to date links to access the API, the documentation, and a brief guide to use the core functions can be found in this link: https:// gee-ccdc-tools.readthedocs.io.

RESULTS
We illustrate the varied usage and multiple applications of the tools presented above with three examples. The time series viewers are demonstrated using a localized change example in Quincy, Massachusetts. The CCDC visualization tool is demonstrated using much larger fire events in northern Argentina. Finally, an example of land cover and land cover change mapping are shown for the Southeastern USA, followed by the estimation of accuracies and change areas. The different areas were selected to demonstrate how the tools can aid in the understanding of land change at different scales, and how they can be used in different geographical areas with varying levels of data coverage and drivers of change.

Time Series Viewers
The Simple Time Series Viewer has been used in university courses (at Boston University) to visualize the history of the landscape in situ through a mobile phone. On field trips, students were able to interact with Landsat time series plots, allowing them to address questions related to the landscape history at their immediate location. One of such locations was in Quincy, Massachusetts. Students accessed the tool through their mobile phones, located themselves in the map using the "User location" button, and clicked on the map to reveal the Landsat time series, as well as the time segments detected by CCDC on the fly. In the example shown in Figure 7, two temporal segments are detected by CCDC: one starting in 1997 and ending in 2004, and another starting in 2005 and ending in 2020. The first segment represents the forested area that was cleared to make space for housing, represented by the second segment. The time range between segments circa 2005 corresponds to the transition between these two land cover classes, when the CCDC algorithm waits for more stable observations to start fitting a new segment. The ability to easily explore historical imagery on mobile phones in-situ to make sense of the land changes over time is a groundbreaking advancement made possible by Google Earth Engine and its Apps platform. More complex temporal dynamics are easier to investigate using the full time series viewer on a larger browser, such as a tablet of laptop. Using the same example presented before, the Advanced time series viewer can be used to investigate other changes in the area in more detail. Figure 8 shows the time series for a pixel located in the golf course visible in the western side of Figure 7. This area, formerly a covered landfill, was transformed into a golf course during the early 2000's using the material extracted from a mega-construction project in the Boston area (commonly referred to as the "Big Dig"). The time series in Figure 8 show three segments representing these land covers and transitions: grass vegetation (first segment), followed by a complete removal and the landscape transformation that took place between 1998 and 2006 (second segment), and then the establishment of permanent, irrigated grass for the golf course around 2007 (third segment).

Spatial Data Viewer
The spatial data viewer tool is demonstrated using an example from northern Argentina. CCDC results for the period between 1985 and 2000 are loaded in the tool, and the "Create synthetic image" section on the right panel is used to create a predicted image for the date 1996-01-01, to use as a snapshot of the landscape. Fire scars are evident in the image, showing nongeometric patterns of change compared to nearby agriculture and roads (Figure 9). After changing the band selection on the left panel to "NBR" (Normalized Burned Ratio) and clicking anywhere in the fire scar, the time series loaded in the bottom panel show negative values circa 1996 that suggest the occurrence of a fire event.
To obtain a map that shows the approximate time of change for the entire area, the "Visualize section" on the right panel can be used to select a period centered in 1996 (e.g., 1994-1997), a spectral band to visualize the magnitude of the change detected by CCDC (i.e., SWIR1) and the desired change layer, in this case, the "Time of maximum magnitude of change." The loaded image shows a clear fire even that happened circa 1994 for this fire scar, as well as the timing for other change events in the surrounding areas.

Classification and Estimation of Area and Accuracy
The results consisted of: (1) Classified time segments between 2000 and 2016 for the study region, equivalent to land cover maps at 30 m resolution for any given date between that period, (2) a single land cover change map showing stable classes and the "New built" change class between 2000 and 206, and (3) estimates of area and accuracy with standard errors for the classes in the land cover change map. Figure 10 highlights the "New Built" class for two subsets in North Carolina and Tennessee. The results are displayed using a custom palette, defined manually in the Visualization Parameters of the Land Cover tool and represented in the map as a legend. The map contains seven stable classes (Water, Snow/Ice, Built, Bare, Forest, Shrub, and Herbaceous) and the change class of interest (New Built). It should be noted that the results in Figure 10 are shown as created "on-the-fly" in the Land cover tool but they must be exported as assets to be used in the AREA2 tools. The runtime to export the classified time segments as an Earth Engine asset from the Classification tool was 5 h, while the runtime to export the change map from the Land cover tool was 10 h, but these times will vary depending on several factors, including the area extent, pixel size, number of input features and time segments, and other factors inherent to GEE.
Area and accuracy estimates produced by AREA2 are shown in Tables 2, 3. The Overall accuracy of the map was 80% and the Producer's and User's accuracies of the "New Built" class were 40 and 15%, respectively. The sample-based estimates revealed 507,000 ha of conversion to Built during the study period. The predominant land cover class at the end of the study period in 2016 was the Tree class, representing 60% of the study area, or 71,142,000 ha.

DISCUSSION
The functionality of the tools presented in this paper is intended to simplify user interaction with time series of Landsat data, and to provide a system to explore and extract valuable information from the data, such as multi-temporal metrics and long term landscape change through the use of the CCDC algorithm. The use of CCDC by the broader community has been historically hampered by storage and preprocessing demands, and computational complexity. With the emergence of cloud computing and public repositories of satellite data, we have an opportunity to expand user access to CCDC outputs. To achieve broader access, we have developed a set of tools with graphical interfaces that enable a range of workflows from simple exploratory analysis to production of more complex products In essence, the tools presented here provide a new way for the public to interact with time series of satellite data with minimal effort. Our hope is that these tools will lead to a broader understanding of the way Earth is changing as well as an appreciation of the value of sustained satellite measurements over long time periods.
The tools presented here facilitate the monitoring and estimation of land change and associated carbon emission by removing the challenges such as the processing of large volumes of Landsat imagery, screening of noise and cloud/cloud shadow, and demands on storage space and computational power. Graphical user interfaces and API functionality, coupled with tutorials and documentation, provide sufficient flexibility to accommodate national and local circumstances and suit users with different levels of technical expertise. For example, users with little or no experience with the CCDC algorithm can use the exploratory tools to get familiarized with the basic concepts surrounding its use, such as the "time segments." Users that are more familiar with the outputs of the algorithm can visualize model coefficients or interact with them without having to process the CCDC results themselves. Similarly, countries without national maps of land cover change can start by exploring the time series data in their region to assess the suitability of the CCDC algorithm to capture change patterns of interest. Countries with higher image-processing capabilities can make use of the API to extract CCDC-derived products and process them with their own custom-made functionality. Furthermore, the availability of tools to create stratifications, assess the accuracy of map products, and estimate areas of change can benefit both new users creating preliminary maps, or experienced users with existing maps, even if created using different algorithms and tools.
An example that illustrates the use of CCDC outputs for monitoring, along with sample-based estimation, of multiple land cover transitions in support of REDD+ reporting was presented by . The study was conducted before the availability of the CCDC algorithm in GEE, and required manual downloading and stacking of more than 5,000 Landsat images. Visualization of time series for a single pixel was possible only after the stacking process was complete, and the CCDC algorithm had to be run for each WRS-2 path and row scene separately; this translated into running the model at least 25 times to cover the Colombian Amazon. Once the time segments were classified, each Landsat footprint was mapped separately, resulting in the need to merge individual maps into a coherent, seamless mosaic. This cumbersome workflow could now be completely reproduced in GEE using the API and tools presented here, without the need for extensive pre-and postprocessing. The gain in efficiency is illustrated in the example presented in Section classification and estimation of area and accuracy for the Southeastern USA. Even though the study area was larger than the area mapped in , the process was completed in just 15 h. Furthermore, the ability to inspect time series and visualize classification results on-the-fly allows for time-efficient quality control that is likely to result in improved maps.
Another example illustrating how the presented tools can accelerate map creation, quality control, and assessment processes is provided by a recent capacity-building collaboration with the Government of Zambia via SilvaCarbon. In this collaboration, we assisted Zambia's forestry department in using a set of existing training datasets to create land cover change maps, and to assess and improve their accuracy. First, we used the tools for time series visualization to determine if the land cover changes of interest were readily captured by the CCDC algorithm. This step was accomplished using the function in the Advanced Spatial Data Viewer that allows visualization of the breaks found by CCDC. Second, we wanted to test which of the various training datasets produced the most accurate map. The use of the CCDC and Classification modules in the API tools enabled iterative production of temporally consistent maps of land cover and land cover change for the entire country using each of the training datasets (Figure 11). Finally, the accuracies of the newly created maps were estimated using AREA2 to determine which of the training datasets produced the most accurate map. An important aspect of the Zambia example is capacity building; all steps were completed in just a few days of work with us the authors, using the tools in the beginning of the workshop but after only a day or two, the participants from the Zambia Forestry Department were producing the maps, inspecting the results and estimating accuracy measures. The example illustrates the potential of the presented tools to enable local practitioners in tropical countries to use advanced environmental remote sensing in support of climate science, and ultimately climate mitigation.
The tools and functionality available in the API are currently being used as part of the mapping and assessment procedures for the Global Land Cover mapping and Estimation (GLanCE) project, a NASA MEaSURES project at Boston University (http:// sites.bu.edu/measures/). GLanCE will provide a data record of twenty-first century global land cover, land use and land cover change at 30 m resolution. The intent is to provide global maps of land cover and land cover change with estimates of area and map accuracy. A similar effort was carried out for the United States by the USGS developed a system for operational monitoring of land change as part of its LCMAP initiative. To implement CCDC at a national scale, they built the computing infrastructure necessary to process the entire Landsat archive for the United States using CCDC, as well as implement similar kinds of tools to those presented in this paper to convert CCDC results to a suite of 10 operational products related to land cover and land cover change FIGURE 10 | Screenshot of the Land Cover Tool used to display the land cover classification at the end of 2016 and conversion to Built land ("New Built"). Insets detailing some of the largest change areas were manually overlayed for display in this figure. (Brown et al., 2020). The resulting global maps produced by the GlanCE project are expected to provide baseline characterization of the ways the land surface is changing and will aid in the characterization of biophysical properties of the Earth's land surface. We also expect the resulting products to contribute to model-based investigations that require land cover and change information, such as carbon balance models. The deliverables are complex and can be thought of as a database of time segments with a set of attributes in space and time. This complexity poses a potential hurdle for users, but we envision that the presented tools will allow the creation of maps of land cover and land cover change for any region in the world where the required data is available.
Once maps have been created, guidance is needed for how to use maps in an inference framework to estimate areas and map accuracy. AREA2 was designed to guide users through the decisions involved in the estimation process and to relieve users of the burden of constructing estimators. The tools in AREA2 have been used for this purpose at capacity building workshops in Costa Rica, Peru, and other countries, and for sample interpretation and area estimation in Guatemala as part of efforts to estimate forest disturbance in Guatemala's protected areas (Bullock et al., 2020b). As it becomes easier to create maps with existing and new tools like the ones presented here, it is crucial to assess and report the accuracy of such products. For this reason, we envision AREA2 as an integrated part of a mapping workflow that uses advanced remote-sensing products, such as those derived from the CCDC algorithm.
Some considerations and limitations must be considered when using the tools and functions presented in this paper. First,   they are provided free of charge "as-is" without any warranty or support, and they are subject to change at any time. We expect to continue improving and adding functionality to the tools, especially to the time series and classification tools, for at least another year or as long as it is feasible. The AREA2 repository will continue to be monitored for potential coding issues. Users are welcome to contribute to the development of new functionality or provide suggestions. Due to the current limitations for multiuser collaboration in a GEE repository, we have created a GitHub repository with a copy of the most recent and stable version of the API and tools and we will update it periodically to ensure the proper archival of new functionality. The repository is located at: https://github.com/parevalo/gee-ccdc-tools. Second, while the Google Earth Engine team has been supportive of the development of these tools, they are not produced or endorsed by Google. Similarly, this paper does not intend to serve as the official description of the GEE implementation of CCDC. Finally, the tools in their current iteration have not been designed or tested with CCDC results from sensors other than Landsat. However, our end goal is to provide a set of tools that are sensor agnostic. Web links to other tools designed to explore Sentinel data will be published on the website referenced previously as soon as they become available. Finally, while a primary objective of the tools presented here is to aid in land cover mapping, various other applications of benefit to users interested in land surface dynamics and their effect on climate are possible. For example, the frequency and magnitude of change events detected by the algorithm could be used to determine the timing of forest loss events or agricultural cycles. The variability in the surface reflectance could be visualized via the RMSE of the time segments. Seasonality information contained in the amplitude and phase of the harmonic terms of the regression could be used to study vegetation phenology as well as mapping of deciduousness. Long-and short-term trends in land cover condition, as captured by the slopes of the time segments, could be analyzed to identify vegetation recovery. Cloud-free synthetic images could benefit studies of the land surface, particularly in areas with low data density or high cloud occurrence. Previously published studies and algorithms that make use of CCDC outputs could be implemented for larger study areas. For example, a carbon bookkeeping model recently implemented to operate at the Landsat pixel scale using CCDC results (Tang et al., 2020) could be implemented and applied to a larger geographical area, such as the pantropics, using the tools we present here. Furthermore, current efforts to use CCDC results and the presented API to improve the mapping of biomass density in the Amazon basin by using the information contained in the temporal and spatial domains are other examples of how the tools could enable the study of land processes and properties directly related to climate. We believe that the tools presented in this paper will make these and other types of analyses more accessible to a broader set of scientists.

CONCLUSION
We have created a set of tools that facilitate the exploration of time series of Landsat observations and the time segments identified by CCDC. We believe these tools will simplify the visualization of time series of Landsat observations, facilitate the exploration of CCDC coefficients and derivatives, and streamline the creation of land cover and land cover change maps, as well as the estimation of areas and accuracies. By removing some of the barriers that prevent users from tapping into the information contained in time series of Landsat observations, we expect that users will extract new useful information that can be directly applied to climate science. Updated documentation and current links to the tools and API are available at https://gee-ccdc-tools.readthedocs.io.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary materials, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
PA and EB implemented the tools. All authors conceived and contributed to the text.

FUNDING
PO and PA were funded by support from the NASA Carbon Monitoring System, grant NNX16AP26G (PI: Pontus Olofsson). PA was also funded through NASA research grant 80NSSC18K0994 (PI: Mark Friedl). EB was funded through NASA research grant 80NSSC18K0994 (PI: Mark Friedl) and SERVIR Applied Science Team grant NNH16AD021 (PI: Sean Healey). CW was funded by support from the USGS/NASA Landsat Science Team. Parts of the code required to visualize the time series and time segments, as well as reprocessing of CCDC outputs were provided by Noel Gorelick and Zhiqiang Yang.