A Modified Self-adaptive Method for Mapping Annual 30-m Land Use/Land Cover Using Google Earth Engine: A Case Study of Yangtze River Delta

Land Use/Land Cover (LULC) data are important in applications such as land management, agricultural monitoring, ecological service research, and climate change assessment (Zhao et al., 2015; Mao et al., 2018; Teluguntla et al., 2018; Gong et al., 2019). With advances in remote sensing technology toward providing satellite images, the corresponding datasets have been effectively applied to classify LULC types at different spatial scales from global to local (Mohajane et al., 2018; Zhang et al., 2019; Li et al., 2020). Among the different remote sensing datasets, Landsat data provides free global coverage of satellite images with a long history in comparison with other open-access remotely sensed data (e.g., MODIS, Sentinel-2), which is advantageous for LULC classification tasks (Li and Gong, 2016; Huang et al., 2017; Liu et al., 2020; Mao et al., 2021).

The conventional large-area long-term LULC classification method uses remote sensing-derived time series with a supervised non-parametric classifier (Zhai et al., 2018; Bullock et al., 2020). In areas with strong landscape heterogeneity, there are limitations in employing this kind of supervised classifier (Ghosh et al., 2014; Müller et al., 2015). Firstly, with strong heterogeneity of landscape, the phenomenon of ‘different objects with the same spectrum, different spectrums within the same object’ often occurs (Long, 2014). Further research is needed to develop high-precision classifiers (Chen et al., 2018). Secondly, there are challenges in obtaining a large amount of long-term training data (Viana et al., 2019). Finally, with a large research area, a large amount of remote sensing data need to be processed, requiring strong computational and storage capabilities (Daldegan et al., 2019; Tamiminia et al., 2020).

The random forest (RF) classifier is a commonly used supervised classifier for LULC mapping (Chakraborty et al., 2016). In general, the use of remote sensing time series data combined with the RF classifier can obtain a LULC map with high classification accuracy in an area with strong landscape heterogeneity (Zeng et al., 2020). This is because the phenological information provided by time series data can partially solve the problem of ‘different objects with the same spectrum, different spectrums within the same object’ (Lopes et al., 2017). However, challenges remain regarding large-area LULC classification based on time series (Ji et al., 2020). In particular, only by using a large amount of high-precision training data can higher classification accuracy be obtained (Ghorbanian et al., 2020).

Over the past decade, it has been demonstrated that it is efficient to use existing LULC maps as a source of training data (Viana et al., 2019). This is advantageous because it: 1) allows classification in an automated manner without the need for interactive manual training data, 2) provides a potentially large and geographically distributed training dataset, and 3) enables satellite data to be mapped with existing LULC maps (Pandey et al., 2021). However, existing LULC maps may have some misclassified results (Mack et al., 2017). Therefore, the use of existing LULC maps as a source should be carefully considered to ensure that the generated training data have a good enough level of accuracy (Li et al., 2017).

Computational power, data storage, management, and processing times have also traditionally been restrictions for using remote sensing time series over large areas (Anchang et al., 2020). Google Earth Engine (GEE) not only provides powerful computing ability for free but can also directly access a variety of open source data (Hird et al., 2017; Capolupo et al., 2020). The workflow designed on GEE not only solves the limitation of computing power and data source but also helps other researchers to conduct similar research (Wu, 2020; Xu et al., 2021).

The aims of the present study are: 1) to classify Landsat time series data supported by GEE to obtain high-precision annual LULC maps of the Yangtze River Delta (YRD) from 1992 to 2015 and 2) to analyze the characteristics of LULC changes in the YRD region. Specifically, we address the following research questions: 1) how to get a high-precision sample dataset generated from multi-source LULC products? 2) How to construct a feature space, so that time series classification can be carried out in a cloudy and rainy area (such as YRD)? 3) How to build a local adaptive classifier to improve the classification accuracy of LULC in an area with strong landscape heterogeneity?

The Yangtze River Delta (YRD) is situated in the eastern China and covers four provinces: Anhui Province, Jiangsu Province, Zhejiang Province, and Shanghai Municipality (Fig. 1), with a surface area of 348 000 km² (Chen et al., 2019). The topography of the YRD is dominated by plains, hillsides, and mountains (Feng et al., 2018). The average elevation is 140.17 m (Shen et al., 2008). This region is under a monsoon climate; the annual average temperature for the growing season (March to October) is 18℃ to 22℃ and the average annual precipitation for the growing season is 800 to 1400 mm (Wan et al., 2020).

Figure 1.  Geographical location of Yangtze River Delta (YRD) in the eastern China

