Integrated Mapping of Spatial Urban Dynamics—A European-Chinese Exploration. Part 1—Methodology for Automatic Land Cover Classification Tailored towards Spatial Allocation of Ecosystem Services Features

Urbanisation processes inherently influence land cover (LC) and have dramatic impacts on the amount, distribution and quality of vegetation cover. The latter are the source of ecosystem services (ES) on which humans depend. However, the temporal and thematical dimensions are not documented in a comparable manner across Europe and China. Three cities in China and three cities in Europe were selected as case study areas to gain a picture of spatial urban dynamics at intercontinental scale. First, we analysed available global and continental thematic LC products as a data pool for sample selection and referencing our own mapping model. With the help of the Google Earth Engine (GEE) platform and earth observation data, an automatic LC mapping method tailored for more detailed ES features was proposed. To do so, differentiated LC categories were quantified. In order to obtain a balance between efficiency and high classification accuracy, we developed an optimal classification model by evaluating the importance of a large number of spectral, texture-based indices and topographical information. The overall classification accuracies range between 73% and 95% for different time slots and cities. To capture ES related LC categories in great detail, deciduous and coniferous forests, cropland, grassland and bare land were effectively identified. To understand inner urban options for potential new ES, dense and dispersed built-up areas were differentiated with good results. In addition, this study focuses on the differences in the characteristics of urban expansion witnessed in China and Europe. Our results reveal that urbanisation has been more intense in the three Chinese cities than in the three European cities, with an 84% increase in the entire built-up area over the last two decades. However, our results also show the results of China’s ecological restoration policies, with a total of 963 km2 of new green and blue LC created in the last two decades. We proved that our automatic mapping can be effectively applied to future studies, and the monitoring results will be useful for consecutive ES analyses aimed at achieving more environmentally friendly cities.


Impact of Urbanisation Processes on Ecosystem Services
Increased urbanisation and the expansion of urban areas cause radical changes to and the eradication and fragmentation of ecosystems and green infrastructure-particularly in mid-1980s), the Urban Atlas (initiated in 2006, with updates in 2012 and 2018) and the Copernicus services. The latter launched its Land Monitoring Service in 2012, and all these platforms provide data on the European environment and various detailed information to users free of charge. So, what we know and what we have at hand are a multitude of publicly accessible products for environmental research containing defined sets of LC categories at defined scales over specific periods of time or time slots. Such readily available spatial information is especially of value for intra-and international comparisons in spatial planning and research tackling urban and regional environmental issues. What is lacking is a standardised set of LC categories across all these systems. They hardly cover similar time spans, and most neither differentiate between dense and dispersed urban fabric, nor between different types of vegetation cover (for details see Table 1). This deficiency may be caused by the specifications under which they were developed, or by the specific scale they refer to. For this reason, these thematic products can be used for general comparisons and as a profound basis for further in-depth analysis.

Further RS Requirements for Capturing Key Elements of Urban LC Categories across Continents
To assess LC in cities around the globe, all spatial information must be provided at the same spatial and temporal scale and with the same detailed categories. If the designs of the above-mentioned publicly available RS catalogues do not meet the requirements for research on LC related to urban ES, then those catalogues can be taken as basic input for training and validation to establish an ambitious classification scheme. Enhanced accuracy would make RS more suitable for monitoring biodiversity loss and ecosystem dynamics, for example, and for other applications that depend on baseline and changing LC [11]. To improve the accuracy of products derived from shared satellite observations of urban areas, a spatial analysis needs to be undertaken developing elaborated methods and tools. When mapping LC elements to determine the extent and condition of natural capital and its biophysical expression over time, then the Google Earth Engine (GEE) serves as an effective analytical platform. GEE is a sophisticated and scalable mapping tool, with which we are enabled to determine the time courses of urban growth rates and quantify the other more natural LC types that are being lost through advancing urbanisation.
The aim of this study is to design an automatic classification method at the intercontinental scale to understand the different urbanisation processes in Europe and China, firstly by exploiting global urban datasets and secondly by a distinct mapping of grey, green and blue features as a basis for further ES analysis. The methodological innovation of our contribution includes an automatic variable selection for a more efficient creation and application of the classification model. To realise this, a well-developed workflow for an automatic samples collection based on the above-mentioned multi-source ground truth and raster-based land cover products is presented. It is a prerequisite to collect spatial information for both continents to obtain training and validation data. For reasons of consistent results, we aim to capture key categories by elaborating a novel classification approach. To make our approach easily reproducible for other cities, we develop this processing workflow based on GEE. This enables us to take into account the existing thematic products for LC, and extract valuable information for our discerning, globally applicable mapping model presented here in Part 1. Because we will perform further ES analyses in urban areas (Part 2), we undertake the effort to differentiate various types of vegetation with distinctions between grassland and cropland, and deciduous and coniferous forests. Up to present, there are no existing products for both China and Europe over decades showing such levels of detail. For this reason, our research tackles this central question: Can we capture a detailed LC mapping that contains highly differentiated categories for selected cities in Asia and Europe by developing a streamlined processing in GEE? This encompasses two secondary questions: Can we make distinctions between deciduous and coniferous forests and between dense and dispersed urban areas? To what extent can we map cropland, grassland, and bare land for each city?  All data from the Landsat series are available at the GEE platform. In this research, we used various datasets from Landsat sensors including the Landsat 5 Thematic Mapper (TM) for 2000 and 2010, and the Landsat 8 Operational Land Imager (OLI)/Thermal Infrared Sensor (TIRS) for 2020. The Landsat series offer the most long-term satellite images at the same spatial scale of 30 m ground resolution and make it possible to extract the same spectral resolution over decades. Their global repetition rate of 16 days means that data is acquired throughout the entire year which makes it easier to find images with no or minimum cloud cover.

