Optical remotely sensed time series data for land cover classiﬁcation: A review

Accurate land cover information is required for science, monitoring, and reporting. Land cover changes naturally over time, as well as a result of anthropogenic activities. Monitoring and mapping of land cover and land cover change in a consistent and robust manner over large areas is made possible with Earth Observation (EO) data. Land cover products satisfying a range of science and policy information needs are currently produced periodically at different spatial and temporal scales. The increased availability of EO data—particularly from the Landsat archive (and soon to be augmented with Sentinel-2 data)—cou-pled with improved computing and storage capacity with novel image compositing approaches, have resulted in the availability of annual, large-area, gap-free, surface reﬂectance data products. In turn, these data products support the development of annual land cover products that can be both informed and con- strained by change detection outputs. The inclusion of time series change in the land cover mapping process provides information on class stability and informs on logical class transitions (both temporally and categorically). In this review, we present the issues and opportunities associated with generating and val-idating time-series informed annual, large-area, land cover products, and identify methods suited to incorporating time series information and other novel inputs for land cover characterization.


Introduction
Land cover relates to the biophysical cover of the Earth's terrestrial surface, identifying vegetation, inland water, bare soil or human infrastructure. Distinct land cover types provide specific habitats and determine the exchange of energy and carbon between terrestrial and atmospheric regions (Houghton et al., 2012;Running et al., 1999). Biophysical processes occurring on the land surface have an impact on the climate system (Mahmood et al., 2014) and also on habitat diversity and the capacity of ecosystems to support human needs. Land cover, recognized as an essential climate variable (GCOS, 2003), is the most important element for description and study of the environment (Herold et al., 2006) and is a key input to models of ecosystem services (Andrew et al., 2014).
Characterizing and mapping land cover is essential for planning and managing natural resources (e.g. development, conservation), modeling environmental variables, and for understanding distribution of habitats. Remote sensing and digital image processing enable observation, identification, mapping, assessment, and monitoring of land cover at a range of spatial, temporal, and thematic scales (Rogan and Chen, 2004). Identification of land cover types provides basic information for generation of other thematic maps, and establishes a baseline for monitoring activities.
Regional and global land cover products have been made possible by the synoptic and periodic observations of satellite Earth Observations (EO) missions, starting in 1972 with Landsat-1 (Belward and Skøien, 2015). Several national land cover mapping programs became operational in the 1990s and the early 2000s (e.g. Vogelmann et al., 2001;Wulder et al., 2008a). National and international agencies have produced regional land cover maps representing single-date conditions (e.g. EOSD 2000 in Canada; CORINE Land Cover 2000 in Europe), with Landsat observations acquired within an interval of the target year (e.g. ±1-3 years). Some programs repeat the land cover mapping periodically (e.g. NLCD 2001NLCD /2006NLCD /2011 in the US; Australian NCAS-LCCP) and even annually for certain land cover types (e.g. Brazilian PRODES) to enable change assessment.
Global maps have also been produced, typically based on cloudfree composites of coarse spatial resolution data from AVHRR, SPOT VEGETATION, MODIS or MERIS data (e.g. IGBP-DISCover, GLC2000, MODIS Collection 4 and 5 Land Cover product, Glob-Cover). Global land cover maps generated from coarse spatial resolution data have typically had low local accuracy (Frey and Smith, 2007;Fritz et al., 2010), particularly in regions with heterogeneous land cover (Herold et al., 2008). Also, due to disparities in classification schemes, spatial resolution, thematic detail, or estimated accuracy, comparison of these maps to inform on change over time is problematic (Bai et al., 2014;Pérez-Hoyos et al., 2012;Pflugmacher et al., 2011) and generally discouraged (Gebhardt et al., 2014;Homer et al., 2007), as inaccuracies in thematic classes compound when change is considered (Fuller et al., 2003). Global land cover maps have also been produced recently with multiple-year data from Landsat (Gong et al., 2013) and singleyear Landsat-like imagery (Chen et al., 2014), with reported overall accuracies ranging from 65-80%.
Land cover maps can be updated by identification and mapping of changed areas, leaving unchanged areas in the original map intact. Some large-area land cover programs currently apply such a change-updating approach, for example European CORINE Land Cover (Büttner, 2014) or the US NLCD (Xian et al., 2009). Remote sensing is well suited to monitoring land cover change, and a myriad of approaches have been developed for the purpose (Coppin et al., 2004;Hussain et al., 2013;Lu et al., 2004;Singh, 1989;Tewkesbury et al., 2015). The change-updating approach for production of frequent land cover maps is conceptually appealing, operationally efficient, and economically convenient; however, a potential downside is the introduction of additional sources of error (i.e. via change detection), particularly in heterogeneous and dynamic regions, with frequent and diverse land cover changes. Alternative and robust methods to generate annual large-area land cover maps, such as from time series of EO imagery, require additional investigation and application development.
Until recently, land cover maps generated from coarse spatial resolution data (e.g. 1 km from AVHRR or MODIS) were deemed adequate for ingestion in global and large-area environmental models (Hansen and Loveland, 2012;Loveland et al., 2000). A confluence of data availability, technological capacity, and increasingly detailed information needs has motivated progression of research and operational efforts for generation of temporally, spatially, and thematically improved land cover databases.
Land cover changes can be related to natural processes (e.g. flooding, wildfire) and anthropogenic activities (e.g. urbanization, agriculture). The rate of change and the nature of land cover transitions can also differ in time and space. Some regions are relatively stable (e.g. permanent forest); whereas, other areas are subject to rapid and persistent transformation (e.g. urban expansion of previously vegetated areas). Increased human population and technological development has been found to accelerate land cover change (Goldewijk, 2001;Ramankutty and Foley, 1999). Assessing the location, extent, type, and frequency of land cover transitions, as well as identifying spatial and temporal patterns of change through interpretation and analysis of frequent land cover maps provide insights into underlying processes and drivers of change . Annual land cover information is valuable to aid in formulation of socio-economic policies (e.g. European Common Agriculture Policy) and data provision for environmental applications (e.g. vulnerability and risk assessment).
Maps that relate change between two dates typically lack information regarding underlying processes and do not enable insights on the nature of the transformations present, such as rate or persistence of change (Gillanders et al., 2008). In order to study changes in land cover, ideally both subtle modifications and rapid changes should be accounted for (Lambin et al., 2003). A temporal series of land cover maps can capture the complexities of Earth's changing surface (Liu and Cai, 2012;Sexton et al., 2013a) and can be used to parameterize biogeochemical models (e.g. Feddema et al., 2005;Running et al., 1994). Furthermore, the level of detail, the definition of land cover categories, and the land cover types selected for generation of a map are often insufficient or inadequate for a broad range of uses, and specific semantics might be misunderstood and wrongly used (Comber et al., 2005). An alternative option is to generate information of essential land cover attributes (e.g. canopy cover) as a versatile approach for generation of widely usable land cover products. Elementary feature layers can support-via rule-based combinations of attributes-multiple classification schemes to suit a range of science communities and policy makers. Demonstrating progress in this direction, estimates of bare ground, woody and herbaceous vegetation proportions are provided globally using MODIS (Hansen et al., 2003), and in the US using Landsat Vegetation Continuous Field products . Likewise, the European CORINE Land Cover project also aims to periodically generate five high resolution attribute layers: imperviousness, forest, grassland, wetland, and water bodies (Blanes Guàrdia et al., 2014).
Time series of medium spatial resolution optical data have demonstrated high capacity for characterization of environmental phenomena, describing trends as well as discrete change events. Time series of medium spatial resolution data have been used to map forest disturbance  and surface water bodies (Tulbure and Broich, 2013), to characterize land cover change (Zhu and Woodcock, 2014a) and to identify the nature of land cover changes (Olthof and Fraser, 2014). Time series enable modeling and estimation of ecosystem structural variables and have been used to estimate aboveground biomass Main-Korn et al., 2013), forest carbon sinks (Gómez et al., 2012), forest degradation (Shimabukuro et al., 2014), and forest disturbance . Strategies have been developed to deal with irregular and spare time series of data , but near-anniversary annual series are most appropriate for information extraction in vegetated ecosystems. Moreover, intra-annual time series have proven of great value to acquire phenological insights for land cover mapping (Melaas et al., 2013).
Herein, we investigate current opportunities and issues for generation and validation of temporally dense, large-area land cover products to meet current science, monitoring, and reporting information needs. Our focus is on data and land cover products with sufficient temporal frequency (i.e. annual) and spatial resolution to support monitoring efforts. Specifically, we focus on time series of medium spatial resolution remotely sensed data (e.g. 10-100 m), primarily Landsat (with full awareness of applicability to the forthcoming Sentinel 2 data stream), with time series being defined as a series of consecutive observations, acquired at regular time intervals. In addition, we focus on products that provide accurate estimates of land cover change (including anthropogenic impacts), enable land cover transitions to be characterized, and that are sufficiently flexible or adaptable to provide versatility for a range of applications. We place emphasis on the identification of data and approaches suited to incorporating time series information and novel inputs for land cover characterization.
2. Developments enabling progress: Data, technology, and analysis 2.1. Medium spatial resolution satellite observations Medium spatial resolution (10-100 m) EO data are captured by an increasing number of satellite missions (Belward and Skøien, 2015). Synergies between operational programs like Landsat and Sentinel-2  will raise the frequency of geometrically and radiometrically compatible acquisitions. Combined, the Landsat OLI and Sentinel-2 dual system are expected to provide global observations with 2-5 day frequency (Drusch et al., 2012;Irons et al., 2012;Wulder et al., 2015). Satellite sensors appropriate for land cover classification provide consistent and repeatable measurements at an appropriate spatial scale (Verbesselt et al., 2010). Landsat data (TM, ETM+, and OLI) are the standard for land cover classification (Cohen and Goward, 2004) and change detection (Wulder et al., 2008a), as a result of their spatial resolution (30 m), revisit period (16 days; 8 days with multiple Landsat sensors in orbit), and wide spatial coverage (185 Â 185 km). Rigorous calibration and consistency in the radiometry of the Landsat sensors, in particular the Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI), makes the Landsat image archive a strong example for the benefits of calibration and data interoperability. Given operational imperatives, other data sources (e.g. SPOT, IRS, ASTER) can be availed upon to offer complementary image coverage (Wulder et al., 2008b).