The YRD is one of China’s most economically developed regions and plays an important role in the social and economic development of the country (Yu et al., 2019). In 2015, the YRD accounted for 16.06% of China’s population, and its GDP accounted for 23.50% of the country’s GDP (Zhang et al., 2020). According to the Statistical Yearbook from China National Knowledge Infrastructure (http://www.cnki.net/), the populations and GDPs of the YRD region increased significantly from 1992 to 2015. This growth in the population and economic development has been accompanied by rapid urbanization and great changes in land use: the loss of farmland and expansion of urban areas, in the region (Xu et al., 2018).

To construct a sample set for LULC classification, we used several datasets of well-recognized LULC products, including 1992–2015 year by year European Space Agency Climate Change Initiative (ESA-CCI 300) LULC products, 2001–2015 annual MCD12Q1 (MODIS LULC cover) products, 2000 and 2010 LULC data developed by the China National Basic Geographic Information Center (GlobeLand30), and the 2015 Finer Resolution Observation and Monitoring of Global Land Cover (From-GLC), developed by Tsinghua University, China (Table 1). The GEE platform can access MCD12Q1 products directly, and these do not require uploading into the application (Gorelick et al., 2017). The other three sets of LULC products need to be uploaded to the GEE platform. All work was performed at a 30 m resolution. The ESA-CCI300 and MCD12Q1 products were resampled on the GEE platform.

Landsat 5/7/8 Surface Reflectance (SR) data from 1992 to 2015 are the main remote sensing images for classification and can be accessed directly in the GEE (Xu et al., 2019). Previous research has shown that the Normalized Difference Vegetation Index (NDVI) is sensitive to vegetation characteristics, the Normalized Difference Water Index (NDWI) can identify bodies of water, and the Normalized Difference Built-up Index (NDBI) can distinguish built-up areas effectively (Simonetti et al., 2015; Bailly et al., 2017; Wagle et al., 2020). In the study, the spectral indices (NDVI, NDWI, and NDBI) calculated by the Landsat SR products in the YRD are also used. The spectral indices were introduced for the calculation of the spectral-temporal metric.

Topographic features can affect regional climatic conditions and vegetation growth, so these are often used as auxiliary data for LULC classification (Adepoju and Adelabu, 2020). Topographic features generated from DEM (Digital Elevation Model) include elevation (affecting temperature and precipitation), slope, and aspect (affecting the results of solar radiation and vegetation growth). To describe terrain features, we use the Shuttle Radar Topography Mission (SRTM) data to generate topographic features (Hurni et al., 2019). The GEE platform has built-in SRTM data, which can be used directly.

We extracted the SR Tier 1 datasets of the Landsat 5/7 (B1–B5 and B7 bands) and Landsat 8 (B2–B7 bands) in the growing season during 1992–2015. The first pre-processing step was to build cloud-free Landsat tile mosaics for each year. For this, we used the CFMASK algorithm provided by the GEE and the pixel quality assessment (pixel_qa) information available in the Landsat collection (Li et al., 2019). The cloud-masked Landsat scenes were combined to produce NDVI, NDWI, and NDBI for the annual growth period of each year. This procedure required optimal spectral contrast and separability amongst the LULC classes.

The annual Landsat mosaics were generated with statistical reducers including median, standard deviation, minimum, and maximum. After building the annual mosaics, the next step was to build the feature space for the self-adaptive RF classifier. For this process, we used the compositional, spectral, and temporal information extracted from the annual image mosaics. The NDVI values of all cloud-free pixels from each year were divided into quartiles; the median values of the higher quartile were considered as the wet-season image. Finally, a total of 39 bands consisting of 4 × 9 temporal feature bands and 3 topographic feature bands (elevation, slope, and aspect) were available for LULC classification.

Supervised classification usually requires a certain number of training samples and verification samples (Zhao et al., 2020). Typically, traditional research uses manual visual interpretation to obtain sample points (Huang et al., 2020). For a study with a large area and long study period, such a method presents considerable practical difficulties (Ghorbanian et al., 2020). This study proposes a new method to achieve highly credible sample points. The specific steps are as follows:

(1) Based on From-GLC (2015), all other LULC products were reclassified, and the LULC cover classes were divided into eight classes including: cropland, forest, shrubland, grassland, wetland, water bodies, built-up land, and unused land.

(2) ESA-CCI (1992–2015), MCD12Q1 (2001–2015), GlobeLand30 (2000, 2010), and the From-GLC (2015) data were overlaid, and pixels with completely consistent LULC types that had not changed from 1992 to 2015 were selected. ESA-CCI300 and MCD12Q1 were reduced to 30 m resolutions in the GEE using the Reduce Resolution function. GEE performs nearest neighbor resampling by default.

(3) From the selected pixels, training sample points were randomly selected; the number of points was slightly adaptive according to the area ratios of the different LULC types.

(4) The created samples still had a high probability of spectral-temporal signature confusion among LULC classes. Therefore, we applied the k-means clustering technique to refine the samples. The clustering process is essentially divided into three distinct processes: 1) classification-trimmed likelihood calculation, 2) cluster computation, and 3) mean discriminant factor value calculation.

