Essential Climate Variables were listed by the Global Climate Observing System as critical information to further understand the climate system and support climate modelling. The European Space Agency launched its Climate Change Initiative (CCI) to provide an adequate response to the set of requirements for long-term satellite-based products for climate.

Within this program, the Land Cover (LC) project aimed at revisiting all algorithms required for the generation of global LC products that are stable and consistent over time. To this end, the LC concept was revisited to deliver a set of consistent 300m global annual LC maps from the 1990s to 2015 describing the land surface into 22 classes using the FAO LC Classification System.

The entire MERIS Full and Reduced Resolution (FR and RR) archive from 2003 to 2012 was first classified into a unique 10-year baseline LC map, which is then back- and up-dated using (i) AVHRR time series from 1992 to 1999, (ii) SPOT-VGT time series from 1998 to 2012 and (iii) PROBA-V time series from 2013 to 2015. This method avoids independent classification from year to year, thus ensuring temporal and spatial consistency between successive maps. The CCI LC project also delivered climatological 7-day time series representing seasonal burned areas, vegetation and snow dynamics of the land surface, 7-day surface reflectance time series for the whole archive of MERIS FR and RR and a global map of open permanent water bodies at 150 m derived from Envisat ASAR images and auxiliary data. All products are delivered along with an aggregation tool, enabling re-projection and re-sampling as well as the translation from LC classes into Plant Functional Types for the different climate models.

The extension of the global annual LC maps series from 2016 to 2019 in a consistent manner will be ensured with PROBA-V within the operational European Commission Copernicus Climate Change Service. Each LC map of a given year will be released the next year, with accuracy assessed in the 3 months following the release.