Reference Data for Sample Points
The selection of reference data was driven by two main factors. On the one hand, spatial coverage was important, as the different urban areas are located across China and Europe. Thus, only global datasets were considered at the beginning. This selection proved to be too narrow, so we expanded our analysis to include European and Chinese datasets. They feature different LC classes in greater detail. The second factor we considered pertains to the temporal coverage, which represents continuity. All the selected thematic mapping products feature different intervals, but they are available for at least three points in time and thereby map the changes in LC over decades. The data used to generate the sample points in this study were sourced from the following ground-truth-based and remotely sensed products.

Datasets for Sample Point Extraction Based on Ground Truth CORINE Land Cover
The Coordination of Information on the Environment land cover, also known as the CORINE Land Cover or CLC, is a European land use and land cover product spanning the years 1990 to currently 2018 in six-to 10-year increments. The classification system is based on 13 main classes and contains 44 classes in total. Its contents are generated through manual and automatic classification of satellite images as well as field surveys. The product is available as a raster dataset with a resolution of 30 m. The overall accuracy is specified as 87% with different classes ranging from 70% to 95%, but several studies have shown this to be inaccurate and, depending on the land cover and land use class, the rate of accuracy ranges from 82.8% to 97.6% for different years [12][13][14]. However, given its good temporal coverage, the product was deemed to be a suitable reference for this study.

Urban Atlas
The Urban Atlas is another European mixed land cover and land use product. It is available for the years 2006, 2012 and 2018 with 21 and 27 classes for the last two versions. The coverage was also increased from 305 to 695 cities in 2012 and then to 788 cities in 2018 [15,16]. Spatial accuracy is presented at 5 m [16]. The accuracy is reported to be above 80% for all classes, whereby the 2012 product has different accuracies for urban and rural areas (98.4% and 96.9% respectively) [15,17]. Agreement between the GlobeLand30 product and the Urban Atlas is reported to be above 85% [18]. The urban setting of this product along with its temporal continuity and high spatial accuracy meant that it was also considered to be a good reference.
The Land Use and Land Cover Area Frame Survey (LUCAS) Sample Points LUCAS is a land cover and land use survey for Europe conducted by the European statistical office [19] on a three-year basis. The dataset includes some 240,000 individual points in the latest version from 2018, all measured in situ. Thus, the dataset not only provides information about the predominant LC, but also about land use. It features 16 thematic classes that are determined either through the interpretation of RS data or by on-site evaluation, along with photo documentation taken in the field. The field work guidelines allow for consistent classification by minimising the influence of the respective surveyor on site [19]. Because of their consistency, the sample points are available over time and match well with the classes used for the mapping. For this reason, they are used to validate the LC mapping of the European cities studied here.

Datasets for Sample Point Extraction Based on RS Products GlobeLand30
GlobeLand30 is a 30-m resolution global LC data product for 2000, 2010 and 2020 that was developed by the Chinese Ministry of Natural Resources [8]. It comprises ten thematic classes including cropland, forest, grassland, shrubland, wetland, water bodies, tundra, artificial surface, bare land, perennial snow, and ice. The multispectral images used to develop and update the LC classification are derived from the sensors TM5, ETM+, OLI of Landsat (USA), as well as HJ-1 and GF-1 (China). Image selection followed the principle that the image is cloudless (or has less cloud), and the multispectral images refer to the growing season of vegetation within ± 2 years of the baseline year when the data was produced or updated. For areas that are difficult to obtain, the image acquisition time could be adapted to ensure the integrity of the global coverage. The total accuracy of GlobeLand30 was 83.5% in 2010 and 85.7% in 2020, their Kappa coefficients were 0.78 and 0.82, respectively. The results were validated by over 230,000 points from the whole dataset using the landscape index sampling model [20]. It can provide us with the latest information on LC in 2020 for validating the current mapping work.

