<figure><img src="image_1.png" /></figure>

HAL
open science

# Application of the symbolic machine learning to copernicus VHR imagery: the european settlement map

Christina Corbane, F. Sabo, V. Syrris, T. Kemper, P. Politis, M. Pesaresi, P.

Soille, K. Ose

► To cite this version:

Christina Corbane, F. Sabo, V. Syrris, T. Kemper, P. Politis, et al. Application of the symbolic machine learning to copernicus VHR imagery: the european settlement map. IEEE Geoscience and Remote Sensing Letters, 2020, 17 (7), pp.1153 - 1157. 10.1109/LGRS.2019.2942131. hal-02610302

HAL Id: hal-02610302

https://hal.inrae.fr/hal-02610302

Submitted on 3 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Application of the Symbolic Machine Learning to Copernicus VHR Imagery: The European Settlement Map

C. Corbane<sup>1,2</sup>, F. Sabo<sup>1,2</sup>, V. Syrris<sup>1,2</sup>, T. Kemper<sup>1,2</sup>, P. Politis, M. Pesaresi<sup>1,2</sup>, P. Soille<sup>1,2</sup>, and K. Osé

**Abstract**— This letter introduces the new European Settlement Map (ESM) production workflow, presents some indicatory results, and discusses the validation process. Unlike the previous ESM versions, the built-up (BU) extraction is realized through supervised learning (not only by means of image filtering and processing techniques) based on textural and morphological features. Input data are the Copernicus very-high-resolution image collection for the reference year 2015 coming from a variety of sensors. The workflow is fully automated, and it does not include any postprocessing. For the first time, a new layer containing nonresidential buildings was derived by using only remote sensing imagery and training data. The produced BU map is delivered at 2-m-pixel resolution (level-1 layer), while the residential/nonresidential layer (level 2) is delivered at 10-m spatial resolution. More than 46000 scenes were processed, and around 6 million km² of the European territory was mapped. The workflow was executed on the JRC Big Data platform. The validation showed a balanced accuracy of 0.81 and 0.91 for level 1 and 2 layers, respectively, and 0.71 for only the nonresidential layer.

**Index Terms**— Building typology, built-up (BU), classification, Global Human Settlement Layer (GHSL), remote sensing.

## I. INTRODUCTION

SUSTAINABLE development of urban areas is one of the major problems of modern civilization, and therefore, it is at the heart of international and European policies, such as the Sustainable Development Goals (Goal 11), the New Urban Agenda for Europe [1], and the European Union (EU) Cohesion Policy 2014–2020. Almost 70% of the EU population lives in an urban area, and these areas generate more than two thirds of the EU’s Gross domestic product (GDP). Projections show that cities, and especially capital cities, will continue to grow with higher levels of migration and natural change, while the rest of the EU will start to lose population [2]. The use of cutting-edge technologies, such as Earth observation to monitor changes in the European landscape, will be essential for informed policy and decision making. In particular, continuous mapping of human settlements along the urban–rural continuum is needed to better capture the dynamics of built-up (BU) areas (and hence, their population) and

understand the drivers and impacts of the human settlements development. Since 2014, the Copernicus Earth Observation program of the EU is providing a set of information products based on remote sensing data for urban growth. Among those, the High Resolution Layer Imperviousness Degree, the Urban Atlas (UA), and Corine Land Cover<fnref n="1" /> have become the *de facto* standard for urban studies in Europe. However, none of these products offers a wall-to-wall consistent and fine-scale representation of human settlements. The UA contains detailed nomenclature for urban areas; however, it considers only 785 cities with 0.25 hectare as a minimum mapping unit. The first product to offer such characteristics was the European Settlement Map (ESM) produced in the framework of the Global Human Settlement Layer (GHSL) project of the European Commission. This layer was extracted automatically from 2.5-m-resolution Satellite Pour l'Observation de la Terre (SPOT)-5/6 images acquired in the context of the pan-European GMES/Copernicus (Core_003) data set for the reference year 2012 (ESM_2012) [3]. The product, released in 2016, demonstrated the added value of finer spatial and thematic granularity for applications that rely on detailed urban maps with a pan-European coverage (e.g., population disaggregation [4] and assessment of human exposure to floods [5]).

