Medium Spatial Resolution Mapping of Global Land Cover and Land Cover Change Across Multiple Decades From Landsat

Land cover maps are essential for characterizing the biophysical properties of the Earth’s land areas. Because land cover information synthesizes a rich array of information related to both the ecological condition of land areas and their exploitation by humans, they are widely used for basic and applied research that requires information related to land surface properties (e.g., terrestrial carbon models, water balance models, weather, and climate models) and are core inputs to models and analyses used by natural resource scientists and land managers. As the Earth’s global population has grown over the last several decades rates of land cover change have increased dramatically, with enormous impacts on ecosystem services (e.g., biodiversity, water supply, carbon sequestration, etc.). Hence, accurate information related to land cover is essential for both managing natural resources and for understanding society’s ecological, biophysical, and resource management footprint. To address the need for high-quality land cover information we are using the global record of Landsat observations to compile annual maps of global land cover from 2001 to 2020 at 30 m spatial resolution. To create these maps we use features derived from time series of Landsat imagery in combination with ancillary geospatial data and a large database of training sites to classify land cover at annual time step. The algorithm that we apply uses temporal segmentation to identify periods with stable land cover that are separated by breakpoints in the time series. Here we provide an overview of the methods and data sets we are using to create global maps of land cover. We describe the algorithms used to create these maps and the core land cover data sets that we are creating through this effort, and we summarize our approach to accuracy assessment. We also present a synthesis of early results and discuss the strengths and weaknesses of our early map products and the challenges that we have encountered in creating global land cover data sets from Landsat. Initial accuracy assessment for North America shows good overall accuracy (77.0 ± 2.0% correctly classified) and 79.8% agreement with the European Space Agency (ESA) WorldCover product. The land cover mapping results we report provide the foundation for robust, repeatable, and accurate mapping of global land cover and land cover change across multiple decades at 30 m spatial resolution from Landsat.

Land cover maps are essential for characterizing the biophysical properties of the Earth's land areas. Because land cover information synthesizes a rich array of information related to both the ecological condition of land areas and their exploitation by humans, they are widely used for basic and applied research that requires information related to land surface properties (e.g., terrestrial carbon models, water balance models, weather, and climate models) and are core inputs to models and analyses used by natural resource scientists and land managers. As the Earth's global population has grown over the last several decades rates of land cover change have increased dramatically, with enormous impacts on ecosystem services (e.g., biodiversity, water supply, carbon sequestration, etc.). Hence, accurate information related to land cover is essential for both managing natural resources and for understanding society's ecological, biophysical, and resource management footprint. To address the need for high-quality land cover information we are using the global record of Landsat observations to compile annual maps of global land cover from 2001 to 2020 at 30 m spatial resolution. To create these maps we use features derived from time series of Landsat imagery in combination with ancillary geospatial data and a large database of training sites to classify land cover at annual time step. The algorithm that we apply uses temporal segmentation to identify periods with stable land cover that are separated by breakpoints in the time series. Here we provide an overview of the methods and data sets we are using to create global maps of land cover. We describe the algorithms used to create these maps and the core land cover data sets that we are creating through this effort, and we summarize our approach to accuracy assessment. We also present a synthesis of early results and discuss the strengths and weaknesses of our early map products and the challenges that we have encountered in creating global land cover data sets from Landsat. Initial accuracy assessment for North America shows good overall accuracy (77.0 ± 2.0% correctly classified) and 79.8% agreement with the 1 INTRODUCTION: CONTEXT, JUSTIFICATION, OUTLINE Land cover plays a key role in the Earth's climate, ecological, and socio-economic systems (Bonan, 2008;Foley et al., 2011), and maps of land cover provide the single most important basis for characterizing the ecological state and biophysical properties of the Earth's land areas (Feddema et al., 2005;Foley et al., 2005). Land cover and land surface properties influence biosphereatmosphere interactions and the role of land cover and land use in the global carbon cycle is both well-documented and profound (Houghton, 2020). Accurate information related to the global distribution of global land cover is therefore required to parameterize land processes in regional-to-global scale Earth system models (Bonan et al., 2002;Fuchs et al., 2018). In addition, as the Earth's population has surged over the last several decades, the global area of land dominated by humans has rapidly expanded, ecologically important regions have been degraded, the area available for new land use (e.g., cities and croplands) has decreased, and land resources have become increasingly scarce (Goldewijk, 2001;Sanderson et al., 2002;Foley et al., 2005;Ellis and Ramankutty, 2008). Hence, accurate and timely information related to land cover and land cover change is essential for managing natural resources and for understanding the ecological, biogeographic, and resource management footprint of society (Song et al., 2018).
Remote sensing has been used to map land cover for over four decades (Strahler, 1980;Loveland et al., 1991;Townshend et al., 1991). Until relatively recently global land cover and land use data were only available at coarse spatial resolution from data sets such as the Global Land Cover Type product (Friedl et al., 2010;Sulla-Menashe et al., 2019) produced by NASA, and the GlobCover (Arino et al., 2008) and Climate Change Initiative Land Cover data sets (Plummer et al., 2017) produced by the European Space Agency (ESA). Because these products are based on coarse spatial resolution imagery (nominally at 500 or 300 m), the representation of land cover they provide is highly generalized. Further, because spatio-temporal variation in most land cover occurs at spatial scales well below the resolution provided by coarse spatial resolution remote sensing, these products do not meet the needs of the large, diverse, and growing community of scientists and applied stakeholders who require high quality and high spatial resolution information related to land cover, and ultimately land cover change, over large areas.
To address this need, global land cover products at medium spatial resolution based on imagery from Landsat and other sensors such as Sentinel 1 and 2 have started to become available. The most widely used of these products tend to be focused on specific themes including urban land use or impervious surface mapping (Liu et al., 2018;Gong et al., 2020;Liu et al., 2021), forest cover (Hansen et al., 2013;Kim et al., 2014;Feng et al., 2016), or surface water (Pekel et al., 2016), and therefore support a focused community of end-users. In parallel, several general-purpose land cover maps that depict a finite set of discrete classes (generally~10) have also been created. The GlobeLand30 data set provides maps of global land cover at 30 m for 10 classes for 2000, 2010, and 2020 (Chen et al., 2015). The iMap World data set provides land cover data and maps at 30 m for the period 1985-2020 based on Landsat imagery . The GLC_FCS30 product provides maps of global land cover for 2015 in three different land cover schemes with 9, 16 and 24 classes (Zhang X. et al., 2021) and, most recently, ESA created the WorldCover data set for 2020, which includes 11 land cover classes and is produced at 10 m spatial resolution using data from Sentinel 1 and 2 (Zanaga et al., 2021). Significantly, these two latter products do not provide information related to land cover change, which is required for many applications and is essential for characterizing and modeling changes in global land cover that have occurred over the last several decades.
To address the need for high-quality long-term records of land cover and land cover change we are using the Landsat archive to create a global record of 21st Century annual land cover and land cover change at 30 m spatial resolution. In this context, the goal of this paper is to describe the methods, data sets, data products, and early results from this project. The paper specifically focuses on land cover mapping results, which provide the foundation and first step towards a broader goal of robust, repeatable, and accurate mapping of global land cover change across multiple decades from Landsat.

