Annual Wetland Mapping in Metropolis by Temporal Sample Migration and Random Forest Classiﬁcation with Time Series Landsat Data and Google Earth Engine

: Wetlands provide various ecosystem services to urban areas, which are crucial for sustainable urban management. With intensiﬁed urbanization, there has been marked loss of urban natural wetland, degradation, and related urban disasters in the past several decades. Rapid and accurate mapping of urban wetland extent and change is thus critical for improving urban planning toward sustainability. Here, we have developed a rapid method for continuous mapping of urban wetlands (MUW) by combining automatic sample migration and the random forest algorithm (SM&RF), the so-called MUW_SM&RF. Using time series Landsat images, annual training samples were generated through spectral angular distance ( SAD ) and time series analysis. Combined with the RF algorithm, annual wetland maps in urban areas were derived. Employing the Google Earth Engine platform (GEE), the MUW_SM&RF was evaluated in four metropolitan areas in different geographical and climatic regions of China from 1990 to 2020, including Tianjin, Hangzhou, Guangzhou, and Wuhan. In all four study areas, the generated annual wetland maps had an overall accuracy of over 87% and a Kappa coefﬁcient above 0.815. Compared with previously published datasets, the urban wetland areas derived using the MUW_SM&RF approach achieved improved accuracy and thus demonstrated its robustness for rapid mapping of urban wetlands. Urban wetlands in all four cities had variable distribution patterns and showed signiﬁcantly decreased trends in the past three decades. The annual urban wetland data product generated by the MUW_SM&RF can provide invaluable information for sustainable urban planning and management, so as for assessment related to the United Nation’s sustainable development goals.


Introduction
With continuous population growth and rapid economic development, urban areas have expanded significantly worldwide [1,2]. Massive urban expansion has brought about a slew of concerns, such as loss of natural habitat, reduction of biodiversity, and degradation of natural landscapes, posing enormous challenges for urban planning and development [3,4]. Wetlands serve urban areas by contributing to various ecosystem services, including water purification, climate mitigation, recreation, and tourism [5,6]. Targets 11.3 and 11.5 of the Sustainable Development Goals (SDGs) emphasize the necessity of sustainable urban development and ecosystem protection. Wetlands serve a valuable ades [22]. All of these cities have populations in excess of 10 million, and their urbanization rates exceed 80% and are continuing to rise [23]. The urban boundary of each study site is defined as the area within the latest metropolitan boundary of the city, separating the urban area from rural lands. These metropolitan boundaries were determined based on geographical knowledge and visual interpretation of the latest Google Earth images.

Imagery and Processing
All available 919 Landsat series images from the annual summer months, from June to September between 1990 and 2020, were obtained in GEE, including 566 Landsat 5, 102 Landsat 7, and 251 Landsat 8. To obtain consistent, high-quality images for annual urban wetland mapping, all images were subjected to a series of processes on the GEE. First, the QA band of the image, which flags the pixels disturbed by clouds and shadows, was used to eliminate bad-quality pixels within the image. Then, terrain shadows were removed

Data and Pre-Processing 2.2.1. Imagery and Processing
All available 919 Landsat series images from the annual summer months, from June to September between 1990 and 2020, were obtained in GEE, including 566 Landsat 5, 102 Landsat 7, and 251 Landsat 8. To obtain consistent, high-quality images for annual urban wetland mapping, all images were subjected to a series of processes on the GEE. First, the QA band of the image, which flags the pixels disturbed by clouds and shadows, was used to eliminate bad-quality pixels within the image. Then, terrain shadows were removed from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes. from each image by utilizing the digital elevation model (DEM) images and terrain algorithm provided by GEE. We used the SLC Gap-Fill algorithm provided by the USGS and implemented it on the GEE to eliminate the scan-line error strips on the Landsat 7 images [20]. Particularly, the temporal aggregation algorithm, which generates a composite image by calculating the statistics of both images over a defined period [24], was used to produce high-quality images annually from 1990 to 2020 in this study.