Recently, as a part of the Copernicus contributing missions, a pan-European coverage (EEA39 countries) of very-high-resolution (VHR) images for the reference year 2015 (VHR_IMAGE_2015) was acquired for the purpose of continuation of local land component applications for hotspots, coastal areas, risk zones, and protected sites. The availability of this new pan-European coverage of VHR optical images for the reference year 2015 offers an interesting opportunity to derive updated settlements characteristics by revisiting the workflow used for the production of ESM_2012 and implementing novel machine learning techniques developed in the frame of the GHSL.

This letter describes the new methodology that builds on the GHSL core technology [6] for automatic extraction of BU areas. It also presents the new updated ESM based on the VHR_IMAGE_2015 collection (ESM_2015) and evaluates its quality and accuracy. The main innovations of the proposed methodology are as follows:

1) the use and adaptation of the symbolic machine learning (SML) method that has been combined with textural and morphological features;

2) the automatic extraction of building typology discriminating between residential and nonresidential BU areas.

This letter is structured in the following way. The input data are described in Section II. Section III gives an overview of the

Manuscript received June 11, 2019; revised August 30, 2019; accepted September 14, 2019. Date of publication December 24, 2019; date of current version June 24, 2020. (Corresponding author: F. Sabo.)

C. Corbane, V. Syrris, T. Kemper, M. Pesaresi, and P. Soille are with the Joint Research Centre, European Commission, 21027 Ispra, Italy.

F. Sabo is with Arhs Developments Italia S.r.l., 20121 Milan, Italy (e-mail: filip.sabo@ext.ec.europa.eu).

P. Politis is with Arhs Developments S.A., 4370 Belvaux, Luxembourg.

K. Osé is with Irstea, UMR TETIS, 34000 Montpellier, France.

Color versions of one or more of the figures in this letter are available online at http://ieeexplore.ieee.org.

Digital Object Identifier 10.1109/LGRS.2019.2942131

<footnote n="1">¹https://land.copernicus.eu/</footnote>
<figure><img src="image_3.png" /><figcaption>Fig. 1. VHR_IMAGE_2015 collection per contributing mission.</figcaption></figure>

ESM_2015 workflow, including a concise description of the SML classifier. Results and validation experiments are summarized in Section IV, followed by a discussion in Section V about the strengths and limitations of the methodology as well as directions for future work.

## II. INPUT DATA

The Copernicus VHR_IMAGE_2015 collection provides an optical VHR multispectral (MS) and panchromatic coverage of the EEA39 countries. Even though it consists of commercial imagery, the collection is free and open for Copernicus users. The data have been acquired within pre-defined time windows corresponding to the vegetation season as regards to the years 2014–2015–2016. The collection is composed of heterogeneous imagery coming from a consortium of image providers. The following satellite missions were selected for this collection: Pleiades (~27000 scenes) from Airbus, Deimos-2 (~4887 scenes) from Deimos Imaging, WorldView-2, WorldView-3, and GeoEye-01 (~7896 scenes) from European Space Imaging (EUSI). Additional scenes from SPOT 6/7 (Airbus) were added in the final gap-filling phase. The data were downloaded and stored on the Joint Research Centre Earth Observation Data and Processing Platform (JEODPP) [7]. The total size of the downloaded data is 127 terabytes. All the data were delivered as a four-band composite (blue, green, red, and infrared) with a separate panchromatic image. The image projection is the ETRS89-Lambert azimuthal equal-area coordinate reference system (EPSG:3035). The VHR_IMAGE_2015 collection per contributing satellite mission is shown in Fig. 1.

Airbus and Deimos data are scene based, while the EUSI data were delivered in tiles. The image size ranges from 4096 × 4096 pixels to more than 15000 × 15000 pixels. The panchromatic spatial resolution is in the range 0.4–1 m, while the MS spatial resolution ranges from 1.4 to 6 m. The Deimos sensor captures the data as 10-bit values, while the Airbus and EUSI missions collect the data in a 12-bit range. All the image data have been resampled to 16-bit radiometric resolution and encoded in Tagged Image File Format (TIF) or JPEG2000 file format.

