Mapping and Quantifying the Human-Environment Interactions in Middle Egypt Using Machine Learning and Satellite Data Fusion Techniques

: Population growth in rural areas of Egypt is rapidly transforming the landscape. New cities are appearing in desert areas while existing cities and villages within the Nile ﬂoodplain are growing and pushing agricultural areas into the desert. To enable control and planning of the urban transformation, these rapid changes need to be mapped with high precision and frequency. Urban detection in rural areas in optical remote sensing is problematic when urban structures are built using the same materials as their surroundings. To overcome this limitation, we propose a multi-temporal classiﬁcation approach based on satellite data fusion and artiﬁcial neural networks. We applied the proposed methodology to data of the Egyptian regions of El-Minya and part of Asyut governorates collected from 1998 until 2015. The produced multi-temporal land cover maps capture the evolution of the area and improve the urban detection of the European Space Agency (ESA) Climate Change Initiative Sentinel-2 Prototype Land Cover 20 m map of Africa and the Global Human Settlements Layer from the Joint Research Center (JRC). The extension of urban and agricultural areas increased over 65 km 2 and 200 km 2 , respectively, during the entire period, with an accelerated increase analysed during the last period (2010–2015). Finally, we identiﬁed the trends in urban population density as well as the relationship between farmed and built-up land.


Introduction
Over the next 25 years, the world's population growth is expected to be concentrated in urban conglomerates within the developing world [1]. The challenges in achieving sustainable urban development will be particularly significant in Africa [2,3]. Over the last decades, urban sprawl resulted in the loss of fertile soil for agricultural food production in the Egyptian Nile Valley. From the time of the 1996 population census onwards, the policy of the Egyptian government has been to avoid new constructions in the Nile floodplain, thus encouraging people to live outside the so-called 'green land' by settling in the arid areas of the eastern and western desert plateau [4]. Despite the restrictions introduced in 1996, inner-city slums have grown and informal settlements have bloomed on the urban fringes [5].  [29]), indicating the interaction area within the yellow polygon.

The Study Area
The study area is located in Middle Egypt ( Figure 2) and includes part of the governorates of El-Minya and Asyut. It stretches from approximately 170 km to 280 km south of Caïro and comprises the regional cities of El-Minya, Mallawi and Dayrut from north to south. The study area comprises the Nile floodplain where almost all the urban areas and agricultural land is situated, but also a stretch of desert to the west and east.
Within the Nile floodplain, two major river channels are running south to north, i.e., the Nile River in the eastern part and the Bahr Youssef River in the western part. The latter branches of the Nile River near Dayrut is the main source of water flowing into the Faiyum depression in the North. The total study area comprised of 6170 km 2 .

Materials and Methods
In order to overcome the limitations of the optical sensors, we employed a satellite data fusion approach, combining Synthetic Aperture Radar (SAR) and optical multi-spectral data, covering the period from 1998 up to 2015 enabling the analysis of the human-environment interactions for almost three decades. Before 1998, there was not enough data available for the study area to apply this data fusion approach. Information on the selected SAR and optical data are provided in Table 1, whereas the footprint is shown in Figure 2.  53 8 To analyse the land cover dynamics in the Middle Egypt region, the methodology employed can be described in the following three steps ( Figure 3): (1) data preparation and extraction of land cover indicators; (2) spatio-temporal land use and land cover (LULC) mapping; (3) computation of urban and agricultural land cover changes, and changes in urban population density. Individual steps are further detailed in the following sections.

Step 1: Data Preparation and Indicators Extraction
Step 1 was implemented within Google Earth Engine (GEE) [30] for the analysis of the Landsat (5 TM, 7 ETM+, and 8 OLITIRS) and Sentinel-1 IW GRD imagery, as it allows direct access and processing of the datasets to export the indicators used as input in step 2.
However, as the historical ERS SAR and Envisat ASAR are not available within GEE, these images were processed on a Virtual Machine provided by the ESA Research and Service Support through their CloudToolbox service [31]. These ERS and Envisat SAR datasets were calibrated using the ESA specifications in NEST [32] and calibrated SAR images were geocoded using the Digital Elevation Model from the Shuttle Radar Topography Mission. The resulting datasets were grouped in three periods: (i) ERS dataset acquired in a time range centered around 1998; (ii) Envisat data acquired around 2004; and (iii) Envisat data acquired in 2010 (see Table 1).
Finally, the images were registered with the corresponding Landsat dataset, using 20 ground control points that were manually selected and uniformly distributed over the study area and represent road intersections on the satellite images. The georeferenced SAR images were resampled using a linear polynomial, obtaining a final Root Mean Square Error (RMSE) lower than 1 pixel (30 m) for all the different cases. Optical images were finally cropped to match the SAR data extent. At the end of the data preparation, all data were unified to the same image resolution, namely 30 m posting, and georeferenced in the Universal Transverse Mercator (UTM) projection (Zone 36N) and World Geodetic System 84 (WGS84) Datum.