Sample Data
The initial sample data were collected in 2020 from a series of field investigations and collaborators. The acquisition dates of these field samples provided by collaborators varied and applied their own described system for classification based on their research topic. Thus, using images from 2020, we checked and modified these sample data to maintain their validity through visual interpretation, and they were relabeled with a uniform format following the definition of urban wetlands in this study, including urban wetlands (waterbody and vegetated wetland) and non-urban wetlands (land cover types other than waterbody and vegetated wetland in urban areas) ( Table 1). In areas where sample data were insufficient, we supplemented the sample by utilizing high-spatial resolution Google images and the visual interpretation approach. Ultimately, a total of 1740 samples were collected across all study sites in 2020, with Tianjin, Hangzhou, Wuhan, and Guangzhou corresponding 530, 309, 532, and 369 samples, respectively ( Figure 1). Based on the initial sample data obtained in 2020, sample data for the each year from 1990 to 2019 were generated at each study site utilizing the sample migration algorithm (see Section 3.2) which was developed for this study. In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and recategorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes.

Other Public Wetland Datasets
In this study, wetland maps extracted from the National Wetland Datasets of China (CAS_Wetlands) [10] and the National Land Cover datasets of China (ChinaCover) [25] in 2015 were used to compare with the generated result. Both the ChinaCover and CAS_Wetlands datasets were derived from Landsat images using an object-oriented classification approach. Wetlands in the ChinaCover were classified as surface water (lake, river, reservoir/pond) and wetland (swamp, marsh), while the CAS_Wetlands dataset involves more detailed other wetland types, e.g., lagoon, estuarine water, shallow marine water, tidal flat, salt pan, and coastal aquaculture pond. We extracted and re-categorized those wetland types from both datasets to keep the mapping categories consistent with our wetland classes.

Methodology
This study developed an annual mapping urban wetlands approach by combining sample migration and RF algorithm, i.e., the MUW_SM&RF. The workflow of this process is shown in Figure 2, and the specific components are as follows:

Methodology
This study developed an annual mapping urban wetlands approach by combining sample migration and RF algorithm, i.e., the MUW_SM&RF. The workflow of this process is shown in Figure 2, and the specific components are as follows:

Data Processing and Building Database of Time Series Feature Images
In this study, many features were extracted to enhance the image spectrum information and improve the accuracy of image classification. Three components of K-T transformations [26], including wetness, brightness, and greenness, were calculated from each image. Among these components, the wetness component can highlight the water content information of the image, which is commonly used in mapping wetlands [27]. In addition, several popular spectral indices were measured, including the Nominalized Difference Vegetation Index (NDVI) [28], Normalized Difference Built-up Index (NDBI) [29], Normalized Difference Water Index (NDWI) [30], and modified Normalized Difference Water Index (mNDWI) [31], which are useful in enhancing the information of major landscapes in urban areas, including vegetation, wetlands, water bodies, and urban areas. These indices are defined as Equations (1)-(4). Several texture feature bands were also derived to enhance the texture information, such as Angular Second Moment (ASM), Contrast, and Entropy. All spectral bands and features mentioned above were combined to generate feature images by the temporal aggregation method described in Section 2.2.1. Finally, the image database of time series features from 1990 to 2020 for each of the four study sites was created.

Data Processing and Building Database of Time Series Feature Images
In this study, many features were extracted to enhance the image spectrum information and improve the accuracy of image classification. Three components of K-T transformations [26], including wetness, brightness, and greenness, were calculated from each image. Among these components, the wetness component can highlight the water content information of the image, which is commonly used in mapping wetlands [27]. In addition, several popular spectral indices were measured, including the Nominalized Difference Vegetation Index (NDVI) [28], Normalized Difference Built-up Index (NDBI) [29], Normalized Difference Water Index (NDWI) [30], and modified Normalized Difference Water Index (mNDWI) [31], which are useful in enhancing the information of major landscapes in urban areas, including vegetation, wetlands, water bodies, and urban areas. These indices are defined as Equations (1)-(4). Several texture feature bands were also derived to enhance the texture information, such as Angular Second Moment (ASM), Contrast, and Entropy. All spectral bands and features mentioned above were combined to generate feature images by the temporal aggregation method described in Section 2.2.1. Finally, the image database of time series features from 1990 to 2020 for each of the four study sites was created. where ρ green , ρ red , ρ nir , and ρ swir are green, red, near-infrared, and short-wave infrared bands of Landsat images, respectively.