## III. ESM PRODUCTION WORKFLOW

Despite the technical specifications and the data set descriptions reporting about 98.5% of cloud-free coverage, residual

<figure><img src="image_4.png" /><figcaption>Fig. 2. ESM_2015 training data, input imagery, and simplified workflow.</figcaption></figure>

clouds were observed in the VHR data. Persistent cloud coverage was noticed in Scandinavian countries, Iceland, and the areas close to the Alps in Switzerland, Austria, and Germany. In addition to this, data gaps over BU areas were noticed in large parts of England and Ireland, smaller parts of Germany, and in Austria, close to the Alps.

An automated cloud screening technique was applied to the input data to evaluate the percentage of cloudy scenes and, hence, exclude them from the VHR collection. The automatic cloud screening is based on the brightness values of the input bands with additional morphological filtering followed by dilation and erosion operations. Around 6500 tiles and scenes were removed due to their cloud contamination, which was higher than the selected 30% threshold. More precisely, 8% of the total Deimos scenes and 3% of the total EUSI tiles were excluded.

The ESM workflow has been designed for extracting BU areas primarily at the finest spatial resolution provided by the input data. However, due to the relevance and the potential usage of the fine-scale maps of BU areas for population disaggregation, the workflow was expanded with additional information extraction modules resulting in a two-level architecture (see Fig. 2).

Unlike the workflows for the ESM_2012 production, which mainly relied on a decision-making process based on auxiliary information about BU presence/absence, the building block of the ESM_2015 workflow is the SML algorithm [8] in combination with textural (Pantex) [9] and morphological features [characteristic saliency and leveling (CSL)] [10]. In similar settings, the performance of the SML was compared in [11] against other available classification techniques, such as maximum likelihood, random forest, and support vector machine. According to those results, the SML outperformed both parametric and nonparametric classifiers in the task of statistical downscaling of uncertain and noisy reference sets. Therefore, the SML was selected for the ESM_2015 production.

### A. SML

The SML method is based on the concept of the empirical multidimensional histogram where associations among class labels and tuples of features derived from the input imagery are specified. A normalized measure scores these tuples according
to the number of occurrences in each class. The analysis of the cumulative distribution of the estimated scores is used at the end in their reclassification [12].

The SML schema is based on two steps: 1) data reduction and 2) association analysis (AAx). In the data reduction step, the distribution of every SML input vector, i.e., of each feature extracted from the MS image, is partitioned in few groups using a uniform quantization approach. The parameter that defines the quantization levels is sensor and landscape specific. Moreover, different quantization levels are necessary to be applied since the sensors collect the data in a bit depth prescribed by their design (10 or 12 bits as regards to the data set used in this study). After data reduction, the unique data sequences are identified together with their frequencies. The sequences are then used in the AAx step where they are scored according to the frequencies of their joint occurrences with positive and negative data instances from the learning set. The output of the AAx processing stage is the evidence-based normalized differential index (ENDI), which is a continuous layer ranging from -1 to +1. Values close to +1 indicate that the data sequence is in a strong association with the class of interest (e.g., BU area).

To achieve a binary classification (i.e., BU versus non-BU surfaces), a cutoff value of ENDI needs to be defined, either manually or automatically; in that way, each data sequence is assigned to a single class. For the data set presented here, the Otsu thresholding approach has been used for automatically dichotomizing the ENDI output. The Otsu method chooses an optimal threshold by minimizing the within-class variance and maximizing the between-class variance [13].

The adaptation of the SML framework for the purpose of extracting BU areas and for describing the building typology is presented in Sections III-B1 and III-B2.

## B. Two-Level Classification

The ESM_VHR_2015 workflow is divided into two levels of processing corresponding to two different output products: level 1 for BU areas extraction and level 2 for the discrimination between residential and nonresidential BU.

In addition to the four spectral bands (radiometric features), the textural (Pantex) and morphological (CSL) spatial features were extracted from the panchromatic image in order to improve, as much as possible, the BU detection in complex urban environments.