The study advances our understanding of climate-human population-tree cover dynamics in sub-Saharan Africa. The tree cover in sub-Saharan Africa is a sink for approximately 25% of global carbon and can act as a buffer against drought at various spatial scales, but remains understudied compared to other regions. Major socio-ecological transitions impacting tree cover have taken place in sub-Saharan Africa over decades. Large declines in tree cover have been detected in humid areas near population centers where demand for agricultural land is high. These areas have widened since the early 1990’s due primarily to rapid population growth. Population growth at drier latitudes is more gradual and has less influence on changes in tree cover. Instead, increased tree cover at drier latitudes has been linked to a wetter climate in the early 1990’s. To monitor tree cover over sub-Saharan Africa, we need long time series before the 1990’s, but we don’t have them. We expand on a method first demonstrated in highly fragmented landscapes of Kenya (Marshall et al., 2017: https://www.earth-syst-dynam.net/8/55/2017/esd-8-55-2017.html) to map climate-human population-tree cover on an annual basis at 5 km spatial resolution for historical (pre-1990) and contemporary (post-1990) periods in sub-Saharan Africa. Five-kilometer spatial resolution is appropriate to establish boundary conditions for regional climate models and related hydrological and geochemical cycle studies. The method uses Random Forest for data reduction and Generalized Additive Models that functionally relate big geospatial data inputs to changes in the proportion of tree cover. Major model inputs include bioclimatic variables (http://www.worldclim.org/bioclim), population density and phenological metrics that were derived from CHIRPS (http://chg.geog.ucsb.edu/data/chirps/), AfriPop (http://www.worldpop.org.uk/) and NASA GIMMS NDVI3g (https://nex.nasa.gov/nex/projects/1349/) from 1982-2015. AfriPop is considered the most accurate population geospatial dataset in sub-Saharan Africa, but is only available for 2010. UNEP GRID and GPW-v4 growth/decay coefficients were used to backcast and forward cast AfriPop to 1982 and 2015, respectively. Twenty additional model inputs were collected and processed that capture topography, soils and other socio-ecological conditions of change. Random Forest and Generalized Additive Models will be trained with 30 m spatial resolution Global Land Cover Facility’s tree cover estimates derived from Landsat 5/7 for 2000, 2005, 2010 and 2015 (http://glcf.umd.edu/data/landsatTreecover/). These data have been processed and brought to a common 5 km with the model inputs. It represents the major milestone of the study. Once the method is trained and cross-validated, independent sources, such as Global Land Cover Reference Data (https://landcover.usgs.gov/glc/) will be used for validation. We expect: 1) population density will be the most important and highly correlated variable with changes in the proportion of tree cover; 2) climate, phenology and other model inputs will be important, but will have significantly less explanatory power; and 3) large declines in tree cover will be estimated near population centers experiencing rapid growth, while modest tree cover growth will be estimated away from population centers experiencing prolonged periods of above-average rainfall.

Land cover has strong implications for the small-scale distribution of tundra snow cover. In particular, shrub height information is required in order to identify snow trapping areas and subsequent impact on ground thermal conditions. Current global and circumpolar maps lack thematic detail and/or spatial resolution to appropriately represent shrubs, their types and height.

The recently launched Sentinel-2 satellite optical data is suitable for land cover monitoring with 10 m spatial resolution. SAR data from Sentinel-1 (10-30m) have also been shown sensitive to vegetation height in tundra.

We have selected a transect in Western Siberia, spanning from the northern tip of the Yamal peninsula (continuous permafrost) to the south into the tundra-taiga transition zone (with discontinuous permafrost) for calibration of vegetation height retrieval. The region represents not only a gradient in vegetation zones but also sites with high heterogeneity of shrubs. Especially central Yamal is characterized by shrubs of up to 1.5 m height within continuous permafrost. Several studies over a CALM (Circum Arctic Active Layer Monitoring) site and additional transects have exemplified the role of the shrubs for active layer thickness in relation to snow using in situ measurements as well as Synthetic Aperture Radar data from satellites in this region. It is one of the key regions of ESA DUE GlobPermafrost. More than 1500 data points are available from surveys in 2017 and 2018.

Data from several sites across the Arctic have been used for validation. This includes a further site on central Yamal, the tundra uplands east of the Mackenzie Delta, the Alaska North Slope and the Yukon delta.
All available cloud-free images of summer 2016, 2017 and 2018 have been combined and several indices derived from Sentinel-2, including Tasseled Cap Wetness and NDVI. Sentinel-1 winter acquisitions (IW VV) have been pre-processed and normalized for comparison. Results show that all tested indices and data types respectively can be used to derive vegetation height up to 160 cm. Specific settings (e.g. aquatic vegetation, moraine dominated landscapes) do however require a combined use of Sentinel-1 and 2.

Spatially explicit accuracy estimates of land cover products are of considerable value for user communities in various fields, from biodiversity to agricultural monitoring. Well-known accuracy values such as kappa, or overall accuracy give a general idea on the thematic accuracy of the map, but don’t provide information on the spatial differences in map quality. Spatial accuracy estimates however, inform users regarding the degree of uncertainty in land cover mapping across space, thus building confidence in using different land cover products. From a user’s perspective, these spatially explicit accuracy estimates help when comparing different land cover products in order to choose the best product for a region of interest. From a developer’s perspective, such an accuracy product can be used to identify hotspots of uncertainty for further map improvements, as well as for monitoring these improvements over time. To date there have been a few publications with different methods on how to map per-pixel accuracies, e.g. kriging or the use of kernel functions. However, it remains unclear regarding how to interpret the results from the different methods, which input parameters should be reported along with the resulting spatial accuracy maps, and is the outcome a true pixel-based accuracy estimate? Finally, how can spatial accuracy maps be compared? On the latter question, users should be aware of the fact that all the existing methods are very sensitive to the input parameters such as the reference data sets, the number of nearest neighbors used for interpolation, spatial relation models, e.g. linear, exponential, and the computational resources. Therefore, one should only compare the spatial accuracy maps of different land cover products when these maps have been produced in the same way.

The aim of this study is to present the results from different methods for mapping spatial accuracy such as kriging and geographically weighted regression, and illustrate the sensitivity of these methods to different input parameters. This comparison has not been done previously. We will provide guidelines on how to produce spatial accuracy maps for different purposes, mainly for highlighting uncertain areas and for error propagation by modelers.

In addition, we will provide the first ever independent spatial accuracy estimates for a few recent land cover products. To undertake this validation, we have used the Copernicus Global Land Service (CGLS) reference data set on land cover. This data set has been developed by a group of experts, trained by staff at IIASA and WUR, using the Geo-Wiki web-application (https://geo-wiki.org/). At each reference location, experts were presented with a 100m x 100m area, subdivided into 10m x 10m cells (with the same geometry as Sentinel-2 imagery) superimposed on high-resolution imagery from Google Maps or Microsoft Bing Maps. Using visual cues acquired from multiple training sessions, the experts interpreted each cell as a land cover class such as trees, shrubs, water, arable land, burnt areas, etc. In total, we have collected circa 150 000 sample sites globally. The reference data set was designed to be used for validating land cover products at resolutions ranging from 10m to 100m, and adjusting to different land cover definitions. The year of reference for the data set is 2015 but there will be updates with changes over the period 2015-2018.

Based on the spatial and temporal characteristics of the reference data set, we have selected the following land cover products for the production of spatial accuracy maps:
• Global land cover maps:
o Copernicus Land Cover Map (https://land.copernicus.eu/global/products/lc) at a 100m resolution for 2015;
o Finer Resolution Observation and Monitoring of Global Land Cover (FROM GLC), at a 30m resolution for 2015;
• Regional land cover maps:
o ESA CCI LC regional maps, e.g. Africa and Mexico;
o EEA Copernicus layers at a 10m resolution;
• Thematic maps:
o Global 30m cropland extent (the Global Food Security-support Analysis Data - GFSAD) for 2015;
o Global forest map by Hansen et al (http://data.globalforestwatch.org/datasets/tree-cover-loss-hansenumdgoogleusgsnasa);
o Urban layers such as GUF+ (https://www.dlr.de/) and GHSL (https://ghsl.jrc.ec.europa.eu/).

The spatial accuracy maps for the listed above products show high variability of maps quality and highlight hotspot areas for improvements. There are areas where none of the maps is performing well and areas where one map performs better than others. We will highlight such areas in our talk.

1 Introduction
══════════════

Since 2016, the French Theia Land Data Centre publishes every year a
land cover map of metropolitan France at 10 m resolution. This map is
produced by CESBIO using Sentinel-2 image time series, and the same
approach has been applied to Landsat data to produce maps over the
same area for the pre-Sentinel years. The nomenclature contains the
following 17 classes: Continuous urban fabric; Discontinuous urban
fabric; Industrial or commercial units; Road surfaces; Annual summer
crops; Annual winter crops; Intensive grassland; Orchards; Vineyards;
Broad-leaved forest; Coniferous forest; Natural grasslands; Woody
moorlands; Beaches, dunes and sand plains; Bare rock; Glaciers and
perpetual snow; Water bodies.

Three years of operational production allows us to propose an
interesting retrospective integrating feedback from both the producer
and the user's point of view. In this talk, we will present the
lessons learned and give an overview of the ongoing and future work
for improving the quality and usefulness of the product.

2 Methodology
══════════════

The approach used is based on supervised classification of Sentinel-2
image time series and has been described in detail in [1]. The
reference data used for training the machine learning models comes
from a variety of existing data source as Corine Land Cover, Urban
Atlas or the French LPIS. Since these databases may be out of date and
have coarser spatial resolutions than the target map, a specific
pre-processing workflow has been designed.

The main difficulty for large area mapping comes from the data volumes
and their heterogeneous quality across the different regions. Indeed,
several orbits are needed to cover the whole area, which implies that
different areas have different acquisition dates. Furthermore, cloud
cover is heterogeneous. All this amounts to having time series with
different number of time stamps and temporal grids. The adopted
approach to use all available data (no image is discarded unless it
has more than 90% cloud cover) and resample them on to a common time
grid. Cloud and cloud shadow masks are used to decide whether a pixel
at a particular date is valid or not. In this way, a cleansed and
homogeneous data cube is produced and fed to the machine learning
step.

In order to take into account climatic variability between regions, a
spatial stratification is performed using eco-climatic areas. A
different classifier is trained for each climatic area allowing to
improve the overall accuracy by an amount close to 5%.

The production workflow is fully automatic and is implemented in the
`iota2' open source processing chain [2].

3 Product validation
═════════════════

The validation of the product is crucial for end-users to integrate
the maps in their workflows. Several levels of validation are
performed. The first one is a classical machine learning validation
which uses a hold-out data set extracted from the reference data and
that is not used for training. This allows to produce detailed
confusion matrices and derived metrics as overall accuracy, kappa
index and per-class metrics as the FScore. This validation allows to
check that the machine learning approach has good performances and
that over-fitting is avoided. However, due to the errors present in
the reference data, disagreements between the map and the reference
data may be difficult to interpret. Furthermore, biases in the
reference data could be present in the produced map and would not be
detected [3].

In order to complement this machine learning validation, independent
smaller data sets are also used. These datasets come from field
surveys, partner institutes and end users. They are small and local by
construction and can have different nomenclatures. Therefore, a
specific validation procedure including nomenclature translation and
class regrouping needs to be used. The results are coherent with the
classical machine learning validation and confirm the good overall
accuracy of the maps and also the main confusions between similar land
cover classes [4].

A third kind of validation, performed by professional
photo-interpreters was also performed for the 2016 map. A fully
independent validation following a sound protocol is costly and needs
skills and expertise that are very specific. SIRS is a company which
is specialised in the production of geographic data from satellite or
aerial images. Among other things, they are the producers of Corine
Land Cover for France and they are also responsible for quality
control and validation of other Copernicus Land products. The analysis
showed that the product meets the usually accepted thematic validation
requirements for operational Copernicus products [5].

4 Main limitations and user feedback
═════════════════════════════

Despite the good conclusions of the different validation approaches,
there are still many issues with the delivered products.

First of all, some individual classes are poorly detected and there
are strong confusions between some classes of the same category. For
instance, continuous and discontinuous urban fabric are sometimes
difficult to split. The same kind of problem occurs with grasslands
and moorlands, mostly in the gradient areas between the two. The
origin of these confusions is twofold:
1. the classes are very similar;
2. the reference data imposes arbitrary boundaries depending on the
area and the context.

The first issue needs better algorithms, the second needs better
reference data.

Second of all, many users ask for a more detailed nomenclature, mainly
in croplands and wetlands. As a consequence, a new nomenclature made
up of 30 classes instead of 17 has been designed and is starting to be
tested. The first results are encouraging, but new algorithmic
developments will be needed.

Finally, producing a new map once a year invites users to analyse land
cover changes. Depending on the scale of the change analysis, our maps
may not be appropriate. While the trends over large areas seem to be
reliable, the land cover evolution analysis at the pixel level is not
possible with the current maps. Indeed, comparing 2 maps which have
90% accuracy does not allow to look for changes which have occurrences
much lower than 10%. In order to help users interpret the changes, a
confidence map is also distributed. This map gives, for every pixel, a
confidence index of the decision of the classifier is provided. This
index can therefore be used to weight the changes detected between 2
maps. Future versions of the map will be accompanied by the
probability of each class, in order to allow more detailed analysis of
the observed changes. Nevertheless, this additional information is not
a solution, and better accuracies should be achieved to grant access
to robust change information.

5 What's next
══════════════

In order to improve the quality of the maps and their usefulness for
end users, several research activities are ongoing. The first one,
already pointed out above, is the extension of the nomenclature to 30
classes to fulfil users' demands. This nomenclature will be presented
in the talk. Added to that, the following algorithmic improvements are
being addressed:

1. In order to reduce confusions between similar classes, a contextual
classification approach needs to be implemented. Indeed, the high
spatial resolution of Sentinel-2 reveals textures and local pixel
distributions which are meaningful for the classification. Deep
learning approaches and particularly Convolutional Neural Networks
seem appropriate for the task. However, their computational cost
may be to high for operational mapping at scale. Other aspects as
the need for dense annotations and nearly error-free reference data
hinders their use in our application. In this context we work on
the development of contextual classification with shallow learning.
The first results show that these approaches are competitive with
Deep learning in terms of accuracy and much cheaper in terms of
computation.
2. Since some of the errors in the maps come from the fact that the
reference data used may be out of date, techniques allowing to
train the classifiers with fewer amounts of reference data for some
classes are interesting. We are developing domain adaptation
techniques which allow to use reference and image data from
previous periods to leverage supervised classification using
outdated reference data.
3. In order to reduce the false alarms in change detection, the use of
imagery spanning several years in the past is of great help, since
it allows to regularise the classification in the temporal
dimension.
4. In order to improve the maps in areas with frequent cloud cover,
our processing chain is able to jointly use Sentinel-1 and
Sentinel-2 time series. However, our experiments up to date do not
show significant increase in accuracy. This may be mainly due to
the increase of dimensionality of the problem. In this regard, we
need to investigate better approaches for the use of multi-modal
imagery.

[1] Inglada, J., Vincent, A., Arias, M., Tardy, B., Morin, D., &
Rodes, I., Operational high resolution land cover map production at
the country scale using satellite image time series, Remote Sensing,
9(1), 95 (2017). [http://dx.doi.org/10.3390/rs9010095]

[2] [https://framagit.org/inglada/iota2]

[3] [http://www.cesbio.ups-tlse.fr/multitemp/?p=12630]

[4] [http://www.cesbio.ups-tlse.fr/multitemp/?p=11835]

[5]
[http://www.cesbio.ups-tlse.fr/multitemp/wp-content/uploads/2018/02/SIRS_CESBIO_OSO_Validation_Report_v1.3.pdf]

Long-term global remote sensing data archives are now available to both the scientific and user communities from many past and current Earth Observation missions [1]. As a consequence, an unprecedented wide range of possibilities opens up for climate change analysis. In the framework of the Climate Change Initiative (CCI), the European Space Agency (ESA) has started the Global Monitoring of Essential Climate Variables (ECVs) to provide an appropriate, comprehensive and timely response to the strict requirements for long-term satellite based products for climate change analysis addressed to Space Agencies through the Global Climate Observing System (GCOS) [2],[3].
In this paper, we present the activities planned for the High Resolution (HR) Land Cover (HRLC) ECV project, launched in the context of the Climate Change Initiative Extension (CCI+) programme. Accurate land cover mapping and land cover change detection are extremely important for a better understanding of the global climate system [4]. Land cover is an essential climate variable that allows the quantification of surface energy, water fluxes and sources of greenhouse gasses, the monitoring of variation in land use as well as the characterization of the impact of extreme events (e.g., floods, heatwaves, droughts or hurricanes) [5]. The main purpose of the HR Land Cover ECV project, which is at its initial stage, is to study in detail the role of the spatial resolution in the mapping of land cover and land cover changes to support climate modelling research. Land cover and related changes are in fact both cause and consequence of climate change either when the change is human-induced or generated from natural events. This has been widely recognized by the scientific community and demonstrated by the previous phase of the CCI programme, focused on the generation of Medium Resolution (MR) Land Cover maps at global scale. Differently from the MR land cover CCI, which provided annual land cover maps at 300m resolution covering the period 1992-2015 [6], the HRLC project aims at producing maps characterized by a spatial resolution of 10m/30m. Moving from 300m to 30m requires the definition of new data analysis and processing methods, reframing the perspective with respect to the MR project both from the theoretical and the operational view point. Although HR potentially increases the capability of analyzing in detail spatial patterns and changes in the land cover, many challenges are introduced with respect to the MR case.

The temporal availability of HR data in the past/current archives is much lower than that of the MR ones and it strongly varies across the years. Differently from the MR case (e.g., SPOT-Vegetation archive), no daily acquisitions are available and only in the very recent years it was possible to get a dense temporal sampling thanks to the Sentinel and Landsat-8 missions. Prior to them, the amount of yearly-based images available in archives dramatically reduces (being Landsat Thematic Mapper the most relevant source), yielding to a much more challenging framework for the development of HRLC products. This scenario leads to a complex process to produce historical time series of products. Moreover, it requires a shift in the processing paradigm that moves from the analysis of many images per year acquired at MR to a few (for some areas and years very few) images characterized by high spatial resolution.
The primary output of the project consists of: (i) a HR land cover map at subcontinental scale at 10m as reference static input to the climate models, (ii) a long-term record of regional HR land cover maps at 30m in the regions identified for the historical analysis every 5 years, and (iii) change information at 30m and yearly scale for supporting the updates of the HR land cover maps. The 10m resolution subcontinental static maps will be produced mainly with Sentinel 1 and Sentinel 2 images. To generate the long term record of HR land cover maps the following processing chains will be defined: (i) an optical processing chain (focusing on Landsat data sets, but capable to integrate other possible multispectral images like SPOT), (ii) a SAR processing chain (focusing on ERS, ASAR), and (iii) a data fusion processing chain which combines the multi-sensor land cover maps. The static HR land cover will be used as reference map that will be back-updated every 5 years by integrating the change information at 30 m derived with a yearly scale.

The Climate User Group involved in this project defined three large regions of particular interest (see CCIplus_areas.png) to study the climate/LC feedbacks in three continents, involving climate (tropical, semi-arid, boreal) and complex surface atmosphere interactions that have significant impact not only on the regional climate but also on large-scale climate structures. These interactions and the associated vegetation-atmosphere feedbacks are still poorly captured by state-of-the art Earth System Models (ESM), not to mention their future evolution which is highly uncertain. Additionally, these regions are critical for the global carbon cycle (CC) through the uptake of carbon by terrestrial ecosystems; they cover major biomes that are vulnerable from the point of view of land carbon stocks (tropical forest, permafrost...).

The first region is located in the Amazon basin. Although Amazonian forests have considerable resilience to moderate annual drought, the interactions between deforestation, fire and drought potentially lead to losses of carbon storage and changes in regional precipitation patterns and river discharge [7] with some signs of a transition to a disturbance-dominated regime.
The second region that draws our attention is the Sahel band in Africa which is a complex climatic region experiencing severe climatic events (droughts in the 70's and 80's and floods more recently) often attributed to climate warming and for which the future predictions (amplitude of the regional warming and rainfall changes) are very uncertain [8].
The third region is situated in the northern high latitudes, for which future climate changes are expected to be particularly strong, a phenomenon known as polar amplification [9]. In Siberia, complex climate feedbacks over land, implicating natural and human factors, may further amplify these changes and make this region as a possible hot spot of future climate changes [10]. Here, changes in LC may impact directly the fate of the carbon stored in permafrost, which in turn will affect long-term terrestrial carbon balance and ultimately climate change [11].

[1] Bruzzone, Lorenzo, Begüm Demir, and Francesca Bovolo. "Recent Advances in Remote Sensing Time Series Image Classification." Pattern Recognition and Big Data. 2017. 713-734.
[2] ESA – European Space Agency: ESA Climate Change Initiative description, EOP-SEP/TN/0030-09/SP, Technical Note – 30 September 2009, 15 pp., 2009.
[3] Hollmann, Rainer, et al. "The ESA climate change initiative: Satellite data records for essential climate variables." Bulletin of the American Meteorological Society 94.10 (2013): 1541-1552.
[4] Bontemps, Sophie, et al. "Revisiting land cover observation to address the needs of the climate modeling community." Biogeosciences 9.6 (2012).
[5] Loveland, Thomas R. "History of land-cover mapping." Remote Sensing of Land Use and Land Cover. CRC Press, 2016. 32-41.
[6] Bontemps, S., et al. "Consistent global land cover maps for climate modelling communities: current achievements of the ESA’s land cover CCI." Proceedings of the ESA Living Planet Symposium. 2013.
[7] Davidson, E. A., de Araújo, A. C., Artaxo, P., Balch, J. K., Brown, I. F., Bustamante, M. D. C. & Munger, J. W. (2012). The Amazon basin in transition. Nature, 481(7381), 321-328.
[8] Roehrig, R., D. Bouniol, F. Guichard, F. Hourdin, and J.-L. Redelsperger (2013), The present and future of the west African monsoon: a process-oriented assessment of cmip5 simulations along the AMMA transect, Clim. Dyn., 26, 6471-6505.
[9] Holland, M. M., & Bitz, C. M. (2003). Polar amplification of climate change in coupled models. Climate Dynamics, 21(3-4), 221-232.
[10] Piao, S., Nan, H., Huntingford, C., Ciais, P., Friedlingstein, P., Sitch, S., ... & Levis, S. (2014). Evidence for a weakening relationship between interannual temperature variability and northern vegetation activity.
[11] Schuur, E. A. G., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., ... & Natali, S. M. (2015). Climate change and the permafrost carbon feedback. Nature, 520(7546), 171-179.