Methodologies for Sample Migrations
3.2.1. Extracting Unchanged Samples from 1990 to 2020 by Temporal Analysis Using the initial samples in this study (described in Section 2.2.2), we selected samples that have not changed significantly in type between 1990 and 2020 by analyzing their standard deviation of wetness (SD_WET) and the change in the NDVI and wetness values. As shown in Figure 3, a lower SD_WET value indicates a lesser change in land cover type during these decades (Figure 3A), whereas a higher SD_WET demonstrates a higher change probability ( Figure 3B,C). Therefore, we can use an SD_WET threshold to divide the initial samples into unchanged and changed parts throughout the whole period. Time series wetness and NDVI were derived from the initial samples between 1990 and 2020 to determine the optimum SD WET threshold. The variation characteristics of time series wetness and NDVI for each sample type were then analyzed. Finally, the SD-WET threshold was determined to obtain unchanged samples for each study site.

Methodologies for Sample Migrations
3.2.1. Extracting Unchanged Samples from 1990 to 2020 by Temporal Analysis Using the initial samples in this study (described in Section 2.2.2), we selected samples that have not changed significantly in type between 1990 and 2020 by analyzing their standard deviation of wetness (SD_WET) and the change in the NDVI and wetness values. As shown in Figure 3, a lower SD_WET value indicates a lesser change in land cover type during these decades (Figure 3A), whereas a higher SD_WET demonstrates a higher change probability ( Figure 3B,C). Therefore, we can use an SD_WET threshold to divide the initial samples into unchanged and changed parts throughout the whole period. Time series wetness and NDVI were derived from the initial samples between 1990 and 2020 to determine the optimum SD WET threshold. The variation characteristics of time series wetness and NDVI for each sample type were then analyzed. Finally, the SD-WET threshold was determined to obtain unchanged samples for each study site.

Extracting Unchanged Samples from 1990 to 2020 by Temporal Analysis
With the unchanged samples identified from the preceding step, we selected all water body and vegetated wetland samples among them to analyze their spectral characteristics, and then generated a reference spectral extent of the water body and vegetated wetland for each site. To enhance the difference in spectral information, besides the spectral bands, spectral indices and the main components of K-T transformations were also used. The spectral characteristics of the water body and vegetated wetland in each study site are shown in Figure 4. Compared to a vegetated wetland, the water body has a higher wetness level and a lower NDVI, NDWI, and mNDWI, making it possible to accurately separate them. According to the reference spectral ranges of vegetated wetlands and water bodies, the average of the ranges was selected as the reference spectral for each type at the four study sites.
ter body and vegetated wetland samples among them to analyze their spectral characteristics, and then generated a reference spectral extent of the water body and vegetated wetland for each site. To enhance the difference in spectral information, besides the spectral bands, spectral indices and the main components of K-T transformations were also used. The spectral characteristics of the water body and vegetated wetland in each study site are shown in Figure 4. Compared to a vegetated wetland, the water body has a higher wetness level and a lower NDVI, NDWI, and mNDWI, making it possible to accurately separate them. According to the reference spectral ranges of vegetated wetlands and water bodies, the average of the ranges was selected as the reference spectral for each type at the four study sites.

Reclassifying Changed Samples by the Reference Spectral and Decision-Tree Model
A decision-tree model was developed to relabel the changed samples identified in the above steps ( Figure 5). Firstly, the changed samples were divided into urban wetland and non-urban wetland samples using the wetness image and Otsu's algorithm (OTSU) [32]. OTSU is a well-known and efficient image segmentation algorithm that has been employed in various remote sensing studies [33]. Using the OTSU algorithm and the wetness images, it is possible to divide the images into two distinct parts based on the difference in wetness and thus separate the samples into urban and non-urban wetland samples.
Using the Spectral Angle Distance (SAD) (Equation (5)), the urban wetland samples obtained following the above steps (comprising both vegetated wetland and water body samples) were further subdivided into vegetated wetland and water body samples. SAD is the distance metric between two spectral vectors, commonly used to estimate the similarity between them, where the greater the similarity, the larger the SAD value [34].