The Pantex is an edge density detector and can be considered as an index indicating the presence of BU. The higher the edge density is, the higher the probability of a pixel is to contain information about the presence of a building in its neighborhood. The window size for the calculation of the Pantex corresponds to the local area of the textural measure (contrast). The window size value is estimated empirically. In this letter, several values for the Pantex window size were tested, ranging from 25 × 25 m to 300 × 300 m. It was demonstrated that increasing gradually the window size leads to better detection of large buildings and industrial facilities.

CSL is a compact representation of the multiscale information that is conveyed by the differential morphological profiles (DMPs) computed over the raster layers. The CSL reduces the dimensionality of the DMPs by mapping three layers, the scale, the saliency, and the level, features that highlight the dominant image components and actually lead to the segmentation of the image. It is used for object classification based on shape, color, and other features.

The unique sequences of all the radiometric and morphological features are then associated with the annotated learning set through the AAx step.

Two different learning sets, available at European scale, were used for guiding the SML training:

1) latest ESM_2012 with 2.5-m-pixel size;

2) the Enhancing Activity and Population Mapping (ENACT) Map [14] with 100-m-pixel size.

The class representing BU areas from the ESM_2012 (class 50) was used as a learning set for BU area detection.

The ENACT map is a thematic (50 classes) and spatial refinement (1-ha minimum mapping unit) of the Corine Land Cover. For defining the learning set for nonresidential areas, the ENACT map was combined with the ESM_2012 class 50. The intersection of the ESM_2012 class 50 with commercial (class 1212), production facilities (class 1211), and port areas (class 1230) from the ENACT map was used to train the SML in the second classification stage.

1) Level 1 Classification: Delineation of BU Surfaces: Delineation of BU area is achieved with the level 1 classification. AAx for level 1 works at the scene level; it considers both the ESM_2012 training set for a specific scene and the unique sequences of radiometric and morphological characteristics extracted from that input scene. Both training and classification were performed per scene at the spatial resolution of the MS bands.

ENDI from level 1 processing was weighted next by the Pantex layer (see Fig. 2) in order to remove the overdetections in the non-BU areas. Combined with the SML, Pantex contributes to the improved BU detection by giving high confidence to BU while removing overdetections, such as roads, bare land, sand dunes, and other open spaces. However, a drawback of using Pantex is that it can remove relevant information when the edge density is low. This can be the case with large industrial facilities or large commercial and public buildings. Multiplying the ENDI with the Pantex in this case and then thresholding can lead to complete or partial removal of the buildings in the binary output layer. In order to overcome this issue, the Pantex feature was calculated for two aggregated panchromatic images considering two different window sizes as follows:

1) panchromatic image at 1 m and window size 100 × 100 m (ptx1);

2) panchromatic image at 2.5 m and window size 25 × 25 m (ptx2).

The ptx1 setup provides good results in dense urban areas; however, in the low-density urban areas, it underdetects isolated buildings. Therefore, the ptx2 setup was used to correct these omissions. The value for the window size for ptx1 setup was set to 100 × 100 m since larger values were saturating the BU detection accuracy. Increasing the window size to 100 × 100 m contributed to better detection of large residential and industrial buildings at the cost of overdetecting roads and railroads as BU areas. This type of commission error (CE) was allowed in order to capture as many BU pixels as possible. Furthermore, in order to increase the probability for the large nonresidential facilities to be captured, ENDI was not weighted by the Pantex layer in the areas where nonresidential buildings can be found as defined by the ENACT training set. These rules are explained with the
following equations using a logical OR (∨) operator:

$$BU_{ptx} = BU_{ptx1} \lor BU_{ptx2} \quad (1)$$

$$BU_{no\_ptx} = BU_{ptx} \lor BU_{org} \quad (2)$$

where BU is the binarized ENDI using the Otsu method, $BU_{ptx1}$ is the BU with ptx1 setup, $BU_{ptx2}$ is the BU with ptx2 setup, and $BU_{org}$ is the BU not weighted by the Pantex. $BU_{no\_ptx}$ is then derived by applying (2) in areas where large buildings can be found according to the ENACT learning set.

2) Level 2 Classification (Discrimination of Residential and Nonresidential BU Areas): At this level of classification, the discrimination of residential and nonresidential BU areas is achieved by combining the ENACT map with the ESM_2012, as explained in Section III-B and with the following logical AND (∧) operation:

$$LRN_{non\_res} = ESM_{2012_{50}} \land ENACT_{non\_res} \quad (3)$$

where $LRN_{non\_res}$ is the learning set used for defining nonresidential BU, $ESM_{2012_{50}}$ are the buildings from the ESM_2012 (class 50), and $ENACT_{non\_res}$ are the nonresidential pixels corresponding to the ENACT map. The training for level 2 has been implemented at the MS pixel size, but the classification has been performed at the output resolution of 10 m. The AAx was applied by considering only the nonresidential training set and unique sequences of radiometric and morphological features from the input imagery. Nonresidential buildings are larger than residential buildings and generally with different shapes and radiometry represented by white, gray, and black color. Morphological features can contribute in discriminating the residential from nonresidential buildings based on their size and shape. On the other hand, the radiometric features can also improve the SML output by giving high confidence to the pixels indicating buildings, which corresponds to a specific reflection in the visible and near infrared response domain. The output of the AAx is the confidence for level 2 layer (ENDI BU nonresidential), see Fig. 2. Generally, the confidence is high for large industrial facilities, port areas, large public buildings, stadiums, and other commercial facilities. It must be noted that at this stage, the Pantex was not used to weigh the ENDI BU nonresidential because it can lower the confidence for very large buildings. The main issue for the level 2 classification was the overdetection of nonresidential BU when using only the ENDI BU nonresidential. To bypass this issue, the ENDI from the level 1 (not weighted by the Pantex) was combined with the ENDI from level 2 by applying the fuzzy union operation, that is, the maximum of the two ENDI memberships. In this way, the Otsu threshold value is automatically increased for each scene while keeping the ENDI for nonresidential BU areas high; this operation allowed the detection of new nonresidential BU areas outside the learning zones. Finally, the level 2 ENDI was aggregated to 10-m output resolution and binarized following Otsu's approach.

## IV. RESULTS AND VALIDATION

ESM_2015 processing was performed using the JEODPP data-intensive computing platform [7]. Specifically, the JEODPP batch service with the HTCondor scheduler for high-throughput computing and the Docker-based containerization was used to deploy the jobs to the cluster nodes. At the time of the processing, the maximum processing

capacity of the platform was ~1500 processing cores with 10-19-GB RAM per core. In addition, JEODPP provides a petabyte (PB)-scale distributed storage system and a plurality of data sets. More than 46680 jobs (1 job/scene) were submitted via the HTCondor job scheduler for the processing of the SML-based workflow. The requirements for a single job in terms of processing resources were 40 GB of RAM and four processors. The total processing time, including reading input and storing output, ranged from 200 to 600 s per scene. Total time to process France was ~36 h by deploying 90 concurrent jobs according to the processing resources that were allocated exclusively for our needs by the JEODPP. Almost entire France is covered with EUSI imagery; therefore, it was the most computationally intensive.

### A. Level 1 Classification Results

The level 1 layer is provided as a discrete data set at 2-m spatial resolution. The BU detection results of level 1 processing are shown in Fig. 3 with red. Increasing the window size of the Pantex contributed to better detection of large residential buildings. The exclusion of the Pantex in the nonresidential areas improved the detection of industrial and commercial facilities with some CEs on parking lots or roads between them.

### B. Level 2 Classification Results

The level 2 layer is also a discrete data set providing residential and nonresidential BU areas at 10-m pixel size. The level 2 layer is a generalized BU area; the purpose of which is to identify nonresidential BU zones rather than single buildings. Therefore, the spatial resolution of the output is at 10 m. The results of the BU detection for the nonresidential class are shown in Fig. 3 with yellow color. In addition, new buildings outside training zones were detected, providing, in that way, not only the spatial refinement but also new nonresidential BU detections. In some cases, entire new industrial zones were detected (see Fig. 3). However, there are examples where large residential buildings were also captured as nonresidential due to their size and radiometry. As it can be seen from Fig. 3, roads, railroads, bare land, and agricultural areas were not detected as BU despite the exclusion of textural features (Pantex) at this level of classification.