Landsat pixel-based image composites
Despite long-term plans for periodic systematic acquisitions (Markham et al., 2004;Williams et al., 2006) and the improved accessibility to Landsat data through global archive consolidation efforts Wulder et al., 2016) the amount of Landsat data available for many locations is less than desirable (Roy et al., 2010a). In persistently cloudy areas, the number of available observations may be inadequate for analysis (due to cloud, shadow and haze), both within and between years. To help address the limitations of scene-level data quality, pixel-based compositing with Landsat data increased in popularity soon after the opening of the USGS Landsat archive in 2008 . Image compositing is an established approach, common to global and regional studies using lower spatial resolution data (as reviewed in White et al., 2014;Hansen and Loveland, 2012). Prior to free and open access to the Landsat archive , compositing with higher resolution data was largely precluded due to data costs  and a lack of image data sources with both large area and time series coverage.
Repeat observations over the same location increase the chance of cloud-free data, even in areas of persistent cloud coverage. Large-area and cloud-free surface reflectance composites can be generated from preferred Landsat observations (Griffiths et al., 2013a;White et al., 2014), selecting the best available observation subject to user defined criteria (e.g. sensor preference, acquisition day of year, free of cloud, etcetera). These best-available-pixel (BAP) composites offer opportunities for monitoring with frequent data in areas where restrictions exist due to cloud persistence. BAP image composites, produced with rule-based criteria and preferences, can be optimized for specific applications , including land cover. Furthermore, the spatiotemporal context of the time series provides a means for improvement of residual missing or noisy pixels (i.e. spectral spikes) in image composites (Hermosilla et al., 2015a;Weiss et al., 2014), with these gap-free, infilled composites referred to as proxy BAP composites.
BAP image compositing approaches can take advantage of the robust radiometric calibration associated with the Landsat program, and automated approaches for generating surface reflectance values (e.g. LEDAPS, Masek et al., 2006). Automated processes can also be implemented for cloud and shadow detection and subsequent masking (e.g. Fmask; Zhu and Woodcock, 2012;Zhu et al., 2015a). For efficient, consistent, and timely large-area land cover mapping, scene-based Landsat inputs are not practical (Hansen et al., 2014), and the feasibility of image composites as foundation for generation of large-area products is a topic of active research (Franklin et al., 2015). To date, Landsat BAP image composites have been used for a range of applications (Table 1).

Towards annual land cover products derived from timeseries data
The availability of analysis-ready data products, increased technological capacity, and robust time series analysis approaches are guiding the emergence of methods for generation of annual thematic products informing current and historical land cover dynamics. The greatest challenge for novel methods is the complete integration of the temporal dimension into modeling. Temporal series approaches have been demonstrated as superior to singledate methods for a range of applications (Broich et al., 2011;Gómez et al., 2014), including land cover (Franklin et al., 2015).