MAPPING PROCESS
The Global Land Cover Estimation (GLanCE) mapping process includes five distinct elements: 1) a global land cover classification legend; 2) a global data set of training data used for classification model estimation; 3) algorithms for fitting models to pixel-scale time series of Landsat reflectance values; 4) classification of land cover at each pixel; and 5) post-processing adjustments applied to annual classification products to correct errors in the raw classification results. In the sections that follow, we provide a concise description for each of these elements.

Classification Legend
The GLanCE land cover classification scheme includes two nested sets of classes. Currently we are only mapping classes included in Level 1 ( Table 1). The classes are mutually exclusive and are similar (but not identical) to the land cover classes used in the United States Geological Survey (USGS) Land Change Monitoring and Assessment Program (LCMAP; Brown et al., 2020). A key attribute of the legend is that it is designed to map land cover, not land use, although it does include a "Developed" class for land areas composed of buildings, roads, and other features associated with the built environment. Another feature of the legend is that it corresponds to the IPCC top-level land categories for greenhouse gas inventory reporting (GFOI, 2020). A variety of additional land cover attributes (e.g., leaf type, phenology, agricultural land use; Figure 1) are being recorded by interpreters as they compile the training site database that will support the GLanCE Level 2 classification scheme. In the long term, we plan to add this level to our mapping system. In the near term, we are focusing on mapping the Level 1 classification.