### C. Validation With LUCAS

Points from Land Use/Cover Area Frame Statistical Survey (LUCAS) 2015 survey database were used in order to validate both the level 1 and level 2 layers as well as the BU area layer of ESM_2012 at European scale. Total of 340 144 LUCAS points are distributed in 28 EU countries. From the total number of points, 4275 belong to the BU area classes. LUCAS points were also used for the nonresidential layer (non-res) validation from level 2 since they contain detailed land use. If the point belongs to the BU class (classes A11, A12, or A13 from LUCAS) and the Land Use column (LU1) corresponds to non-residential buildings (industrial, commercial, and so on), then the point was considered for the validation of nonresidential BU versus all other classes (including residential BU). A buffer (circle with a radius of 2 m) was applied around LUCAS point in order to follow the theoretical definition of the points and to adapt it to the pixel size of the BU layers.

The balanced accuracy (BA) was used as the main performance metric because there is a high share of non-BU
<figure><img src="image_5.png" /><figcaption>Fig. 3. Level 1 BU detection (red) and nonresidential BU from level 2 (yellow). (a) BU detection in a dense urban area in Milan, Italy. (b) BU detection in a low-density urban area in Belgrade, Serbia. (c) Additional detection of nonresidential BU near Milan (ENACT training shown with transparent purple). (d) New nonresidential BU detections in Stavanger, Norway.</figcaption></figure>

<table><caption>TABLE I
PERFORMANCE METRICS CALCULATED FROM THE CONFUSION MATRIX ELEMENTS</caption><thead><tr><th></th><th>Balanced Accuracy</th><th>Omission error</th><th>Commission error</th></tr></thead><tbody><tr><td>ESM_2012</td><td>0.82</td><td>0.34</td><td>0.63</td></tr><tr><td>ESM_2015 Level 1</td><td>0.81</td><td>0.38</td><td>0.60</td></tr><tr><td>ESM_2015 Level 2</td><td>0.91</td><td>0.13</td><td>0.76</td></tr><tr><td>ESM_2015 Level 2 non-res</td><td>0.71</td><td>0.58</td><td>0.70</td></tr></tbody></table>

areas compared to BU. According to the performance results given in Table I, the BA of the level 2 layer is the highest (0.91). However, the CE increased to 0.76 due to the aggregation from 2 to 10 m. Consequently, the omission error (OE) dropped to 0.13. Many of the false positives (FPs) for all layers actually belong to the non-BU area or linear features in the LUCAS data set. These classes correspond to roads, railroads, and parking lots in residential areas. Despite the fully automated process and its independence from external ancillary sources, such as OSM, ESM_2015 performs very well. The ESM_2015 BU layer (level 1) has a performance similar to that of ESM_2012.

Level 2 nonresidential validation resulted in 0.71 BA, 0.58 OE, and 0.70 CE. Similar to the level 1 layer, the FPs of the nonresidential layer also belong to the class roads as identified by the LUCAS set. In addition, the CE can be attributed to the confusion between the nonresidential and residential buildings in high BU density areas.

## V. CONCLUSION

In this letter, we presented the new ESM_2015 workflow, results, and validation at the pan-European scale. The workflow is fully automatic, and it does not include any postprocessing. The workflow was split into two levels: level 1 for BU extraction and level 2 for the discrimination between residential and nonresidential BU's. The novelty of the study is twofold: the application of the SML to the Copernicus VHR data at large scale for BU extraction at 2 m and, for the first time, the production of residential/nonresidential buildings map at 10-m spatial resolution. Level 2 nonresidential layer was produced by using only remote sensing imagery, SML, and learning sets. The SML framework can provide continuous updates of

the residential/nonresidential layer without the ENACT map updates because of its incremental learning characteristics. The validation against LUCAS gave a BA of 0.81 and 0.91 for level 1 and 2 layers, respectively. A decrease in the BA was observed when only validating the nonresidential layer. The release of the new ESM presents the next step toward an operational production of VHR BU information. The limitations of the method are mainly due to the propagation of errors provided by the learning set. Ongoing experiments are also exploring different quantization approaches, other than uniform, which could improve the results of the SML. Future work will also explore the added value of convolutional neural networks for the multiclass segmentation of the input imagery [15]. The acquisition of VHR imagery every three years by the Copernicus program provides an opportunity for regular monitoring and change detection starting with this layer.

