Earth Science Data Records of Global Forest Cover and Change

tree to 60% in areas tree cover 0.35. and omission of were to tree cover. The omission was 45% and error was < 70% in areas with 0.3 - 0.6 tree cover but > 50% in high or low tree cover.


Approach
Global, spatially and temporally comprehensive forest-cover change Earth Science Data Records were inferred from high-(30-m) and moderate-(250-m) resolution satellite data. At 30-m spatial resolution, forest cover and changes in and between 1990, 2000, and 2005 were mapped using enhanced Global Land Survey (GLS+) data sets, supplemented with additional images where and when the GLS data were incomplete or inadequate for analysis (Tucker et al. 2004, Gutman et al. 2008. This effort also included production of surface reflectance ESDRs at 30-m resolution for 1990, 2000, and 2005, as well as fragmentation products based on the FCC records. (Note that the years 1990, 2000, and 2005 for all fine-resolution data sets refer to nominal years throughout this proposal, but the actual acquisition year of the GLS+ data set varies from place to place due to cloud cover and image availability.) The fine-resolution ESDRs were produced using algorithms that have been implemented or are now implemented in the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS), which was developed through previous NASA projects and includes algorithms for geometric orthorectification, radiometric normalization, and data quality screening. Atmospherically corrected surface reflectance, which is the basis for many other ESDRs and analyses, was generated as an intermediate product.
For each year from 2000 to 2005, an enhanced moderate-resolution change product was generated as a secondary record of forest-cover change. We also generated products to quantify and monitor fragmentation.
Efforts were restricted to mapping per-pixel gains and losses of forest cover between the epochs at fine spatial resolution and between years for moderate spatial resolution. Also, we restricted our definition of FCC exclusively to changes in forest cover and not to any change in the type of forest land use (cf. FRA 2000). Even within forest cover per se, there are many other types of changes-e.g., selective logging-that are also important for many science and land-management applications (e.g. Muchoney & Haack 1994, Olsson 1994Asner et al. 2005), but a global analysis of these is not yet feasible.
Like any ESDR, the data produced contain uncertainty, but this 15-year record represents a major advance in our understanding of Earth's changing forest cover. In processing the fineand moderate-resolution data sets, we ensured that the data provide coverage of the greatest extent possible and are internally consistent and that errors and uncertainty are thoroughly characterized. Global Land Cover Facility

Significance
These Earth Science Data records provide the first and only consistent, global record of forest cover changes documenting the period from 1990 to 2005, and they enable the first comprehensive assessment of Earth's forest cover at a scale appropriate to recent changes. The data also provide the basis for understanding impacts of forest change on the Earth system, including carbon budgets and the hydrological cycle. The fine-resolution and global extent of the fragmentation products support habitat analyses and other ecological studies at scales ranging from local to global, which is particularly valuable to natural resource managers, especially those responsible for conserving biodiversity (Dudley et al. 2005;Hilli & Kuitunen 2005). The protected-area subsets of the forest change and fragmentation records allow assessment of local conservation effects as well as the broader effectiveness of international environmental and biodiversity agreements. The moderate-resolution products are of particular value to various modeling communities, especially those concerned with regional to global carbon modeling (Ojima & Galvin 1994, DeFries et al. 1999) and regional hydrological modeling (Band 1993, Sahin & Hall 1996, Bounoua et al., 2002. Completion of this project satisfies key components of the GOFC-GOLD requirements for fine-resolution products (Skole et al. 1997, Townshend et al. 2004) and forms a contributory activity to GOFC-GOLD through its Land Cover Implementation Team.