Reclassifying Changed Samples by the Reference Spectral and Decision-Tree Model
A decision-tree model was developed to relabel the changed samples identified in the above steps ( Figure 5). Firstly, the changed samples were divided into urban wetland and non-urban wetland samples using the wetness image and Otsu's algorithm (OTSU) [32]. OTSU is a well-known and efficient image segmentation algorithm that has been employed in various remote sensing studies [33]. Using the OTSU algorithm and the wetness images, it is possible to divide the images into two distinct parts based on the difference in wetness and thus separate the samples into urban and non-urban wetland samples.
Using the Spectral Angle Distance (SAD) (Equation (5)), the urban wetland samples obtained following the above steps (comprising both vegetated wetland and water body samples) were further subdivided into vegetated wetland and water body samples. SAD is the distance metric between two spectral vectors, commonly used to estimate the similarity between them, where the greater the similarity, the larger the SAD value [34].
where θ represents the spectral angle and X i and Y i represent the two spectral vectors to be measured, respectively. The variable i represents the index number of the spectral bands, which varies from 1 to the total number of bands (N). SAD is the distance of spectral angle between the spectral vectors X i and Y i , which ranges from 0 to 1. A sample will be classified as a water body if the SAD between it and the reference water body spectrum is greater (Equation (6)). Conversely, when the SAD between it and the reference vegetated wetland spectrum is larger, it will be classified as a vegetated wetland. Using this rule, urban wetland samples were further subdivided into the water body and vegetated wetland. where TAR T1 represents the spectra of the given sample at time T1 and the Re f _Water and Re f _Wetland are the water body and vegetated wetland reference spectra, respectively. SAD(Tar T1 , Re f _Water) and SAD(Tar T1 , Re f _Wetland) represent the SAD value for a given sample with the water body and vegetated wetland reference spectra, respectively. Finally, by combining the unchanged samples and reclassified samples following the above processing, time series urban wetland sample datasets from 1990 to 2019 were generated for each study site. We divided the samples into two groups equally for training and testing.
A sample will be classified as a water body if the SAD between it and the reference water body spectrum is greater (Equation (6)). Conversely, when the SAD between it and the reference vegetated wetland spectrum is larger, it will be classified as a vegetated wetland. Using this rule, urban wetland samples were further subdivided into the water body and vegetated wetland.
Where 1 represents the spectra of the given sample at time 1 and the _ and _ are the water body and vegetated wetland reference spectra, respectively.
( 1 , _ ) and ( 1 , _ ) represent the value for a given sample with the water body and vegetated wetland reference spectra, respectively.
Finally, by combining the unchanged samples and reclassified samples following the above processing, time series urban wetland sample datasets from 1990 to 2019 were generated for each study site. We divided the samples into two groups equally for training and testing.

Mapping Urban Wetlands and Accuracy Assessment
The RF algorithm is the most commonly used supervised machine learning method in various remote sensing studies, which combines a number of classification trees to obtain an optimal solution to complex problems [35]. It has been widely discussed in the literature and has been proven to be a suitable method for classifying land cover using satellite images at medium and high resolution [14,36]. The GEE platform provides a complete random forest algorithm, which also includes sample collection, feature extraction, and detailed accuracy evaluation parts [37]. Based on RF, annual urban wetland maps from 1990 to 2020 for the four study sites were generated by using time series feature images and training samples.

Mapping Urban Wetlands and Accuracy Assessment
The RF algorithm is the most commonly used supervised machine learning method in various remote sensing studies, which combines a number of classification trees to obtain an optimal solution to complex problems [35]. It has been widely discussed in the literature and has been proven to be a suitable method for classifying land cover using satellite images at medium and high resolution [14,36]. The GEE platform provides a complete random forest algorithm, which also includes sample collection, feature extraction, and detailed accuracy evaluation parts [37]. Based on RF, annual urban wetland maps from 1990 to 2020 for the four study sites were generated by using time series feature images and training samples.
An accuracy assessment was carried out using all of the testing samples in this study. Specifically, we used a confusion matrix to evaluate the generated classification results, containing four indicators, namely overall accuracy, the Kappa coefficient, and user's and producer's accuracy.

