Monitoring Urban Areas with Sentinel-2A Data: Application to the Update of the Copernicus High Resolution Layer Imperviousness Degree

Abstract

Monitoring with high resolution land cover and especially of urban areas is a key task that is more and more required in a number of applications (urban planning, health monitoring, ecology, etc.). At the moment, some operational products, such as the “Copernicus High Resolution Imperviousness Layer”, are available to assess this information, but the frequency of updates is still limited despite the fact that more and more very high resolution data are acquired. In particular, the recent launch of the Sentinel-2A satellite in June 2015 makes available data with a minimum spatial resolution of 10 m, 13 spectral bands, wide acquisition coverage and short time revisits, which opens a large scale of new applications. In this work, we propose to exploit the benefit of Sentinel-2 images to monitor urban areas and to update Copernicus Land services, in particular the High Resolution Layer imperviousness. The approach relies on independent image classification (using already available Landsat images and new Sentinel-2 images) that are fused using the Dempster–Shafer theory. Experiments are performed on two urban areas: a large European city, Prague, in the Czech Republic, and a mid-sized one, Rennes, in France. Results, validated with a Kappa index over 0.9, illustrate the great interest of Sentinel-2 in operational projects, such as Copernicus products, and since such an approach can be conducted on very large areas, such as the European or global scale. Though classification and data fusion are not new, our process is original in the way it optimally combines uncertainties issued from classifications to generate more confident and precise imperviousness maps. The choice of imperviousness comes from the fact that it is a typical application where research meets the needs of an operational production. Moreover, the methodology presented in this paper can be used in any other land cover classification task using regular acquisitions issued, for example, from Sentinel-2.

1. Introduction

1.1. Monitoring Urban Areas with Remote Sensing: A Copernicus Land Mission

Understanding the urban growth phenomenon is among the major issues that public services have to deal with. Today, more than half of the world population lives in urban areas, and it is estimated that this will reach up to two-thirds by 2025 [1]. In the context of globalization and climate change, not only the urbanization process is fast, but the consequences due to the increase of imperviousness poses serious challenges related to the indication of risks (floods, heat wave events, pollution, etc.).

In Europe, the European Commission conducts projects to monitor urban areas in order to serve regional funding and policy making. Among those, the High Resolution Layer Imperviousness Degree (HRL IMD), as well as the Urban Atlas, which are part of the Copernicus Land missions, are based on remote sensing data. One of the objectives of the Copernicus program developed by the European Commission is to provide information products for monitoring the growth of urban areas in order to assess whether policies to reduce urban sprawl are effective.