Overall, the samples were distributed uniformly and randomly throughout the study area; however, for some LULC types with relatively small and patchy areas (such as unused lands) a relatively dense distribution occurred. Therefore, we exported the refined sample point data to ArcGIS and carried out reselection at 1500 m resolution to reduce the correlation between samples. Through this process, we obtained 20 859 reference sample points. Finally, we randomly selected 10 430 points from the sample library to be used for classification training. To ensure the authenticity and reliability of the verification samples, we selected them through visual inspection. In the end, we obtained 2857 reliable sample points in the study area to verify the LULC classification results.

The classification process considered a spatial and temporal stratification, individually training and classifying all 9368 Landsat tiles that cover the YRD. The LULC was divided into eight classes which were exactly the same as used for sample points (details see section 3.2). Considering that LULC classes are highly susceptible to climatic conditions and topography throughout the YRD territory, this approach allowed the classification models to better identify all LULC classes. To minimize the impact of spatial stratification, part of the training set was shared among different classification models. Thereby, the 3 × 3 scenes sample points were used to train the self-adaptive RF algorithm (Fig. 2). Equally important, the temporal stratification was also intended to minimize the impact of spectral and radiometric differences among the Landsat sensors in the classification results, since, for each year, a classifier was trained considering only images obtained by a single satellite.

Figure 2.  Flowchart of self-adaptive RF (random forest) LULC (Land Use/Land Cover) classification

The self-adaptive RF models are constructed by the smileRandomForest function in the GEE (Kakooei and Baleghi, 2020). In order to use the smileRandomForest function, two parameters need to be set: the number of decision trees to create per class (numberoftrees) and the minimum size of a terminal (minleaf) (Gumma et al., 2020). Finally, through repeated comparative experiments, we set numberoftrees to 100 and minleaf to 10 in all models.

To better homogenize these classification results, we applied a spatial filter and a temporal filter, capable of minimizing abrupt and sometimes unrealistic variations, simultaneously considering these two dimensions. First, we used a spatial filter in the edge area of each scene. There are some pixels with inconsistent classification results in these edge areas. These pixels are merged in accordance with the majority agreement rule by the mode function in GEE (Gorelick et al., 2017). All of the classified scenes were merged into a LULC map collection year by year. Next, we applied a temporal filter to the LULC map collection. We detected whether the LULC classes in year N − 1 and year N + 1 are consistent pixel by pixel. If so, the LULC classes in year N must be consistent with those LULC classes. Except for the first and last years, we iterated the above process year by year and obtained final LULC classification results for the YRD.

We used several measures to assess and statistically compare the accuracy of our classifications. First, we calculated the classification error matrices for each year. From these matrices, we then quantified overall accuracy (OA). Second, we calculated producer’s accuracy (PA), user’s accuracy (UA), and F₁-score, which is the harmonic mean between user’s accuracy and producer’s accuracy. PA, UA, and F₁-score can be calculated for each class i as follows:

where r is the number of classes and n_ij is the element of the confusion matrix in row j and column i, that is, the count of elements of class j classified as class i. PA_i, UA_i, and (F₁)_i stand for PA, UA, and F₁-score for class i, respectively.

The F₁-score is particularly useful for class-level accuracy assessment, as it gives equal importance to both PA and UA by combining PA and UA into a single measure that can be compared across confusion matrices. Finally, we estimated the unbiased area (using the sample weight obtained with our reference sample dataset) following a standard good practice protocol.

To evaluate the quality of the 30-m LULC results for the YRD using multi-source LULC products as training samples, an accuracy assessment of the classification was conducted. A confusion matrix, which is the primary tool used in remote sensing for accuracy assessment, was used to evaluate the accuracy. A total of 2857 verification samples were collected for all LULC classes. The evaluation was carried out based on the error matrix produced and analyzed the accuracy of the 30 m LULC classification results.