Training Data
The GLanCE training site database is derived from four sources. First, a team of research assistants at Boston University is using a dichotomous key and online tools for efficiently interpreting high resolution imagery developed to support this project (Arévalo et al., 2020a) to collect training data. As a starting point, we adapted the training site database developed in support of the MODIS Collection 6 Land Cover Type Product (MCD12Q1), which includes 3,395,500 m MODIS pixels selected to provide a representative sample of land cover in all major ecoregions of the world (Sulla-Menashe et al., 2019). To incorporate these data into the GLanCE training database, 10 Landsat pixels located within each MCD12Q1 training pixel were randomly selected, visually interpreted using high resolution imagery, and labeled by research assistants using the dichotomous key shown in Figure 1. Second, to augment these sites, we developed a procedure that performs unsupervised clustering on spectraltemporal features estimated from Landsat image time series using the Continuous Change Detection and Classification (CCDC) algorithm (Zhu and Woodcock, 2014). Clustering was performed at the scale of World Wildlife Fund (WWF) ecoregions (Olson et al., 2001). 14 to 482 pixels were selected in each ecoregion (median = 327 sites/ecoregion) based on this clustering (i.e., using TABLE 1 | Land cover classes included in the GLanCE Level 1 classification scheme. Note that in addition to these layers, the product includes a variety of complementary information (see Table 3).

Name Description
Water Areas covered with water throughout the year: streams, canals, lakes, reservoirs, and oceans Developed Areas of intensive use; land covered with structures, including any land functionally related to developed/built-up activity Barren/Sparsely Vegetated Land comprised of natural occurrences of soils, sand, or rocks where less than 10% of the area is vegetated Tree Cover Land where the tree cover is greater than 30%. Note that cleared trees (i.e., clear-cuts) are mapped according to current cover (e.g., barren/sparsely vegetated, shrubs, or grasses) Herbaceous Land covered by herbaceous cover. Total vegetation cover exceeds 10%, tree cover is less than 30%, and shrubs comprise less than 10% of the area Shrublands Land with less than 30% tree cover, where total vegetation cover exceeds 10% and shrub cover is greater than 10% Ice/Snow Land areas where snow and ice cover is greater than 50% throughout the year FIGURE 1 | GLanCE dichotomous training data key, legend, and attributes used to assign land cover to training sites. Green boxes show the Level 1 land cover classes that are being mapped for the GLanCE data product; Yellow boxes show additional attributes that are recorded to support Level 2 classification. BL, broadleaf; NL, needleleaf.
Frontiers in Remote Sensing | www.frontiersin.org June 2022 | Volume 3 | Article 894571 the clusters to stratify the spectral-temporal feature space), which were then interpreted by analysts to create training sites. In doing so, this procedure provides a sub-sampling of the feature space from Landsat imagery that ensures each ecoregion is adequately sampled while also ensuring that training sites are allocated efficiently. For both the MODIS-and cluster-derived sites, interpretations were performed on-screen using highresolution imagery from openly available sources such as Google Earth in combination with time series of Landsat surface reflectance values and vegetation indices using tools developed on Google Earth Engine (GEE). Multiple analysts reviewed all samples to ensure high-quality training data. The third source of training data is data provided by trusted collaborators. This process is non-systematic, opportunistic, and depends on the availability of data from collaborators. However, because good quality training data is essential to estimating highquality classification models and acquisition of training data is difficult and labor-intensive at global scale, training data provided by collaborators who have regional expertise in land cover and land cover change processes is extremely valuable. In North America, for example, we incorporated a data set that includes 17,476 randomly sampled reference sites distributed across the conterminous United States that were used to assess the accuracy of map products from the USGS LCMap program (Stehman et al., 2021), along with 9,073 sites from Canada and Alaska that were compiled as part of a separate project funded by NASA's Arctic-Boreal Vulnerability Experiment (ABoVE) . Similarly, in South America, we are leveraging training sites distributed across the Northern half of South America (Brazil, Bolivia, Peru, Ecuador, Columbia Venezuela, French Guyana, and Suriname) provided by colleagues involved in the MapBiomas project (Souza et al., 2020). For all such cases, integration of training data provided by collaborators into our database requires collaboration with the provider to crosswalk the training site attribute data to the GLanCE classification legend. In addition, to avoid classifier bias (e.g., Zhu et al., 2016), when a contributed data set is large we sub-sample the data to enforce pseudouniform training data density across space so that regions with data contributed by collaborators are not over-sampled relative to other regions. For example, the MapBiomass data set includes 464,129 pixel-years of data from 2000-2019. For GLanCE, we use a relatively small random sample (n = 20,119) of this larger data set.
The fourth source of training data that we incorporate is derived from the GeoWiki Global Land Cover and Land Use reference data set (Fritz et al., 2017). This data set includes over 150,000,100 m grid cells with land cover information collected via crowd-sourcing. For GLanCE, the GeoWiki data was processed and subset using four steps: 1) we excluded training samples that did not have a year assigned to the attributes; 2) for samples that were assigned multiple land cover types in a given year, we used the mode of the assigned land cover labels; 3) we cross-walked the GeoWiki legend to the GLanCE legend and excluded samples that could not be clearly assigned to a class in the GLanCE Level 1 legend ( Table 1); and 4) we only retained training samples where the land cover label agreed with the ESA WorldCover map aggregated to 30 m. After applying these steps, we retained 16,543 training GEoWiki training samples, all of which were added to the GLanCE global training site database.
Our training data workflow collects data on a continent-bycontinent basis. Once an initial training data set is compiled for a continent, we perform four addition steps to create the final data set that we use for mapping. First, because human error is inherent to all manual interpretation and training data compilation activities (Elmes et al., 2020), we adapted the basic approach described by Brodley and Friedl (1999) that uses machine learning to identify and remove misrepresentative or incorrectly labeled training samples. Second, we inspect the provisional map results and strategically augment the training set in each continent to correct regional under-or over-representation of land cover classes. Third, to augment the training set in regions where training data is sparse we used the algorithm described by Zhang and Roy (2017), which uses the MODIS Land Cover Type product (Sulla-Menashe et al., 2019) to automatically select candidate training samples at Landsat resolution, and removed candidate samples where the label disagreed with the land cover label in the ESA WorldCover data set. We then used the ESA WorldCover data set, cross-walked to the GLanCE land cover scheme, to estimate the expected frequency distribution of GLanCE land cover classes within 1,000 × 1,000 km grid cells, and then sampled the candidate training sites to approximately match this distribution. Note that we only applied this procedure to cells that had less than 250 training sites. These last two steps ensure that the training data has a minimum sampling density and reduces bias in classification results arising from misrepresentative class distributions in training data.