Among those, the HRL IMD, as well as the Urban Atlas (http://land.copernicus.eu/) are key products based on remote sensing data. The strength of these Copernicus Land products relies on the use of high resolution (spatial/temporal) remote sensing data, which provide homogeneity, repetitiveness and objectivity over the whole of Europe, whereas compiling available local and/or national datasets throughout Europe in order to compare cities would be difficult to achieve since the data would originate from different materials and methods.

The HRL IMD consists of a $20\times 20$ m grid covering all of Europe, including Turkey, and aims at representing the degree of imperviousness within each grid cell from 0% to 100%. This layer had already been produced in 2006, 2009 and 2012 and is key for monitoring urban sprawl. The 2015 update will soon be produced and will coincide with an unprecedented coverage of remote sensing data that can be used as input for its production. Previous updates were primarily based on a semi-automated approach (mainly due to the limited availability of input data), and the increased availability of remotely-sensed data can increase the level of automation to produce and update such a large dataset. In addition, because of cloud coverage in a large part of Europe, this amount of images can ensure cloud-free and uniform results at a spatial resolution less than or equal to 20 m to fit the HRL IMD specification. Large datasets are also critical to characterize vegetation phenology false positives to reduce potential confusion between induced agricultural bare soil and artificial surfaces. All of these contributions expected with the use of many data require the design of innovative and efficient methods to process such an amount of data in order to maintain Copernicus Land products’ data production costs down.

The substantial increase in remotely-sensed data availability is primarily due to the launch of Sentinel-2A in 2015, which was designed specifically to meet the needs of the Copernicus program. The full Sentinel-2 constellation (including Sentinel-2A and -2B) will have a revisit cycle of less than five days globally (whereas Landsat-8 is around 16 days) and will provide 13 bands from visible to Short Wave Infrared (SWIR) associated with three spatial resolutions from 10 to 60 m. The multi-temporal resolution ensures a better monitoring of land use and cover with better opportunities to obtain cloudless mosaics; the wide spectral resolution facilitates the thematic identification of land cover [2], while the high spatial resolution allows for the identification of small objects, such as individual houses or landscapes structures [3].

However, because the launch of Sentinel 2A is very recent, few studies exist on how it could be used for monitoring urban areas [4]. Therefore, in this paper, a method based on exploiting Sentinel-2’s unique properties is proposed to update the Copernicus High Resolution Imperviousness Layer. However, the large amount of data generated raises methodological questions to deal with reliable, uncertain and contradictory information. These issues are discussed below.

1.2. Updating Copernicus Products: Change Detection and Data Fusion

Two main approaches are possible to detect structural changes in remote sensing data: directly detect changes from a pair of raw data or performing single classifications that are compared to highlight altered regions. The first approach is often preferred since a single process has the advantage of preventing errors from being accumulated through the two classification steps. However, in the context of the “online” updating of urban areas where classified regions in the first image are already available, the comparison with a new classification issued from the new image appears more rational. This is the process chosen in this paper. Here, in order to prevent from the accumulation or errors issued from independent classifications, we propose to accurately perform this step using specific data fusion rules. Before describing our approach, we first review briefly the main change detection approaches:

• Image to image comparison: A naive way to detect changes is to rely on differences, either computed from the input images or on features derived from the images. Relying on the simple difference of raw data is limited (because of known problems related to differences in terms of acquisition), and a substantial number of techniques has been proposed to provide reliable measurements [5,6,7,8], for example statistical hypothesis tests [9] or statistics on the spatial relations between points [10,11]. From the feature differences, the detection of changes can also be performed using supervised learning techniques [12]. The reader can refer to [13] for more details on associated approaches.
• Post-classification comparison: Classification plays a key role in change detection, as it often enables one to refine the results. Strategies range from the classification of the various objects in the images [14,15,16] (in this case, the change is immediate by comparison of the two classified data) to post-classifications or to classification of the change features [17,18,19]. As for the first situation, the quality of the result is strongly related to the segmentation. In [20], the authors have found that post-classification usually underestimates the areas of land cover change, but where the change was detected, it was overestimated in magnitude. Therefore, classifying change features seems a more appealing choice. In that context, one can mention the work of [21], where the authors use an unsupervised classification mixed with a visual interpretation of aerial photographs to detect land cover changes. We refer the readers to the complete reviews in [22,23,24,25] for an overview of the change detection issue.

Though well performing techniques exist, only a few of them specifically take into account the problem of handling accurately the errors issued from independent classifications. In [26], the authors propose to fuse the classification results in order to minimize the omission or the commission errors. The approach introduced in this paper relies on the same idea. More precisely, we use the Dempster–Shafer fusion theory, since this framework is especially adapted to deal with uncertainty that can be derived according to previous classification performances (this theory is presented in Section 3.4). In addition, our process deals with cloud coverage, enabling a continuous growth of the area covered by our approach.

Classifications of single remote sensing images generally contain errors due to constraints originating from both the data and the algorithms used. This inaccuracy (or uncertainty associated with classification results) depends on:

• the characteristics of the sensor:

  • its spatial resolution (pixels contain different ratios of pure and mixed materials)
  • its spectral resolution and the acquisition date;
• the efficiency of the classification. Though in general a simple algorithm such as the maximum likelihood can be expected to be less efficient than a more complex one, such as Support Vector Machine (SVM) or Random Forest (RF), no general rule can be extracted;
• the accuracy of the labeled data used to train and test the classifiers.

By taking into account all images issued from existing satellites, a large number of data can be acquired on a given area. However, because of the different nature of associated images, their classification may provide results associated with various levels of quality. Although selecting the best result among all available classifications would seem a rational approach, combining them by taking into account their qualities should make it possible to reach a higher level of accuracy. This is the idea behind the concept of data fusion. A large number of techniques is available to fuse data, and grouping them into families is a difficult task. Nevertheless, one can roughly distinguish the two main groups of techniques based on:

• probability theory, which aims at determining the state of a variable given various observations. For these families, we mainly find the Kalman filter and other data assimilation techniques depending on the presence of models for sensors (see, for example, [27,28]);
• evidence theory, where each decision is represented with a belief function associated with uncertainties (or in some approaches, a paradox). In these families, we find Dempster–Shafer Theory (DST) and its variants [29] or the Dezert–Smarandache paradoxical approach [30].

In a remote sensing context, the Dempster–Shafer Theory has mainly been applied to facilitate the decision process for the classification in urban environments [31,32], but also in other various contexts (sea breeze front identification, forestry monitoring or more generally, land cover [33,34,35,36]). Though widely used in remote sensing data, the DST has only rarely been applied in the context of large dataset classification (such as a country or continent) where information availability issues are critical (in terms of repetitiveness, cloud cover, acquisition window, etc). As will be shown later, this idea is explored in this paper where, by fusing individual classifications acquired at different dates, urban classifications are generated combining both an accuracy greater than the individual ones and cloud-free coverage. It is also proposed to combine Sentinel-2 with Landsat-8 images and to evaluate the potential of our approach as a multi-source fusion.

1.3. Highlights of the Paper

In this paper, it is proposed to exploit Sentinel-2 data to design a fully-automatic technique that updates imperviousness over Europe by ensuring:

• a more precise quality than existing products;
• an associated degree of confidence;
• less sensitivity to cloud coverage.

As will be shown in this paper, the gain in terms of quality is due to the use of Dempster’s fusion rule that accurately deals with uncertainties issued from all classifications. It also enables one to give a degree of confidence related to all estimations. As for the last point, because of the cloud coverage all over member states’ territories, the production of cloud-free data requires a huge amount of images to ensure uniform results at a spatial resolution less than or equal to 20 m to fit the specification of Copernicus products.

2. Study Areas and Data

2.1. Study Areas

We choose to validate our approaches on two cities of different sizes with different urban planning policies: Prague (over one million inhabitants, 6900 km²), the capital of the Czech Republic, and Rennes (over 200,000 inhabitants, 2500 km²), a mid-sized city in France. The areas are delimited by the contour of the Copernicus Former Urban Areas (FUA) of the Urban Atlas dataset and are visible in Figure 1. In practice, their internal organization is different: while Rennes is mainly composed of medium-sized buildings, individual houses and industrial/commercial areas, Prague is larger and embeds a large historical center connected with big avenues together with old and high level apartments and some modern skyscrapers in its periphery. Therefore, these characteristics are complementary and illustrate the variety of situations one encounters in practice.

2.2. Data

We took a pair of Sentinel-2 images for each study site, the time shift between the two acquisitions being about a week. This is too short to observe significant vegetation changes, but these two images remain interesting to overcome gap filling issues (mainly because of clouds).

To simulate more complete series of Sentinel-2, we also added some Landsat-8 images, since this sensor is at the moment the most similar in terms of specifications (spectral and spatial resolution). Landsat-8 has 11 bands from coastal blue to middle infra-red and has spatial resolution varying from 15 m to 60 m. Images acquired in spring (or early summer) and autumn were selected. All characteristics of the data are visible in

Table 1. List of input images.

#	Prague	#	Rennes
1	12/10/2015- Landsat 8	1	09/09/2015 - Landsat 8
2	30/08/2015 - Sentinel-2A	2	28/08/2015 - Sentinel-2A
3	13/08/2015 - Sentinel-2A	3	21/08/2015 - Sentinel-2A
4	17/07/2015 - Landsat 8	4	11/05/2015 - Landsat 8

. This time period is more suitable to take into account different seasonal agricultural practices and therefore to evaluate the ability of our approach to detect urban changes while being insensible to agricultural ones. These Landsat images are first used as a comparison with the Sentinel-2 dataset. In a second step, they are included in the classification fusion process to assess the potentialities of Sentinel-2 in a multi-sensor approach.

3. Method

The overall approach is composed of six steps:

• Data preparation: extracting the cloud mask and texture descriptors of all images
• Supervised data learning: sample data from all images and build classification models
• Map production: classify each datum from models of the previous stage
• Classification fusion: combine all classifications
• Change detection: detect changes on temporal data
• Accuracy assessment: validate our products

The workflow is presented in Figure 2, and each step is detailed in the followings sections.

3.1. Data Preparation

The data preparation step consists mainly of two operations: cloud mask detection and computation of the texture descriptor for each image. As they are classified independently, no specific spectral processing, nor radiometric corrections are required on each datum.

3.1.1. Cloud Mask Detection

We perform cloud detection in a semi-automatic way through k-means classification on Band 1 (coastal blue), Band 8A (near-infrared) and Band 10 (water vapor), which are the most interesting to highlight cloud. Because of the differences in these bands, clouds and shadows fall in distinct clusters that are manually labeled as clouds or shadows. Here, we perform a manual labeling since we have few images, but in an operational context, automatic multi-temporal approaches could be used, such as MTCD (Multi-Temporal Cloud Detection [37]). Their final shapes are delineated after a morphological opening in order to prevent from border effects.

This method has the advantages of being simple, not requiring atmospheric correction and providing efficient results. An example in Figure 3 presents the result of the cloud mask delineation over Prague. Even though the cloud opacity is low, we observe in this image that clouds are accurately detected.

Concerning Landsat-8 images, cloud masks were directly downloaded from the USGS servers. These cloud masks rely on the Fmask (Function of mask) method, which is dedicated to Landsat missions and extracts both clouds and shadows [38].

3.1.2. Computation of the Texture Descriptor

To characterize the content of high or very high resolution remote sensing images, texture descriptors are powerful, especially in urban areas, as proven by numerous studies [39,40,41,42,43,44]. Among the existing techniques, the so-called PANTEXmethod [44], which relies on the analysis of the Grey Levels’ Coocurrence Matrices (GLCM) and the computation of associated indexes [45], has the properties of being simple and efficient. To be invariant by rotation and less sensitive to small edges that appear in rural or natural landscapes, PANTEX keeps the minimum value of the energy computed on GLCM in various directions (10 in practice [44]). This enables one to highlight local contrast variation, which is a particularly crucial criteria in an urban area. PANTEX has already proven its efficiency in urban areas at a global level [46], and an example of PANTEX criteria is visible in Figure 4. In practice, PANTEX criteria are computed in the blue band (10-m resolution) in the case of Sentinel-2 and in the panchromatic one (15 m resolution) with Landsat-8.

3.2. Supervised Data Learning

This supervised data learning step aims at sampling data and building a classification model per image that will be applied in the further step (Section 3.3). Two classes, urban areas and the rest of the landscape, are foreseen in the classification output. To extract learning points inside urban areas, we use the Copernicus High Resolution Layer Imperviousness Degree 2012 (HRL IMD 2012) as a reference. It is indeed assumed that urban areas in 2012 are still labeled as urban regions in the future. Indeed, urban areas grow slowly and on some small areas compared to the entire arable land available. More quantitatively, given that urban areas represent just over 1% of the overall land in Europe and they grew by about 0.6% per year between 1990 and 2000 [47], the probability of sampling changes remains below 2%. Changes selected in a random sampling can then be considered as noise in the data source and should be neglected. This enables one to provide a large number of samples over the entire area on which a classification model can be designed.

In practice, a maximum sampling of 1,000,000 pixels over the image is performed. The number of samples is not equal for each class, but depends on their proportion in the image. For example, given that the urban areas represent about 4% of the image, the distribution of samples is about 4% in the ‘urban area’ class and 96% for the ‘non-urban’ one. This is more adapted to represent the different spectral characteristics of the rural landscape’s classes (vegetation, bare soil, water, forest).

The classification algorithm used for this application is bagging trees [48], which is not sensitive to a large number of training data. Bagging is part of the ensemble methods, which also include random forest and boosting. Ensemble methods have already proven their efficiency for the processing of remote sensing images, in particular in urban environments [44,49]. They constitute a set of techniques for improving the performance of simple decision trees. The global idea is to perform multiple random samplings and to construct separate tree models. Each tree provides a criterion related to the fact that the input data belong to the urban class or not. Then, a choice procedure on the basis of each result is applied (average in general).

Though many other classification techniques could be used, bagging remains a suitable choice to deal with a large number of observations (as well as neural networks, for example), since their calibration is fast. Moreover, with such an approach, it is also possible to easily evaluate the discrimination power of each input variable by evaluating their impact on each node, which is of prime importance to assess the discriminative power of Sentinel-2 bands [2].

All spectral bands are kept in the feature vector on which we add PANTEX criteria. Features are then resampled by nearest-neighbors at ($20\times 20$ m) with respect to the specifications of the Copernicus IMD HRL product.

3.3. Map Production

Each image has been classified using its own bagging tree model and validated by cross-validation. Therefore, in each pixel $(x,y)$, a series of criteria $\mathcal{C}(x,y,t)$ (corresponding to the concatenation of the results of all models for each time t) that belongs to the urban thematic class is available. A global quality criterion can then be estimated to evaluate its accuracy and provide an overall uncertainty for each date. Based on the kappa index, the global uncertainty ${\Delta }_{global}\left(t\right)$ at time t is computed as:

$${\Delta }_{global}\left(t\right)=1-k\left(t\right)$$

(1)

where $k\left(t\right)$ is the kappa index of the classification at time t. With this criterion, a good classification will generate low global uncertainty and vice versa. This ensemble of criteria associated with their global uncertainty needs now to be fused in each pixel to generate a global map. This is the scope of the next section.

3.4. Classification Fusion

In this step, each pixel contains a value related to its similarity to the urban class associated with a global uncertainty. We now need to synthesize this information to provide a unique map. To this end, we rely on evidence theory and in particular on the Dempster–Shafer Theory of evidence (DST). The DST is based on a Bayesian approach and fuses a set of mass functions issued from various sources of observations associated with the belief on some hypotheses. A key advantage is that uncertainty (union of hypotheses) is accurately managed with the Dempster’s fusion rule. Applied to our classification problem, let $\mathcal{U}$ and $n\mathcal{U}$ be the hypotheses that a given pixel belongs to the urban class or not. The uncertainty is formalized as the union $\mathcal{U}\cup n\mathcal{U}$. For each time step t, if one computes three belief functions $m_t^{\mathcal{U}}$, $m_t^{n\mathcal{U}}$ and $m_t^{\mathcal{U}\cup n\mathcal{U}}$, such that $m_t^{\mathcal{U}}+m_t^n\mathcal{U}+m_t^{\mathcal{U}\cup n\mathcal{U}}=1$ and:

• $m_t^{\mathcal{U}}$ is the belief that the corresponding pixel belongs to the urban class;
• $m_t^{n\mathcal{U}}$ is the belief that the corresponding pixel does not belong to the urban class;
• $m_t^{\mathcal{U}\cup n\mathcal{U}}$ is overall uncertainty;

it is possible to combine two sources of information at times $t_1$ and $t_2$ (associated with mass functions $\{m_{t_1}^{\mathcal{U}},m_{t_1}^{n\mathcal{U}},m_{t_1}^{\mathcal{U}\cup n\mathcal{U}}\}$ and $\{m_{t_2}^{\mathcal{U}},m_{t_2}^{n\mathcal{U}},m_{t_2}^{\mathcal{U}\cup n\mathcal{U}}\}$, respectively) using Dempster’s fusion rule. This is given, by assumption $\mathcal{H}\in \Theta$ of the powerset $\Theta =\{\mathcal{U},n\mathcal{U},\mathcal{U}\cup n\mathcal{U}\}$, by:

$$m\left(\mathcal{H}\right)=\left[m_{t_1}\oplus m_{t_2}\right]\left(\mathcal{H}\right)=\frac{{\displaystyle \sum _{\mathcal{A},\mathcal{B}\in {\Theta }^2,s.t.\mathcal{A}\cap \mathcal{B}=\mathcal{H}}}{m}_{t_1}\left(\mathcal{A}\right)m_{t_2}\left(\mathcal{B}\right)}{\underset{\mathcal{K}}{{\displaystyle \underbrace{1-\sum _{\mathcal{A},\mathcal{B}\in {\Theta }^2,s.t.\mathcal{A}\cap \mathcal{B}=\varnothing }{m}_{t_1}\left(\mathcal{A}\right)m_{t_2}\left(\mathcal{B}\right)}}}}$$

(2)

Roughly, the numerator combines the4 information of the two sources $m_{t_1}$ and $m_{t_2}$ in hypothesis $\mathcal{H}$, while the denominator $\mathcal{K}=1-{\displaystyle \sum _{\mathcal{A},\mathcal{B}\in {\Theta }^2,s.t.\mathcal{A}\cap \mathcal{B}=\varnothing }}{m}_{t_1}\left(\mathcal{A}\right)m_{t_2}\left(\mathcal{B}\right)$ is related to the conflict (its second term is 0 if sources are in accordance and grow with conflict). As this fusion rule is associative, to combine N sources of information, we start by fusing the two first ones, then fusing the result with the third one, and so on (this reads: $m_f\left(\mathcal{H}\right)=\left[\left[\left[m_{t_1}\oplus m_{t_2}\right]\oplus m_{t_3}\right]...\oplus m_{t_N}\right]\left(\mathcal{H}\right)$, and only pairwise fusions defined in Equation (2) are involved.).

To apply this fusion rule to the classification of urban patterns, we need to define the mass functions $m_t^{\mathcal{U}}$, $m_t^{n\mathcal{U}}$ and $m_t^{\mathcal{U}\cup n\mathcal{U}}$ at each time t. Let us remind that, after the classification step, each pixel contains in each time t a value $\mathcal{C}(x,y,t)$ related to its similarity to the urban class. In addition to this value, a global accuracy ${\Delta }_{global}\left(t\right)$ (see Equation (1)) is also available in each time step.

The way we get $m_t^{\mathcal{U}}(x,y)$, $m_t^{n\mathcal{U}}(x,y)$ and $m_t^{\mathcal{U}\cup n\mathcal{U}}(x,y)$ in $(x,y)$ is therefore a normalization step:

$$\mathcal{S}\left(t\right)=1+{\Delta }_{global}\left(t\right)\mathrm{and}\left\{\begin{array}{c}{m}_t^{\mathcal{U}}(x,y)=\mathcal{C}(x,y,t)/\mathcal{S}\left(t\right)\hfill \\ m_t^{n\mathcal{U}}(x,y)=\left(1-\mathcal{C}(x,y,t)\right)/\mathcal{S}\left(t\right)\hfill \\ m_t^{\mathcal{U}\cup n\mathcal{U}}(x,y)={\Delta }_{global}\left(t\right)/\mathcal{S}\left(t\right)\hfill \end{array}\right.$$

(3)

In practice, to accurately deal with clouds, each pixel identified in Section 3.1 as a cloud has uncertainty $m_t^{\mathcal{U}\cup n\mathcal{U}}$ equal to 1.

Finally, once fused, the final decision of the classification between urban/non-urban areas is chosen as the maximum between $m_f^{\mathcal{U}}$ and $m_f^{n\mathcal{U}}$ in each pixel $(x,y)$ (with $m_f$ the fused mass functions).

3.5. Change Detection and Imperviousness Degree Computation

Each time a new datum is available, a change detection can be performed by the comparison of classifications. This enables one to easily update databases, as will be shown in the experimental part on the Copernicus HRL IMD 2012, which has been validated by the EEA (European Environment Agency).

The derivation of imperviousness degree (1% to 100%) is produced using an automatic algorithm based on a calibrated Normalized Difference Vegetation Index (NDVI). A description of the Copernicus Imperviousness Layer methodology was described by [50] for the 2009 update and by [51] for the 2012 update. Similar methods were also applied in the USA for the development of the National Land Cover database [52].

3.6. Change Detection Validation

Change detection validation is made from a stratified sampling on the change stratum and unchanged stratum.

A density threshold of 30% was used to derive the built-up layer from the imperviousness layer [53]. This was not intended to be a separate product, but instead was calculated for the verification process only, because density products cannot be verified.

Based on the assumption that urbanized areas spread more than they appear randomly in the landscape, sampling unchanged areas was done in a 250-m buffer zone along urbanized areas.

In practice, we follow the same process as the Copernicus products to validate our results [54]: 200 samples were selected including 100 samples in the change stratum and 100 samples in the unchanged stratum. Each sample represents a $100\times 100$ m polygon in which a manual interpretation was performed. The assigned class matches the major land use surface inside a polygon.

4. Results and Discussion

The overall processing chain presented in this paper has been applied on the Sentinel-2 and Landsat-8 images for the two study areas. Independent classifications have been performed in each image, and the fusion has been applied in two different ways:

• Multitemporal fusion: classifications issued from images of the same sensor have been fused. This yields two results: fusion with Sentinel-2 and fusion with Landsat-8;
• Multi-source fusion: all classifications have been fused to provide a single classification.

Results are depicted in Figure 5 and Figure 6 for Prague and Rennes, respectively, and quantitative criteria are presented in

Table 2. Intermediate result for Prague.

Classification Sensor (date)	Overall Accuracy	Uncertainty (1-Kappa)	Cloud Coverage
Landsat-8 (17/07/2015)	97.5%	0.25	21%
Sentinel-2 (13/08/2015)	97.3%	0.27	17%
Sentinel-2 (30/08/2015)	97.2%	0.28	7%
Landsat-8 (12/10/2015)	97.3%	0.28	24%
Fusion of Landsat-8 classifications	97.6%	0.24	0%
Fusion of Sentinel-2 classifications	97.3%	0.23	0%
Fusion of Sentinel-2 and Landsat-8 classifications	97.8%	0.22	0%

and

Table 3. Intermediate result for Rennes.

Classification Sensor (date)	Overall Accuracy	Uncertainty (1-Kappa)	Cloud Coverage
Landsat-8 (11/05/2015)	97.7%	0.24	9%
Sentinel-2 (21/08/2015)	97.3%	0.27	0%
Sentinel-2 (28/08/2015)	97.1%	0.30	44%
Landsat-8 (09/09/2015)	97.7%	0.24	3%
Fusion of Landsat-8 classifications	97.8%	0.23	0%
Fusion of Sentinel-2 classifications	97.3%	0.26	0%
Fusion of Sentinel-2 and Landsat-8 classifications	97.9%	0.21	0%

, both for individual classifications and fusions of them. Because of cloud coverage, results cannot be provided in each pixel, and we also show in these tables the percentage of data for which estimations are available.

From these table, it is worth noting that overall accuracies and kappa indexes are good in all situations. This first reveals the ability of our texture index based on PANTEX to accurately extract urban areas. However, looking at individual classifications, we observe that overall urban areas are not completely covered because of image properties (spatial coordinates) and cloud coverage (third line of Figure 5 and Figure 6). When fusing classifications issued from images of the same sensor (fourth line), one observes more consistent maps covering the entire region. The fusion process enables one not only to provide more complete maps, but as shown on Table 2 and Table 3, uncertainty is also reduced by combining in an efficient way all classifications, yielding overall accuracies and kappa index that are higher. The same observation holds when fusing images of various sources (last line). This demonstrates the ability of Dempster–Shafer theory to optimally combine various classifications.

In the following sections, we go into more detail, and we show the benefits of the approach to classify Sentinel-2 images and its utility to update large datasets, such as Copernicus HRL IMD.

4.1. Variable Importance of Single Image Classifications

In Figure 7, we have plotted the contribution of the different bands (all spectral values and PANTEX criteria) in the classification models of each image. The impact of each variable is evaluated during the tree construction through the following principle: a given node associated with a specific variable is split (using a threshold on this variable) if it is considered as heterogeneous. Once split, if the children nodes are not heterogeneous, therefore the contribution of this variable is significant, since it enables one to separate inhomogeneous groups into homogeneous ones. More precisely, each importance is estimated based on the differences of impurity between parent and children nodes. More details are given in [55].

Several remarks can be made from Figure 7: as one can observe with the significant contribution of PANTEX, texture plays a key role in each classification model. This is all the more true for Sentinel-2 images whose spatial resolution is higher than Landsat-8, and therefore, PANTEX is more sensitive in urban areas (since it relies on local variations of spectral values). To illustrate this, Figure 8 presents PANTEX maps extracted on Sentinel-2 and Landsat-8, and we observe sharper delineation of urban areas from Sentinel-2’s PANTEX band.

To explain the relative better efficiency of the City of Prague compared to Rennes, let us remind that this region of France is characterized by a very fragmented agricultural landscape with parcels of small sizes, unlike the Prague area, whose agricultural crops are larger. Therefore, some rural areas are more similar in terms of texture than urban ones, and this reduces the rule of PANTEX. Nevertheless, the overall accuracy is still very good and even larger for the city of Rennes.

As for the spectral bands, Figure 7 shows that the near-infrared (B8, B8A Sentinel-2 and B5 for Landsat-8) and blue (B01, B02) bands are particularly valuable to extract urban areas. These bands are then followed by the Short Wave Infrared (SWIR) bands (B10, B11, B12 for Sentinel-2 and B6 for Landsat). These most important Sentinel-2 bands are in accordance with [2]. Man-made features have low reflectance in the near-infrared and can be separated from the vegetation. They also have strong reflectance in the blue and SWIR that can separate them from the rest of the agriculture or semi-natural features.

4.2. Benefits of Classification Fusion

The fusion process enables one to provide a cloud-free classification and also to improve the overall accuracy, as shown in Table 2 and Table 3. This illustrates the complementarity of available sources:

• Multitemporal data (one of the main pros of Sentinel-2) enable one to consider different phenological stages of the agricultural landscape and to limit commission errors, such as bare soil and other low chlorophyll activity land covers. As observed in Figure 9, the first image shows a set of bare soils in which urban membership criteria are high. On the contrary, the second image shows the set of parcels with high chlorophyll activity and low membership values. The fusion of both criteria results in low membership values. These differences are highlighted by a high value in the conflict map $\mathcal{K}$ of Equation (2) (Figure 9f). A comparison with the mean of classifications shows that urban estimates are higher and easily separable (Figure 9g).
• Multi-source fusion (combination of Landsat-8 and Sentinel-2 classifications) enhances also the geometric accuracy from 30 m to 10 m. As presented in Figure 10, the delineation of the urban areas is more sharp on Sentinel-2 image. Fusion of both classifications provides a sharp delineation with high membership values. The conflict map shows also that differences between images were mainly located in the edges of the urban areas. Figure 10g also shows that classification fusion provides sharper edges than the mean of classifications.

4.3. Application to Change Detection and Copernicus Products’ Updates

The resulting classifications were then compared to the Copernicus HRL IMD 2012 to detect changes between 2012 and 2015. The change maps were validated with an overall accuracy over 94.75% and a kappa index over 0.90 (

Table 4. Change detection validation for Prague (OA = 95%; K = 0.90).

		Reference
		No Change	Change	Total
Classification	No Change	194	14	208
	Change	6	186	192
Total		200	200	400

and

Table 5. Change detection validation for Rennes (OA = 94.75%; K = 0.90).

		Reference
		No Change	Change	Total
Classification	No Change	183	4	187
	Change	17	166	213
Total		200	200	400

). Commission errors are mainly due to the detection of water (lake, reservoir) or long-term bare soil. On the contrary, omission errors correspond to some mixed urban and vegetation patterns, such as residential subdivisions, where the vegetation is particularly dense.

We illustrate in Figure 11 and Figure 12 the changes observed between our product (issued from multi-temporal multi-source fusion of Sentinel-2 and Landsat-8) and the Copernicus HRL IMD 2012 derived from it. As one can observe, changes are mainly identified as pixel blocks along the urban areas, illustrating the urban growth.

To visually analyze these changes, we have depicted in Figure 13a, Figure 14a and Figure 15a Google Earth image close to 2012 (the date of the generation of the Copernicus HRL IMD; Figure 13a, Figure 14a and Figure 15a), the latest Sentinel image in 2015 (Figure 13c, Figure 14c and Figure 15c) and the detected changes (Figure 13b, Figure 14b and Figure 15b). From these figures, it is interesting to observe that either new commercial or industrial areas (Figure 13b) or thin transport infrastructures (Figure 14b) and residential subdivisions (Figure 15b) are accurately detected, illustrating the efficiency of our process to detect various kinds of urban areas.

5. Conclusions

In this work, we have proposed an exploitation of Sentinel-2 images to monitor urban areas and investigate their utility to update Copernicus Land products. The proposed approach takes advantage of the regular availability of new images provided by Sentinel-2 to reduce the uncertainty in the detection of urban areas and the estimation of imperviousness. It relies on several separate image classifications that are fused using the Dempster–Shafer theory. The regular acquisition of images feeds our processing chain to enhance the classification accuracy and provide a cloud cover-free product. In practice, in addition to all spectral bands, a texture index (PANTEX) has been added as a feature descriptor to accurately represent urban areas.

The approach has been evaluated on two large urban areas (Prague in the Czech Republic and Rennes in France) following the Copernicus process. Results exhibit very good accuracies with a kappa index over $0.9$. This illustrates the ability of Sentinel-2, both in terms of spectral and spatial (through the texture index) resolutions, to detect urban areas.

We also have successfully combined Sentinel-2 and Landsat-8 data to show that Sentinel-2 images increase the geometric accuracy of Landsat-8 classifications, and conversely, Landsat-8 increases the Sentinel-2 repetitiveness and its thematic accuracy.

Therefore, the presented approach is applicable both to heterogeneous landscapes and heterogeneous data. It is also able to reach the needs of operational projects, such as the Copernicus High Resolution Layer Imperviousness, and can be conducted on very large areas, such as Europe. Though applied to urban areas, the methodology can be used in any other land cover classification tasks.