Operational strategies for generating annual land cover maps
Methods for generating large-area land cover products must be robust, consistent, and repeatable (Franklin and Wulder, 2002). Although automation can minimize the rate of random error and reduce economic costs, operational experiences have uncovered the difficulties for automated mapping of large areas with acceptable accuracy (Bontemps et al., 2012;Gong et al., 2013). Fundamentally, two approaches can be followed for generation of annual land cover products supported by time-series data. First, independent land cover maps can be generated for each time step (i.e. for each year), using pre-and post-change information derived from the time series (Fig. 1A). This approach assumes that the models used to generate the land cover information are portable, and that land cover signature extension is enabled through the use of surface reflectance values (Song et al., 2001;Gray and Song, 2013). MODIS Collection 5 Global Land Cover product  is an example of an annual series of land cover maps produced independently using a time-series of intra-annual spectral data. A series of annual maps derived from inter-annual spectral time series, such as that proposed in Fig. 1A, was demonstrated by Franklin et al. (2015). In this study, the use of timeseries metrics in the classification approach improved the overall accuracy by 6.38% relative to single-date results. A second approach for generating annual land cover products from timeseries data would be to generate a base map representing reference date conditions for a single year, which could then be updated and backdated annually with change information obtained from the spectral time series (Fig. 1B). In this case annual spectral data provides information to identify the nature of change (e.g. subtle versus drastic), including the no-change option. As an example, Pouliot et al. (2014) mapped the forest area of Canada (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) annually with MODIS data and a change-updating approach. All efforts to generate time series informed land cover maps benefit from expert knowledge of ecologically feasible land cover transitions, as this contributes in building confidence, evaluating accuracy, and improving classification with an iterative strategy (Liu and Cai, 2012).
The complexity of the land cover characterization problem, together with increasingly detailed information needs for science, monitoring, and reporting, have advanced methods for mapping of single land cover classes (i.e. forest, water) (Frazier and Page, 2000;Schneider et al., 2010), in contrast to the more common multiple-class characterization required for simultaneous and spatially exhaustive mapping. In single-class mapping, binary classifiers segregate the target from all other cover types, thereby saving sampling resources . Building upon the single-class identification approach, cascade or hierarchical strategies progressively identify land cover types one by one (Chen et al., Multi-year (target ± 2 years) Six epochs: 1985/1990/1995/2000/2005/2010 Mapping forest disturbance and change in forest type Intra-annual and inter-annual spectral series offer different opportunities for a wide range of classification approaches. Intraannual series inform on the average phenology of individual land cover types, enabling identification of subtle differences and variations over multiple years. Inter-annual series provide medium-to long-term unique spectral profiles of land cover types and inform subtle modifications or drastic changes. Plus, inter-annual series provide evidence in support of possible changes to land cover. Several single date land cover maps have been generated with intraannual time-series data over large areas (Table 3). Also, some series of annual land cover maps have been produced over large areas (Gebhardt et al., 2014;Friedl et al., 2010;Pouliot et al., 2014).

Variables derived from spectrotemporal feature space
Spectral signatures of distinct land cover types vary with soil moisture, sun elevation, topography, atmospheric conditions, and view angle (sensor dependent). Spectral variability and signature overlap increase for land cover over large areas and in complex environments (e.g. mixed forests versus coniferous forests). Much of the spectral variability of land cover types can be captured with repeat observations over time, and thus the value of time-series data for discriminating between land cover types. Intra-annual variability is particularly helpful for identification of land cover types dominated by vegetation of marked phenology (e.g. deciduous forest, annual pasture). Phenological differences within a single-year or over consecutive years has helped in discriminating between cover types such as herbaceous crops and savannah (Senf et al., 2015) or crops and pasture (Müller et al., 2014). Inter-annual variations captured with anniversary measures inform processes of change or stability (Gómez et al., 2015;Lehmann et al., 2013); some processes are only feasible in certain cover types (e.g. stability of rock outcrop), while other processes indicate transitions from one land cover to another (Olthof and Fraser, 2014).
Multi-temporal spectral data encapsulates rich information ) that requires specialized tools. Numerous approaches and perspectives are possible for disentangling this multidimensional feature space (Table 2). Multi-temporal spectral data captured at various dates of the same or consecutive years can be summarized by statistical metrics (e.g. average, variability), which in turn can function as descriptive or predictive variables. Spectrotemporal statistical metrics have shown to be a viable means for discrimination of land cover types (DeFries et al., 1995;Gebhardt et al., 2014;Petitjean et al., 2012) and for quantification of forest cover change (Hansen et al., 2014).
A multi-year series of radiometrically calibrated anniversary spectral data is known as a temporal trajectory. Temporal trajectories provide insights of the spectral behavior of corresponding land cover elements over time. A number of methods have been developed to unveil and interpret the rich information provided by temporal trajectories. One strategy partitions the temporal trajectory into linear segments and interprets individual components (Land-Trendr; Kennedy et al., 2010). LandTrendr and similar algorithms that segment the temporal trajectory (e.g. C2C, Hermosilla et al., 2015a; DBEST, Jamali et al., 2015) derive descriptive measures of individual segments known as spectrotemporal change metrics (e.g. magnitude, duration). Spectrotemporal change metrics have been used successfully for characterization of land cover (Franklin et al., 2015), as well as forest disturbance and recovery (Hermosilla et al., 2015a;Kennedy et al., 2012).
Temporal trajectories can also be treated with a continuous approach as stationary (e.g. periodic) or non-stationary shape variables (i.e. referring to the shape of the temporal trajectory). Deviations from modeled stationary patterns, derived from multiple observations per year, have been used for classification of land cover types (Zhu and Woodcock, 2014b). Idealized non-stationary shapes have been used for identification of land cover transitions like stand-replacing forest disturbance (Kennedy et al., 2007) and for characterizing the nature of land cover change (Olthof and Fraser, 2014). Temporal trajectories capturing processes of natural succession (e.g. growth) constitute a spectrotemporal signature of specific land cover types. The length or duration of the shape variables necessary to identify significant patterns is presumably variable, as in the estimation of biomass Pflugmacher et al., 2014), and dependent on land cover type. Different strategies for analysis of time series are the decomposition into trend, seasonal, and break elements (BFAST, Verbesselt et al., 2010) enabling analyses of components in relation to landscape processes of change and land cover (DeVries et al., 2015;Verbesselt et al., 2012), or the simple analysis of trend (Eastman et al., 2009;Parmentier and Eastman, 2014). The existing range and variety of approaches and more than 30 years of radiometrically and spatially consistent Landsat TM/ETM+ observations will enable retrospective land cover mapping based on the spectrotemporal feature space.