Determination of SD_WET Threshold and Sample Migration Results
To determine reasonable SD_WET thresholds for extracting samples of the water body and vegetated wetland not changing in type between 1990 and 2020, we examined the changes in NDWI and wetness for all initial samples (2020) at different SD_WET thresholds. To illustrate the detailed process of determining this threshold, we have chosen Site C as an example. Firstly, the threshold selection can be narrowed down to 0.01-0.05. Within this range, samples of all types maintain their spectral characteristics related to NDVI and wetness that they should have, i.e., water body and vegetated wetland samples have higher wetness values ( Figure 6A), and NDVI values for water body samples are the lowest ( Figure 6B). Accordingly, samples in this range have high purity and have barely undergone a change in type between 1990 and 2020. Then, by analyzing the number of samples of each type within the above range (0.01-0.05), we can determine that the most reasonable threshold is 0.05. Utilizing this threshold will ensure that sufficient samples of each type are obtained ( Figure 6C). Using the same approach, the optimal SD_WET thresholds were determined to be 0.04, 0.06, and 0.05 for Sites A, B, and D, respectively. The sample numbers for each type after migration in the four study sites are shown in Figure 7. In 2020, the number of initial samples of each type was more evenly distributed across the four study site, while the number of samples generated by migration during 1990-2019 has changed, indicating that samples have changed differently over the three decades. Moreover, the amounts of each sample type in the four study sites have various change trends characterized by different SD values from 1990 to 2019. The sample numbers for each type after migration in the four study sites are shown in Figure 7. In 2020, the number of initial samples of each type was more evenly distributed across the four study site, while the number of samples generated by migration during 1990-2019 has changed, indicating that samples have changed differently over the three decades. Moreover, the amounts of each sample type in the four study sites have various change trends characterized by different SD values from 1990 to 2019.

Accuracy Assessment of Mapping Urban Wetlands
The confusion matrix showed that our annual urban wetland maps from 1990 to 2020 achieved good classification results in the four study sites. As shown in Figure 8, the OA for each year is greater than 87% in all study sites, and the mean values from 1990 to 2020 in four sites are greater than 90.31%. For all study sites, the Kappa coefficient is greater than 0.82 from 1990 to 2020, indicating that the mapping results and ground-truth exhibit a high level of agreement. Moreover, as shown in Figure 9, both PA and UA of both urban wetland types have high values in all study sites, with the mean value of PA and UA for each site being higher than 84.5%.

Accuracy Assessment of Mapping Urban Wetlands
The confusion matrix showed that our annual urban wetland maps from 1990 to 2020 achieved good classification results in the four study sites. As shown in Figure 8, the OA for each year is greater than 87% in all study sites, and the mean values from 1990 to 2020 in four sites are greater than 90.31%. For all study sites, the Kappa coefficient is greater than 0.82 from 1990 to 2020, indicating that the mapping results and ground-truth exhibit a high level of agreement. Moreover, as shown in Figure 9, both PA and UA of both urban wetland types have high values in all study sites, with the mean value of PA and UA for each site being higher than 84.5%.

Spatial Patterns and Temporal Trends of Studied Urban Wetlands
As shown in Figure 10, urban wetlands dominated by water bodies were widely distributed over the four study sites. The majority of vegetated wetlands are found adjacent to open water bodies. The vegetated wetlands and water bodies in small areas experienced evident changes and thus have lower occurrence frequency. The subset views ( Figure 10A-D) show the most representative change of urban wetlands for each study site. In Tianjin (Site A), small water bodies suffered the most severe change ( Figure 10A). In Hangzhou (Site B), changes of vegetated wetland have notable spatial heterogeneity ( Figure 10B). In Wuhan (Site C), urban wetlands with small areas spreading throughout the city experienced significant changes ( Figure 10C). In Guangzhou (Site D), vegetated wetlands distributed around a larger water body suffered the most severe losses ( Figure 10D).