Continuous Change Detection and Classification
The core algorithm being used in this project to perform mapping is the Continuous Change Detection and Classification (CCDC) algorithm (Zhu and Woodcock, 2014). CCDC assumes that noise is ephemeral and occurs at the time scale of individual images (e.g., as a result of undetected clouds), that land cover change is persistent, and uses all available Landsat observations to map land cover and identify the timing of land cover change at each pixel. To accomplish this, the algorithm includes two core steps: (i) Identification of change points and modeling of stable time segments. Time series of surface reflectance in each band at each pixel are first filtered to remove observations affected by clouds, cloud shadows, and snow. The resulting time series is then modeled as a Fourier series with three harmonics. To fit the time series at each pixel, CCDC estimates an initial model based on the first 3 years of observations at each pixel.
Subsequent reflectance values at each pixel are then compared against fitted values, and change points are identified based on persistent large residuals (across time) between new observations and model fits based on Landsat visible, near-infrared, shortwave infrared reflectance values. When no change is detected after 3 years, new data are appended to the time series and the model is re-fit. This process is applied recursively until a land cover change (if any) is detected. Change points are identified using a change vector metric that integrates differences between observed and predicted reflectance in each spectral band, weighted by the root mean squared error of the model fit for each spectral band (Zhu and Woodcock, 2014). The concept of a time segment-the period between two consecutive land cover changes (if any), bounded by the start and end of the time series-is central to the CCDC algorithm. To illustrate, Figure 2 shows a time series of shortwave infrared reflectance values for a single Landsat pixel, where the red dot corresponds to a change point identified by CCDC. In the example shown, CCDC identifies two distinct time segments, one beginning in 2004 and ending 2012 (when a change is detected), and the second extending from 2012 to the end of the time series. (ii) Assignment of class labels to time segments. Spectraltemporal features for each time segment are used to assign a land cover label to each pixel for each segment. This approach has two important advantages. First, data from the entire time segment contributes to the classification. This provides stronger statistical support for classification at each pixel than is provided using data from individual dates with no temporal information. Stated another way, the phenology or seasonal variation in reflectance at each pixel is a key feature that is diagnostic of land cover in each time segment at each pixel. Second, in the classification stage, instead of using time series of surface reflectance as input to the classifier, CCDC uses model parameters estimated from harmonic models fit to time series of surface reflectance in each Landsat band, which has been shown to be highly effective (Zhu and Woodcock, 2014). For this project we adopt a strategy that was developed as part of the USGS LCMAP project, which is also using CCDC, and estimate classifications at annual time steps. To perform classification, we use the Random Forest classification algorithm (Breiman, 2001) estimated using training data described in Section 2.2. We estimate a single classifier for each continent, incorporating a variety of ancillary data layers to account for sub-continental variation in the distribution of land cover types (see next section).

Applying CCDC at Global Scale
CCDC has been successfully and widely used to map land cover and land cover change at regional scale (Arévalo et al., 2020b;Brown et al., 2020;Wang et al., 2020). However, given the data volumes and challenges involved in mapping land cover and land cover change at global scale, implementation of CCDC required two key modifications relative to how it has been applied in the past at local or regional scales: • First, to accommodate the computational requirements involved in mapping land cover and land cover change across multiple decades at Landsat spatial resolution for all global land areas (excluding Antarctica), we are using an implementation of CCDC on Google Earth Engine (GEE) (Gorelick et al., 2017). The cloud-based infrastructure provided by GEE supports rapid iteration and refinement of map products over very large areas in ways that are both essential for this project and impractical on highperformance compute facilities. Note that an added advantage of using GEE is that it allows us to easily share the data sets and tools that we create with a large (and growing) community of GEE users (Arévalo et al., 2020). • Second, a variety of ancillary data sets and land cover map products are used to complement information from Landsat and improve classification results generated by CCDC ( Table 2). We use the World Settlement Footprint and Worldpop Global Population data sets as features in the classification process to help distinguish the Developed land cover class from the Bare/Sparsely Vegetated class. The Global Surface Water and Topography classes (from Landsat and SRTM/ALOS, respectively) are included as features to help distinguish water bodies from other land cover classes and reduce spurious mapping of water in shadowed areas located in regions with complex topography. Finally, we include mean annual temperature, precipitation, and aridity (the ratio of total evapotranspiration to actual evapotranspiration) from the TerraClim data set as features, which helps the classifier label sub-classes in the Level I legend that occur at continental scale (e.g., woody and herbaceous cover at low versus high latitudes). Note that these ancillary data sets are provided as static features, and hence do not change.
Year-to-year change is captured by variation in the Landsat imagery.