Classification methods incorporating time series data
The choice of classification algorithm requires considering multiple aspects of the problem: type of data, statistical distribution of classes, target accuracy, ease of use, speed, scalability, and interpretability of the classifier. Direct tradeoffs of some of these factors exist, and it is also important to establish a balance between acceptable accuracy and optimal use of resources. Algorithms that cluster elements by similarity of attributes without a priori human intervention (i.e. unsupervised classification) are typically used when scarce knowledge of the land cover types is available (Chen and Gong, 2013;Eva et al., 2004). Clustering algorithms (e.g. kmeans, ISODATA) run iteratively until convergence of an optimal set of clusters. Because clusters produced automatically do not necessarily correspond with the land cover types (Loveland et al., 2000), post-classification refinement techniques (e.g. merging and splitting clusters) are necessary before labelling (e.g. Wulder et al., 2008a). Also, to avoid classes with high internal variance (e.g. water, bare soil, snow) dominating the clustering (Loveland et al., 1991), prior stratification and masking are common practices (Furby et al., 2008;Wulder et al., 2008a). Despite the attractiveness of the automatic character of clustering algorithms, they become time consuming when the data dimension is high or the data volume large (Chen and Gong, 2013), and interpreting clusters properly is a challenging and intensive process. Alternatively, supervised classification approaches assimilate the data to a number of reference land cover samples labelled a priori. Selecting an adequate number of good quality training samples is crucial (Bruzzone and Demir, 2014;Foody and Mathur, 2006;Shao and Lunetta, 2012), a time consuming task typically done manually, although accumulated experience and improved databases enable semiautomatic selection in certain conditions (e.g. Radoux et al., 2014). Selecting and labeling samples is not error-free and potentially the cause of poor and biased classification performance Pal and Mather, 2006). Supervised methods require that the training data completely represent the classification problem, the classifier is incapable of identifying what is unknown to the training sample (Foody, 2000(Foody, , 2004. The state of practice for large-area land cover mapping has evolved during the last decade, from a preponderance of unsupervised techniques (Franklin and Wulder, 2002) to an increased use of supervised techniques (Khatami et al., 2016), in part due to increasing ancillary data that facilitates sample collection for training data (e.g. Colditz et al., 2011;Homer et al., 2007;Radoux et al., 2014). Other approaches involve various classifiers used in parallel or in succession (e.g. Bauer et al., 1994;Bontemps et al., 2011;Chen et al., 2014;Tateishi et al., 2011) that can be both supervised and unsupervised. Table 3 illustrates how several of the common classification algorithms have been applied for large-area land cover mapping using time series of optical EO data. Table 4 synthesizes the strengths and weaknesses of algorithms included in Table 3, as identified in the literature.
Partitioning is the only clustering category (Han and Kamber, 2001) widely used for land cover mapping with remotely sensed data, although sporadic efforts of hierarchical clustering exist (Shenthilnath et al., 2011). For large datasets, k-means and ISO-DATA algorithms are preferred, as they are less time consuming than others. Parametric supervised classifiers (e.g. maximum likelihood, minimum distance, discriminant analysis) are difficult to use with multi-temporal data of many spectral features and multimodal distributions (Glanz et al., 2014). As a general rule, their reduced flexibility in decision boundaries makes parametric classifiers unsatisfactory for characterization of land cover in large areas and complex environments (Hubert-Moy et al., 2001). On the other hand, non-parametric classifiers (e.g. k-Nearest Neighbor, kNN; decision trees, DT; neural networks, NN; and Support Vector Machines, SVM) impose boundaries of arbitrary geometries, and provide higher flexibility at the expense of computationally intense iterative processes. In general, non-parametric classifiers that focus decision rules on class boundaries are proficient when the statistics and distribution of land cover types are unknown (Foody and Mathur, 2006;Hansen, 2012)-a common scenario for large areas-thereby making non-parametric classifiers more appropriate than parametric classifiers, which focus on central tendency statistics.
Significant effort has been dedicated to evaluating the performance of land cover classification algorithms, and identifying their relative strengths and shortcomings (Huang et al., 2002;Mather, 2003, 2005;Rodríguez Galiano et al., 2012). A few efforts have specifically compared algorithms using time-series data (Jia   Schneider, 2012). The strengths of an algorithm may be general (e.g. easiness of application and interpretation) or they may manifest only in certain scenarios (e.g. capacity to handle missing data) (Table 4). DT, which are based on recursive binary partitions complying with a set of optimized rules (Breiman et al., 1984), are an attractive option for large-area land cover classification for a number of reasons, primarily their ease of application and interpretation, and their capacity to handle data measured on different scales, non-linear relationships, and missing data . DT can be trained quickly and perform rapidly (Pal and Mather, 2003); however, DT have been shown to perform less well than algorithms such as SVM and NN in feature spaces with high dimensionality (Hansen, 2012;Pal and Mather, 2003), and are also sensitive to noisy observations and over-fitting (Ghimire et al., 2012). Random Forest (RF) (Breiman, 1996) is an improved implementation of a DT, which votes for the best tree out of a forest of trees created by recursive modification of the sampling data (Belgiu and Drǎguț, 2016). RF provides higher classification accuracy than other forms of DT and avoids over-fitting, at the cost of increased computational intensity and a black box nature that obfuscates decision rules (Rodríguez . SVM algorithms identify one or more hyperplanes separating target groups in a multi-dimensional space. SVM is less sensitive to the Hughes phenomenon (Hughes, 1968)-whereby classification accuracy decreases when the number of input features to the classifier is over a given threshold-than other algorithms, and performs well with a large number of variables relative to training size (Pal and Mather, 2005;Shao and Lunetta, 2012), a significant benefit when abundant remotely sensed data exists and little ground truth is available. Neural networks are accurate classifiers, but tend to overfit the data (Rodríguez  and remain a black box for interpretation. Both NN and SVM are computationally intense and require a number of parameters to be adjusted (Table 4).
Combining the strengths of various algorithms into an ensemble classifier tends to increase accuracy (Bauer and Kohavi, 1999) and provides information on classification uncertainty  or confidence (McIver and Friedl, 2001), but it also tends to reduce interpretability and increase computational complexity and overhead (Table 4). Ensemble learning algorithms use the same or a combination of base classifiers to produce multiple classifications of the same data (e.g. Random Forest, bagging, boosting) (Rodríguez . Ensemble approaches can be dependent, whereby the output of one classifier is used to inform the next classifier, or independent, whereby each classifier is run independently and the outputs are combined using some weighting or voting mechanism (Rokach, 2010). Bagging methods resample the original training set with replacement each time, while boosting methods resample or reweigh the training data by emphasizing more those instances that are misclassified by previous classifiers. Boosting has been found to be a convenient technique for large classification problems .
The incorporation of complex temporal data into land cover classification has not been paired with novel classification algorithms and in most cases the same rules are applied as for classifying single-date data. A notable exception by Liu and Cai (2012) demonstrated the utility of considering spectrotemporal context in the classifier. It has been noted that the inclusion of timeseries data can have a more positive impact on results (i.e. accuracy) than the choice of classification algorithm (Jia et al., 2014b). Novel algorithms are needed to leverage the predictive potential of time-series data; consider the case of hyperspectral and multispectral datasets, whereby techniques that function well with single-date data, may not offer the same performance when using time-series data (Bruzzone and Demir, 2014).

Novel inputs for land cover classification
Progress in the development of robust methods for the generation of repeat land cover products over large areas is made possible with contributions from novel inputs that were not necessary or possible before. Knowledge of ecological succession can facilitate refinement of a time series of land cover maps, while information on land cover dynamics can be used to impose model restrictions on short-and medium-term land cover transitions, contributing to the temporal consistency of land cover products. Variables with an inherent temporal dimension and temporal trajectories of multiple vegetation indices provide complementary insights regarding the status and change of essential biophysical attributes associated with land cover types. Improved quality of training samples, and combinations of multi-scale and multi-sensor data for enhancement of temporal and spatial resolution has direct impacts on the accuracy of land cover characterization.

Land cover transitions and temporal stabilization of classification
The likely character of land cover transitions depends on the time interval considered and on the environmental context. Likely successional transitions (e.g. open to dense forest) may take just a few years to manifest in tropical areas but longer in boreal or Mediterranean regions. Possible land cover transitions by loss of vegetation (e.g. shrub to bare ground) could happen rapidly anywhere following a sudden disturbance. Furthermore, transitions related to subtle, long-term processes occur in different directions (e.g. canopy closure by natural growth, loss of density by diseaserelated mortality, change in dominant species). Understanding ecological processes is essential for production of reliable time series land cover maps, regardless of the mapping strategy applied. Insights of possible and likely transitions at different time intervals can be used to impose logical restrictions for generation of consecutive maps (Liu and Cai, 2012;Pouliot and Latifovic, 2013) and also for refinement and verification. Land cover transitions can be identified by modifications in spectral seasonal curves over several years (Gutiérrez-Vélez and DeFries, 2013;Reed et al., 1994) and by specific patterns or trends in inter-annual spectral trajectories that can inform on change types (Olthof and Fraser, 2014).
Developing accurate and consistent annual land cover maps requires dedicated techniques for temporal stabilization of the classification results . Furthermore, for some change assessment applications, temporal consistency is more important than the overall accuracy of individual maps (Radoux et al., 2014). Markov Random Field (MRF) models are particularly attractive techniques in temporal pattern recognition, as a result of their ability to identify relationships of temporally dependent data, and have been applied in land cover map series to identify and correct illogical class transitions (i.e. transitions that are ecologically impossible) (Wehmann and Liu, 2015), and to incorporate spatial context into classifications (Li et al., 2014). Liu and Cai (2012)   a spatial-temporal consistency model, improving the classification of single date maps produced with multi-year MODIS data. Other approaches incorporating transitions in the process of classification have also demonstrated an improvement to the temporal consistency of land cover maps series. For example, working with annual series of land cover class probabilities derived from MODIS via Random Forest classifiers, Yin et al. (2014) identified land cover transitions by the change in probability of occurrence for the five land cover types of their classification system. Postclassification approaches that apply temporal filters with transition rules (Clark et al., 2010) or change transition matrices (Pouliot et al., 2014) have also been shown to be useful for identification and amendment of illogical transitions. Utilizing the legend relating the land cover present over the forested area of Canada (after Wulder et al., 2008a) a class transition matrix is generated and presented (Fig. 2). The above logic is incorporated to create scores in the transition matrix from likely, through probable and possible, to not likely. Such a matrix can be used to reconcile the consistency across annual land cover maps and to incorporate and inform with knowledge of change (such as from spectral trends or change detection). Class stability in a non-changing pixel-series can be enforced through rules. More sophisticated approaches to use of the class transition matrix may incorporate disturbance and rules informed by knowledge of successional trends and spectral evidence.

Multi-scale and multi-sensor combination of spectral data
Bringing together medium spatial and high temporal data from complementary sensors like Landsat and MODIS or MERIS has been shown to be successful for land cover applications (Amorós-López et al., 2013;Gong et al., 2013). Combining these data sources is increasingly frequent, either via concurrent use or by fusing data sets to generate synthetic images. Supported by the spatial detail of Landsat, the high temporal data from MODIS/MERIS informs phenological traits and helps identify difficult land cover types (Jia et al., 2014a(Jia et al., , 2014c. In their Landsat-based global map, Gong et al. (2013) supported training data selection with MODIS Enhanced Vegetation Index (EVI) time series from 2010, aiming to improve spectral separation between cultivated bare land and natural barren lands, and to distinguish sparsely vegetated land from more persistent barren land.
Among data fusion algorithms, the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM)  with improvements (Hilker et al., 2009;Zhu et al., 2010) stands out for the algorithm's ability to generate Landsat-like images, although unmixing-based algorithms and the Spatial and Temporal Reflectance Unmixing Model (STRUM) are preferable when few Landsat data are available (Gevaert and García-Haro, 2015). In a land cover change updating approach context, Chen et al. (2015) employed downscaled intra-annual MODIS NDVI time series for Fig. 2. Class transition matrix for the land cover types present over the forested area of Canada (after Wulder et al., 2008a). For a single year period four kinds of transition are considered according to their likelihood to occur: from ''likely" (stability or no change) and ''probable" (natural succession transition), through ''possible" (transition that could happen) to ''not likely" (ecologically illogical transition). ''Not likely" transitions point to potential error in the prior or/and later land cover classification. identification of real rather than spurious change, improving the overall land cover classification accuracy from 0.76 to 0.89. Jia et al. (2014a) fused Landsat and MODIS data via STARFM into 23 16-day composites and used statistical metrics (maximum, minimum, mean, and standard deviation) for supervised classification improving single date classification accuracies by $4%. In this case no land cover change is assumed during the period of blended data.
Despite successful efforts blending MODIS and Landsat data, an important restriction for this approach is the limited availability of MODIS data, going back just to year 2000. Also, Senf et al. (2015) demonstrated this type of synthetic seasonal data improves single date land cover classification, but the authors show that the results from the synthetic data are inferior to those from the original Landsat data. Synthetic Landsat images free of clouds, shadows, snow and SLC-off gaps, and without phenological and sun orientation issues can also be generated from Landsat models in regions where abundant observations are already available (Zhu et al., 2015b). The capacity of the synthetic Landsat data developed by Zhu et al. (2015b) for land cover mapping has yet to be demonstrated. Another approach for improving land cover characterization is to combine data from different temporal scales acquired by a single sensor. For instance, both an annual composite and a time series of bimonthly MERIS composites were used to improve classification results of the GlobCover 2005 and 2009 products (Bontemps et al., 2012). The advantages of using data from a single sensor are in the consistency of the spatial registration and radiometric calibration.

Multiple temporal spectral variables
For accurate characterization of land cover types and optimal detection of land cover transitions multiple spectral variables are needed (Zhu and Woodcock, 2014a). The temporal trajectory of vegetation indices like the Normalized Burnt Ratio (NBR) is particularly useful for identification of land cover transitions associated with drastic loss of vegetation cover (e.g. forest to bare soil), because values are markedly changed after years of relative stability (Hermosilla et al., 2015b). Other spectral variables provide information of change causality (Hais et al., 2009;Schroeder et al., 2011) informing possible land cover transitions. Furthermore, subtle processes responsible of transitions between similar cover types (e.g. open forest to dense forest by recruitment or growth) can only be identified by temporal spectral trajectories robust to noise and sensitive to variations in vegetation density (e.g. Tasseled Cap Angle).
Smoothing algorithms have also been applied to reduce noise inherent in the time series data (e.g. Shao et al., 2016). For generation of rich and detailed land cover dynamics databases over large areas, multiple spectral variables are required. Parmentier and Eastman (2014) combined time series of various vegetation indices and successfully identified land cover transitions, although the contribution of each variable remained unsolved.
Including multiple spectrotemporal variables in the classification problem increases the feature space, hampering some classifiers' performance. Dimensionality reduction by linear or nonlinear techniques (Yan and Roy, 2015) facilitates data compression while retaining relevant information. Principal Component Analysis (PCA) is a well understood dimensionality reduction technique, but its components can be difficult to interpret (Homer et al., 2004). Vegetation indices are an efficient way to reduce the amount of data and some of them provide the additional advantage of being linked to vegetation's physical characteristics (e.g. Tasseled Cap Transformation).

Optimized training sample data
Reference data should represent all land cover types within the area of interest, and should be defined in the context of the classi-fication approach to be applied. Classifiers based on statistical distributions (e.g. Maximum likelihood) perform better with pure samples (Lillesand and Kiefer, 2008), whereas classifiers with rules targeting boundary regions between classes (e.g. SVM, DT) benefit from critical border samples (Foody, 2004;Hansen, 2012). Nonparametric methods such as kNN require the reference data to represent the full range of variability in the population, as these methods cannot extrapolate beyond the range of the training data. The adequate size of training sample is also algorithm dependent (Foody and Mathur, 2006;Pal and Mather, 2003) and although no guidelines exist, a general assumption is that the size of the training sample should be proportional to the independent space size (i.e. the number of variables in the model). Classification algorithms frequently suffer from inability to handle high-dimensional data. Therefore, supervised classifications aim to target an adequate sample with a minimum number of elements that maximize accuracy. To improve the accuracy/sample size relation, dimensionality reduction via feature transformation (e.g. PCA; Tasseled Cap Transformation, TCT) and feature selection (Elghazel and Aussem, 2015) are useful tools.
Acquiring reference data is expensive and time consuming, but also critical for the accuracy of classification results (Shao and Lunetta, 2012). The use of time-series data for land cover classification presents additional challenges for reference data sampling. Signature extension by generalization of a single set of class spectral signatures across time and space might yield substantial interclass mixing. Distributing samples in space and time with the objective of capturing spectral variability (Sexton et al., 2013a(Sexton et al., , 2013b) relies on good radiometric normalization between image dates. Alternatively, automatic and adaptive generalization of spectral signatures in temporally irregular (non-anniversary) time series (Gray and Song, 2013) is possible. Annual series of BAP image composites represent a seamless radiometrically corrected dataset that enables internally consistent and repeatable land cover characterization (Hansen et al., 2013;Pengra et al., 2015). Furthermore, datasets with a temporal dimension (e.g. interannual and intra-annual spectral series) provide opportunities to augment the sample size, assimilating patterns of identical spectrotemporal trajectories across space and time to a good reference (Xue et al., 2014) with reliable measures of time series similarity (e.g. Dynamic Time Warping or variations; Maus et al., in press). Active learning (semi-supervised) methods may also optimize the choice of training samples, proposing training samples to the operator in an iterative form until satisfactory classification is achieved (Tuia et al., 2011). Active learning can be beneficial at the initial phases of the mapping project, reducing the expensive phase of redundant sample labelling (Bruzzone and Demir, 2014), but can be difficult to implement over large areas.

Validation and accuracy assessment
Validation of land cover products is essential to demonstrate the quality of the products for operational applications (Cihlar, 2000;Foody, 2002). The accuracy of classified data has to be assessed and reported with informative metrics adequate for the user community. Although there are no standard methods to assess the accuracy of thematic maps, best practices have emerged in the literature, indicating the inclusion of an error or confusion matrix (Foody, 2002), which can aid in identifying confusion between classes, as well as potential sources of error. Furthermore, quantitative metrics that are derived from the confusion matrix, such as overall accuracy, as well as area-weighted metrics and confidence intervals are useful measures (Olofsson et al., 2014). Validation of large-area land cover products is challenging and requires appropriate sampling strategies for statistical assessment of the product's accuracy. Moreover, the reference data must have sufficient accuracy to enable robust validation (Wulder et al., 2006). The reliability of large-area land cover maps is frequently tested against reference maps produced at a higher spatial resolution assuming good, although not perfect accuracy.
Validation of two or more consecutive land cover maps with reference datasets is a difficult task, because acquiring reference data for multiple years over large areas can be unfeasible, even if auxiliary maps exist at certain dates. Adequate reference data is particularly elusive if land cover products span a historical period. In such cases, manual interpretation of the Landsat time-series data is becoming an accepted approach for generating the required reference data . In the change-based updating approach (Fig. 1B), if the area changed is proportionally small, the accuracy of land cover maps in a time series is often assumed to be close to that of the base map (Pouliot and Latifovic, 2013). In theory, validation of a time series of land cover and change detection results can be achieved with independent error matrices (Mertens and Lambin, 2000;Yuan et al., 2005). However, to ensure robustness, this approach requires no-change samples to be acquired as a component of the reference data (Olofsson et al., 2014). Post-classification comparison of a time series of land cover maps can aid in identifying illogical land cover transitions in space and time (Liu and Cai, 2012), which, in well-registered maps of the same spatial resolution, can be indicative of classification error (Townsend et al., 2009). Moreover, explicit information of change incorporated in land cover products provides a powerful tool for self-assessment and validation. The temporal consistency between consecutive maps can, to some extent, be evaluated against land cover changes. Acceptable land cover transitions in consecutive maps should conform to the time interval separating them, that is, short intervals impose tighter restrictions and provide more reliable judgement. For example, in an annual series of land cover maps, grassland can transit to pasture but not to open forest, but the same transitions would be more uncertain in a five year period. Permitted transitions between land cover classes provide a means for self-assessment, but only partially, because even ecologically logical class transitions could be incorrect.
Two main difficulties arise with validation via logical land cover transitions. First, transitions associated with subtle or progressive spectral variations are difficult to assess, because enough evidence or historical data are seldom available. However, a short delay in detection of such transitions should be considered a minor error compared with, for example, missing a forest to bare ground disturbance-related transition. Second, validation of transitions becomes even more challenging with a higher number or very similar classes (Olofsson et al., 2013). Likewise, confusion between discrete classes of a continuum (e.g. between herb and shrub) does not have the same implications as other confusion types (e.g. between herb and water). Adapted or new measures of error are also needed to report the accuracy of time series land cover maps (e.g. Tsutsumida and Comber, 2015).

Synthesis: Issues and opportunities
Despite progress in technology and data availability, generation of large-area land cover products is challenging. While some of these challenges relate to the migration of algorithms from a research to an operational phase (e.g. handling and processing massive amounts of data), we have identified some challenges specific to the use of time-series data, such as annual image BAP composites, for land cover mapping. We have also identified opportunities afforded by novel data analysis techniques developed by the remote sensing community and other scientific disci-plines. These techniques merit further exploration for the purpose of land cover characterization with time-series data.

Time-series data and temporal variables
Generating Landsat annual BAP image composites (Roy et al., 2010b;White et al., 2014) and improved annual series of proxy image composites (Hermosilla et al., 2015a) over large areas has been demonstrated. Until now, single-date land cover maps have been produced with multiple intra-annual composites (Yan and Roy, 2014). When generating a series of annual land cover maps, consideration of phenological information from seasonal image composites (Roy et al., 2010b)-in areas where there is sufficient data to do so-would presumably facilitate better discrimination of land cover types. However, the irregular availability of medium spatial resolution data over large areas and across many years, as well as in high latitudes with persistent snow coverage, may jeopardize the robustness and spatial-temporal consistency of the land cover product. MODIS data can support Landsat data via combination or fusion, and inform on phenological traits for dates after 1999 only (Chen et al., 2015;Jia et al., 2014a).
Among the common approaches applied thus far for characterizing the spectrotemporal domain of Landsat time series (e.g. change metrics, trends), some are better suited for certain applications. Statistical metrics have demonstrated efficiencies for synthesizing temporal data to distinguish among various land cover types, but their specificity is not demonstrated: various land cover types could be represented by the same statistical metrics over a certain period if the metrics are not properly defined. Change metrics describe the evolution over time of biophysical parameters correlated with the spectral features observed. Change metrics are intuitive and easy to characterize and interpret. The choice of an optimal set of change metrics, including single or multiple spectral features, as well as the interval prior and after the target date, are crucial and application specific. Shape variables are visually appealing and easy to interpret, but characterizing and modeling the domain of shape variables is not an easy task. The relevance of slight pattern variations, which might be difficult to capture with parametric or functional approaches, is also application dependent; for example, a shortly delayed maximum in a density indicator does not modify biomass estimation  but is a determinant of land cover change types (Olthof and Fraser, 2014). Parameterized trends are not as intuitive as other variables, but have been effective in some land cover problems (Eastman et al., 2009;Parmentier and Eastman, 2014). The capacity of these approaches to provide temporal variables that are resilient to inter-annual reflectance variability and capable of robustly and consistently discriminating land cover types in annual series of large-area land cover maps needs to be demonstrated. The spatio-temporal context provided by an annual series of proxy images and insights of ecologically feasible land cover transitions facilitate development of meaningful rules with which to define spatially and temporally consistent maps. Hermosilla et al. (2015b) present an example for using Random Forests to label disturbance types combining change metrics, object-level spectral information, and the geometric characteristics of the objects. Disturbance year is assigned using a breakpoint analysis of the Landsat time series, with 90% accuracy found for the disturbance types labeled over the 30-year analysis period.

Emerging options for automatic classification of land cover
Despite the widespread incorporation of temporal data into approaches for characterizing the status and dynamics of ecosystems, and recognition of the superiority of models that integrate temporal data, development of novel land cover classification methods that take full advantage of this information remains to be realized. Temporal data are operationally analyzed disregarding the full capacity of the time dimension and presumably missing the opportunity to unveil more detailed information. We mention here some techniques that are novel or are being used in an exploratory phase, and that could benefit the land cover mapping community, some of which incorporate the temporal dimension.
Unsupervised clustering is currently only marginally used in large-area land cover mapping (e.g. GlobCover 2009; Wulder et al., 2008a), despite the attractiveness of its automatic character and independence from field samples. In the context of time series informed land cover, clustering deserves renewed attention because the strengths of novel techniques could be decisive for the land cover problem. Traditional clustering focuses on discovering a single grouping of the data, which does not necessarily correspond with the land cover types, particularly when dealing with multidimensional feature domains. Some novel approaches define multiple sets of clusters that describe alternative aspects and characterize the data in different ways. Three novel clustering paradigms are attracting considerable research effort (Müller et al., 2015): alternative clustering, clustering ensemble, and subspace clustering. Alternative clustering and clustering ensemble search for an alternative to a given solution and a consensus of various clustering options respectively, and neither of them seem practical for the large-area land cover problem. In contrast, subspace clustering focuses on detecting multiple clusters in arbitrary, sometimes overlapping, subspace projections of the high dimensional data (Parsons et al., 2004). For land cover characterization with time-series features, subspace clustering is an attractive tool for exploration and identification of relevant subsets of variables that characterize essential land cover attributes individually. Subspace clustering also provides a means for interpretation of complex data spaces.
Clustering of time series differs from clustering of static data mainly in how the similarity between elements is computed (Liao, 2005). Dynamic Time Warping, Mahalanobis, or Triangle distances are measures generally preferred over Euclidean distance, but the adequacy of a time series similarity measure is application dependent (Gorecki, 2014;Lhermitte et al., 2008). Time-series spectral data can be directly clustered as raw data, or indirectly clustered via features (e.g. statistics, change metrics) or modeled series (e.g. shape, trends) (Liao, 2005). The second approach (i.e. features) facilitates reduction of noise and data size, as well as classification of sequences of unequal length (e.g. missing data). Techniques for clustering continuous series are used in medical (Aach and Church, 2001;Wismüller et al., 2002), pattern recognition (Liao, 2005;Maharaj, 1999), and finance applications, accommodating single and sometimes multiple variables (Ramoni et al., 2000), and with hidden Markov models frequently underlying time series models for clustering (Zhang et al., 2011).
In supervised classification, machine learning and ensemble algorithms are more accurate and efficient than conventional parametric classifiers when faced with large dimensional and complex data spaces, and are currently used for large-area land cover mapping. However, supervised classification methods can benefit from semi-automated techniques for speeding up training data collection (Schneider, 2012). Semi-supervised learning (SSL) (also known as deep learning and active learning) methods for classification exploit both training data and unlabeled samples in the learning phase of the classifier, in order to obtain a general decision function that considers the training set information and the structure of all data in the feature space (Bruzzone and Demir, 2014). SSL can therefore outperform other classifiers when few training samples are available. SSL deals with the intrinsic variability of spectral signatures in land cover classes and has already been used for land cover characterization with hyperspectral data (Camps-Valls et al., 2007;Jun and Ghosh, 2013).
Another relatively new classification technique with potential application in large-area land cover mapping with time-series data is the Relevance Vector Machine (RVM), a Bayesian extension of the SVM (Tipping, 2001). As a kernel method like SVM, RVM classifies large dimensional data with few training samples, achieving high accuracy. Unlike SVM, RVM provides an estimate of the posterior probability of class membership, which may be used to illustrate the uncertainty of the class allocations (Foody, 2008). RVM provides more flexibility for the choice of kernel type than SVM, and its implementation does not require the penalty and kernel parameters, saving considerable computational resources (Gómez-Chova et al., 2011). RVM is a technique frequently used in pattern recognition (Caesarendra et al., 2010;Wang et al., 2009), and to date, just a few studies can be found in the land cover literature, using optical (Foody, 2008) and hyperspectral (Camps-Valls and Bruzzone, 2005;Demir and Ertürk, 2007) data.

Summary
Characterization of large-area land cover dynamics is necessary to inform monitoring, reporting, science, and policy. The increased availability of EO data at a range of spatial and temporal resolutions, coupled with increasingly complex information needs for science and policy making, and technological improvements in image processing capacity and storage, guide the advance of methods to generate frequent and accurate land cover products. Timeseries spectral data provides better characterization of land cover types, with intra-annual series informing on phenological differences and inter-annual series informing on land cover dynamics.
Techniques for data analysis and interpretation that fully incorporate the temporal dimension remain an area of intense research and represent an important challenge for operational mapping. For land cover characterization, intra-annual series of coarse spatial data have long informed on phenological traits and differences among land cover types. Medium spatial resolution intra-annual series with less frequency (i.e. Landsat) are feasible in certain areas, but not a global reality, due to variation in data availability. Interannual time series, on the other hand, are facilitated by BAP image composites, which compile the best possible observation at each pixel location, subject to user criteria (e.g. day of acquisition, sensor preference). Annual series of gap-free proxy image composites represent a valuable source of data for generation of time series land cover maps as well as other applications.
To date, information from Landsat time-series data has taken the form of statistical metrics, change metrics, or pattern/trend components. These metrics or components are subsequently used in applications, such as land cover mapping. Despite the temporal character of the data employed, classifiers used in operational mapping are those typically used for single date spectral data classification; however, the special traits of time-series data are quickly driving evolution of novel techniques, such as RVM. Improvement of existing approaches, as well as inclusion of novel techniques, imported and adapted from other disciplines, is important to make the most of data available and to produce annual land cover maps that meet a wide range of needs. Subspace clustering, time series clustering, and semi-supervised learning are techniques with potential to facilitate annual large-area land cover maps production.
Validation of annual series of land cover maps, and particularly historical maps, requires methods not fully dependent on field or ancillary reference data. Ecological insights of land cover transitions over different time intervals provide opportunities for assessment of land cover products. Incorporating spatiotemporal context into classification and validation improves the accuracy and temporal consistency of land cover time series. In this context, the quality of annual series of proxy image composites, without data gaps or anomalous values, represents a real opportunity for generation of accurate and consistent annual land cover maps.
Through this review we have demonstrated the need and related utility of annual land cover products to support increasingly complex information needs for monitoring, reporting, science, and policy. Time-series data offers new opportunities and challenges for land cover mapping. There is currently a dearth of land cover classification methods that can take full advantage of the rich information available from time-series data. This review presented a range of options for source image data (via image archives, image compositing, synthetic data generation) as well as a number of classification routines that avail upon both the spatial and temporal information present. The production of land cover maps was among the first applications of satellite multispectral remotely sensed data. Over four decades later, new and improved classification algorithms continue to emerge. Furthermore, sole reliance on spectral information for classification purposes is not required with supplemental datasets available to inform on conditions that can further aid in differentiation between classes. Incorporation of the temporal element into classification approaches is already benefiting from sophisticated spatial and spectral algorithms, pointing to the potential for further improvements in land cover classification outcomes. It is unlikely that there will be a single, or optimal, method for land cover mapping incorporating time-series data. Rather, approaches that incorporate spatial, spectral, and temporal data, and the knowledge of the underlying processes (e.g. phenology, disturbance, and succession) of the classes of interest are poised to meet a wide range of information needs.