Fig. 3 reports the classification accuracies that quantify the level of agreement between the classification and the sample data. The overall accuracies of the classification results are greater than 84.50% for every year, and the average overall accuracy is 86.33%. This shows that the overall classification level of agreement is high. The F₁-score annual accuracies were close to 90% for the following main classes: cropland, forest, and water bodies. The F₁-scores of built-up areas ranged between 75% and 85%, mainly due to the misclassification of many small villages in the northern part of the study area. Because shrubland, grassland, wetland, and unused land account for a very small proportion in the study area, we analyzed these LULC classes as a whole. The F₁-scores were < 80% for all these classes for each year except for 2014 and 2015. The poor classification accuracy of these classes can be explained on the one hand because their proportions are low and on the other hand because the sample points of these LULC classes are small. In general, these results indicate quite reasonable classification accuracies. The results show that using multi-source LULC products to generate training data is a viable and effective option.

Figure 3.  Classification accuracy for the LULC (Land Use/Land Cover) in the Yangtze River Delta, China from 1992 to 2015

By using the Landsat time series data with the self-adaptive RF method, we developed annual LULC datasets for the YRD from 1992 to 2015 (Fig. 4). The LULC classification results showed extensive change across the study area. Specifically, cropland was the most extensive LULC type but continuously decreased from 1992 to 2015. In contrast, built-up land continuously and significantly expanded due to rapid urbanization. We observed that the expanded built-up areas were mainly concentrated near the peri-urban region, possibly due to the topographic conditions. The expansion of built-up land was mainly centered on cropland around the city, which led to the continuous decrease in cropland area. Fig. 5 also indicates that areas of forest and water bodies remained largely unchanged during the study period.

Figure 4.  Annual LULC (Land Use/Land Cover) maps for the YRD from 1992 to 2015. Maps for four selected years (1992, 2000, 2008, and 2015) are enlarged

Figure 5.  Area percentage of different LULC (Land Use/Land Cover) classes for the YRD from 1992 to 2015

Fig. 5 illustrates how the primary LULC types of the YRD changed from 1992 to 2015. It can be seen that the greatest change in the YRD occurred in the case of cropland, whose proportion decreased from 44.6% in 1992 to 42.5% in 2015. This was followed by forest, whose proportion decreased from 37.5% to 36.2%. In contrast, built-up land’s proportion increased from 0.8% to 2.1%. Bodies of water increased by 0.2%. The smallest changes were in grassland, whose proportion decreased by just 0.1%. Shrubland, wetland, and bare land occupy relatively small proportions of the study area. Compared with the primary LULC types (cropland, forest, water bodies, and built-up land) they did not change noticeably. In addition, due to their small proportions in the study area, they were more affected by classification uncertainty. Therefore, changes in these LULC types were excluded from the analysis in this paper.

It can be seen from Figs. 6–8 that our results are similar to other LULC products in terms of geographic distribution, and the overall layout of cropland, forest land, and water bodies is similar. Overall, our results are more accurate in identifying small ground objects than other LULC products.

Figure 6.  Comparison with other common LULC (Land Use/Land Cover) products for the year 2015. The first column shows Landsat-8 OLI images with 30 m resolution, the second is our result with 30 m resolution, the third is ESA-CCI 300 with 300 m resolution, and the fourth is MCD12Q1 with 500 m resolution. Row (a) is a selected typical urban area in plain, row (b) is a selected typical transitional area from mountain to plain, and row (c) is a selected typical mix area

Figure 7.  Comparison of our results with GlobeLand30 for the year 2010 in the YRD, China

Figure 8.  Comparison of our results with From-GLC for the year 2015 in the YRD, China

Compared with ESA-CCI300 products, the main difference lies in the built-up lands in the plain areas and the forest in the mountain areas (Fig. 6). In plain areas, more small built-up land areas were identified in our results, which may be related to the higher spatial resolution of our results compared to those of the other LULC products. In mountainous areas, our results showed more forest area compared to results from the other LULC products, with ESA-CCI300 classifying more as croplands. Compared with MCD12Q1 products, the main difference is that MCD12Q1 misclassified many shrublands while our method correctly classified them as croplands. In plain areas, large tracts of land are misclassified as built-up areas in MCD12Q1, whereas our results are more refined. We detected these inconsistent areas visually through Landsat-8 OLI images and found that our results are more accurate among the results from all the LULC products.

In Fig. 7, we can see that our results are generally consistent with GlobeLand30 in mountainous areas, and the distribution of forest and water bodies in plain areas are also similar. The main difference is in the distribution of built-up lands and croplands in plain areas. GlobeLand30 showed many small built-up lands, and our results classify these areas as croplands. This is mainly because, in addition to Landsat data, GlobeLand30 also uses other multispectral images with higher spatial resolution and processes them manually, so that it can classify more built-up lands in the mixed pixels of rural areas (Fig. 7 left). Due to the introduction of other multispectral data, however, whatever the percentage of built-up land is in the mixed pixel, it will still be classified as built-up land (Fig. 7 right). In our results, due to the influence of MODIS and ESA data on the training points generated, the number of built-up land samples in rural areas was small. As a result, only contiguous built-up lands can be accurately identified in the classification results, while small built-up lands in rural areas showed a poor level of accuracy.