ACCURACY ASSESSMENT AND REFERENCE DATA
For land cover and land cover change maps to be scientifically valid and useful, their accuracy and uncertainty need to be quantified and documented using high quality observations of reference land cover conditions collected via probability sampling (Stehman, 2000;Olofsson et al., 2014). To provide this, we are compiling reference data sets and using them to estimate: 1) the overall accuracy of the land cover and land change maps; and 2) the global user's and producer's accuracy of each map category, along with estimates of associated 95% confidence intervals. Two sets of random samples are being collected to support characterization of error and uncertainty in both land cover and land cover change. The first sample is designed to characterize the overall accuracy of the land cover data sets in any given year. This reference sample is drawn from all of the Earth's land areas by simple random sampling. The second reference sample is designed to characterize the accuracy and uncertainty of land cover change. Because land-cover change often impacts only a small proportional of the landscape, this latter data set is selected by stratified random sampling with strata corresponding to map classes identified as stable versus changed by CCDC. This latter data set will be compiled over the next year once the GLanCE land cover change products are produced.
For both reference data sets, the design-based inference framework that we are employing requires high quality reference observations for the assessment statistics to be reliable and robust. To ensure this, our data collection process includes reviews of each reference site by multiple analysts with deep experience in land cover mapping. Interpretations for sites that remain uncertain after multiple reviews are removed from the sample and new units are selected using the original sampling design.
A key objective of our accuracy and uncertainty assessment is to quantify uncertainty in the estimated area of land cover types and change areas. Because classification errors introduce area bias and smaller or more complex categories such as deforested lands are especially susceptible to bias and large errors, the area of individual map categories cannot be estimated by simply counting pixels (McRoberts, 2011). To account for this, we apply the stratified estimators provided in Olofsson et al. (2013) to the sample data for area estimation. Note that because our data products are still in development and the final strata are not yet defined, the results we present below are preliminary and do not include an assessment of GLanCE land cover change results and we only provide results for North America, where our map products and reference data are more mature. We include the above description of our overall approach for completeness and to emphasize that rigorous error and uncertainty assessment are key elements that will be provided with the final data products when they are released.

DATA PRODUCTS
Once complete, the GLanCE data record will cover the period from 2001-2020 at annual time steps. The science data sets (SDSs) included are designed to provide three main types of information for each year in the data record (Table 3). First, five SDSs provide information related to: 1) the land cover class, 2) the estimated quality associated with the assigned label, 3) whether the assigned label reflects a change in land cover from the previous year, and for those pixels where change has occurred, 4) an index that quantifies the magnitude of change, and 5) the approximate date of change during the year. Second, three SDSs are included that characterize the overall greenness, the annual amplitude in greenness, and the trend in greenness (if present) at  Huete et al., 2002). Third, two SDSs are included that indicate the inferred leaf type and phenology for pixels classified as trees. Note that because we are developing the GLanCE data sets on Google Earth Engine, all the data sets listed above will be made available via this platform as well as the LP DAAC.

PROGRESS AND INITIAL RESULTS
To create the data products identified in Table 3, our mapping process includes three core steps, each performed independently at continental scale: 1) creation of training data; 2) classification and mapping; and 3) assessment using reference data. Currently, we have completed a first version of our global training site database, created a first version of our global maps for 2001-2020 (hereafter, V1.0), and we are in the process of refining and assessing our mapping process and results, focusing on North America. Below we summarize our initial results from each of these core steps.

Training Data Collection
Excluding regional training data provided collaborators, the global training database includes 87,333 sites, each composed of a single Landsat pixel visually interpreted for a specific year in the time series (Figure 3). The number of sites in each continent is roughly proportional to the area of each continent-Eurasia, which is by far the largest continent, has 35,490 sites and Oceania (the smallest) has 4,631 sites. Because we allocate sites by sampling within ecoregions, the density of training sites is relatively uniform across space, with modestly lower site density in larger and more uniform ecoregions, and vice versa. Note that the density of points in Brazil in Figure 3 is less dense than elsewhere because the figure only shows training data collected through this effort, which does not include the densely sampled training data provided from the MapBiomas project (Souza et al., 2020). Also, as we briefly describe in Section 2.2, we use machine learning to filter the training set in each continent, which removes roughly 10% of the total training data. The distribution of classes included in the training data set is variable across continents and reflects the biogeography of vegetation, modified to varying degrees, by land use in each continent. Herbaceous and Tree Cover are uniformly the most frequent classes in the training site database across all continents, followed by Shrubland and Barren/Sparsely Vegetated sites. The remaining three classes (Developed, Water, Snow/Ice) comprise a much lower proportion of the training data, which reflects the relative frequency of these classes in global land areas.