Spatial Patterns and Temporal Trends of Studied Urban Wetlands
As shown in Figure 10, urban wetlands dominated by water bodies were widely distributed over the four study sites. The majority of vegetated wetlands are found adjacent

Spatial Patterns and Temporal Trends of Studied Urban Wetlands
As shown in Figure 10, urban wetlands dominated by water bodies were widely distributed over the four study sites. The majority of vegetated wetlands are found adjacent In general, urban wetlands in each city experienced a significant areal decline from 1990 to 2020 ( Figure 11). The most remarkable areal reduction of urban wetlands occurred in Tianjin (Site A), with a rate of 67.71% and an average loss rate of about 2.82 km 2 per year, followed by Guangzhou (35.51% and 5.53 km 2 ) and Wuhan (28.85% and 5.74 km 2 ). The declining trend of urban wetlands in Hangzhou (Site B) is relatively slow between 1990 and 2020, with a rate of 18.93% and the smallest reduction of about 0.57 km 2 per year. In general, urban wetlands in each city experienced a significant areal decline from 1990 to 2020 ( Figure 11). The most remarkable areal reduction of urban wetlands occurred in Tianjin (Site A), with a rate of 67.71% and an average loss rate of about 2.82 km 2 per year, followed by Guangzhou (35.51% and 5.53 km 2 ) and Wuhan (28.85% and 5.74 km 2 ). The declining trend of urban wetlands in Hangzhou (Site B) is relatively slow between 1990 and 2020, with a rate of 18.93% and the smallest reduction of about 0.57 km 2 per year.  (Site(A-D)). The value between 1 and 31 indicates the number of urban wetland occurrences from 1990 to 2020. The subset views (A-D) present the detailed changes of urban wetland for each study site with a 10-year interval from 1990 to 2020, respectively. The backgrounds in the views are the Landsat images.  Figure 12 overlays the products from the ChinaCover, the CAS_Wetlands, and the MUW_SM&RF. The ChinaCover and the CAS_Wetlands had a larger extent of water bodies than the result in this study. This difference can be attributed to the images and method. In this study, the images were generated by the temporal aggregation method rather than at a specific time. Thus the extent of extracted water varied. In addition, vegetated wetlands in this study are more detailed than the results of both the CAS_Wetlands and ChinaCover. Both the CAS_Wetlands and ChinaCover used the object-oriented classification method to identify vegetated wetlands. With a spatial resolution of 30 m, it is difficult to segment small vegetated wetland objects in urban areas. In this study, we generated reliable urban wetland samples by sample migration and used a pixel-based RF algorithm for urban wetland extraction. As a result, the range of vegetated wetlands obtained in this study is more refined.  Figure 12 overlays the products from the ChinaCover, the CAS_Wetlands, and the MUW_SM&RF. The ChinaCover and the CAS_Wetlands had a larger extent of water bodies than the result in this study. This difference can be attributed to the images and method. In this study, the images were generated by the temporal aggregation method rather than at a specific time. Thus the extent of extracted water varied. In addition, vegetated wetlands in this study are more detailed than the results of both the CAS_Wetlands and ChinaCover. Both the CAS_Wetlands and ChinaCover used the object-oriented classification method to identify vegetated wetlands. With a spatial resolution of 30 m, it is difficult to segment small vegetated wetland objects in urban areas. In this study, we generated reliable urban wetland samples by sample migration and used a pixel-based RF algorithm for urban wetland extraction. As a result, the range of vegetated wetlands obtained in this study is more refined.

Performance of the MUW_SM&RF
In this study, MUW_SM&RF was effective at mapping annual urban wetlands over the long term. By developing an effective sample migration approach, MUW_SM&RF can produce sufficient and high-quality training samples for the classifier and generate annual urban wetland maps. The generated maps allow us to understand the effects of urbanization on wetlands, which is significant not only for the conservation of wetland ecosystems, but also for sustainable urban planning. Additionally, the open-accessed Landsat time series imagery and the computing capability provided by the GEE platform enabled MUW_SM&RF to produce annual and time series urban wetland maps with uninterrupted data in a moderate 30 m spatial resolution.