Primary data inputs 2.1 Landsat images 2.1.1 Enhanced Global Land Survey
The primary data sources for generating the fine-resolution ESDRs were the GLS Landsat image datasets centered around 1990, 2000, and 2005. The GLS is a partnership between USGS and NASA, in support of the U.S. Climate Change Science Program and the NASA Land-Cover and Land-use Change (LCLUC) Program. Building on the existing GeoCover dataset developed for the 1970s, 19901970s, , and 20001970s, (Tucker et al. 2004, the GLS was selected to provide wall-to-wall, orthorectified, cloud free Landsat coverage of Earth's land area at 30meter resolution in nominal "epochs " of 1990, 2000(Franks et al. 2009, Gutman et al. 2008. The GLS was intended to provide one clear-view image acquired during the peak growing season of each epoch for each World Reference System (WRS) scene. The 1990 epoch ranges from 1984 to 1997 and is composed of 7,375 Landsat-5 Thematic Mapper (TM) images from 1984 to 1997. The GLS 2000 is composed of 8,756 Landsat-7 Enhanced Thematic Mapper Plus (ETM+) images from 1999 to 2002. The GLS 2005 is composed of 7,284 gap-filled Landsat-7 images and 2,424 Landsat-5 TM images acquired between 2003 and 2008. In many cases, however, images had to be selected with a date outside this range, mostly due to lack of cloud-free images during the growing season (Franks et al. 2009, Gutman et al. 2008. Because images have been selected from somewhat different dates, there are variations in phenology which account for the patchiness of image mosaics in some locations (Kim et al. 2011;Townshend et al. 2012). Global Land Cover Facility The original GLS data set did not fully cover Earth's terrestrial surface in all epochs; gaps were filled to the degree possible with newly available images (Figure 1). A major hole over northern South America in 1975 was filled using Landsat images from the Brazilian National Institute for Space Research (INPE) orthorectified using our own modules. However, no data exist to fill an expansive coverage gap over central and eastern Siberia in the 1990 epoch. Smaller, isolated holes also persist where coverage is missing in one or several adjacent WRS tiles for individual epochs; we obtained the best available Landsat images to fill these gaps. Finally, GLS images acquired near or during the leaf-off season, which are not suitable for forest cover change analysis, were replaced with images acquired during the local "leaf-on" growing season to use in our forest cover change analysis, pending availability (Kim et al. 2011.

Phenological selection
A challenge in using GLS data sets for analysis is that many of the GLS images were acquired near or during leaf-off seasons. Because the spectral differences between leaf-on and leafoff deciduous forests can be great, automated FCC analysis based on leaf-off images can result in widespread, erroneous changes. Prior to classification and forest change analysis, each Landsat image was evaluated to determine its phenological suitability for forest-coverchange analysis. We used the NDVI temporal profiles calculated using the GIMMS AVHRR and MODIS data record (Tucker et al., 2005) to determine whether an image was acquired near or during leaf-off seasons. The GLS 1990, 2000 images were evaluated using the GIMMS record directly.

Orthorectification
Many non-GLS Landsat images were needed to supplement the GLS dataset to produce the fine-resolution ESDR products. Many of these non-GLS images were not orthorectified and might therefore have contained significant geolocation errors. We developed and implemented an orthorectification algorithm in the LEDAPS software that automatically orthorectifies a Landsat image to match the GLS data set (Gao et al. 2009). The module was used to orthorectify over 500 images in North America, ~30 images in Madagascar, and ~20 images in Africa, as well several SLC-off ETM+ images. During extensive validation, residual misregistration errors in the orthorectified products were found to be less than 1 pixel.

Digital Elevation Model: ASTER GDEM (v2.0)
We used the Global Digital Elevation Model, version 2.0 (GDEM v2.0) as an ancillary layer in many analyses. Produced from images acquired by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) the GDEM dataset was jointly released by the Ministry of Economy, Trade, and Industry (METI) of Japan and NASA. The first and second versions of the ASTER GDEM were released in June 2009 and October 2011 respectively. The 30-meter resolution ASTER GDEM was generated using stereo-pair images collected by the ASTER instrument onboard the Terra satellite. The dataset is distributed in GeoTIFF format, spanning from 83S to 83N.

MODIS VCF Tree Cover Layer
The MODerate-resolution Imaging Spectroradiometer (MODIS) Vegetation Continuous Fields (VCF) Tree Cover dataset, Version 5, was produced at 250-m resolution globally from 2000 to 2010(DiMiceli et al 2011. In contrast to methods based on linear mixture models (e.g., DeFries et al. 1999, Asner et al. 2005, the MODIS VCF is based on a flexible regression tree algorithm, which is more capable of incorporating empirical information to improve correlation of estimates to measured tree cover. Although the MODIS Tree Cover VCF has been used for a wide range of continental-to global-scale assessments, many land cover changes occur in patches beneath its 250-m resolution (Townshend and Justice 1988). Higher-resolution continuous-field datasets had been generated for limited areas based on Landsat data (e.g., Homer et al. 2004), but there were currently no global datasets representing tree cover at resolutions finer than that of the MODIS sensor.
The spatial and thematic scale of the MODIS VCF and other continuous-field datasets (e.g., Asner et al. 2005 have made reference data difficult to acquire, and so quantitative error estimates of these datasets are quite limited. Hansen et al. (2002) provided the first de facto-although not independent-estimates of MODIS VCF accuracy by comparing an experimental version of the dataset to the Landsat data used to train the generating model. Later, White et al. (2005) compared the MODIS VCF Version 1 to independently gathered field data across the arid southwestern US, and Montesano et al. (2009) validated the Version-4 MODIS VCF against independent reference data derived from photo-interpreted high-resolution images across the boreal-taiga ecotone. Also, Heiskanen et al. (2008) and Song et al. (2011) compared the MODIS VCF to other remotely sensed global datasets. Global Land Cover Facility Across all biomes and types of reference data, these independent assessments found that saturation of the optical signal, phenological noise, and confusion with dense herbaceous vegetation led to errors in the MODIS VCF between 10-31% Root-Mean-Squared Error (RMSE), over-estimation in areas of low cover, and under-estimation in areas of high cover.

Introduction
Reflectance is defined as the fraction of incident radiance within a specified interval, or band, of the electromagnetic spectrum that is reflected (i.e., neither absorbed nor from a target. Directional surface reflectance is further specified as the ratio of the radiance reflected from a surface to the incident radiance incoming from a direct source of illumination in a given infinitesimal solid angle. Estimated by atmospheric correction of satellite images, directional reflectance ideally decouples the surface properties from the atmospheric signal, thus representing the value that would be measured by an ideal sensor held just above the Earth's surface at a given solar and viewing geometry and without any atmospheric effects. Directional surface reflectance is the most basic remotely sensed surface parameter in the solar-reflective wavelengths and therefore provides the primary input for essentially all higher-level surface geophysical parameters, including vegetation indices, albedo, Leaf Area Index (LAI), Fraction of absorbed Photosynthetically Active Radiation (FPAR), burned area, land cover and land cover change. Directional surface reflectance is also directly used in various applications to visually or quantitatively detect and monitor changes on the Earth's surface. Because they enable other comparisons among data imaged under various illumination and atmospheric conditions, reflectance data products have value independently of their utility for monitoring forest cover change. For example, the ESDR Community White Paper on Surface Reflectance  notes that validation of global reflectance data sets from AVHRR, MODIS, and VIIRS will need to rely on reflectance products derived from high-resolution sensors.
Nearly half of the original GLS-1990 dataset did not have correct radiometric gain and bias coefficients at the time of data acquisition; thus atmospheric correction and conversion to surface reflectance were not possible (Chander et al. 2003(Chander et al. , 2009Townshend et al. 2012). These un-calibrated GLS images were replaced after the original GLS compilation with substitutes from the updated USGS archive within the epoch wherever possible (Figure 1). To perform the selection of replacement imagery while minimizing phenological or atmospheric noise, a tool was constructed to query the USGS Global Visualization Viewer (GloVis) database (glovis.usgs.gov/) for appropriate images based on phenological time series of Normalized Difference Vegetation Index (NDVI) from the MODerate-resolution Spectroradiometer (MODIS) (Kim et al. 2011;Townshend et al. 2012). Global Land Cover Facility Each image of this enhanced GLS dataset was then atmospherically corrected to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al. 2006b). Atmospheric inputs and parameterization of LEDAPS are described by Feng et al. (2013). The surface reflectance data set from the enhanced version of GLS-1990 is available from the Global Land Cover Facility and use of these data is strongly recommended for studies based on the GLS-1990 data. Clouds were identified in a spectral-temperature space (Huang et al. 2010) and removed from subsequent analysis. This "aggressive" cloud-detection algorithm's low rate of omission error makes it suitable for masking pixels from forest-cover change analysis. Cloud shadows were identified by projecting cloud masks onto a digital elevation model through solar geometry at the time of image acquisition (Huang et al. 2010) and were also removed from analysis.

Algorithm
3.1.1.2.1 Radiometric calibration and estimation of top-of-atmosphere reflectance The Landsat-7 ETM+ instrument has been carefully calibrated and monitored since launch in 1999, and the calibration has been stable since shortly after launch . The Landsat-5 calibration history has recently been updated (Chander & Markham 2003, Chander et al. 2009) and is compatible with subsequent Landsat-7 ETM+ data. LEDAPS uses updated calibration histories to convert 8-bit quantized Landsat data to at-sensor radiance and then to top-of-atmosphere (TOA) reflectance using solar geometry and instrument band pass.
3.1.1.2.2 Atmospheric correction to estimate surface reflectance Atmospheric correction seeks to estimate surface reflectance by compensating for the scattering and absorption of radiance by atmospheric constituents. In practice, atmospheric correction is typically achieved by inverting a highly parameterized model of atmospheric radiative transfer coupled to a surface reflectance model. For speed and simplicity, the reflecting surface is often assumed to be Lambertian. Atmospheric radiative transfer modeling is relatively mature, and so several methods may be used to model the surface/atmosphere interaction (e.g., Successive Order of Scattering, Doubling adding). The main challenge to the operational implementation of these models lies in the assignment of the atmospheric parameters and the a priori knowledge of the surface BRDF -strictly necessary for a full inversion. Approaches to operationally retrieving the atmospheric parameters have advanced considerably in the last 10 years as remote sensing instruments capable of retrieving atmospheric properties (aerosol, ozone, water vapor, etc..) have been put into operation. In the absence of operational retrievals, atmospheric climatology or forecasted values can be applied, although product accuracy degrades considerably without coincident atmospheric measurements. The determination of surface BRDF is currently practical operationally only for satellite sensors with single-pass multi-angular capability, such as MISR or POLDER. Thus, the uncertainty introduced by surface BRDF was assumed to be constant inter-annually and to not have significant influence on analyses at this temporal scale. Global Land Cover Facility The atmospheric perturbation of the directional surface reflectance signal depends on the type and characteristics of atmospheric particles interacting with the radiation. Atmospheric gas molecules (N2, O2, O3, H2O, CO2, etc.) scatter radiation according to Rayleigh's theory (i.e., molecular scattering) and absorb radiation over the spectrum varying by species. These specific scattering effects are governed by atmospheric pressure and the vertical temperature profile. Aerosols (i.e., suspended particles ranging from about 10 -3 μm to about 20μm) scatter and absorb radiation according to the Mie and Geometric Optics theories; the former applies to aerosols with diameters on the order of the radiation's wavelength, and the latter idealizes particles larger than the wavelength of radiation as individual spheres with given real and imaginary refractive indices.
Atmospheric correction removes or reduces the effects of these atmospheric perturbations.
In an idealized case of a Lambertian surface (i.e., with angularly uniform reflectance) and in a narrow spectral band (here referred to with the index i ) outside of the main absorption feature of water vapor, the top-of-atmosphere signal can be written as (Vermote et al. 1997): ρTOA is the reflectance at the top of the atmosphere; Tg is the gaseous transmission by a gas species (g), e.g., water vapor (TgH 2 O), ozone (TgO 3 ), or other gases, TgOG (e.g. CO2…); ρatm is the atmosphere's intrinsic reflectance; Tratm is the total atmospheric transmission (downward and upward); Satm is the atmosphere's spherical albedo; A is the atmospheric pressure, which influences the number of molecules and the concentration of absorbing gases in the radiation's path; τA, ω0 and PA describe the aerosol properties and are spectrally dependent: τa is the aerosol optical thickness; ω0 is the aerosol single scattering albedo; PA is the aerosol phase function; UH 2 O is the integrated water vapor content; UO 3 is the integrated ozone content; m is the air-mass, computed as 1/cos(θs)+1/cos(θv); and ρS is the surface reflectance to be retrieved. The geometrical conditions are described by the solar zenith angle (θs), the viewing zenith angle (θv), and by Φ, the difference between θs and θv. The effect of water vapor on the intrinsic atmospheric reflectance is approximated as: where ρR represents the reflectance of the atmosphere due to Rayleigh scattering and ρR+Aer represents the reflectance of the mixing molecules and aerosols. Accounting correctly for mixing and coupling effects is important for achieving high accuracy in modeling the atmospheric effect. Eqn.
(2) conserves the correct computation of the coupling and assumes that the water vapor is mixed with aerosol and that the molecular scattering is not affected by water vapor absorption.
The transmission, intrinsic reflectance, and spherical albedo terms are computed using the vector version of the 6S radiative transfer code (Kotchenova et al. 2006). Since the cost of running 6S for each pixel would be prohibitive, 6S was run early in the process to generate a look up table (LUT) accounting for pressure, water vapor, ozone, and geometrical conditions over the whole scene for a range of aerosol optical thicknesses. The LUT was created for each TM band and was used both in the aerosol retrieval process as well as in the correction step at the end.
Ozone concentrations were derived from Total Ozone Mapping Spectrometer (TOMS) data aboard the Nimbus-7, Meteor-3, and Earth Probe platforms. The gridded TOMS ozone products are available since 1978 at a resolution of 1.25º longitude and 1.00º latitude from the NASA GSFC Data Active Archive Center (DAAC). In cases where TOMS data were not available (e.g., 1994-1996), NOAA's Tiros Operational Vertical Sounder (TOVS) ozone data were used. Column water vapor was taken from NOAA National Centers for Environmental Prediction (NCEP) reanalysis data available at a resolution of 2.5 by 2.5 degrees over the Landsat era. Digital topography (1 km GTopo30) and NCEP sea-level surface pressure data were used to adjust Rayleigh scattering to local conditions.
Like other atmospheric correction schemes for MODIS and Landsat, the Dark, Dense Vegetation (DDV) method Remer et al. 2005) was used to infer aerosol optical thickness (AOT) from each image. Based on the correlation between chlorophyll absorption and bound water absorption, this method postulates a linear relation between surface reflectance in the atmospherically insensitive shortwave-infrared (2.2 μm) and surface reflectance in the affected visible bands. The method then uses this relation to calculate surface reflectance for the visible bands and estimate aerosol optical thickness by comparing the result to the TOA reflectance. For LEDAPS AOT estimation, each image was averaged to 1-km resolution to suppress local heterogeneity, and candidate "dark targets" of TOA reflectance were selected. For these targets, correlation was assumed only between the blue (0.45-0.52) and SWIR (2.2μm) bands, such that water targets were excluded. The Global Land Cover Facility specific relation was derived from an analysis of data from Aerosol Robotic Network (AERONET) sites where AOT is measured directly. The calculated AOT in the blue wavelengths was propagated across the spectrum using a continental aerosol model. A "sanity check" for the aerosol was performed by analyzing the surface reflectance derived in the red band for each 30-m pixel contained in the 1-km grid cell; if too many "unphysical" values were found, the aerosol retrieval at this 1-km location was rejected. The valid aerosol optical thicknesses at 1 km were interpolated spatially between the dark targets using a spline algorithm. The interpolated AOT, ozone, atmospheric pressure, and water vapor were supplied to the 6S radiative transfer algorithm, which then inverts TOA reflectance for surface reflectance for each 30-m pixel.
As noted above, water targets were excluded from the aerosol retrieval. However, interpolation of valid (i.e., land) aerosol targets occurs across the entire scene. Thus, the surface reflectance of small lakes surrounded by land was likely to be reasonable, while the reflectance of open ocean water (far from any valid aerosol target) was likely to be problematic.

Cloud-and shadow-masking
Removing pixels contaminated by clouds and their shadows was necessary to avoid erroneous retrieval of surface reflectance and false detection of forest-cover change. LEDAPS implemented two cloud masks -a version of the Landsat Automated Cloud Cover Assessment (ACCA) algorithm (Irish, 2000) and a more aggressive mask based on MODIS spectral tests (Ackerman et al. 1998). Shadows were located from the latter using solar geometry and an estimate of cloud height based on the temperature difference between known cloudy pixels and NCEP surface temperature. A third cloud-masking algorithm has been developed by Dr. Vermote through his USGS-funded Landsat Science Team project -"A Surface Reflectance Standard Product for LDCM and Supporting Activities". Quality Assessment codes for this algorithm are listed in Table 1. Finally, an automated cloud and shadow masking algorithm has also been developed by Huang et al. (2010) as part of the TDA-SVM algorithm.

Validation
Landsat surface reflectance products were validated in two ways. Internal aerosol optical thickness (AOT) estimates retrieved by LEDAPS have been compared to measurements taken at Aerosol Robotic Network (AERONET) observations (Masek et al. 2006b), and surface reflectance was compared to simultaneously acquired MODIS daily reflectance and Nadirand BRDF-Adjusted Reflectance (NBAR) images (MOD09 and MOD43, respectively) (Feng et al 2013). These paired validations provide an internal check on a driving parameter of the LEDAPS algorithm (AOT), as well as a consistency check against the thoroughly calibrated and validated MODIS product. Internal cloud mask (1=cloudy, 0=clear) 9 Cloud shadow 10 Snow mask 11 Land/water mask based on spectral test 12 Adjacent cloud 13-15 unused 3.1.1.3.1 Comparison of retrieved AOT to AERONET measurements Aerosol Robotic Network (AERONET) sites measure and record aerosol properties across the globe, with records at some sites extending back to the early 1990's (Holben et al. 1998). Aerosol optical thickness estimates from pixels processed through LEDAPS were compared to coincident measurements from 21 of these AERONET sites (Table 2, Figure 2). All AOT values reported are for the blue wavelengths. Results suggest reasonable agreement with AERONET observations, and the discrepancies between LEDAPS and MODIS reflectance products were generally within the uncertainty of the MODIS products themselves-the greater of 0.5% absolute reflectance or 5% of the retrieved reflectance value. Spatial patterns for the sites suggested that land cover type may influence the aerosol retrievals ( Figure 3), although this artifact was slight in comparison to the direct effect of reflectance itself and therefore appears to have little impact on the retrieved surface reflectance values.   (Table 3), the MODIS sensor aboard the Terra platform follows the same orbit and crosses the equator roughly 30 minutes behind Landsat 7. MODIS surface reflectance data products (MOD09) have been calibrated and validated comprehensively (Vermote et al. 2002, Kotchenova et al. 2006, Vermote and Kotchenova 2008) and may be used as a reference to validate Landsat surface reflectance products ).
We developed an online tool for validating Landsat surface reflectance estimates against coincident MODIS estimates and used it to validate the 2000-and 2005-epoch SR products. Initial tests for WRS-2 scenes over eastern Africa showed strong agreement between Landsat-7 ETM+ and MODIS surface reflectance products, with the majority of R 2 values above 0.9 ( Figure 4). Landsat scenes with R 2 values below 0.8 were inspected individually, revealing explanations for the discrepancy. Of the poor quality images, one was corrupted and others were either cloudy or dominated by ocean.

Tree Cover (2000 and 2005)
Spatio-temporal estimates of tree-canopy (or simply "tree") cover provide a biophysically relevant, sensible, and consistent basis for monitoring forest cover and change (Sexton et al. 2016). The following algorithm and its results have been peer-reviewed and are described by Sexton et al. (2013b).

Model
Tree cover (C) was estimated as a piecewise-linear function of surface reflectance and temperature: where X is a vector of surface reflectance and temperature estimates, ε is error in the estimates produced by f() applied to X, subscript i denotes the pixel's location in space, indexed by pixel, and t refers to its location in time, indexed by year. Continuous measurements, such as percent cover and surface reflectance, are robust to changes in resolution (Gao et al. 2006, Feng et al. 2013; although the data were derived from Landsat; the model makes no specification of scale and thus may be calibrated and applied at arbitrary, even different, resolutions between those of Landsat (30 m) and MODIS (250 m).
To estimate tree cover at 30-m resolution in 2000 and 2005, MODIS-based, 250-m tree cover estimates were overlaid on rescaled Landsat surface reflectance layers in each year, and a joint sample of cover and reflectance variables was drawn to generate a training dataset for each Landsat scene in each epoch ( Figure 5). (Throughout this section, we refer to the data used to estimate model parameters as "training" data, and we refer to data whose accuracy is assumed as "reference" data.) The model was thus fit locally to each scene of the Landsat tiling system of WRS-2 in each epoch. The model was fit using the Cubist™ regression tree algorithm and applied using CubistSAM, an open-source parser for Cubist (Quinlan 1993). Except for an allowance for extrapolation within the range [0,100], our application of regression trees was standard (i.e., neither sample boosting or bagging nor ensemble "random forests" or "committee models" were employed). Cubist -as well as regression trees in general -has been found to provide accurate estimates of percent-scale land cover attributes in numerous studies (e.g., Sexton et al. 2006Sexton et al. , 2013a. Because regression trees can over-fit the data and there are often few data points at the extremes of the range of the response variable (e.g., tree cover), Cubist gives an option for either estimating within the range of the response variable at each node (the default) or extrapolating within a specified range. To avoid over-fitting to the sometimes small samples at terminal nodes with extreme cover values, we allowed for extrapolation within the range of 0-100% tree cover. The fitted model was then applied to Global Land Cover Facility the original, 30-m Landsat data in order to estimate tree cover at the Landsat spatial resolution. Random errors (i.e., those which were not systematic, e.g., bias) were minimized by using the six-year median of cover for each pixel. Land-cover changes between 2000 and 2005 were removed by calculating Global Land Cover Facility the standard deviation of annual tree cover estimates for each pixel over that interval and removing pixels in the top 10% of the distribution of standard deviations of each Landsat scene. Because only six years of MODIS VCF data were available, we used the median, which is a better representation of central tendency than the mean in small samples such as the six values of cover from 2000-2005.
Pure (i.e., 0% or 100%) and near-pure pixels are rare in the MODIS data, and tree cover tends to be over-estimated in areas of low cover, especially agricultural fields. To ameliorate under-representation of low tree-cover in the training sample, we augmented the MODISderived reference data with information from the Training Data Automation and Support Vector Machines (TDA-SVM) automated classification algorithm (Huang et al. 2008) and the MODIS Cropland Probability Layer (Pittman et al. 2010). Cropland Probability and Tree Cover images were overlaid within each Landsat scene, and Landsat pixels with crop probability > 0.5 and tree cover < 50% were selected. This selection comprised Landsat pixels with either crop or sparse vegetation cover. Within the selection, Landsat pixels identified by TDA-SVM as "non-forest" in both 2000 and 2005 were assumed to be sparsely vegetated and were labeled as 0% tree cover. The remaining (i.e., crop) pixels in the selection were ranked by their NDVI values and divided into three sub-strata: high, medium, and low NDVI. Pixels from each of these sub-strata were randomly sampled such that the maximum proportion of Landsat "crop" pixels was the proportion of MODIS pixels within the scene whose crop probability was > 60%. All of the sparsely vegetated pixels and the sample of crop pixels were then pooled with the MODIS-based reference data to form an ensemble training sample of tree cover and reflectance.

Post-processing
3.1.2.2.1 Water mask Surface-water bodies were masked from the tree and forest-cover & change data, and the surface-water layer is a useful input to many other applications. The following algorithm is described by Feng et al. (2015).
Water cover was defined as a state of the landcover domain c ϵ C, and its probability of occurrence in each pixel was modeled as a function of reflectance and topographic covariates (X): where f is a binary decision tree fit by the See5™ algorithm (Quinlan 1986(Quinlan , 1993. The topographic covariates were elevation and slope derived from the ASTER GDEM (Tachikawa 2011), reflectance covariates were Landsat Band-5 (SWIR) surface reflectance, the Normalized-Difference Water index (NDWI) (McFeeters 1996): to distinguish water from other cover types, as well as the Normalized Difference Vegetation Index (NDVI) (Tucker et al. 2005) to distinguishes water from vegetation specifically. The optimal threshold of each index for separating water varies regionally and over time due to mixing and local similarities with other cover types (Ji et al. 2009;Jiang et al. 2014).
Water was detected in each 30-m Landsat pixel with a classification-tree model (Quinlan 1986) parameterized through an automated, two-stage procedure. An initial, deductive stage identified reference water pixels of varying certainty by comparing multi-spectral water and topographic indices to coarse-resolution (MODIS) water estimates. This stage leveraged prior knowledge with multiple sources of independent information to stratify the decision space into regions of possible water with varying degrees of certainty. An inductive stage then optimizes rules based on high-resolution estimates of surface reflectance, brightness temperature, and terrain elevation.
The first stage of classification generates local reference data with varying levels of certainty. The pixels, identified as water by multi-spectral indices, were compared with a priori water pixels resampled from the 250-m resolution MODIS water mask to the spatial resolution and extent of each Landsat image. This comparison resulted in four possible levels of certainty, through which weights were assigned to each reference datum (Table 4).
Topographic, spectral, and brightness temperature variables were first stratified into generic cover types: water, land, snow and ice, and cloud. A loose and a strict threshold-equaling -0.1 and 0.1-were applied each to NDWI and MNDWI to distinguish water with low and high certainty. Terrain shadows were identified as pixels with hill-shade value <150 (on a scale from 0 to 255) and slope >20 degrees, as discussed in section 2.1.2.
Snow and ice show high reflectance values in the visible and NIR bands and low reflectance in SWIR bands, leading to high MNDWI but low to moderate NDWI. A strict difference threshold (0.7) was used to reduce confusion of water with snow and ice, and a criterion of brightness temperature <1.5 ℃ was also included to further improve the discrimination: MNDWI > NDWI + 0.7 and 6 < 1.5 ℃ .
(8) Global Land Cover Facility Redundancy among multiple images was leveraged to maximize certainty in each location within each epoch. This is accomplished by a "best-pixel" compositing rule, taking the treecover estimate with the lowest estimated uncertainty from those available.
For each location (x,y) in each epoch (t), there could be any number, k ϵ K= [1,2…n], of cover and uncertainty estimate-pairs ( , ) , , , . The "best-pixel" approach takes the least uncertain estimate of cover, as well as its corresponding estimate of uncertainty: For a static, continuous variable (e.g., tree cover), is quantified as the root-mean-square error of the estimate. For a static, categorical variable (e.g., forest cover), is quantified as the complement of the probability of class membership--i.e., 1-p( ).
This selection was applied up to twice for each (x,y,t) pixel: first if there were multiple Landsat images available for a scene ("within-scene compositing") and second if multiple WRS-2 scenes overlapped at that pixel location ("sidelap compositing"). Missing estimates in any of the contributing images (due, e.g., to clouds, cloud-or terrain-shadows, or scan-line gaps) were treated as having maximum uncertainty so that pixels were filled with clear-view estimates wherever available.

Methods
The uncertainty of the tree-cover estimate in every pixel was assessed relative to the training data by ten-fold cross-validation. Pixel-level uncertainty was quantified at each terminal node of the regression tree and assigned to pixels identified with that node. Because these pixel-level uncertainties were assessed only relative to their training data, errors between the reference data and actual cover were not included at the pixel level. As described in a later section, training (MODIS) and output estimates were compared to Global Land Cover Facility i=1 =1 approximately coincident measurements derived from small-footprint lidar measurements in order to assess their accuracy relative to more direct measurements of actual cover. (We use the term "measurement" to refer to lidar-derived values of cover -which are calculated without statistical inference -and the more general "estimate" to refer to values derived statistically from MODIS and Landsat images.) All comparisons were made at 250-m resolution, using MODIS estimates from 2005 and Landsat estimates from the 2005 epoch. Preliminary analyses comparing Landsat estimates to lidar measurements at 30-m resolution were not appreciably different than those reported here, although there was a small reduction of correlation believed to be due to spatial misregistration of Landsat data. Uncertainty metrics were based on average differences between paired model and reference (or training) values (Willmott, 1982), quantified by Mean Bias Error (MBE), Mean Absolute Error (MAE), and Root-Mean-Squared Error (RMSE): where Mi and Ri are estimated and reference tree cover values at a location i in a sample of size n.
After modeling the relationship between M and R by linear regression, their (squared) difference was disaggregated into systematic error (MSES) and unsystematic error (MSEU) based on the modeled linear relationship (Willmott 1982): where is the cover value predicted by the modeled relationship (Willmott 1982). Accuracy is thus quantified by the difference between the trend of model over reference cover, and precision is quantified by the variation surrounding that trend. MSES and MSEU sum to Mean-Squared Error (MSE), and therefore: (Willmott 1982). To maintain consistency, we report the square roots of MSES and MSEU, i.e., RMSES and RMSEU, in units of percent cover.  The Costa Rica site is dominated by tropical moist broadleaf evergreen forest surrounded by livestock pastures. The Utah site is an ecotone of temperate evergreen needle-leaf conifer forest, deciduous broadleaved shrubland, and annual grasses. The California site is dominated by tall, mixed-species temperate evergreen conifer forests of varying cover. The Wisconsin site is dominated by a mixture of temperate deciduous broadleaf hardwood and coniferous needle-leaf tree species with significant coverage of herbaceous agriculture, including corn. All lidar measurements were acquired during the growing season of each respective site, with mean point densities > 1 return/m 2 . The Costa Rica dataset, collected in 2006, is described by Kellner et al. (2009), and the Wisconsin dataset is described by Cook et Global Land Cover Facility al. (2009). Figure 7 shows an example of the 3-dimensional distribution of lidar measurements in the California site. All sites were assessed visually for obvious changes in cover between data acquisitions; in the WI dataset, obvious cover changes due to forest harvesting between Landsat and lidar acquisitions (totaling 21 pixels) were delineated manually and removed. Tree cover (C) was calculated from lidar returns by dividing the number of returns above a criterion height by the total number of returns within a 10-m radius:

Reference data
where n is the number of returns and nh is the number of returns above the specified height (h) (Korhonen et al. 2011). In accordance with the International Geosphere-Biosphere definition of forests, we specified the criterion nh = 5 meters. Following calculation of tree cover at 10-m resolution, rasters were aggregated to 250-m resolution by averaging the values within the extent of each 250-m pixel. In pixels with steep underlying terrain (as might be likely especially in CA and UT), the varying ground elevation in large pixels can cause spurious detection of tree cover as lidar returns above 5-m height; first computing cover in small, 10-m pixels and then aggregating to 250-m pixels avoided this possibility. Also note that Relative Height (i.e., RH100) and other waveform-based metrics (Hyde et al. 2005, Dubayah et al. 2010 were not used; only height of the (discrete-return) lidar posts was used to calculate canopy height. Global Land Cover Facility

Consistency of Landsat-and MODIS-based (VCF) tree cover estimates
The relationship between Landsat estimates of tree cover and the MODIS data on which they were based was very strongly linear, near parity, and consistent among biomes ( Figure  8, Table 5, Table 6). Relative to the MODIS estimates, Landsat estimates exhibited MBE of -6%, MAE of 8%, and RMSE of 10% cover (   Table 5) in the biome samples of 2005 data. The modeled linear relationship explained 88% of the variation between the two datasets, and RMSE was equally partitioned between systematic and random components, with both RMSES and RMSEU equaling approximately 7% cover (Table 7). Although significantly different from zero, the intercept of the linear relationship was relatively small (4.5%).
The    RMSEu is mean "unsystematic", or "residual" error between original and calibrated measurements, and RMSEs is the "systematic" error remaining between calibrated and reference measurements (see text for full explanation). Unless otherwise noted, all coefficients are significant at Pr(>|t|) < 0.01 2 Accuracy of Landsat-based tree cover estimates relative to lidar reference data Across the four sampled biomes, the correspondence of Landsat-based estimates of tree cover to reference lidar measurements was similar to the relationship between MODISbased estimates and lidar-based measurements (Figure 10). Across the biomes, RMSE of Landsat estimates relative to lidar-measured cover was 17%, with MAE of 15% and MBE of -11% cover (Table 6). However, the overall linear relationship between Landsat estimates and lidar measurements was stronger (R 2 = 0.81) than that of MODIS estimates relative to lidar measurements (R 2 = 0.71). This strong linear trend resulted in a greater dominance of systematic (RMSEs = 15%) over unsystematic, or random noise (RMSEU =9%) in the Landsat estimates compared to MODIS, suggesting a greater potential for empirical calibration of Landsat estimates than is possible for the MODIS dataset. Although still present, saturation of Landsat estimates relative to lidar measurements was reduced slightly compared to the saturation seen in MODIS-based estimates. Landsat estimates reproduced the spatial pattern of tree cover in most sites with greater fidelity than did MODIS estimates ( Figure 10). The exception to this was the UT site, where there was no clear correspondence between either Landsat or MODIS estimates and the lidar measurements. Another artifact shared in both the Landsat and MODIS data was the slight compression of the actual frequency distribution of values, such that there were more intermediate values and correspondingly fewer values near the extremes of cover (i.e., 0 and 100%). It should be stressed, however, that even considering the minor artifacts, Landsat estimates resolved greater spatial variation in tree cover than did the relatively coarse MODIS estimates.

Definitions
Based on the global land cover classification system developed by the International Geosphere Biosphere Programme (IGBP) (Belward & Loveland 1996), forest is defined as a minimum area of land of 0.27 hectares with ≥30% tree cover-i.e., as land cover, as opposed to land use (Sexton et al., 2016;Townshend et al., 2012). Only net forest/non-forest typeconversion changes were included in the fine-resolution FCC ESDR products. We defined "forest gain" as categorical change from non-forest to forest and "forest loss" as change from forest to non-forest. Stasis of forest or non-forest classification in a pixel over a period was defined respectively as "persistent forest" and "persistent non-forest".

Forest cover and change from 2000-2005
The following algorithm and its results have been peer-reviewed and are described by Sexton et al. (2015).

Defining forest cover in terms of tree cover
"Forest" is defined as a class of land cover wherein tree (-canopy) cover, c, exceeds a predefined threshold value, c * . The probability of belonging to "forest", p(F), is therefore the probability of c exceeding the threshold c * (Figure 11)-i.e., the integral of the density function of c above c * : Complementarily, the probability of membership in non-forest is simply 1-p(F).
In any location i, tree cover ci is estimated by a model f of remotely sensed variables X (Hansen et al. 2003, Homer et al. 2004, Sexton et al. 2013b): Global Land Cover Facility where β is a set of empirically estimated parameters, and ε is residual error.

Figure 11. Estimation uncertainty of tree and forest cover within a pixel, modeled as a normal probability density function of tree cover, with probability of forest (shaded) and non-forest (unshaded) defined relative to a threshold of tree cover, c*.
Given a joint sample of locations i = [1,2,…n] with coincident true and estimated values of a continuous variable such as tree cover (ci, i), error may be quantified as the Root-Mean-Square Error (RMSE), which for large samples approximates the standard deviation of estimates of the true value of cover: −1 Thus, given ci, and an estimator (e.g., linear regression) producing estimate and rootmean-square error σi = σ, a Normal probability distribution of possible values of ci may be assumed (Snedecor and Cochran 1989, Hastie et al. 2001, Clark 2007:

Global Land Cover Facility
Given paired estimates of cover and its RMSE, this model provides a probability density function of tree cover p(c) (Eqn. 13) and therefore the probability of identifying forest for each pixel i (Eqn. 10).

Change detection based on bi-temporal class probabilities
Given the probability of detecting forest in a location i = (x,y) at each of two times t, four dynamic classes (D) are possible: stable forest (FF), stable non-forest (NN), forest gain (NF), and forest loss (FN). Calculating the probability of each of these dynamics at that location simply requires calculating the following joint probabilities: where subscripts denote observation times ( Figure 12). In practice, the model of error is approximate, and so carets (^) denote that the resulting values are estimates. These joint probabilities sum to unity at each location i, and because they are merely transformations of the original cover and error values in every pixel, they may be mapped geographically without gain or loss of information from those estimates. In order to produce a categorical map of change classes, each pixel may be assigned either the most probable class at i, or some other criterion of probability may be set (e.g., p ≥ 0.9) to filter detection based on certainty of the tree-cover and derived forest-cover and -change estimates. Global Land Cover Facility

Forest cover and change from 1990-2000
The following algorithm and its results have been peer-reviewed and are described by Kim et al. (2014).

Forest-cover retrieval using stable pixels
We inferred forest cover in 1990 and change from 1990 to 2000 using a signature-extension approach based on stable pixels hindcast from 2000 and 2005 epochs ( Figure 13). For the purpose of large-area mapping, extrapolation of models beyond the immediate temporal and spatial domain in which they were trained has been explored by many researchers (Sexton et al. 2013b;Gray and Song 2013). Termed as "generalization" or "signature extension", this approach has been successfully applied for the classification of forest cover ) and change  using Landsat data. This approach also has been implemented by deriving training data from one date and using it to train a classifier on a different image from the same path/row scene but different acquisition date ). Complementary to the traditional signature extension method, Gray and Song (2013) combined a procedure to identify stable pixels to deal with irregular time-series images. This approach has been found to be effective for the automated classification of large areas, especially when there are actual changes in class spectral signatures from phenological variability, atmospheric differences, or land cover changes (Fortier et al. 2011, Gray andSong 2013). Global Land Cover Facility

Reference forest/non-forest data
Persistent forest (F) and non-forest pixels (N) were sampled from forest-cover change maps between 2000 and 2005 GLS epochs and then filtered so that only "stable" pixels-i.e., those whose class did not change between 1990 and 2000 epochs-were retained for analysis. The details of the filtering process are presented below.
For each WRS-2 scene, an annual rate of forest-cover (F) change, , and an annual rate of non-forest-cover (N) change, , were calculated as: where F and N are the percentage of forest and non-forest pixels, respectively, and t1 and t2 were respectively the acquisition years of the Landsat images for 2000 and 2005 GLS epochs.
The spectral difference (∆SR) -quantified as the Euclidean distance between two pixels over time in the spectral domain-was calculated for 1990-2000 (ΔSR1) and 2000-2005 (ΔSR2). To minimize impact from accelerating or decelerating rates of forest-cover change between two periods, a parameter α was defined as the ratio of the sums of spectral difference of all persistent pixels and was calculated as: Given the large number of available pixels within the overlapping portion of two Landsat images within the same WRS-2 scene, α was doubled to increase the selectivity of filtering Global Land Cover Facility for stable pixels. A percentage of forest equaling α x 2 x 100 × and non-forest pixels equaling α x 2 x 100 × were thus removed per year of difference between 1990-and 2000-epoch images in the order of spectral difference (∆SR). Limiting the sample to pixels that were stable from 2000 to 2005 minimized inclusion of erroneous data, and filtering the most spectrally different pixels from 1990 to the later epochs removed the pixels most likely to have changed over that period.

Forest-cover classification
Using the sample of stable-pixel locations, a forest/non-forest reference sample was extracted from forest-cover maps in 2000 and 2005. This sample was then filtered to maximize certainty and minimize change between observation periods ( Figure 13).
Forest cover in circa-1990 was retrieved by a classification-tree algorithm. The probability of forest cover, p(F), in each pixel i at time t ≈ 1990 was estimated by a conditional relationship (g) to remotely sensed covariates ( ): where is a vector of surface reflectance and temperature estimates; subscripts i and t denote the pixel's location in space, indexed by pixel, and time indexed by year. The relation g was parameterized using the C 5.0 ™ classification-tree software (Quinlan 1986), trained 9 on a sample of pixels within each Landsat image; the model was thus fit locally within each Landsat World Reference System 2 (WRS-2) scene. Reflectance and temperature covariates were acquired from the 1990-epoch Global Land Survey collection of Landsat images (Gutman et al. 2008) and other Landsat images selected from the USGS archive, each of which was atmospherically corrected to surface reflectance and converted to radiant temperature by the LEDAPS implementation of the 6S radiative transfer algorithm (Masek et al. 2006b). Whereas retrievals from within the period of overlap between the Landsat-5, Landsat-7, and MODIS eras may be based on general-even global-models based on phenological metrics that require dense image samples within each year (e.g., Hansen et al. 2013), this local fitting instead maximizes use of the single-image coverage characteristic of much of the history of Earth observation. Use of atmospherically corrected surface reflectance fulfills the conditions for signature extension in space .
Decision trees and other empirical classifiers are sensitive to bias in training samples relative to class proportions within their population of inference (Borak 1999, Carpenter et al. 1999, Sexton et al., 2013c and to uncertainty in the training data set (McIver 2002, Strahler 1980. To minimize these effects, we maintained a large sample with representative class proportions by removing a small, but equal fraction of the least stable pixels from each class while maintaining the class proportions from reference epoch to training sample. Further, we weighted each pixel's contribution to the classifier's parameterization based on the pixel's classification certainty in the reference data. A weight w was adopted for each pixel as the classification probability of the estimate (pmax) of forestor non-forest cover (C) from the 2000-epoch dataset: Global Land Cover Facility = ( ).

(29)
The weights were then applied to adjust the objective (i.e., purity) function maximized by the iterative binary recursion algorithm employed by C5.0™ (Quinlan 1986).

Forest-cover change
Classification trees estimate the probability p(C) of each class in each pixel as a conditional relative frequency. Given C = "F" (i.e., "forest"), each pixel was labeled either "forest" or "non-forest" based on p(F): Forest-cover change between 1990 and 2000 epochs was detected given the joint probabilities in 1990 and 2000 epochs : That is, given the probability of forest P(F) vs. non-forest P(N) in a pixel i in the 1990-epoch (t1) and 2000-epoch (t2), four classes were derived: stable forest (FF), stable non-forest (NN), forest gain (NF), and forest loss (FN). A categorical map of change classes was then produced by assigning each pixel the class with the highest probability.

Hedge rule
In the forest cover change products, the forest dynamics (i.e., forest loss and forest gain) between two periods were determined by checking the joint probabilities of forest and nonforest estimated for each of the dates (Kim et al., 2014;Sexton et al., 2015). Dynamic classes are more difficult to detect than stable classes, and a criterion is applied to filter the detected change estimates (Kim et al., 2014;Sexton et al., 2015). As an example, Figure 14 presents the accuracies for forest loss and gain between 2000 and 2005 calculated with criterions increased from 0.1 to 0.9 at 0.1 intervals. The global overall accuracy of the forest cover change data is insensitive to the changing criterions because areas of forest change are small fractions of global land area, but the criterion has a significant impact to the accuracy of the estimates of the forest dynamics. Higher criterion rules out less certain changes but also leads to high omission errors. On the contrary, lower criterion reduces omission errors, but introduces higher chance of commission errors.
According to the investigation presented in Figure 14, commission and omission errors reached the closest point when the criterion is near 0.6 for both forest loss and forest gain. Global Land Cover Facility Legend UA PA OA The threshold 0.6 was, therefore, applied to the production of the forest cover change datasets to produce unbalanced estimations of the global forest dynamics.

Minimum Mapping Unit
A minimum mapping unit (MMU) was applied to comply with the forest definition and also to minimize erroneous detection of change due to spatial misregistration of Landsat images. Raster polygons smaller than the threshold MMU (0.27 hectare, or 3 pixels) were replaced by the class of the largest neighboring polygon. An eight-neighbor rule was used to delineate patches, which includes diagonally connected neighbors.

Sampling design
Accuracy assessment employed a two-stage, stratified sampling design (Cochran, 1977;Sannier et al., 2014;Särndal et al., 1992;Stehman, 1999;Stehman & Czaplewski, 1998). To increase the representation of rare classes, reference data were sampled across the global land area in two stages, first selecting Landsat WRS-2 tiles within predefined global strata and then sampling pixels within each selected tile. The spatial location of sample points was held constant for all time periods.

Biome definition
Biome-level stratification was based on the 16 major habitat types delineated by the Nature Conservancy (TNC) Terrestrial Ecoregions of the World dataset (TNC, 2012). Excluding deserts and xeric shrublands, inland water, and rock & ice, we merged the major habitat types into eight forest and non-forest biomes (Table 8). Among the 7,277 WRS-2 tiles in the 8 biomes, the 5,294 tiles completely contained within any biome were assigned to their respective biomes, and tiles spanning biome boundaries (including land/ocean Global Land Cover Facility boundaries) were excluded. This reduced the land area for each of the 8 biomes available for sampling by 18.7 -58.2% of each biome (

.2 Tile selection
Sampling within biomes focused on WRS-2 tiles exhibiting high rates of vegetation change, detected using the Training Data Automation and Support Vector Machines (TDA-SVM) change-detection algorithm (Huang et al., 2008). The median vegetation-change rate for each biome was then used as the threshold for discriminating high-and low-change strata for that biome. Within each biome, eight tiles were then randomly selected in the Table 8. Reclassification of TNC major habitat types (TNC, 2012) into biome strata. The land area for each biome is reported in "Land area (km 2 )" column, and the percentage of that area reduced by excluding tiles spanning boundaries is reported in "Spanning biome WRS-2 tiles (%)" column. The percentage of the remained area after the "spanning biome" exclusion that further reduced by excluding edge pixels is reported in the "Edge pixels (%)" column. Global Land Cover Facility  (36) where is the desired number of sampled tiles within the population of the stratum ( ); was set to 4 and 8 for low-and high-change strata, respectively. A random number * was assigned to each tile, and tiles with * < ( | ) were selected as the sample tiles. Globally, 89 tiles were selected out of the intended 96 because only one tile met the criterion for the "high-change" stratum in the boreal non-forest biome.

Point selection
Following biome-level sampling, each selected tile was divided into 8 strata representing forest/non-forest status in each of the two periods, 1990-2000 and 2000-2005. This preliminary forest/non-forest discrimination was again performed by TDA-SVM. All pixels identified as cloud, shadow, water, or no-data, as well as pixels located at the edge of two classes, were excluded from the population. This exclusion reduced the available land area for each of the 8 biomes by 3.8 -13.2% (Table 8).
The inclusion probability for each stratum was calculated as: where the probability ( | ) is the ratio of the desired number of pixels ( ) to the total number of pixels in the stratum ( ). As recommended by Congalton (1991) and Olofsson et al. (2014), was set to 50 for each stratum (S). A random number * was assigned to each Global Land Cover Facility pixel, and pixels with * < ( | ) were selected as the sample points. A total of 27,988 points were thus collected across the globe. Figure 15 shows the selected points in WRS-2 tile p224r078, located at the boundary of Paraguay, Argentina, and Brazil.

Response design
Forest or non-forest cover in each pixel and each epoch was visually identified by experienced image analysts using a web-based tool presenting the GLS Landsat image(s) covering each location, as well as auxiliary information, including: Normalized Difference Vegetation Index (NDVI) phenology from MODIS, high-resolution satellite imagery and maps from Google Maps, and geotagged ground photos ( Figure 16) (Feng et al., 2012b). The Landsat images were presented in multiple 3-band combinations-e.g., near infrared (NIR)red (R)-green (G), R-G-blue (B), and shortwave infrared (SWIR)-NIR-R. The extent of each selected 30-m Landsat pixel was extracted in the Universal Transverse Mercator (UTM) coordinate system and delineated in both the Landsat image and in Google Maps to facilitate visual comparison. The NDVI profile was derived from the 8-day composited Global Land Cover Facility surface reflectance data (MOD09A1; Vermote & Kotchenova, 2008;Vermote et al., 2002) with nearest-neighbor interpolation, excluding data labeled as cloud or shadow in the MOD09A1 Quality Assurance (QA) layer (Feng et al., 2012b).
The selected points were randomly distributed among 12 experts for interpretation (Table 9). Experts visually checked the information provided by the tool and labeled each point either "forest" or "non-forest" for each of the 3 epochs individually. Points with Landsat pixels contaminated with cloud or shadow were labeled as "cloud" and "shadow" respectively. If an expert was unable to identify the cover of a pixel, he or she was instructed to label it as "unknown" for further investigation.  Over 1,000 collected points were located in each decile of tree cover, with nearly uniform sample size across the range of tree cover > 10% cover (Figure 17). Of these points, > 90% were labeled as forest or non-forest by visual interpretation of TM or ETM+ images in the 1990, 2000, and 2005 epochs, with only 6 % of the points remaining as "unknown". Less than 1 % of the points across all epochs were interpreted as "cloud" or "shadow". The distribution of the unknown points in the 2000 epoch revealed that these difficult points were rare (< 4 %) in areas of low or high tree-canopy cover but were much more frequent in areas with 5 -35 % tree cover ( Figure 18).

Validation metrics
Based on the independent reference sample, the labeled points were used to quantify the accuracy of the global forest-cover and -change layers using validation metrics weighted by area (Card, 1982;Congalton, 1991;Stehman & Czaplewski, 1998;Stehman, 2014). For each reference datum, i, the agreement between estimated and reference cover or change, y, was defined: (38) 0 if ≠ Weights were applied to the data to remove the effect of disproportional sampling, by standardizing the inclusion probability of each observation proportional to the area of each stratum (Sexton et al., 2013b). Each point's weight, , was calculated as the inverse of the joint standardized probability of its selection at the tile-and pixel-sampling stages: where ( | ) is the inclusion probability of the desired number of pixels ( ) to be randomly selected from the number of pixels in the Landsat scene ( ), and ( | ) is the probability of the desired number of Landsat tiles ( ) selected from the total number of Landsat scenes ( ) located inside the corresponding biome. Adjusting the weight by the cosine of the pixel's latitude (ϕ) corrects the sampling bias due to the increasing density of WRS-2 tiles with latitude.
Overall accuracy (OA) was calculated as the weighted number of points showing agreement between the estimated and the reference (i.e., human-interpreted) class-i.e., Percent of " unknown" points (%) Global Land Cover Facility ∑ a y elements of the diagonal of the confusion matrix-divided by the weighed total number of points ( ): The conditional probability of the estimate given the reference (i.e., human-interpreted) class, P(c|ĉ) (i.e., User's Accuracy, UA) and the conditional probability of the reference class given the estimate P(ĉ|c) (i.e., Producer's accuracy, PA) were calculated as: where n were the points identified as type c (e.g., forest, non-forest, forest gain, or forest loss) by the GLCF layers, and n c were the points identified as type c by the reference (Stehman, 2014). The inverse of P(c|ĉ) and P(ĉ|c) were interpreted as errors of commission and omission respectively.

Validation metrics
The variance of the accuracy metrics is described below. The points in each forest/nonforest status stratum were randomly selected. Hence, the variance of the OA for the stratum and the UA and PA of class c (i.e., forest and non-forest for forest cover; FF, FN, NF, and NN for forest-cover change) in the stratum were calculated following ( where was the number of points in the error matrix at cell (i, j), and + and + were respectively the summaries of row (i) and columns (j) in the matrix.
The estimated variances ( ( )) for the accuracy metrics (i.e., OA, UA, and PA) of the globe and each biome were calculated following (Cochran, 1977): where is the central latitude of tile (j). A tile (j) consisted of forest status strata, and the accuracy for the tile ( ) was estimated: where was the weight for a forest status stratum (i) within tile (j): ∑ where was the number of pixels in stratum (i) of tile (j). The mean ( ) of accuracy ( ) for tile (j) was calculated: The standard error (SE) of each accuracy metric was calculated as the square root of its variance:

Accuracies of forest-cover layers
Accuracy of forest-cover detection was consistently high across all biomes and epochs, with OA = 91% (SE≈1%) in each of the 1990, 2000, and 2005 layers ( Figure 11, Table  10). Commission errors (CE = 1 -P(c|ĉ)) and omission errors (OE = 1 -P(ĉ|c)) were < 10% for both forest and non-forest classes in all epochs, for which SE < 2.3%. The original, unadjusted estimates showed a bias toward detection of non-forest, with the forest class having a higher rate of omission errors (<21%) than commission errors (<3%) and the nonforest class having a higher rate of commission errors (<13%) than omission errors (<2%) in all epochs and biomes (Table 11).
These estimates of accuracy are likely conservative, given our exclusion of treeless biomes and the uncertainty of reference data generated by identifying forest cover by visual interpretation of satellite images (Montesano et al., 2009;Sexton et al., 2015a). Montesano et al. (2009) found that human experts achieved 18.7% RMSE in visual estimation of tree cover in high-resolution imagery, and Sexton et al. (2015a) found that visual confusion was greatest near the threshold of tree cover used to define forests, especially when interpreting change. To investigate the relation between accuracy and tree cover, OA of forest/nonforest cover in 2000 was plotted over the range of coincident tree cover estimated by the NASA GFC tree-cover dataset (Sexton et al., 2013a). A distinct concavity was evident in the relation, which reached its minimum near the 30% tree-cover threshold used to define forests ( Figure 19). The OA was large (> 80%) where tree cover was < 0.1 or > 0.35. Commission and omission errors were also investigated in relation with tree cover ( Figure  20). Commission error of the forest class was < 10% except among pixels with tree cover < 0.35, where the commission error was < 20%. Omission error of forest was < 20% in areas with > 0.4 tree cover but increased in areas of sparse tree cover.

Accuracies of forest-change layers
Globally, overall accuracy (OA) of the 1990-2000 forest-change layer equaled 88.1% (SE = 1.19%), and OA = 90.2% (SE = 1.1%) for the 2000-2005 forest-change layer (Table 12). In each period and biome, OA ≥ 78.7% (SE < 5%) ( Table 13). The global accuracies and standard errors of stable forest (FF) and stable non-forest (NN) classes were similar respectively to those of the stable forest and non-forest classes in the 1990, 2000, and 2005 Global Land Cover Facility layers, but the change classes-i.e., forest loss (FN) and forest gain (NF)-had larger error rates than the static classes in the respective epochs.
Commission and omission errors for forest loss were between 45% and 62% globally, with SE between 1.72% and 23.48%. Forest-loss was detected most accurately, with errors dominated by commission, in temperate and tropical evergreen forest biomes (PA ≥ 71.7%; UA ≥ 49.6%). This was likely due to relatively minimal impact of vegetation phenology on canopy reflectance in evergreen forests. Whether in temperate or tropical regions, detection of forest loss was more accurate in evergreen forests than in their deciduous counterparts (30% ≤ PA < 39%; 36.1% ≤ UA ≤ 50.1%). In non-forest biomes, accuracy of forest-loss detection was very low and dominated by omissions, but the rarity of forests and their loss in these biomes made the impact of these errors on overall accuracy small.
Forest gain was consistently the most difficult dynamic to detect, with OE and CE each > 60% in all epochs (SE < 17%). This was likely due to the long traversal of intermediate tree cover during canopy recovery from disturbance, compounded by the uncertainty of human identification of change (Sexton et al. 2015a). Producer's accuracies tended to be largest in tropical evergreen forests (24.9% ≤ PA ≤ 75.7%), where canopy recovery following disturbance is fastest, and smallest in non-forest biomes (PA < 19%; UA < 17%), where recovery is slower and locations spend more time in intermediate ranges of canopy cover.
The effect of tree cover on accuracy was investigated using the 2000-2005 forestchange layer (Figure 21). Similar to that of the 2000 forest-cover layer, a distinct concavity was evident in the relationship between overall forest-change accuracy and tree cover, and accuracy was lowest between 0.2-0.3 tree cover. Commission and omission errors of stable forest and non-forest in relation to tree cover were similar to those of forest and non-forest in the static layers ( Figure 22). The commission and omission errors were large in areas with tree cover < 0.35 and decreased to < 60% in areas with tree cover > 0.35. Commission and omission errors of forest gain were both correlated to tree cover. The omission error was < 45% and commission error was < 70% in areas with 0.3 -0.6 tree cover but > 50% in high or low tree cover. Global Land Cover Facility

(B) NN (A) FF
Accuracy Accuracy Overall Accuracy (OA) Accuracy Accuracy

Forest Edge
Forest edge was mapped as the Euclidean distance from each forest pixel to its nearest non-forest pixel at 90-m resolution. For each 90-m "forest" pixel , the Euclidean distance was calculated to the nearest "non-forest" pixel ′: where x and y are meters of longitude and latitude in Lambert Azimuth Equal Area projection, respectively. The calculation was not performed in non-forest pixels, and so the resulting map represents only forest-edge effects. Non-forest pixels were coded with null values. Because this process takes an (area-weighted) average of all forest pixels within the extent of each 1-km pixel, even 1-km pixels with only one 90-m forest pixel show a distance value. Pixels-especially those with small edgedistances-should not be interpreted as fully forested. Histograms were tallied from these data for each continent and summed globally before coarsening to 1-km resolution via bilinear interpolation.

Forest-patch area
Forest patches were defined by applying an 8-neighbor rule and a 1-ha minimum mapping unit to the binary forest/nonforest map at 90-m resolution. Each forest patch was then labeled with a unique value, i, and area calculated as the sum of all (forest) pixel's area within i: where dx = dy = 90 m.
To enable computation, the calculation was performed for each continent individually and the results merged; buffers were used to avoid truncation of patch size near continental borders.

Forest-patch isolation
Isolation of forest patches was calculated as the least edge-to-edge distance from each forest patch to the nearest forest patch: Ii was calculated based on the (x,y) location of pixel centers, and so the metric has a minimum value of 180 m. To enable computation, Eqn (54) was calculated for a random 20% of patches i, but against all patches i'. Quantitative analyses were performed using these data before coarsening to 5-km resolution.

Data-product access and computation
All of the NASA GFC datasets have been made available via the Global Land Cover Facility, via the GLCF Earth Science Data Interface (ESDI) and File Transmission Protocol (FTP). Developed with support from the NASA REaSON program, ESDI is a web-based tool for users to search and download data from GLCF's archive using spatial and non-spatial queries. FTP is used by those who are more familiar with the structure of the GLCF archive, those who want to automate data downloading using scripts, and for those who use GLCF as a read-only "cloud" storage solution. ESDI's mapping interface uses Java Server Pages (JSP) coupled with MapServer. JSP handles the user clicks for selecting data and selecting the type of query and passes the attributes to MapServer for displaying the data coverage on the map (Figure 23). This is helpful for users to know if their area of interest has data coverage.

Landsat Global Land Survey
Although the entire Landsat archive is now available through the USGS EROS Data Center, a large portion of users prefer optimized data collections such as the Global Land Survey (GLS). The GLS is the result of a partnership between the USGS and NASA in support of the U.S. Climate Change Science Program (CCSP) and the NASA Land-cover and Land-use Change (LCLUC) Program. Building on the existing GeoCover dataset (Tucker, et al, 2004) developed for the 1970's, 1990, and 2000, the GLS was selected to provide wall-to-wall, orthorectified, cloud-free Landsat coverage of Earth's land area at 30meter resolution in nominal "epochs" of 1975, 1990, 2000, 2005 and 2010. The GLCF currently houses and distributes the GLS Landsat dataset for 1975, 1990, 2000, 2005 and 2010 epochs. Depending on the epoch, approximately 7-10,000 Landsat scenes have been compiled to cover the global land area (Gutman et al. 2013;Feng et al. 2013). Over 534 terabytes of the GLS collection have been distributed during the NASA MEaSUREs project. This number does not include the additional Landsat data  that was downloaded to improve the characterization of tree cover.

Surface Reflectance
The GLCF built the LEDAPS modules on their cluster and started processing the GLS2000 collection in 2009. In processing the GLS collection we identified issues in the data and notified USGS. Upon correction of the data, we downloaded the data, and reprocessed it to SR. On our cluster we were able to process the GLS2000 data in about two weeks. As we added more nodes to the cluster we could ultimately process it in about 4 days. In June of 2011, we launched the first global Landsat based SR product upon the submission of the peer reviewed paper  to Remote Sensing of Environment, which was accepted in 2012.  The product was a success and significantly added to the total volume of data distributed via the GLCF. Roughly 15 terabytes of the GLS-2000 SR dataset were distributed during the first month that the product was available. The subsequent spikes in data distribution were due to either new versions or additional epochs of data being added to the online archive. The large spike in June 2015 likely was due to USGS LPDAAC copying the entire GLCF entire SR archive as part of the MEaSUREs data-archiving process. A total of roughly 192 terabytes of data have been distributed to date.

Data Formats and Values
The SR data product was distributed in GeoTIFF file format. Each bands with in a folder were individually compressed and separately made available via FTP.  Land mask (0= water, 1=land) Bit 6-15 Unused

Tree Cover
The tree cover data product, though initially an ancillary layer for generating the forest and forest cover change product, quickly became an important source of information for the user community. The team started generating the product late 2012 and stated distributing the data in early 2013. We have processed GLS 2000GLS , 2005GLS and 2010 to tree cover for this project. Since there is no MODIS tree cover product for the 1990s, we have not been able to generate tree cover before the year 2000. Further improvement of the tree cover product is ongoing with funding support from NASA Carbon Cycle Science and Land Use Land Cover Change programs.
Landsat Tree Cover Downloaded in GB

Data Formats and Values
The derived tree cover product was tiled using the WRS-2 two tiling scheme and kept the native projection information from the Landsat tile. Each tree cover data folder has 6 files associated with it; a browse file, preview file, data file, a per pixel uncertainty layer, an index file, and a text file. See the example below: p015r033_TC_2000: The tree cover data folder is named using the following convention: p stands for path, followed by three digits which represent the WRS-2 path, then r which stands for row followed by the three digits which represent the WRS2-row. Between the underscore are two letters (TC) which is the short name for the tree cover product, followed by four digits which represents the year for which the dataset was generated.

Forest Cover and Change
Results of the world's first global forest cover and change product (version 0) were presented at the NASA LCLUC meeting in 2012 and were subsequently published . We made the beta release of the forest cover and change in May 2013 to select user group to assess the data get feedback. Once we received the feedback we improved the product and released version 1 of the product in May of 2014. We are currently distributing the forest cover and change product from 1990 to 2000 and 2000 to 2005.

Data Formats and Values
The derived forest cover product was tiled using the WRS-2 tiling scheme and kept the native resolution information from the tree cover product that was used to generate the forest cover and change product. Each forest cover folder has 4 files associated with it; a browse file, a preview file, the change map file and the change probability file. See example below: p015r033_FCC_19902000: The forest cover and change data folder is named using the following convention: p stands for path, followed by three digits which represent the WRS-2 path, then r which stands for row followed by the three digits which represent the WRS2-row. Between the underscore