Initial Global and North American Mapping Results
Our mapping process proceeds on a continent-by-continent basis. V1.0 of the GLanCE global land cover data set includes annual maps at 30 m spatial resolution for the period 2001-2020 based on the Level 1 legend in Table 1. At global scale, the initial results realistically capture the expected geographic distribution of land cover (Figure 4). Consistent with the distribution of training data, Herbaceous and Tree Cover are the most extensive classes, Developed and Water are the least common, with Shrublands and Barren/ Sparsely Vegetated land cover occupying a substantial, but variable, proportion of the landscape across continents.
In North America, Herbaceous land cover is the most extensive land cover type, followed by Tree Cover, Shrubland, and Barren/Sparsely vegetated land areas ( Figure 5). Developed land cover occupies 0.8% of North America, with isolated alpine areas covered by Snow/Ice. As we previously alluded to, classes that include mixtures of land cover and plant functional types were particularly challenging to map. At high latitudes, for example, two factors made it difficult to distinguish between Herbaceous cover and Shrublands. First, relative to other regions, there is relatively little high-resolution imagery available at high latitudes. Hence, compiling training data in boreal and Arctic areas is challenging and both the density and quality of training data in boreal and Arctic North America is lower than in other biomes. Second, much of the landscape at high latitudes is composed of mixtures of lichens and mosses (neither of which technically constitute Herbaceous cover or Shrublands), low stature trees, herbaceous cover, and dwarf shrubs. Hence, the GLanCE land cover legend is not well suited for this region and a large proportion of 30 m Landsat pixels at high latitudes are mixed in terms of their land cover composition (Beamish et al., 2020). Consequently, much of the landscape north of the tree line where lichens, moss, and dwarf shrubs are present is mapped as Herbaceous cover, and Shrubland cover is probably under predicted in non-forested boreal and Arctic regions (cf., Myers-Smith et al., 2011, 2020Myers-Smith et al., 2019).  We encountered similar challenges mapping land cover in arid and semi-arid regions at lower latitudes. Specifically, distinguishing among Barren/Sparsely Vegetated areas, Herbaceous cover, and Shrublands is difficult because (similar to issues encountered at high latitudes), most pixels in these regions include mixtures of at least two cover types. This issue is well documented at coarse spatial resolution such as MODIS (e.g., Friedl et al., 2010;Sulla-Menashe et al., 2019). The higher spatial resolution of Landsat reduces, but does not resolve, the impact of this issue. Hence, it is important to note that while we assign a unitary label to each pixel, relatively few pixels are actually uniform in regard to their land cover composition, even at 30 m spatial resolution.
To provide a qualitative assessment and comparison at local scale, Figure 6 compares land cover from 2019 from GLanCE with results from the ESA WorldCover map for 2020, crosswalked to the GLanCE legend. The first column (A) compares ã 400 km 2 area in Louisiana where the GLanCE map shows significantly more Developed and Herbaceous land cover relative to WorldCover, which reflects the different threshold for tree cover that is used in each map (10 vs. 30%). The second column shows a~900 km 2 area in southern Amazonia where, aside from minor differences, the two maps show strong agreement. The third column shows a~100 km 2 area in Southern China that includes a complex mix of Trees, Water, Herbaceous, and Developed land cover. The two maps show significant disagreement at this location, particularly regarding the extent of water, but also regarding the amount of Bare and Sparsely Vegetated land cover, which occupies a substantial portion of the landscape in the WorldCover map but is not present in the GLanCE map. This location is a coastal area with a high concentration of industrial fish farms, and visual inspection of these results suggest that the WorldCover map provides a more realistic representation of land cover in this location than the GLanCE map.