## REFERENCES

[1] P. Neyens. (Feb. 3, 2017). What is the Urban Agenda for the EU? *FUTURIUM*—European Commission. Accessed: Jun. 3, 2019. [Online]. Available: https://ec.europa.eu/futurium/en/urban-agenda-eu/what-urban-agenda-eu

[2] (2016). *State of European Cities Report*. Accessed: Apr. 17, 2019. [Online]. Available: https://ec.europa.eu/regional_policy/en/information/publications/reports/2016/state-of-european-cities-report-2016

[3] A. J. Florczyk et al., “A new European settlement map from optical remotely sensed data,” *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 9, no. 5, pp. 1978–1992, May 2016.

[4] S. Freire, A. Florczyk, D. Ehrlich, and M. Pesaresi, “Remote sensing derived continental high resolution built-up and population geoinformation for crisis management,” in *Proc. IEEE Int. Geosci. Remote Sens. Symp. (IGARSS)*, Jul. 2015, pp. 2677–2679.

[5] EU Science Hub. (May 12, 2017). *Atlas of the Human Planet 2017: Global Exposure to Natural Hazards*—EU Science Hub—European Commission. Accessed: Sep. 14, 2017. [Online]. Available: https://ec.europa.eu/jrc/en/publication/eur-scientific-and-technical-research-reports/atlas-human-planet-2017-global-exposure-natural-hazards

[6] C. Corbane et al., “Automated global delineation of human settlements from 40 years of Landsat satellite data archives,” *Big Earth Data*, vol. 3, no. 2, pp. 140–169, Apr. 2019.

[7] P. Soille et al., “A versatile data-intensive computing platform for information retrieval from big geospatial data,” *Future Gener. Comput. Syst.*, vol. 81, pp. 30–40, Apr. 2018.

[8] M. Pesaresi, C. Corbane, A. Julea, A. Florczyk, V. Syrris, and P. Soille, “Assessment of the added-value of sentinel-2 for detecting built-up areas,” *Remote Sens.*, vol. 8, no. 4, p. 299, Apr. 2016.

[9] M. Pesaresi, A. Gerhardinger, and F. Kayitakire, “A robust built-up area presence index by anisotropic rotation-invariant textural measure,” *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.*, vol. 1, no. 3, pp. 180–192, Sep. 2008.

[10] M. H. F. Wilkinson, M. Pesaresi, and G. K. Ouzounis, “An efficient parallel algorithm for multi-scale analysis of connected components in gigapixel images,” *ISPRS Int. J. Geo-Inf.*, vol. 5, no. 3, p. 22, Mar. 2016.

[11] M. Pesaresi, V. Syrris, and A. Julea, “Benchmarking of the symbolic machine learning classifier with state of the art image classification methods—Application to remote sensing imagery,” Publications Office Eur. Union, Ispra, Italy, JRC Tech. Rep. EUR 27518, 2015.

[12] M. Pesaresi, V. Syrris, and A. Julea, “A new method for Earth observation data analytics based on symbolic machine learning,” *Remote Sens.*, vol. 8, no. 5, p. 399, May 2016.

[13] N. Otsu, “A threshold selection method from gray-level histograms,” *IEEE Trans. Syst., Man, Cybern.*, vol. SMC-9, no. 1, pp. 62–66, Jan. 1979.

[14] K. Rosina, F. B. e Silva, P. Vizcaino, M. M. Herrera, S. Freire, and M. Schiavina, “Increasing the detail of European land use/cover data by combining heterogeneous data sets,” *Int. J. Digit. Earth*, pp. 1–25, Dec. 2018.

[15] V. Syrris, P. Hasenohr, B. Delipetrev, A. Kotsev, P. Kempeneers, and P. Soille, “Evaluation of the potential of convolutional neural networks and random forests for multi-class segmentation of sentinel-2 imagery,” *Remote Sens.*, vol. 11, p. 907, Apr. 2019.