Performance of the MUW_SM&RF
In this study, MUW_SM&RF was effective at mapping annual urban wetlands over the long term. By developing an effective sample migration approach, MUW_SM&RF can produce sufficient and high-quality training samples for the classifier and generate annual urban wetland maps. The generated maps allow us to understand the effects of urbanization on wetlands, which is significant not only for the conservation of wetland ecosystems, but also for sustainable urban planning. Additionally, the open-accessed Landsat time series imagery and the computing capability provided by the GEE platform enabled MUW_SM&RF to produce annual and time series urban wetland maps with uninterrupted data in a moderate 30 m spatial resolution.
The MUW_SM&RF performed well in mapping annual urban wetlands in the four study sites and can be applied to other cities worldwide. SDG targets 11.3 and 11.5 focus on sustainable urban planning and management and emphasize the importance of urban ecosystem conservation. Annual urban wetlands data generated by MUW_SM&RF offer valuable insights into sustainable urban programming, as well as an important reference for assessing these goals. Additionally, the Wetland City accreditation scheme (Resolution The MUW_SM&RF performed well in mapping annual urban wetlands in the four study sites and can be applied to other cities worldwide. SDG targets 11.3 and 11.5 focus on sustainable urban planning and management and emphasize the importance of urban ecosystem conservation. Annual urban wetlands data generated by MUW_SM&RF offer valuable insights into sustainable urban programming, as well as an important reference for assessing these goals. Additionally, the Wetland City accreditation scheme (Resolution XII.10), sponsored by the Ramsar Convention on Wetlands, encourages cities that value the conservation and proper application of wetlands, as well as highlighting the importance of urban wetlands. MUW_SM&RF can provide an effective method to evaluate the conservation and utilization status of urban wetlands for this scheme.

Inadequacies of the MUW_SM&RF
The accuracies of the annual urban wetland maps were inevitably affected by some factors. Even though the temporal aggregation algorithm was adopted and most of the pixels affected by the cloud and shadows were removed from each image, it remains difficult to obtain high-quality images for low-latitude cloud-prone areas (i.e., Site D). Additionally, the MUW_SM&RF may have some shortcomings when it is applied to produce appropriate amounts of samples. If extreme weather events (drought or flood) occur in the study area, the landscapes will be drastically altered, making it difficult for MUW_SM&RF to generate appropriate samples. It will directly affect the accuracy of urban wetland extraction. For example, in 1998 and 2016, the extreme floods in Wuhan (Site C) have greatly changed the pattern of urban wetlands. As a result, the samples of urban wetland in 1998 and 2016 for Site C produced by the MUW_SM&RF were much more uneven than those in other years (Figure 7). To overcome the above shortcoming, we can shorten the detection period, divide the entire time period into several shorter periods, and produce training samples in each shorter period using the MUW_SM&RF to mitigate the effect of extreme weather on change detection.

Conclusions
In this study, continuous urban wetland mapping of four metropolises of China, namely Tianjin, Hangzhou, Guangzhou, and Wuhan, was conducted using Landsat time series from 1990 to 2020 by applying an effective approach (MUW_SM&RF) with the support of the GEE cloud platform. The MUW_SM&RF can automatically produce sufficient and high-quality training samples through spectral angular distance (SAD) and time series analysis and generate annual urban wetland maps using the random forest classification. The generated annual wetlands maps had an overall accuracy of over 87%. Compared to previously published datasets, the urban wetlands delineated using the MUW_SM&RF achieved better mapping accuracy, thus demonstrating its robustness for rapid mapping of urban wetlands. The generated annual urban wetland dataset revealed significant declining trends for each city. The MUW SM&RF can be applied to metropolises worldwide, and the annual urban wetland data it generates can offer valuable insights into sustainable urban programming and facilitate the evaluation of related SDG goals.