Preliminary Accuracy Assessment
At this point in time, our maps are not finalized, and we have not completed compiling our global reference data set. However, to provide a preliminary assessment of our initial mapping results, we used a simple random sample of 1,630 reference observations from North America that have been collected and undergone quality assurance to assess the overall and class-specific accuracies of GLanCE land cover mapping results (Figure 7). Results from this comparison are shown in Table 4.
The overall accuracy (percentage correctly classified) of the V1.0 GLanCE land cover map for North America was 77.0 ± 2.0% (95% confidence interval). Producer accuracies ranged from 40.7 ± 8.9% for the Barren/Sparsely vegetated class to 96.8 ± 3.5% for the Water class. User accuracies ranged from 68.9 ± 9.3% for the Herbaceous class to 95.8 ± 4.0% for the Water class. Tree Cover, which is geographically extensive in North America, shows high User accuracy (83.9 ± 7.4%), which is encouraging. At the same time, the Shrublands and Herbaceous classes show the lowest user accuracies (69.7 ± 9.2% and 68.9 ± 9.3%, respectively) that are below our target of 75%, which reflects the challenges involved in mapping these classes that we alluded to in Section 5.2.
It is important to note that while the sample we used for this analysis is designed to provide an estimate of the overall accuracy with a margin of error of ±5%, which was achieved, the sample size fore reference data in several classes is low. Specifically, the Developed, Ice and Snow, and Barren/Sparsely Vegetated classes have 25, 27, and 118 reference observations, respectively. Hence, the uncertainty in estimated accuracy for these classes is high.
As an additional comparison, to complement the accuracy assessment presented in Table 4, we used the reference locations shown in Figure 7 to compare results for North America from GLanCE against land cover labels from the ESA WorldCover  (Table 5). To perform this comparison, we used GLanCE results from 2019, which is the most recent year we have mapped at the time of this writing. To crosswalk the WorldCover land cover legend to the GLanCE legend, we merged the Grassland, Cropland, Herbaceous wetland, and Moss and Lichens classes into a single class (Herbaceous  (Table 4), and to compare the 2019 GLanCE land cover map against the ESA WorldCover 2020 data set ( Table 5). cover), and the Mangroves class was merged with Tree cover. We then upscaled the 10 m WorldCover data to 30 m using the mode of the 9 WorldCover land cover labels located in each GLanCE pixel. Overall agreement between the two maps is 79.8%, which is slightly higher than the reported accuracy of each map (77.0 and 74.4% 1 for GLanCE and WorldCover, respectively). Disagreement is highest among the Barren/Sparsely Vegetated, Shrublands, and Herbaceous cover types, which is consistent with the accuracies presented in Table 5 and accuracies reported for the WorldCover data set.

Land Cover Mapping Challenges
The qualitative and quantitative results for the V1.0 GLanCE product provide confidence that global land cover can be effectively and accurately mapped at 30 m spatial resolution from Landsat using the training data sets we've compiled in association with strategically selected ancillary feature sets. At the same time, our initial results point to some key challenges that need to be addressed. Most significantly, as we previously alluded to, sub-pixel mixtures of land cover, even at the 30 m spatial resolution of Landsat, continues to be a significant challenge for land cover mapping, especially in regions with sparse or low stature vegetation such as temperate drylands and Arctic ecosystems.
In Arctic and boreal regions, where land cover and plant functional types vary across short length scales (e.g., a few meters) because of microtopography and hydrology, assignment of land cover labels is challenging, both for analysts compiling training data and for machine learning-based classification algorithms. Indeed, sub-pixel mixing of land cover in many high latitude landscapes, even at the resolution of Landsat, effectively means that assignment of a pixel-scale land cover label as Shrublands or Herbaceous cover is often challenging, and the resulting labels have systematically lower confidence.
It is also important to reiterate in the specific context of Arctic regions that the land cover legend we are using for GLanCE is not well-suited for regions where lichen and moss cover is common and extensive. Ideally, we would add a class similar to the Lichen and Moss class included in the ESA WorldCover data set. However, short growing seasons, cloud cover, and large solar zenith angles make acquisition and interpretation of high-quality training sites that distinguish herbaceous life forms from lichens and mosses challenging. Hence, even if this class was added to the GLanCE legend, it may not be possible to populate our training data set with sufficient samples to support mapping it. At lower latitudes, we encountered similar challenges distinguishing between Shrublands and Herbaceous Cover in arid and semiarid systems. However, greater availability of high-quality highresolution imagery at lower latitudes in North America substantially reduces this challenge relative to the challenges at high latitudes described above.

Challenges Mapping Global Land Cover Related to the Landsat Archive
A key requirement for the land cover and land cover change mapping approach that we have developed for GLanCE is that relatively dense time series of Landsat observations are available. Specifically, previous work has shown that seven clear sky observations per year are required to fit robust CCDC models . Unfortunately, the number of images in the Landsat archive is quite variable both in space and time (Wulder et al., 2016) and this requirement is not met in some regions. Generally, more observations are available during periods when two Landsat sensors are operating (i.e., after the launch of Landsat 7, which overlaps with both Landsat 5 and Landsat 8). Since 1999, the USGS has been responsible for Landsat operations and implemented an acquisition strategy designed to provide systemic global imagery. This is the primary reason why we restrict the GLanCE product to 2001-present. However, the number of available cloud-free images at any location is still quite variable globally and depends on which Landsat sensors were in operation at any given time. Additional variability is introduced by the presence of clouds, cloud shadows, snow, sensor saturation, hazy observations (based on atmospheric opacity), and lack of aerosol optical depth information (Zhang et al., 2022).
The net result is that while some regions (e.g., the conterminous US, Canada, Australia) have deep and dense time series, the number of observations in some parts of world is insufficient to support CCDC-based mapping. Key examples of such regions include coastal West Africa, the Pacific coast of Colombia and Ecuador, and much of Siberia. Our current results are based on Landsat Collection 1 and preliminary analyses suggest that for most of the world the number of available images is higher in Landsat Collection 2, especially in the Landsat-8 era. Future versions of the GLanCE product will be based on Collection 2, which will help to mitigate-but not eliminate-regional challenges related to data density. For some regions (e.g., Siberia, West Africa), however, data density will continue to be too sparse to use CCDC for mapping. Hence, an alternative strategy will be required to map land cover and land cover change in these regions. To address this, we are currently investigating mapping approaches that have previously demonstrated the ability to map land cover and land cover change using comparatively sparse time series (e.g., LandTrendr) (Kennedy et al., 2010). For regions where data density is insufficient to support the use of CCDC, we will create maps using this "back-up" algorithm and identify pixels where this approach was applied via quality assurance layers.