Forest Maps
In order to further distinguish between different forest types (deciduous and coniferous) in our mapping model, the 2010 China Forest Type Classification and Mapping Products were used to select sample points for the three Chinese cities [21]. The 30 m resolution forest map for China was established using Landsat TM, by the Moderate Resolution Imaging Spectrometer (MODIS) time-series at 250 m ground resolution, by the Enhanced Vegetation Index (EVI) and other auxiliary datasets. The classification of forested areas is comprised of six types: coniferous/deciduous broadleaf, coniferous/deciduous coniferous, mixed forests, and bamboos. To assess the accuracy of the forest/non-forest classification, 2195 test sample units were used independently of the training sample. The quality indicates a producer's accuracy (PA) of 92.0%, user's accuracy (UA) of 95.7% and an overall accuracy (OA) of 72.7%.

The Chinese Academy of Sciences (CAS) LC Dataset
To classify the Chinese cities more efficiently with regard to the distinction of bare land, waterbodies and grassland, we used the 30 m resolution land use and land cover products derived from Resource and Environment Science Data Center, CAS (www.resdc.cn. Note: This source of land use and land cover products is not entirely public. The products with 1-km resolution are open access after free registration, while the 30-m products should be requested and purchased on the website. All information on the website is in Chinese only). Here, we extracted the maps for 2000, 2005, 2010 and 2015 that are based on multisources RS images from Landsat 5 TM, Landsat 7 ETM, and Landsat 8 OLI. This sequence of datasets was validated by sample points through field investigation with Kappa coefficients of more than 0.8 and accuracy of bare land, waterbodies and grassland of more than 85%, 90%, and 90% respectively [22].
The Peking University UrbanScape Essential Dataset (PKU-USED) for Beijing To distinguish between dense and dispersed built-up areas, we made use of the USED maps to create sample points to validate our current mapping for Beijing. We extracted the dataset using the functional zones mapping tools [23,24] which are based on the Chinese satellites ZY-3, GF-6, and GEE imagery repository [25]. The urban functional zones are divided into twelve categories and have a spatial resolution of 2.4 m. The Thematic Mapping Product GAIA As a high-resolution global artificial impervious area product, GAIA acquired Landsat Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operational Land Imager (OLI) data in different years to map the artificial impervious areas. It features a 30 m resolution and contains more than 30 years of records (from 1985 to 2018) with a mean overall accuracy of above 90% [9]. In this study, GAIA was used to help generate samples for built-up areas in the six cities in 2000, 2010, and 2020. The GAIA dataset can be freely downloaded from http://data.ess.tsinghua.edu.cn and is also available on GEE (Tsinghua/FROM-GLC/GAIA/v10).

Input Data for Samples Generation
We extracted all the thematic products at the global and continental scales to generate samples. This sample selection is essential for training the mapping model and for validation purposes. These data served as our input data source pool, out of which we extracted the best fit for our own mapping model. The details of the LC products we used are listed above, the methodology of how we used the samples selection is explained in Section 3.3.2.