In Fig. 8, it can be seen that From-GLC 2015 misclassifies part of the cropland as forest in the plain area. Our results do not show such a misclassification. The main reason for this type of misclassification is that there are some spectral differences between different tiles of Landsat data. When we employed the self-adaptive RF classification model, the classification features are mainly phenological features, and the training data in 3 × 3 adjacent tiles are used for model training; this effectively reduces the impact of spectral differences. The self-adaptive RF model can be effectively applied because we have a sufficient amount of training data.

LULC classification over large areas and a long study period is challenging due to the large volume of data pre-processing required and the cost and difficulty of collecting representative training data that enable classification models to be both globally consistent and locally reliable (Millard and Richardson, 2015; Mahdianpari et al., 2020; Liu et al., 2021). With ready-to-use data, such as those from Landsat products, and the employment of non-parametric classifiers, the major challenge is training data collection (Gumma et al., 2020).

The detailed results illustrated in Fig.s 6–8 are representative of the results across the study area. In general, our classification results are similar to other LULC products, which indicated that our results are reasonable. The high overall accuracies and generally high F₁-score in the main classes (that quantify the level of agreement between the classification and the training pool data) underscore the utility of the self-adaptive RF classification approach. However, the effectiveness of the self-adaptive RF classification method depends upon the availability of sufficient local training data. In this study, this was not an issue as the training pool data derived from multi-sources LULC products were geographically well distributed.

In this study, a novel method was used to classify Landsat data using high quality training data derived from multi-source LULC products. A training data pool was extracted from the multi-source LULC products by judicious quality and k-means clustering filtering. In addition, the training data selection was undertaken in a geographically systematic manner while ensuring that the selected class ratios were the same as the YRD ground object’s ratios. This is advantageous as it: 1) enables the classification to be undertaken in an automated manner without the need for manual training data collection, 2) provides a large geographically distributed training data set, and 3) results in the generation of a 30 m Landsat LULC product with the same classification legend as the multi-source LULC products.

Although our method of generating training data has many advantages, it also has certain shortcomings. Firstly, due to the limitation of multi-source LULC products, the training data generated a large number of accurate sample points for contiguous LULC types but fewer for the small LULC types (such as built-up lands in rural areas). Secondly, different LULC products have different classification systems, resulting in fewer LULC classes in the training data generated in our study. Finally, for obtaining multi-source LULC products of high time consistency and data quality, we do not use 30 m LULC products after 2015 because there is no free 30 m LULC product at a global scale after 2015 when we finished the experiment. Therefore, the period of our research is selected from 1992 to 2015. In future work, we will consider whether the Landsat data has changed and generate sample data with more LULC classes and a longer period.

Based on the four types of LULC products (ESA-CCI300, MOD12Q1, Globland30, and From-GLC), a data training pool was extracted by judicious quality and k-means clustering filtering. A self-adaptive RF classification was employed to classify 9368 Landsat image tiles. Using these methods, we constructed a map set of LULC in the YRD region from 1992 to 2015 and analyzed the spatiotemporal changes in the main LULC types in the YRD during this period. The primary conclusions of this study are as follows:

(1) The results indicated that the self-adaptive RF classifier can be used successfully to classify Landsat time series into LULC data, with an average overall classification accuracy of 86.33%. LULC classification results with 30 resolution of YRD appeared geographically reasonable and were similar in comparison to other LULC products. Therefore, the approach is suitable for analyzing changes in LULC.

(2) The training sample generation approach from multi-source LULC products described in this study can provide a large number of sample points in an acceptable quality level. This is because k-means clustering filtering can effectively eliminate non-homogeneous sample points, so that the generated sample points have higher accuracy.

(3) The construction of feature space is also an important advantage of this study. The feature space can eliminate the negative impacts of clouds and rain to the time series data. Moreover, it can reduce the difference between different Landsat titles. Therefore, the construction of feature space is essential for time series classification in a large cloudy and rainy area (such as the YRD).

In future work, we will continue to study the use of multi-source LULC products to generate sample points. We will consider applying some spectral unmixed models to the sample point generation process, to obtain more sample points of non-dominant LULC types (e.g., built-up land in rural areas). In the meanwhile, we consider extending our research workflow from level 1 to level 2, so that we can determine whether our workflow performs the same way with more LULC types.