Land Cover Change Detection and Mapping
Moving forward, as we mature our mapping process and the quality of our land cover maps improve, our efforts are heavily focused on land cover change mapping. To illustrate the change detection and mapping process implemented based on CCDC, Figure 8 shows two time series of GLanCE classification results for areas undergoing land cover change, along with time series of surface reflectance from the SWIR1 band of Landsat. The top panel shows a GLanCE time series that captures expansion of Developed land cover associated with construction of an airport terminal in Fort Meyers, Florida. The second row shows the SWIR1 time series for a single pixel (identified by the arrow in the image) that illustrates spectral changes as the land cover transitions from Tree Cover to Herbaceous cover, and finally to Developed cover as the construction proceeds. The third row shows a time series of GLanCE land cover for an area that was flooded for a dam in Minas Gerais, Brazil in 2007, and the fourth row shows the SWIR1 surface reflectance time series for a pixel at the center of these maps (also identified by an arrow) that was initially herbaceous, was then flooded, but returns to land as shoreline after water levels in the dam recede.
More generally, qualitative inspection of the V1.0 GLanCE results regarding land cover change suggests three key conclusions. First, for most locations, especially in humid regions outside of the tropics where land cover change is not extensive, the V1.0 land cover maps are quite stable and robust. Second, in locations where land cover change occurs and is permanent (or at least semi-permanent, e.g., Figure 8), the GLanCE map products generally capture changes with good realism. Third, the V1.0 GLanCE results tend to over-estimate land cover change in semi-arid regions that include a mix of Shrublands, Herbaceous Cover, and to varying degrees, Barren/ Sparsely Vegetated land cover. This issue is at least partly caused by the challenges we discuss above, i.e., the GLanCE data products show spurious changes that are a direct outcome of assigning unitary class labels to pixels that are mixed. This challenge is further complicated by year-to-year variation in phenology and greenness in semi-arid systems, which can be significant. Previous work has used hidden Markov models (Abercrombie and Friedl, 2016; Hermosilla et al., 2018) and Bayesian methods (Cardille and Fortin, 2016) that post-process classification results to reduce spurious interannual variation in map labels arising from classification uncertainty or interannual variation in phenology and greenness. We are currently exploring solutions based on these approaches. Preliminary results suggest that the combination of these post-processing methods, in combination with the use of Landsat Collection 2 imagery, which improves both the density and geolocation of available Landsat observations at each pixel, helps to reduce overestimation of change in regions where spurious change events are most frequently mapped. Because provision of information related to land cover change is a core goal of the GLanCE data product, resolving these issues and creating data sets with high-quality representation of change is a core priority of the project moving forward.

CONCLUSION
Maps of land cover provide essential information related to the ecological state of land areas and how they have been modified by humans. Hence, accurate information related to land cover is essential for both managing natural resources and for understanding society's ecological, biophysical, and resource management footprint. In this paper, we describe methods, data sets, and early results from a 5-year effort with the goal of using the global record of Landsat observations to create annual maps of global land cover from 2001 to 2020 at 30 m spatial resolution. Initial results indicate that the data sets and methods we are using to generate land cover maps are effective, but that the maps have somewhat lower accuracy for specific classes. In particular, and unsurprisingly, regions and classes at high latitudes and in semi-arid systems that tend to have mixtures of land cover below the 30 m spatial resolution of Landsat are more challenging to map and show modestly lower accuracy based on reference data than humid regions with uniform land cover. Moving forward, our efforts will focus on improving the quality of our land cover maps, expanding our reference data collection and accuracy assessment outside of North America, characterizing land cover change with well-defined accuracy, and, in the near future, releasing the data sets we are creating for use by the broader scientific community.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation. Data will also be made available via NASA's Land Processes DAAC and Google Earth Engine by the end of 2022.

AUTHOR CONTRIBUTIONS
MF, CW, PO, TL, and ZZ conceived the research. PA, RS, EB, K-TH, YZ, KRT, KT, KM, and JW developed data sets, implemented the mapping process, and performed the analysis. NG provided invaluable support porting CCDC to Google Earth Engine, and CSJr, provided the MapBiomas data. MF wrote the paper, but all authors reviewed the manuscript and contributed edits and text to the paper.

FUNDING
This research was supported by the NASA Making Earth System Data Records for Use in Research Environments (MEaSUREs) program, grant number 80NSSC18K0994.