Study Area
The study area is composed of three cities in Europe (Paris Region, Aarhus and Velika Gorica) and three in China (Beijing, Shanghai, Ningbo). Figure 1 illustrates their locations. These urban areas vary in their spatial extent, population number and density, and the velocity of their urbanisation processes. Their diversity is valuable for testing existing mapping products at an intercontinental scale and demonstrating the applicability of our mapping model. By choosing a diverse variety of cities, we can ensure the transferability of our GEE approach. The common interest among the six cities is their commitment to foster green infrastructure and their related ES in the respective urban areas. Paris Region is the capital region of France with a land area of 12,012 km 2 , home to 12 million inhabitants or 18.2% of the population of France. 21% of the territory is urbanised (16% totally impermeable), 23% covered by forests and 47% of the territory is composed of cultivated, mainly intensive open field crops. Biodiversity has greatly depleted over the last decade, which is especially visible in agricultural areas and in urban parks and gardens. Green space per inhabitant is very scarce in the four dé partements within the inner ring of the Paris Region-it varies between 2 and 10 m 2 depending on administrative unit. The region experiences substantial land take with more than 900 hectares of rural area Paris Region is the capital region of France with a land area of 12,012 km 2 , home to 12 million inhabitants or 18.2% of the population of France. 21% of the territory is urbanised (16% totally impermeable), 23% covered by forests and 47% of the territory is composed of cultivated, mainly intensive open field crops. Biodiversity has greatly depleted over the last decade, which is especially visible in agricultural areas and in urban parks and gardens. Green space per inhabitant is very scarce in the four départements within the inner ring of the Paris Region-it varies between 2 and 10 m 2 depending on administrative unit. The region experiences substantial land take with more than 900 hectares of rural area consumed each year by urbanisation.
City of Aarhus is the second largest city in Denmark with 273,077 inhabitants with a land area of 468 km 2 and a total municipal population of 341,000. The urban area of Aarhus is 91 km 2 and its population density is around 3000 inhabitants per km 2 , compared to an average of 272 people per km 2 for the entire municipality. Over 90% of the inhabitants have access to green space within 500 m and the authorities aim to sustain or increase green area per inhabitant through green/blue structure planning despite densification. The city anticipates densification with an additional 75,000 inhabitants by 2030.
City of Velika Gorica is the sixth largest city in Croatia with 63,517 inhabitants and is located in close vicinity to the City of Zagreb, the capital of Croatia. The municipality covers 328 km 2 and has 193 inhabitants per km 2 . In contrast, the urban area covers only around 13 km 2 (ca. 4%) and experienced intensive construction during the 1970s and 1980s. Here, 32,000 inhabitants live with a population density of around 2500 people per km 2 . In the near vicinity, lowland, humid wooded areas host the highest woodland biodiversity nationally.
Beijing is the capital city of China with a land area of 16,410 km 2 . Beijing has gone through a rapid urbanisation process in the past four decades. Beijing's permanent population increased from 13.8 million in 2000 to 21.7 million in 2016, while the built-up area expanded from 488 km 2 to 1400 km 2 in the same period. Urban green space accounts for about 821 km 2 , corresponding to 21 m 2 per inhabitant. By the end of 2015 there were 339 registered municipal parks.
Shanghai is the largest city in China. It is located on the banks of the Yangtze River Delta in Eastern China and covers a total area of 6340 km 2 . Shanghai has the world's largest port and a population of 24.2 million, including a transient population of more than two million. The city owns various types of natural wetlands which account for 23.5% of its total area, but these have suffered from a loss of biodiversity and vegetation decline. As the administration aims to transform the city into an "ecological city", construction and protection of green spaces, forests, and wetlands have been launched across the city.
Ningbo is an important port city in China. Located 220 km south of Shanghai, it belongs to the densely populated Yangtze River Delta region. Having undergone heavy urbanisation processes including land reclamation that involved dams, bridges, and port constructions, it has expanded its surface area. Presently, it covers an area of 9670 km 2 and is home to about 5.7 million inhabitants, representing a population density of 584 per km 2 . Although there is an emphasis on environmental sustainability with high investment in green and blue infrastructures, there is a lack of design and management at the urban level.
To gain an understanding of the historical urbanisation developments of the study areas, we illustrate the urban LC product GAIA. In Appendix A, Figure A1a-f illustrates the dynamic urban expansion of the six cities in their various expressions. All the information, originally mapped at a 30 m ground resolution, is aggregated to a spatial resolution of 100 m to ease visualisation.
Regarding the European cities, it is striking that there is only some development in Aarhus ( Figure A1a), mainly in the western and northern parts. Radial urban dynamics following the lines of transport infrastructure are discernible in the Paris Region ( Figure A1b) that reach far into the surrounding rural area. In Velika Gorica ( Figure A1c), the urban LC is expanding the urban core area, extending further towards Zagreb Airport in the north, and along the main transport lines.
The monitoring result for the Chinese cities shows quite a different picture that is much more intense in terms of velocity of change and coverage of land take: Beijing is expanding intensively from the core area (baseline 1985). Due to its geotopography with mountainous ranges from the west to the north and northeast, further LC expansion shows a dispersal pattern severely occupying the plains towards the southwest and southern region ( Figure A1d). The most extreme urbanisation process is witnessed in Shanghai ( Figure A1e), with the highest density development patterns now reaching the distant Chongming Island. Ningbo ( Figure A1f) has extended its urbanised area by growing like a network, changing the urban LC by forging connections with neighbouring cities.