Satellite-Derived Indicators for Land-Cover Classification
Different sets of indicators, feeding the land-cover supervised classifier, are defined for SAR and multi-spectral data. Regarding the SAR data, we selected indicators (see Table 2) with significant temporal statistical information per pixel [18]. Different land cover classes are expected to exhibit different statistical backscatter properties. Among those, urban areas have the lowest temporal standard deviation as opposed to crop fields or water, which exhibit larger fluctuation of the signal measured over time.
Regarding the multi-spectral data, instead of using single spectral bands, often a linear combination of bands is used on normalised differences of different bands [28,33,34]. Traditionally, many indexes (e.g., NDVI) have been defined as the normalised difference of two bands with the aim of highlighting specific features of the land, such as vegetation, water surfaces and built-up areas [34][35][36][37]. The advantages of this method are manifold. Firstly, data are globally consistent. Secondly, the information provided by each index gives the opportunity of analysing the contribution of different features in mixed urban areas. Finally, the risk of ambiguities is minimised since the different indexes complement each other in distinguishing the land cover classes. This approach has been found to be effective by several studies [34,[38][39][40] and has been selected as the starting point of this work.
As shown in Table 2, for SAR data, the following indicators were investigated: (i) temporal average calibrated SAR backscatter expressed in decibel scale (dB); (ii) temporal standard deviation of calibrated SAR backscatter signal per pixel in dB; (iii) temporal coefficient of variation of calibrated SAR signal per pixel in dB and; (iv) coefficient of dispersion in linear scale. Table 2. Formulas of the four selected indicators of SAR derived information employed as input in the land use/land cover (LULC) classifier. Each index is calculated within the temporal data series of each pixel i and K images, with K ∈ [1, N], being N the total number of images.

Indicators Formula
Backscatter average [dB] µ a i = 10log  Table 3, for multi-spectral data the following normalised ratios were computed: (i) Normalised Difference Built-up Index (NDBI); (ii) Normalised Difference Water Index (NDWI); (iii) Normalised Difference Vegetation Index (NDVI) and; (iv) Normalised Difference between SWIR-1 and SWIR-2. Of these indicators, the 95th percentiles, respectively, were computed from the temporal distribution (N images) of each pixel (i) and selected as the final indicators. Using this 95th percentile, we avoid artefacts in the data such as the ones produced by the cloud covered data. Table 3. Formulas of the four normalised ratios derived from each of the multi-spectral datasets. Each index is calculated within the temporal data series of each pixel i and K images, with K ∈ [1, N], being N the total number of images.

Indicators Formula
Normalised Difference Built-up Index Step 2: Land Cover Mapping In this work, we have chosen a pixel-based feature fusion approach, using the aforementioned pixel-based indicators as an input of the supervised classifier.

Supervised Classifiers
We made use of an open-source C++ suite of utilities for remote sensing (RS) image processing [41]. This suite implements support vector machine (SVM) models that are machine learning supervised learning models with associated algorithms. In particular, the suite implements a supervised classification SVM model (C-SVM) based on the library libSVM and it uses a radial basis function (RBF) kernel [42,43]. Several steps are needed in order to perform a classification: (i) the input indicators composing the dataset to be classified were scaled in order to avoid values spanning greater numeric ranges dominate those with smaller numeric ranges (e.g., normalised between −1 and 1); (ii) a training dataset (a shapefile of points) was prepared for each target class using Google Earth where each point corresponds to an array containing the correspondent values of the input indicators for that pixel; (iii) using the RBF kernel, the optimal parameters C (penalty parameter for the wrong classification) and g (transformation parameter in the kernel) were obtained maximising the accuracy in classifying test data. These C and g parameters were obtained using the pkopt function provided within the pktools package, which iteratively explores in the user-defined range of C and g values until finding the local minimum that minimises the errors produced by the classifier. At this stage, the classification engine was ready to process and classify unseen data.

Land Cover Classes
We have defined five classes to be employed to create our classification maps. Four generic classes: (i) built-up area, hereafter named 'urban'; (ii) sandy and rocky desert combined into a 'desert' class; (iii) vegetation, crops, garden, grass, and agricultural fields named as 'field' and; (iv) 'water' class mainly formed by the Nile River and smaller water bodies present in the scene. Finally, our fifth class named 'dunes' identifiable thanks to the data fusion approach, not only inside the desert but also inside the floodplain. This class is expected to be characterised by having high temporal variability of the radar signal, with high optical radiance values characteristic of a sandy desert. 'Dunes' and 'desert' behave differently from the SAR response of each of them.

Land Cover Classification
Land cover classification was performed by labelling land cover maps using the aforementioned five classes. We investigated the temporal continuity of land classes in subsequent LULC maps, with special focus on the transitions: (i) 'urban' → 'field'; (ii) 'field' → 'desert'; and (iii) 'dunes' → 'crop'.
The training dataset retrieved from the labelled classes using Google Earth, which was split into 11 sets, of which 10 were used to train as many classifiers, and 1 used as ground truth dataset for final validation. Subsequently, each classifier was applied to the datasets referring to the various years studied, thus returning 10 intermediate classifications per dataset. Unique final classification maps were obtained by feeding the intermediate results to a majority voting condition. The majority voting strategy assigns a pixel to class 'X' if it was classified as 'X' in 6 or more classification maps obtained for the same period.

Validation Approaches
The proposed method is validated by evaluating the goodness of the resulting individual classification maps against the ground truth datasets by calculating the standard metrics typically employed in classification purposes, such as overall accuracy and Cohen's Kappa index [44], as detailed in Section 3.2.4.1, as well as qualitatively comparing our resulting LULC maps with state-of-the-art available datasets, as detailed in Section 3.2.4.2.

Validation against Ground Truth Datasets
The ground truth datasets were selected from the high-resolution (<3 m) historical imagery copyrighted by DigitalGlobe and available in Google Earth as of 2016. Different point sets were selected for the time periods 2004, 2010, and 2015. No data referring to 1998 was available when this work was performed. We computed the accuracy of the final classification maps using the overall accuracy and Cohen's Kappa (K) index, as performed by Congalton [45] and Mather [46], in order to keep into account agreement occurring by chance.

Comparison with State-of-the-Art Datasets
We employed independent state-of-the-art datasets to compare our results, namely: (i) the Global Urban Footprint (GUF) [47,48], processed by the German Aerospace Centre (DLR); (ii) the Prototype of Sentinel-2 Land Cover Map of Africa at 20 m [49] released by the European Space Agency Climate Change Initiative (ESA-CCI); (iii) Global Human Settlements Layer using Landsat data for 2014 [50] (GHSL-L8), created by the Joint Research Center (JRC) and; (iv) Global Human Settlements Layer using Sentinel-1 data for 2016 [51] (GHSL-S1) with a 19 m spatial resolution [52], also created by the JRC. Note that the GUF, ESA-CCI and GSHL-S1 have a higher resolution than our results (30 m) and none of them totally matches the observation time of the data analysed here.
Due to the heterogeneity of the different datasets (data sensor employed, spatial resolution, dates), in order to compare these, a harmonisation process is needed. All datasets were resampled with 30 m resolution and re-projected in the same projection system of our dataset.
Finally, we compared the urban class across the different datasets: GUF, GHSL-S1 and GHSL-L8, ESA-CCI and our classification maps for 2010 and 2015.

Step 3: Multi-Temporal Evolution Analysis
After obtaining the five-class LULC maps for the four time periods, we identified pixels that remained stable and those that underwent changes from one land-use class to another between 1998 and 2015. Non-realistic changes were also detected and labelled as miss-classified pixels, based on our "multi-temporal consistency rules", defined as follows: (1) once a pixel is urban, it will stay urban in the following periods; (2) water cannot turn into urban, field or desert; (3) fields can only become urban, and not desert or water. In our study area the water class mainly corresponds to the Nile River, and was stable during the considered time period and; (4) desert can change to urban land or crop fields, as well as to dunes. In the latter case, it implies that dunes have been migrating into rocky deserts or across desert pavement surfaces.
In the Middle Egypt region, rule 2 is a simplification as there has effectively been a change in the size of the islands in the Nile channel which may have impacted the area classified as water or agriculture. However, as this does not have an impact on the urbanised areas, being the focus of this work, we believe that this simplification is justified.
With the temporal evolution maps for the urban and field classes, together with the population data from Central Agency for Public Mobilization and Statistics of Egypt (CAPMAS) [53], we have computed the trends of urban population density and field area per person. Since these data were not available for the entire study area, this analysis was limited to El-Minya governorate which is centrally located in the study area, as shown in Figure 4.   The accuracy of the land use and land cover maps is expressed in the confusion matrix, see Table 4. Overall, the agreement between the constructed LULC maps and ground truth data is very high as is illustrated by a Kappa index above 0.98 for the different time periods, which is considered an almost perfect agreement [54]. Note that for the period of 1998, historical data in Google Earth over our study area was not available.    The temporal evolution of the different land use classes is extremely interlinked (field → urban, field → desert/dune, desert → dune, desert → field, desert → urban, dune → field, dune → urban), since an increase in one class is directly translated into the reduction of others. This is the case, for example, for the urban increase on the first periods analysed (1998-2010), when most of the newly built-up area was located on the edges of the cities and villages (see Figure 6), and the increase in fields was due to the "land reclamation" phenomena mainly linked to the reduction of desert classes by introducing crop fields with irrigation systems inside the arid environment (see Figure 7).  The urban extent and its evolution per district in the El-Minya governorate are summarised in Table 6. Growth rates vary between 0.01 and 0.61 km 2 /yr. Over the entire period, urban areas have increased from 33 to 888% by 2015 compared to 1998. The lowest percentage increase is for El-Minya city itself, whereas a more prominent increase is observed for the new cities in the desert (New Minya and unorganized El-Minya district). Table 6. Urban spatial extent (km 2 ) within the different analysed districts in Middle Egypt for four time periods, and urban growth rates (% and km 2 /yr) for the period 1998-2015 (total size of the study area equals 6170 km 2 ).  Figure 8 shows the evolution in population density for the built-up environment in El-Minya governorate. These values were calculated by dividing the population census data by the urban spatial extent as shown in Table 6. Overall, there is a decrease from 46 thousand people/km 2 urban area in 1998 to 35 thousand people/km 2 urban area in 2015. Figure 8. also shows how with an increasing population but a decrease in the class 'fields', the average size of agricultural land available per inhabitant is decreasing from 4.7 ha/person in 1998 to 3.6 in 2015. Finally, we compared the urban features detected using our approach with the state-of-the-art urban global datasets, such as GUF and GHSL, as well as the urban class from the ESA-CCI Land Cover Prototype over Africa. Additionally, we computed the urban extent detected on El-Minya governorate, in order to understand whether it is possible or not to use these datasets together for urban expansion analysis. Details of the urban extent measured per dataset are shown in Table 7. Additionally, we put together all the datasets shown in Table 7. to visually compare the urban extent over two smaller areas contained within our study area, one over El-Minya City and surroundings ( Figure 9) and the second one over Mallawi and surroundings (Figure 10), from where it was possible to obtain additional information regarding old urban extent and modern cemeteries from Willems et al. [55]. The main differences are highlighted within colour ellipses in these figures and described in the discussion section. It is important to note the nonrealistic changes in the temporal evolution of detected urban extent for different datasets as seen in Table 7. and Figure 9. These inconsistencies reveal the limitations of performing a blind analysis using data obtained from different producers with their particular data sources and methods.   [55]. In colour circles different areas highlighting different classification performances. The blue rectangle highlights our results. Table 7 and Figures 9 and 10 compare the results of urban expansion obtained in this study with the data fusion approach to urban global datasets, as well as with a high-resolution visual Google Earth image. A direct and accurate comparison, however, is complicated as the timing of the imagery is different, such that discrepancies between the various urban maps are not solely related to differences in methodology but partly also to real changes in urban extent between the periods covered by each map. Nevertheless, some differences emerge from this comparison. The Global Urban Footprint dataset [47,48] does not detect flat urban areas, such as the Minya airport tracks (see Figure 9-yellow oval), as the method employed is optimised for the detection of man-made building structures with a vertical component. Both Sentinel products (ESA-CCI and GHSL-S1) are also having difficulties in mapping the spatial extent of the airport correctly. Only the data fusion approach (this study) and the GHSL-L8 dataset correctly identifies these flat urbanised areas. On the other hand, narrower strips of urban land, such as road surfaces, are in some cases mapped by the GUF and ESA-CCL dataset but not by the data fusion approach (Figure 10 -blue oval). This could be related to the resolution of the imagery used to create these maps. Indeed, the GUF and ESA-CCL maps are using 12 and 20 m resolution data, respectively, whilst in this study, we used 30-m resolution data. A 30-m resolution image could be too low to identify narrow roads accurately. The newly constructed Nile bridge and access road (Figure 10-red oval) is identified relatively well on most images, except the GSHL-S1 dataset. Note that road construction started in 2008-2009 and that the bridge and road were only finished in 2013-2014. Hence, the differences in mapping quality between the various datasets likely correspond to the building progress. The GHSL-L8 dataset performs very poorly when it comes to detecting the new urban expansions of New Minya (red oval in Figure 9). The GUF and GSHLS1 data seem to overestimate the urban density of this new settlement and suggest a similar density as the ancient city center. However, this is not observed, i.e.,the new town has more open spaces and parts are not yet developed. Here, the data fusion approach correctly predicts the expansion of the new town between 2010 and 2014 but at a lower density than the El-Minya city center itself. Furthermore, the ESA-CCI classifies many city centers as bare soil, thus showing a lower built-up density (Figure 9 El-Minya city and Figure 10 Mallawi city center). We attribute this to the fact that many buildings have similar spectral properties as bare soil: many buildings have been constructed from soil material (mud bricks) and deposition of dust on rooftops and against buildings walls in this extremely dry environment may also contribute to this effect. Finally, some large cemeteries can be found along the Eastern Desert margin (Figure 10-yellow). These cemeteries are mapped as built-up on the GUF map but could be erroneously understood as urban.

Quality of the Data Fusion Approach Compared to Single Platform Approaches
Overall, the detailed comparisons of the various urban datasets show that in general, the data fusion approach outperforms the quality of the global urban datasets that are constructed from a single data platform, whether it is optical or radar data. In particular, the data fusion of radar and optical imagery is able to correctly capture the extent of urban areas s.s., i.e., buildings. Errors in the classification of roads are probably related to the lower resolution. Our results also show that comparing data on the urban extent mapped using different approaches is not an accurate method to map and quantify rates of urban expansion through time. Table 7 could incorrectly suggest urban expansion took place between 2010 and 2014 (from 109 to 207 km 2 ), followed by a significant reduction in urban expansion in 2015 (again 148 km 2 ). In the next section, we will, therefore, discuss changes in urban and agricultural land areas using a single methodology, i.e., the data fusion approach applied to 1998, 2004, 2010 and 2015 imagery.

Urban and Agricultural Land Dynamics
Urbanisation follows a trend that accelerated during the last period analysed (see Table 6). Whereas urbanisation rates were already high in the first two periods, 2.5 and 1.8 km 2 /year for the periods 1998-2004 and 2004-2010, respectively, it is mainly in the period 2010-2015 that new urban land was created at a rate of 6.1 km 2 /year. Urban expansion is mainly linked to a reduction in arable land ("fields" class) within the floodplain (compare Figures 6 and 7). This type of urban expansion typically involves a diffusion process whereby cities and villages grow along their edges. This is illustrated in the area around the city of Mallawi in the El-Minya governorate ( Figure 11). No less than 60% of the urban growth near Mallawi between 1998 and 2015 took place in the last five years. It involved mainly new building blocks to the north and south of the city of Mallawi. In more recent time periods, however, urbanisation is no longer restricted to the floodplain. Table 5 shows how between 2010 and 2015, a significant part of the urban expansion took place in the desert. This is strongly linked to government-led urban planning whereby new cities are built in desert areas, such as New Minya, located within the eastern desert (see Figure 9, red circle). Around Caïro, desert cities such as 6th of October City, 10th of Ramadan City and New Cairo, amongst others, have been developed since the 1970s and 1980s [56]. This process has also gradually impacted Middle and Upper Egypt from the 1990s and 2000s onwards, with new towns being planned near Minya, Assiut and Qena, amongst others [57].
Urban expansion goes hand in hand with a reduction in land used for agricultural activities in the Nile floodplain: 53 km 2 of arable land has been irreversibly converted into urban land between 1998 and 2015 in the entire study area (Table 5). However, this does not imply that agricultural land is decreasing. In the same period, no less than 260 km 2 of barren desert or former dunes has been reclaimed for agricultural purposes. This includes areas in the western and eastern desert where large-scale irrigation schemes have been set up, as well as dune-levelling in the interaction area [14]. Figure 7 shows that the levelling of dunes in the interaction area already started in the late 1990s, whereas land reclamation in the desert is a more recent phenomenon. It is mainly the latter process which is responsible for the increasing rate at which new agricultural land is created. These two phenomena are also illustrated in Figures 12 and 13 respectively.
The aim of the Egyptian government with the new-town policy has been to release pressure on the fertile agricultural land with a rising population and to improve the living conditions in highly populated ancient urban centers [56]. However, whilst our data indeed shows that urban expansion in the desert takes place at increasing rates, it is clear that urban expansion in the Nile Valley is far more important in absolute numbers. Contrary to the planned desert cities, this type of urban expansion is not planned and often informal [58]. Between 1998 and 2015, approximately 1500 ha of fertile land has still been converted into urban land each year. Despite the overall increase in arable land through land reclamation in the desert, it cannot keep pace with increasing population numbers. In the El-Minya governorate, the average extent of agricultural land per inhabitant has decreased from 4.7 to 3.7 ha/person, i.e., a reduction of 21% over 17 years. Hence, we can conclude that at least for Middle Egypt, the new town policy did not halt the reduction in the availability of fertile soils for crop cultivation. Newly reclaimed areas are also more vulnerable as water has to be pumped from the Nile Valley or from groundwater resources, which is not only very costly but it may also lead to failures in water supply in the case of malfunctioning systems, thus questioning the sustainability of this type of agriculture (e.g., [59,60]). On the other hand, urban population density decreased from 46,000 to 35,400 inhabitants/km 2 built-up area in the same period ( Figure 8). Despite this reduction, these average urban population densities of more than 350 people per hectare remain very high to global standards [61]. Furthermore, population densities in the new desert towns are often lower as more open spaces are provided and not all towns are yet successful. For New Minya, the population is currently estimated at 45,000, whilst it should reach 638,000 in 2050 [62].This means that population densities in the towns and villages in the Nile Valley so far remain densely populated, similar to many residential areas in Greater Cairo [63]. Furthermore, contrary to the planned desert cities, this type of urban expansion is not planned and often informal [58]. Hence, living conditions remain in general very poor.

Conclusions
In this study, we applied a data-fusion approach to map land-use dynamics in Middle Egypt for the period 1998-2015. Our data-fusion approach proved to be an effective method to discriminate urban areas where other approaches fail, especially when urban areas are built with mud bricks or constructed within desert areas. The combination of optical imagery (Landsat 5 TM, 7 ETM+ and 8 OLITIRS) with radar imagery (ERS-2 SAR, Envisat ASAR and Sentinel-1 IW) resulted in four multi-temporal land cover maps at a resolution of 30 m. Comparison with high-resolution optical images (Google Earth) showed very high accuracy, stressing the efficacy of the data-fusion approach in mapping land cover dynamics. Furthermore, comparison with globally available urban datasets shows that the data fusion approach is performing better than the single platform-based approaches when it comes to identifying new settlements in the desert or to discriminate urban buildings from bare arable land. Our study also suggests that multi-temporal studies on the land cover that combine results from various producers obtained with different methodologies could lead to erroneous nonrealistic interpretations.
Our data show a rapidly increasing trend in urbanisation in Middle Egypt (65 km 2 ), in particular, along the margins of existing towns and villages in the Nile Valley. As a result, a continuous decrease of fertile Nile floodplain soils for agricultural practices can be seen. This loss is compensated by land reclamation processes whereby former dunes in the Nile Valley are levelled, and by irrigation practices in the desert (over 200 km 2 ). Additionally, the results show that both the urban population density and the amount of agricultural land per person have decreased by more than 20% since 1998. Finally, we have demonstrated that the proposed data-fusion approach is a viable tool to continuously monitor future land cover changes and can be used to update the management and planning of urban areas.