Defining Essential LC Mapping Categories to Achieve Research Objectives
We developed a profound understanding of the existing thematic mapping products that visualise the spatial dynamics of the entire study area. When analysing the definitions of their quantified categories, we discovered gaps in the knowledge about intra-urban development, i.e., a differentiation between dense and dispersed urban areas in the global datasets. Such an inner-urban distinction would make it easier to monitor densification processes and reveal needs for green spaces inside cities. Above all, we discovered a lack of distinctions between different types of vegetation cover, which is critical for measuring the loss and gain of LC features towards spatial allocation of ES over time. For this reason, we aimed to extract distinct categories such as cropland, grassland and bare land, and differentiated between deciduous and coniferous forests. Table 1 provides an overview of the LC products we used and shows how it was necessary to produce the essential LC categories in this study.

GEE Mapping Procedure
In order to develop a consistent, automated workflow to map distinct LC features as the basis for further ES analysis, our workflow was developed using full temporal Landsat images, Shuttle Radar Topography Mission (SRTM) data, multiple open-source LC products, and the GEE platform. Our workflow is divided into three parts: (1) preprocessing of satellite images, (2) reference samples selection, (3) mapping of grey, green and blue infrastructures ( Figure 2).

Preprocessing of Satellite Images
As GEE provides free access to all Landsat surface reflectance images, we used the Landsat 5 TM surface reflectance Tier 2 product from GEE for 2000 and 2010 LC mapping. For 2020 we opted for the Landsat 8 OLI/TIRS surface reflectance Tier 2 product. All the cloud and cloud shadow pixels were removed using C implementation of the Function of Mask (CF Mask) algorithm [26].
To test which variable was important for mapping urban LC categories, we prepared a total of 32 variables. The specific variables are as follows: 1.
Three texture variables from the Grey-Level Co-occurrence Matrix (GLCM) measurement 3.
Annual Four seasonal indices derived from multi-temporal spectral indices 5.

Reference Samples Selection
As a novelty in the reference samples selection, we elaborated a sophisticated workflow for an automatic samples collection based on multi-source ground truth and rasterbased LC products. As shown in Section 2.2, both ground-truth samples and remotely sensed datasets were used to generate reference samples. Firstly, the vector datasets, such as LUCAS sample points, were filtered and matched for each category. In addition, readymade raster products such as GlobeLand30 were exploited. We collected sample points in three steps: (1) generating random points; (2) extracting values from raster products; (3) rebuilding classification code for each category. We depict the training and validation sample points in Table 3.

Assessment of Variable Importance
Ideally, more variables would lead to a better representation of the features of LC and improve the mapping accuracy. However, previous studies have shown that the accuracy and efficiency of the classification is affected by the high correlation of the data, the noise in the data collection and correction process, and the increase in computational complexity [48,49].
As a measurement for node impurity, the value of GINI can indicate how often a random instance will be misclassified. Therefore, GINI can be used to evaluate the importance of each variable [50,51]. In the RF model, GINI is calculated over all trees as the averaged reduction of node impurity on one splitting variable. Since significant variables can substantially reduce the impurity of the sample, the higher GINI values represent greater importance of the variable. This index was used as a criterion to evaluate the importance of the variables [49]. Technically, the GINI index can be calculated in GEE directly.

Mapping of Grey, Green and Blue Infrastructures Random Forest Classifier
To date, Random Forest is the most widely used classifier for LC classifications [52][53][54][55]. Previous studies have shown that, compared with support vector machine (SVM), k-nearest neighbour (KNN), maximum likelihood classifier (MLC) and other models, RF has a higher classification accuracy, makes more effective use of high-dimensional data, and has higher model efficiency. Therefore, RF was used as the classifier for LC mapping based on the training samples from 3.2.2 and filtered variables. To benefit most of RF, we designed an automatic variables selection for a more efficient creation of a classification model. Based on the knowledge from previous studies, good results can be obtained using number of trees (Ntree) as 100, and higher values do not significantly increase the accuracy of the classification [56][57][58][59][60][61]. In addition, previous studies suggested that the RF method is not very sensitive to the number of features randomly chosen to split each node (Mty) [62]. Therefore, the Ntree was set as 100 and the default parameter of Mtry was utilised in our study.

Validation
In this study, confusion matrixes were calculated based on reference validation samples to evaluate the accuracy of the classification result [63]. Several evaluation indicators were extracted, including overall accuracy (OA), producer's accuracy (PA) and Kappa value. We used 70% of collected samples of each category as training samples, while the remaining samples (30%) were utilised for accuracy evaluation. The code to reproduce the validation in GEE is https://code.earthengine.google.com/888298b27c9b084154088bee968953d2.

Post-Classification of Dense and Dispersed Built-Up Areas
The dense built-up category refers to built-up areas with dense building stock, including informal settlements and asphalt roads. The dispersed built-up category represents impervious surface areas with ventilation space, open spaces for ventilation and green areas, which are typical of scattered residential areas, high-rise building areas, the urban fringe, and the rural hinterland [64].
To differentiate between the built-up urban areas on the basis of density, the built-up kernel density was generated applying a 17 × 17 pixels moving window according to the equation below. The size of the moving window follows the rule of the CORINE LC definition. Thus, the coverage of the built-up area was further reclassified into dense and dispersed built-up areas. The equation is as follows: where x ij represents the value of pixel i j. The value is 1 when it belongs to the built-up and 0 when it is the non-built-up LULC, n is the size of the moving window, which is set to 17 (with reference to CORINE land cover). ρ b signifies the dense built-up area, if ρ b is at least 50% or greater in moving window. If ρ b is less than 50%, then the value is assigned to dispersed built-up area. The code for this step in GEE is https://code.earthengine.google. com/18c309ec9665c9895850390fe18f37fd.

Feature Selection and Ranking
In Figure 3, we have summarised the results of scores to determine the importance of the variables for our RF mapping model. The results show that the seasonal index was one of the most significant variables for prediction, suggesting the value of temporal information. The seasonality enabled us to distinguish between cropland and grassland as well as between the two forest types. By doing so, we proved the seasonal index by Reschke and Hüttich [45] to highlight the differences among various types of vegetation. High-ranked variables also comprised TCP-based greenness, MNDWI, Red band, CIgreen index and seasonal index VI-MNDWI. However, several variables including four kinds of topographic variables and NIR, BT, wetness, SVVI and TDVI, were lower ranked. In hilly areas, for instance in Beijing and Ningbo, the topographic variables such as slope and DEM ranked highly. accuracy as the categorical variables increased. The results show that the initial variable input clearly improves the classification accuracy, but when the number of variables reaches a certain level, the classification accuracy can no longer improve. Instead, some classification accuracies actually decrease as the variables increase, such as for Aarhus, Beijing and Shanghai. In order to find a balance between the classification efficiency and accuracy, the study chose the peak accuracy as the node to determine the number of input variables for further classification.

Accuracy Assessment
Among all the LC products of the study sites in 2000, 2010 and 2020, the results produced high accuracies (OA ranges from 73% to 95%). By calculating the average producer's accuracy (PA) it was found that the highest classification accuracy was achieved for cropland, where PA could reach 92%, followed by the extraction accuracy for water bodies, which reached 87%, while the extraction accuracy for grassland and bare land were the worst with 64% and 73%, respectively (Table 4). Those two also showed the low- In addition, Figure 3 reflects the change that took place in the overall classification accuracy as the categorical variables increased. The results show that the initial variable input clearly improves the classification accuracy, but when the number of variables reaches a certain level, the classification accuracy can no longer improve. Instead, some classification accuracies actually decrease as the variables increase, such as for Aarhus, Beijing and Shanghai. In order to find a balance between the classification efficiency and accuracy, the study chose the peak accuracy as the node to determine the number of input variables for further classification.

Accuracy Assessment
Among all the LC products of the study sites in 2000, 2010 and 2020, the results produced high accuracies (OA ranges from 73% to 95%). By calculating the average producer's accuracy (PA) it was found that the highest classification accuracy was achieved for cropland, where PA could reach 92%, followed by the extraction accuracy for water bodies, which reached 87%, while the extraction accuracy for grassland and bare land were the worst with 64% and 73%, respectively (Table 4). Those two also showed the lowest classification accuracy. It is meaningful to note that the accuracy evaluation was conducted with the country as the scale. A detailed classification accuracy evaluation table and confusion matrix are available in the Appendix A.2, Tables A1-A7.

Quantitative Analysis of LC and Its Changes
In order to understand the dynamics of LC change in each city over the past two decades, the area of each type of LC was summarised and illustrated in Tables 5 and A2,  Tables A3-A7, and Figure 4a,b. Figure 5a-f shows the area of various LC types in different stages. These differences in China and Europe can be described as follows: (1) Built-up area ratio By calculating the percentage of dense and dispersed built-up areas to total built-up areas, we can assess the degree of built versus more natural LC in the six cities. Based on the results for 2020, this is how the three European cities ranked in terms of share: Paris Region (75.6%), Velika Gorica (69.2%) and Aarhus (66.9%). The three Chinese cities with the highest share of built-up areas were Shanghai (73.2%), followed by Ningbo (69.08%), and finally Beijing (65.2%).
(2) Spatial urbanisation processes We compared the space covered by urbanisation processes in European and Chinese cities in 2020, then estimated the rate, i.e., the ratio of built-up area to total city area, and came up with the following ranking: Shanghai (40.6%), Ningbo (20.8%), Aarhus (18.1%), Beijing (17.4%), Paris Region (15.6%), and Velika Gorica (5.2%). In addition, Chinese cities have experienced dramatic urban expansion over the past 20 years. Specifically, the builtup area of the three cities has grown from 4077.6 km 2 to 7508.7 km 2 . In the two decades between 2000 and 2020, the built-up space increased by 84%. In contrast, the built-up area of the three European cities has been rather stable over the last two decades, with the total area slightly increasing from 1780.5 km 2 in 2000 to 1986.3 km 2 in 2020, an increase of only 11.6%. (

3) Reduction of green and blue spaces in connection with urban expansion
In order to understand the characteristics of LC transformation in each city from 2000 to 2020, the study further designed an LC area transfer matrix as shown in Table 5. It shows that the rapid urbanisation process had a direct impact on the cities' green and blue LC. As a direct result of the urban expansion of the Paris Region, 646.1 km 2 of cropland, 50.5 km 2 of deciduous forest and 10.9 km 2 of grassland have been replaced by built-up areas. A similar situation occurred in the three Chinese cities, where the urbanisation process in Beijing and Shanghai resulted in the loss of 1950.2 km 2 and 1920.3 km 2 of cropland each. As coastal cities, Shanghai and Ningbo have respectively converted 25.8 km 2 and 100.5 km 2 of water bodies into built-up surfaces due to human land reclamation.
Remote Sens. 2021, 13, x FOR PEER REVIEW 17 of 35 ecological restoration driven by ecological restoration policies, or what we know as the turning green effect [65]. Specifically, in the last 20 years, 451.2 km 2 of built-up areas have been restored to green and blue LC in the three European cities and 963.0 km 2 were restored in the three Chinese cities.
(a)   (4) Changes in agriculture In the last two decades, cropland slightly increased in the three European cities from 7486.6 km 2 to 7576.0 km 2 , an increase of merely 1.2%, while the total cropland in the three Chinese cities shrank significantly from 15,097.2 km 2 in 2000 to 11,466.0 km 2 in 2020, a decrease of 24%. This result indicates that the intense urbanisation process in China has led to a great loss of agriculture.

(5) The ecological restoration effect in the selected European and Chinese cities
Although the dramatic urban expansion of the past two decades has caused damage to most ecosystem features, at the same time some areas showed the effects of significant ecological restoration driven by ecological restoration policies, or what we know as the turning green effect [65]. Specifically, in the last 20 years, 451.2 km 2 of built-up areas have been restored to green and blue LC in the three European cities and 963.0 km 2 were restored in the three Chinese cities.

Information on Urban LC in Global Thematic Products
This study aimed to analyse LC changes to understand urban growth patterns at an intercontinental scale, as well as to provide a basis for urban ES analyses. In this analysis, we successfully tracked urbanisation processes in a historical context by interpreting global datasets such as GAIA [9] and GlobeLand30 [8,20]. Each global dataset underlies different specifications and shows the variety in their individual categorical qualities. Each dataset follows specific goals such as mapping the artificial impervious land surface cover at a global scale, or to quantify the ten most important LC categories globally. What they have in common is that they provide insights into their specific LC categories over a long period of time and spatial precision due to their 30-m resolution.
These thematic mapping products are valuable for global analyses and help us to understand how and where LC changes take place, which LC patterns exist, and in which direction loss and gain of these categories have taken place [10]. They also reveal the speed, direction and extent at which this transformation occurs [66]. These LC maps were used for visualisation purposes to gain an initial understanding of the problems that ecosystems face as a consequence of urbanisation processes, and to capture how urban growth occurs in different places.
Harmonising all the information is a challenge, as the products cover different periods and feature different levels of categorical and spatial detail. Urbanisation has a great impact on ES that relate to LC, especially to differentiated vegetation cover [67]. At the European level, highly demanding thematic mapping products have been developed that not only provide long-term information on LC, but also serve as information platforms for land use categories and their respective changes [16]. Although they are high quality and cover a long period of time, their disadvantage is that they are restricted to the European continent [18].
Regarding the Chinese cities, several thematic mapping products provide precise information for analysing LC and correspondent ES changes at fine scale with annual time steps [22]. We showed some lacks that hinder comparability of LC quantifications in Chinese and European cities. This is because of the differences in their categorical definitions and usages. In order to do justice to our targeted research, that is to use the mapping product for further ES analysis, none of the existing products were sufficient in terms of spatial extent or thematic resolution.
Hence, we overcame this challenge and managed to profit from both the global information and the details from specific European and Chinese information layers which we exploited for our distinct classification system. Another major advantage of the thematic mapping products was that we extracted all the information, layers, and categories in order to collect detailed sample points which we then used to build the map model and validate the mapping results. We take into account the inherent inaccuracies of the reference products and are aware of the potential error propagation.

Integrated Mapping Model for Relevant Urban LC Categories Tailored towards ES Feature Allocation
By using existing LC products to integrate and extract classification sample points, we could solve the problem of their difficult acquisition that was depicted in previous studies [68]). As a novelty, we propose an algorithm for automatic sample points' acquisition. This method of classification feature screening can now be applied to other studies. By evaluating and screening multiple feature variables with the help of GINI coefficients on model evaluation, we obtained a balance of good classification accuracy as well as a high classification efficiency. We achieved this discerning classification scheme at the intercontinental level. However, we understand the restriction of our approach. As dominant artificial and natural features and thus LC may vary over different countries and cultures, they might not be represented optimally by our categories.
Our methodology ideally exploits all the spatial LC information, and therefore we were able to design a well-developed mapping procedure. With our elaborated GEE workflow, no additional infrastructure is needed. Our processing is easily reproducible in other cities contrasting other approaches [69]. The performance revealed that our design is scalable while it simultaneously only uses a viable computational time. There are no existing products for both China and Europe in the past decades showing this level of detail which we now present, especially regarding the differentiations in vegetation cover (grassland and cropland, and different types of forests). Our study conceptualised a set of automated land classification processes based on the GEE platform and geospatial big data. It is important to discuss the challenges involved in mapping LC categories at urban scale in great detail. The value of precise knowledge about ecosystems in the urbanised environment is important for gaining a deeper understanding of environmental pressures at the urban scale in a global context. Our results contribute to a better capture of LC and our approach maximises the integrated efficiency of the relevant calculations. Furthermore, the extraction of different densities in the urbanisation processes can now be tackled at intercontinental scale. In the course of urbanisation, it is of little surprise that our results show a substantial growth in artificial surfaces that subsequently caused reductions in cultivated land. We are able to quantify the loss of vegetation at the scale of the urban area. That makes it possible to monitor the direction of changes and variations of this decline at a higher level of differentiation. Another important contribution is the successful spectral Remote Sens. 2021, 13, 1744 20 of 32 separability of forest areas. Deciduous forests store more carbon than coniferous forests, so the spectral separability of forest is essential [70].

Conclusions
It is worth investigating the definitions and criteria of existing thematic mapping systems in order to take full advantage of their capabilities when pulling information together and extracting the appropriate data for one's own research.
Our approach demonstrates that a requisite number of thematic classes can be derived for the mapping of specific LC categories that relate to the allocation of ES. We conclude that for reproduction purposes, a good balance between the necessary computational resources and the number of indices in use is as essential as deploying ground truth information to guarantee high quality products. Yet, there is growing imbalance between the increasing amount of satellite images and ground observations [70]. It is a challenge to gather a plentiful of samples for training and validation across continents. We hence illustrate the potential of samples collection from existing local and continental databases. We argue to keep collating ground truthing information in the light of global satellite image programs. Paying tribute to the high heterogeneity of urban areas, the more recent Sentinel data products that have been available since 2016 provide timely and continuous satellite images for LC monitoring. For less historical change detection research, the Sentinel series support mapping and monitoring at higher resolution than the Landsat series [71,72].
In order to improve the universality and model accuracy of this research method, the following aspects can be optimised in future research: First, more land-cover products as well as citizen science datasets (such as Open Street Map, OSM) can be further considered as the source of samples collection. Second, change detection methods can be employed as ways to expand the year for samples collection. Last, the accuracy evaluation system of sample points needs to be further constructed. We can, for example, use the purity and connectivity of LC to evaluate the reliability of sample points, or combine multiple product voting systems to improve the reliability of sample points, to build a higher accuracy classification model. This level of mapping detail in urban LC facilitates for the next stage of scientific analysis to quantify the benefits and losses of related ES features. Furthermore, the outcomes of this study may also help to assess the changes in ES at an intercontinental scale, and to support the recognition of ES's benefits and co-benefits for a more sustainable urban growth. This will be the object of research in part 2. Hence, our next task will be to estimate the equity in ES provision for this study area by further investigating the quantified spatial patterns of the diverse urbanisation processes.    Figure A1.
(a-f) Visualising the dynamic urban LC expansion of the study area by means of the thematic GAIA product.