Mapping Tidal Flats with Landsat 8 Images and Google Earth Engine : A Case Study of the China ’ s Eastern Coastal Zone circa 2015

Accurate and up-to-date tidal flat mapping is of much importance to learning how coastal ecosystems work in a time of anthropogenic disturbances and rising sea levels, which will provide scientific instruction for sustainable management and ecological assessments. For large-scale and high spatial-resolution mapping of tidal flats, it is difficult to obtain accurate tidal flat maps without multi-temporal observation data. In this study, we aim to investigate the potential and advantages of the freely accessible Landsat 8 Operational Land Imager (OLI) imagery archive and Google Earth Engine (GEE) for accurate tidal flats mapping. A novel approach was proposed, including multi-temporal feature extraction, machine learning classification using GEE and morphological post-processing. The 50 km buffer of the coastline from Hangzhou Bay to Yalu River in China’s eastern coastal zone was taken as the study area. From the perspective of natural attributes and unexploited status of tidal flats, we delineated a broader extent comprising intertidal flats, supratidal barren flats and vegetated flats, since intertidal flats are major component of the tidal flats. The overall accuracy of the resultant map was about 94.4% from a confusion matrix for accuracy assessment. The results showed that the use of time-series images can greatly eliminate the effects of tidal level, and improve the mapping accuracy. This study also proved the potential and advantage of combining the GEE platform with time-series Landsat images, due to its powerful cloud computing platform, especially for large scale and longtime tidal flats mapping.


Introduction
Tidal flats, often defined as sandy and muddy flats, are the important parts of coastal zones and highly productive, providing numerous minerals as well as biological and oceanic resources for human beings [1,2].They are essential not only for marine animals and migratory birds, but for the protection of coastal zones against maritime inundation including storm surges [3,4].A fairly large number of tidal flats are distributed along the coastline of China, and they are important land resources from a perspective of economic construction [5].Intertidal flats are major components of tidal flats, which topographically feature tidal inundation [2].Coastal vegetation is likely to colonize the upper margins of intertidal flats (supratidal flats) with years of vegetation succession [6,7], which can be defined as supratidal vegetated flats (hereinafter referred to as vegetated flats) [8], or salt marshes in terms of the presence of halophytic vegetation (Figure 1) [9,10].Due to their high development potential and wide extent, intertidal flats and supratidal flats are generally reclaimed for coastal development and thus significantly shrinking [4,11].Therefore, there is an urgent need to investigate the distribution and status of unexploited tidal flats to allow the implementation of extra conservation strategies and adjustment of land-use planning.We mapped the extent of tidal flats in this study comprising intertidal flats, supratidal barren flats and vegetated flats.The upper limit of supratidal flats is artificial borders (roads, seawalls and offshore buildings) [10], and thus the tidal flats that were situated inside the artificial borders and reclaimed for coastal development were explicitly excluded.
Traditional field surveys are generally employed to investigate tidal flats but hard to be carried out due to their difficult accessibility and physical constraints over a large extent of tidal flats (>1000 km 2 ) [12].Light Detection and Ranging (LiDAR) and Synthetic Aperture Radar (SAR) techniques have been applied for tidal flats mapping with high precision at regional scale [13,14], however, they are not efficient and economical because of low performance at national scale and high cost for airborne acquisitions [12].
An optical remote-sensing technique, as a macro-scale and reliable alternative to ground-survey methods in remote or inaccessible regions [15,16], has been increasingly utilized to delineate tidal flats.The final maps are predominantly generated through visual interpretation [17,18], supervised classification methods [19][20][21] and object-based approaches [22].These methods merely produced accurate tidal flat maps at a specific time, since tidal flats are only exposed fully during ebb tide and inundated during flood tide.Some studies extracted waterlines to demarcate exposed tidal flats from water in the temporal observations-from threshold segmentation [23][24][25][26] to manual digitization [27].The knowledge of accurate tidal levels at the time of image acquisition is the common requirement for waterline extraction methods, which hinders the upscaling analyses for larger scales.
Satellite sensors can capture land cover when it repeatedly passes by, and thus the periodic change information has been delineated by time-series images.Thanks to recording tidal inundation and vegetation phenology, time-series remotely sensed images have been applied for mapping intertidal zones [28], mangroves [29], paddy rice [30], cropland [31] and rubber plantations [32].We consider it possible to map tidal flats more accurately with serial remotely sensed images, but computational capacity needs to be handled due to the sheer volume of data processing.
Recently, some cloud computation platforms for geospatial data processing have become available with high-performance computational power and big data-processing tools [30], including Google Earth Engine (GEE), Amazon Web Service (AWS) and National Aeronautics and Space Administration (NASA) Earth Exchange (NEX).They possess abundant imagery archives and data products, and also can be easily carried out for thematic mapping as well as spatiotemporal analyses, with the support of parallel-processing computation and advanced machine learning algorithms [33].The advent of cloud computation platforms has changed the way of storing, managing, processing and analyzing of massive amounts of large-scale geospatial data [34].
This study explicitly aimed at accurate tidal flats mapping over large area, using time-series Landsat 8 Operational Land Imager (OLI) imagery.The GEE platform was employed for high-performance data processing, which also provided freely accessible remote sensing imagery.The coastal area from Hangzhou Bay to Yalu River in China's eastern coastal zone was taken as the study area.Firstly, we introduced time-series images to eliminate the effect of tidal levels at different observation times, and analyzed the characteristics of tidal flats with time-series remote-sensing images; secondly, 48 spectral features were calculated using the GEE platform and used for classification tasks; thirdly, 1633 samples were used to train a random forest (RF) classifier and generate an initial classification map; finally, a morphological post-processing approach was proposed to obtain the final map.The study result may be helpful in monitoring the coastal zone and developing coastal conservation measures at nominal national scale.

Study Area
The coastal wetlands in China are geographically divided into two different kinds by Hangzhou Bay, on the north of which sandy and muddy flats predominate while on the south rocky beach accounts for a large proportion [35].The study area covers the region from Hangzhou Bay to Yalu River in China's eastern coastal zone (29° 25′-41° 6′N, 117° 25′E-124° 16′E), including Zhejiang (partly), Jiangsu, Shandong, Hebei and Liaoning Province as well as several economically developed cities-Shanghai, Tianjin, Qingdao and Dalian.The study area is defined on the basis of a costal buffer (Figure 2) that consists of 10 km buffer landward and 40 km buffer seaward along the coastline, following the regulation of the National Multipurpose Investigation of the Coastal Zone and Tidal Wetland Resources [36].According to the figures in the China Second Wetlands Survey [37], the area of wetlands of coasts and seashores covering our study area was 39,445 km 2 .

Satellite Imagery and Auxiliary Data
There are 33 tiles of Landsat Worldwide Reference System 2(WRS-2) paths/rows covering the coastal buffer completely (Figure 2).A total of 1802 Landsat 8 OLI surface reflectance images in Tier 1 were acquired between 1 January 2014 and 31 December 2016, which were used to extend the data range circa 2015, filling the data gap because of bad observations (e.g.cloud and cloud shadow) removal.Landsat 8 surface reflectance data in Tier 1 were available on GEE as image collections.These images were orthorectified and generated based on Level 1 Precision Terrain (L1TP) corrected images as well as processed through Landsat Surface Reflectance Code (LaSRC) [38], which are considered of high quality and ready-to-use [29].Cloud and cloud shadow affecting observation quality were identified by CFMask in surface reflectance data [39].Therefore, clear observations without clouds and cloud shadow were employed in this study, and they include visible (Blue, 0.452-0.512μm; Green, 0.533-0.590μm; Red, 0.636-0.673μm), near infrared (NIR, 0.851-0.879μm) and shortwave infrared (SWIR1, 1.566-1.651μm; SWIR2, 2.107-2.294μm) bands.
For defining the administrative division of study area, the Global Administrative Areas dataset (GADM) [40] was used as auxiliary data, which includes administrative boundaries of national, provincial and even city level, and helped to calculate the detected area of tidal flats in terms of administrative division.Additionally, the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) that comprises World Vector Shorelines (WVS) [41] was applied to extract shorelines of the study area.We downloaded the latest shoreline files of version 2.3.7 released on 15 June 2017 and then generated a 50 km coastal buffer using ArcGIS software.All available

Study Area
The coastal wetlands in China are geographically divided into two different kinds by Hangzhou Bay, on the north of which sandy and muddy flats predominate while on the south rocky beach accounts for a large proportion [35].The study area covers the region from Hangzhou Bay to Yalu River in China's eastern coastal zone (29 • 25 -41 • 6 N, 117 • 25 E-124 • 16 E), including Zhejiang (partly), Jiangsu, Shandong, Hebei and Liaoning Province as well as several economically developed cities-Shanghai, Tianjin, Qingdao and Dalian.The study area is defined on the basis of a costal buffer (Figure 2) that consists of 10 km buffer landward and 40 km buffer seaward along the coastline, following the regulation of the National Multipurpose Investigation of the Coastal Zone and Tidal Wetland Resources [36].According to the figures in the China Second Wetlands Survey [37], the area of wetlands of coasts and seashores covering our study area was 39,445 km 2 .

Satellite Imagery and Auxiliary Data
There are 33 tiles of Landsat Worldwide Reference System 2(WRS-2) paths/rows covering the coastal buffer completely (Figure 2).A total of 1802 Landsat 8 OLI surface reflectance images in Tier 1 were acquired between 1 January 2014 and 31 December 2016, which were used to extend the data range circa 2015, filling the data gap because of bad observations (e.g.cloud and cloud shadow) removal.Landsat 8 surface reflectance data in Tier 1 were available on GEE as image collections.These images were orthorectified and generated based on Level 1 Precision Terrain (L1TP) corrected images as well as processed through Landsat Surface Reflectance Code (LaSRC) [38], which are considered of high quality and ready-to-use [29].Cloud and cloud shadow affecting observation quality were identified by CFMask in surface reflectance data [39].Therefore, clear observations without clouds and cloud shadow were employed in this study, and they include visible (Blue, 0.452-0.512µm; Green, 0.533-0.590µm; Red, 0.636-0.673µm), near infrared (NIR, 0.851-0.879µm) and shortwave infrared (SWIR1, 1.566-1.651µm; SWIR2, 2.107-2.294µm) bands.
For defining the administrative division of study area, the Global Administrative Areas dataset (GADM) [40] was used as auxiliary data, which includes administrative boundaries of national, provincial and even city level, and helped to calculate the detected area of tidal flats in terms of administrative division.Additionally, the Global Self-consistent, Hierarchical, High-resolution Geography Database (GSHHG) that comprises World Vector Shorelines (WVS) [41] was applied to extract shorelines of the study area.We downloaded the latest shoreline files of version 2.3.7 released on 15 June 2017 and then generated a 50 km coastal buffer using ArcGIS software.All available datasets, i.e., Landsat imagery and auxiliary data, were projected into the World Geodetic System 1984 (WGS84).

Methods
Aiming to develop a consistent, automated and provincial scale method for mapping tidal flat extent of the coastal zone from Hangzhou Bay to Yalu River, the workflow was deployed by combining GEE platform with Google Earth Pro and ArcGIS software, and it is divided into four parts: (1) field investigation and reference sample selection, 2) statistics-based time series image processing, 3) RF machine learning algorithm via GEE, and 4) morphological post-processing for tidal flats (Figure 3).

Methods
Aiming to develop a consistent, automated and provincial scale method for mapping tidal flat extent of the coastal zone from Hangzhou Bay to Yalu River, the workflow was deployed by combining GEE platform with Google Earth Pro and ArcGIS software, and it is divided into four parts: (1) field investigation and reference sample selection, (2) statistics-based time series image processing, (3) RF machine learning algorithm via GEE, and (4) morphological post-processing for tidal flats (Figure 3).

Field Investigation and Reference Sample Selection
Considering previous studies in coastal land use/cover [42][43][44] and tidal flats mapping [9,21,45,46], non-tidal flats classes were categorized into land, water and offshore ponds, and their detailed descriptions are listed in Table 1.Field investigation work was conducted during July, 2018.After our investigation, we found that all these categories can be obviously recognized from highresolution imagery provided by Google Earth Pro (GEP) software.
Following the imagery interpretation rules developed for coastal land use in Tianjin, China [17] and coastal wetlands in land cover project of CORINE (Coordination of information on the environment) program [47], random stratified sampling was adopted to obtain training and validation data within 50 km costal buffer for algorithm development and accuracy assessment by visually interpreting high-resolution images circa 2015 from GEP. Intertidal flat samples were selected in the inundated area lying outside seawalls, and the samples of supratidal barren flats and vegetated flats were extracted from barren expanses of little inundation and vegetated area, respectively, which are geographically adjacent to the intertidal flats (Figure 4 a-c).Reclaimed tidal flats, which are distinguishable by parcellation and dike systems, were not included.In total, 2076 sample points in six classes were obtained by manual interpretation (Figure 1), 80% of which were used for training and the remaining 20% for accuracy assessment (Table 1).

Field Investigation and Reference Sample Selection
Considering previous studies in coastal land use/cover [42][43][44] and tidal flats mapping [9,21,45,46], non-tidal flats classes were categorized into land, water and offshore ponds, and their detailed descriptions are listed in Table 1.Field investigation work was conducted during July, 2018.After our investigation, we found that all these categories can be obviously recognized from high-resolution imagery provided by Google Earth Pro (GEP) software.
Following the imagery interpretation rules developed for coastal land use in Tianjin, China [17] and coastal wetlands in land cover project of CORINE (Coordination of information on the environment) program [47], random stratified sampling was adopted to obtain training and validation data within 50 km costal buffer for algorithm development and accuracy assessment by visually interpreting high-resolution images circa 2015 from GEP. Intertidal flat samples were selected in the inundated area lying outside seawalls, and the samples of supratidal barren flats and vegetated flats were extracted from barren expanses of little inundation and vegetated area, respectively, which are geographically adjacent to the intertidal flats (Figure 4a-c).Reclaimed tidal flats, which are distinguishable by parcellation and dike systems, were not included.In total, 2076 sample points in six classes were obtained by manual interpretation (Figure 1), 80% of which were used for training and the remaining 20% for accuracy assessment (Table 1).

Statistics-Based Time-Series Images Processing
Tidal flats in this study are characterized by the expanse of mud in terms of tidal inundation and vegetation coverage.The reflectance of tidal flats varies along with the wavelength of the Landsat satellite sensor [23].
The basic elements of the muddy coastal area can be summarized as soil, water and vegetation.Therefore, many studies considered spectrum-derived indices for coastal mapping.Several studies employed some vegetation and water indices to map tidal flats and coastal vegetation [6,12,24], such as Modified Normalized Difference Water Index (mNDWI) [48], Normalized Difference Vegetation Index (NDVI) [49], Land Surface Water Index (LSWI) [50], Modified Soil-Adjusted Vegetation Index (mSAVI) [51] and Enhanced Vegetation Index (EVI) [52].Therefore, we introduced the spectral indices mentioned above as sensitive features for tidal flats along with six visual and infrared bands of the Landsat OLI sensor (Section 2.2).In addition, Soil Brightness from Tasseled Cap Transformation [53] was also included as sensitive feature because of its high sensitivity for barren soil surface.The definitions of indices mentioned above are formularized as:

Statistics-Based Time-Series Images Processing
Tidal flats in this study are characterized by the expanse of mud in terms of tidal inundation and vegetation coverage.The reflectance of tidal flats varies along with the wavelength of the Landsat satellite sensor [23].
The basic elements of the muddy coastal area can be summarized as soil, water and vegetation.Therefore, many studies considered spectrum-derived indices for coastal mapping.Several studies employed some vegetation and water indices to map tidal flats and coastal vegetation [6,12,24], such as Modified Normalized Difference Water Index (mNDWI) [48], Normalized Difference Vegetation Index (NDVI) [49], Land Surface Water Index (LSWI) [50], Modified Soil-Adjusted Vegetation Index (mSAVI) [51] and Enhanced Vegetation Index (EVI) [52].Therefore, we introduced the spectral indices mentioned above as sensitive features for tidal flats along with six visual and infrared bands of the Landsat OLI sensor (Section 2.2).In addition, Soil Brightness from Tasseled Cap Transformation [53] was also included as sensitive feature because of its high sensitivity for barren soil surface.The definitions of indices mentioned above are formularized as: where ρ BLUE , ρ GREEN , ρ RED , ρ NIR , ρ SWIR1 and ρ SWIR2 are the surface reflectance in blue (B2), green (B3), red (B4), near-infrared (B5) and shortwave infrared (B6 and B7) bands of Landsat OLI imagery, respectively.We can calculate the spectral indices for each pixel using time-series images.However, the number of covered images for different pixels might be quite different.For example, in some areas, there are more than 80 available images, while in some areas only 20 images exist.Therefore, the feature dimensions of different pixels are not the same, bringing a huge challenge for subsequent processing.Normalizing high-dimensional features and retaining useful information is the key problem of this work.Moreover, the reduction of the feature dimension could significantly reduce the computational cost.
Due to diurnal inequality of tide, the water levels at different observations (once per 16 days) are not the same.Assuming that with long-term and periodic remote-sensing observations, the observed lowest tidal level could be approximate to the true value, we can thus obtain a more accurate tidal map than by using only one image.
From the perspective of feature statistics, the features fluctuate at different observations.Figure 5 illustrates the time-series variation of spectral indices in intertidal, supratidal barren and vegetated flats in Yellow River Delta (Figure 5a). Figure 5b    To figure out the spectral change pattern in a three-year observation, we computed the minimum, maximum, mean and standard deviation of original surface reflectance and spectral indices for all sample points.Intertidal flats commonly have higher LSWI value than other land use when it comes to the maximum of LSWI during the observation period (Figure 6b), but always with minus values when the NDVI minimum is counted (Figure 6a).Hence, based on such a statistical result, Landsat time-series images were reduced to form an image with a total of 48 bands, including the minimums, maximums, means and standard deviations of 6 spectral indexes and 6 original spectral channels.Intertidal flat is susceptible to flood and ebb tides, and its values of LSWI and mNDWI markedly fluctuated during the observation period (Figure 5c), with higher maximum and mean values due to high soil moisture content.The NDVI values of intertidal and supratidal barren flat remained steady at below 0 and around 0.1, respectively, while vegetated flat had a great fluctuation around 0.25.The value change of EVI and mSAVI shared the same pattern with that of NDVI, and soil brightness did not explicitly indicate three different flats.
To figure out the spectral change pattern in a three-year observation, we computed the minimum, maximum, mean and standard deviation of original surface reflectance and spectral indices for all sample points.Intertidal flats commonly have higher LSWI value than other land use when it comes to the maximum of LSWI during the observation period (Figure 6b), but always with minus values when the NDVI minimum is counted (Figure 6a).Hence, based on such a statistical result, Landsat time-series images were reduced to form an image with a total of 48 bands, including the minimums, maximums, means and standard deviations of 6 spectral indexes and 6 original spectral channels. .RF classifier consists of a number of decision trees, involving voting mechanism for improving accuracy [54], with high cost in computation when the number of trees is set very large.As the composite image with 48 bands was used as target image, the number of decision trees was limited to 100, keeping a trade-off between accuracy and efficiency.

Morphological Post-Processing for Tidal Flats
Through the steps described above, the preliminary tidal flat map was obtained.Considering that intertidal flats are the base of tidal flats and adjacent topographically to supratidal barren flats and vegetated flats (Figure 1), the map was exported from GEE to ArcGIS software for further postprocessing.A morphological method has been used for burned area mapping in some studies [55,56], and it is hugely useful and necessary to merge topographically adjacent and homogeneous components together as well as to omit unlikely parts.In this study, the morphological method is possible, because artificial maritime borders are almost large enough to be presented in 30m satellite RF classifier consists of a number of decision trees, involving voting mechanism for improving accuracy [54], with high cost in computation when the number of trees is set very large.As the composite image with 48 bands was used as target image, the number of decision trees was limited to 100, keeping a trade-off between accuracy and efficiency.

Morphological Post-Processing for Tidal Flats
Through the steps described above, the preliminary tidal flat map was obtained.Considering that intertidal flats are the base of tidal flats and adjacent topographically to supratidal barren flats and vegetated flats (Figure 1), the map was exported from GEE to ArcGIS software for further post-processing.A morphological method has been used for burned area mapping in some studies [55,56], and it is hugely useful and necessary to merge topographically adjacent and homogeneous components together as well as to omit unlikely parts.In this study, the morphological method is possible, because artificial maritime borders are almost large enough to be presented in 30m satellite images (Figure 4a,c,d) and they are obvious demarcation lines to separate tidal flats from land or reclaimed ones.Hence, 30m spatial solution was applied to pick up eligible tidal flats' pixels.
The connected pixels labeled as intertidal flats were firstly merged into patches with the kernel of eight-connected neighbors, and smaller patches with pixel count less than 81 [57] were removed.An iterative procedure of aggregating eight connected neighbors was then performed around each qualified pixel patch of intertidal flats, and the pixels labeled as supratidal barren flats or vegetated flats eight-connected with confirmed patches were selected.Finally, the resultant pixel patches located in unreasonable places, which make up less than 1% (26.4km 2 ) of final tidal flats area, were manually removed in a couple of minutes.For instance, the patches are far from the coastline and have actually been reclaimed.Figure 7 illustrates an example of morphological process for tidal flats over Yellow River Delta area (Figure 5a).Most of the pixels were well classified (Figure 7a), however, some errors still occurred to supratidal barren flats and vegetated flats, which unreasonably existed in a place far away from coastline.Thus, to exclude these false detections, the pixels identified as intertidal flats after connected aggregation (Figure 7b) and smaller patches removal (Figure 7c) were considered as a foundation to aggregate the connected pixels labeled as supratidal barren flats or vegetated flats (Figure 7d) into the final tidal flat map.
images (Figure 4a, c, d) and they are obvious demarcation lines to separate tidal flats from land or reclaimed ones.Hence, 30m spatial solution was applied to pick up eligible tidal flats' pixels.
The connected pixels labeled as intertidal flats were firstly merged into patches with the kernel of eight-connected neighbors, and smaller patches with pixel count less than 81 [57] were removed.An iterative procedure of aggregating eight connected neighbors was then performed around each qualified pixel patch of intertidal flats, and the pixels labeled as supratidal barren flats or vegetated flats eight-connected with confirmed patches were selected.Finally, the resultant pixel patches located in unreasonable places, which make up less than 1% (26.4km 2 ) of final tidal flats area, were manually removed in a couple of minutes.For instance, the patches are far from the coastline and have actually been reclaimed.Figure 7 illustrates an example of morphological process for tidal flats over Yellow River Delta area (Figure 5a).Most of the pixels were well classified (Figure 7a), however, some errors still occurred to supratidal barren flats and vegetated flats, which unreasonably existed in a place far away from coastline.Thus, to exclude these false detections, the pixels identified as intertidal flats after connected aggregation (Figure 7b) and smaller patches removal (Figure 7c) were considered as a foundation to aggregate the connected pixels labeled as supratidal barren flats or vegetated flats (Figure 7d) into the final tidal flat map.

Accuracy Assessment
Accuracy assessment is necessary for map generation from remote-sensing data.We derived validation data from the remaining 20% reference sample points, which were mentioned in Section 2.3.1.The samples of land, water and offshore ponds were integrated as non-tidal flats samples and named as Others, but three categories of tidal flats were retained.Then, we worked out a confusion matrix to evaluate the accuracy of the resultant map (Table 2).

Accuracy Assessment
Accuracy assessment is necessary for map generation from remote-sensing data.We derived validation data from the remaining 20% reference sample points, which were mentioned in Section 2.3.1.The samples of land, water and offshore ponds were integrated as non-tidal flats samples and named as Others, but three categories of tidal flats were retained.Then, we worked out a confusion matrix to evaluate the accuracy of the resultant map (Table 2).

Tidal Flat Map
A 30 m tidal flat map circa 2015 (Figure 8) covering the coastal zone from Hangzhou Bay to Yalu River was obtained by using RF machine learning algorithm with Landsat time-series images via GEE.The proposed approach obtained a high overall accuracy of 94.4% in collaboration with a small number of training sample points.Both producer's and user's accuracy were basically above 90% (Table 2).But vegetated flats had a lower producer's accuracy at 78.3%, which was mainly caused by the similar phenology between coastal vegetation and inland plantation.The tidal flats from Hangzhou Bay to Yalu River in China are mostly distributed along the coastlines of Jiangsu, Shandong and Liaoning Province, with a total area of 4629.67 km 2 , including 3959.90 km 2 intertidal flats, 293.53 km 2 supratidal barren flats and 376.24 km 2 vegetated flats (Table 3).The tidal flats in Jiangsu Province account for the largest proportion of about 42%, most of which are intertidal flats (up to 1767.56 km 2 ).Liaoning and Shandong Province hold the nearly same area of intertidal flats of around 620 km 2 , which is almost twice as large as that of Shanghai.But the area of vegetated flats in Shanghai is significantly larger, representing 124.74 km 2 nearly as large as that of Jiangsu Province.As for supratidal barren flats, Shandong Province owns the largest area of 259.52 km 2 , although most supratidal flats of other places have been reclaimed for aquaculture ponds or construction land [5].

Extensive Tidal Flats towards Land
Some researches merely considered intertidal flats as tidal flats and excluded supratidal vegetated flats that could be defined as salt marsh from the extent of tidal flats [9,21].However, from the perspective of tidal flats development, supratidal flats should be seen as a part of tidal flats for exploitation and management [2,8].
Tidal flats develop continuously towards the sea, leading to a broad and low-sloping extent [2].The intertidal flat is the major component of tidal flats and almost without any significant vegetation.

Extensive Tidal Flats towards Land
Some researches merely considered intertidal flats as tidal flats and excluded supratidal vegetated flats that could be defined as salt marsh from the extent of tidal flats [9,21].However, from the perspective of tidal flats development, supratidal flats should be seen as a part of tidal flats for exploitation and management [2,8].
Tidal flats develop continuously towards the sea, leading to a broad and low-sloping extent [2].The intertidal flat is the major component of tidal flats and almost without any significant vegetation.Landward lies supratidal flat and below the mean low tide line is the subtidal flat [7].All three components are topographically adjacent [7].That is, it is feasible to recognize supratidal flats from satellite images if intertidal flats are confirmed, while subtidal flats that are concealed by water are hardly captured by spectral satellite sensors.In this study, we selected reference samples of supratidal and intertidal flats in terms of the study of Li et al. [17], and utilized a morphological post-process to identify the extent of supratidal flats as intertidal flats were confirmed.The upper limit of supratidal flats is the artificial borders that are explicitly recognizable from 30 m pixels (shown in Figure 4).Overall, we delineated a broader extent of tidal flats towards land, which has been not reclaimed.

Sensitive Testing of Temporal Window and Training Samples
In this study, we delineated the annual extent of tidal flats in 2015 through a 3-year temporal window.We randomly set 80% reference points as training data and the remaining 20% for validation.It was exactly the considered allocation after we tested how the number of training points and temporal windows affected the overall accuracy (Figure 9).
Initially, we tested the classification accuracy under the condition where 18-month temporal window was set from 1 October 2014 to 31 March 2016 and the ratio of training samples to validation samples was 4:1.Then, a temporal window was added with a six-month interval; that is, the first day and last day of the previous testing period have been advanced and extended by three months, respectively, but the ratio did not change.We have carried out each temporal window 30 times through randomly selecting training and validation samples in the ratio of 4:1, and finally worked out the mean of testing result-from 18-to 36-month temporal windows.When examining the performance of 36-month plus temporal windows, GEE faced that computation timed out.Thus, how three-year plus temporal windows affected the overall accuracy was unknown, but explicitly a 36-month temporal window was ideal investigating period with 80% reference samples points for training to obtain reasonable overall accuracy (Figure 9a).
In order to understand how the proportion of reference samples for training affects the overall accuracy under the condition of 36-month observation, we randomly selected training samples from reference samples 30 times in different proportions, and worked out the mean and standard deviation of overall accuracy (Figure 9b).Obviously, a higher proportion leads to higher overall accuracy.But considering the problem of overfitting, 80% of reference samples for training is a trade-off, through which satisfactory accuracy can also be yielded.
Landward lies supratidal flat and below the mean low tide line is the subtidal flat [7].All three components are topographically adjacent [7].That is, it is feasible to recognize supratidal flats from satellite images if intertidal flats are confirmed, while subtidal flats that are concealed by water are hardly captured by spectral satellite sensors.In this study, we selected reference samples of supratidal and intertidal flats in terms of the study of Li et al. [17], and utilized a morphological postprocess to identify the extent of supratidal flats as intertidal flats were confirmed.The upper limit of supratidal flats is the artificial borders that are explicitly recognizable from 30 m pixels (shown in Figure 4).Overall, we delineated a broader extent of tidal flats towards land, which has been not reclaimed.

Sensitive Testing of Temporal Window and Training Samples
In this study, we delineated the annual extent of tidal flats in 2015 through a 3-year temporal window.We randomly set 80% reference points as training data and the remaining 20% for validation.It was exactly the considered allocation after we tested how the number of training points and temporal windows affected the overall accuracy (Figure 9).
Initially, we tested the classification accuracy under the condition where 18-month temporal window was set from 1 October 2014 to 31 March 2016 and the ratio of training samples to validation samples was 4:1.Then, a temporal window was added with a six-month interval; that is, the first day and last day of the previous testing period have been advanced and extended by three months, respectively, but the ratio did not change.We have carried out each temporal window 30 times through randomly selecting training and validation samples in the ratio of 4:1, and finally worked out the mean of testing result-from 18-to 36-month temporal windows.When examining the performance of 36-month plus temporal windows, GEE faced that computation timed out.Thus, how three-year plus temporal windows affected the overall accuracy was unknown, but explicitly a 36month temporal window was ideal investigating period with 80% reference samples points for training to obtain reasonable overall accuracy (Figure 9a).
In order to understand how the proportion of reference samples for training affects the overall accuracy under the condition of 36-month observation, we randomly selected training samples from reference samples 30 times in different proportions, and worked out the mean and standard deviation of overall accuracy (Figure 9b).Obviously, a higher proportion leads to higher overall accuracy.But considering the problem of overfitting, 80% of reference samples for training is a trade-off, through which satisfactory accuracy can also be yielded.

Comparison with Other Available Tidal Flat Maps
The resultant tidal flat map was compared with other available maps with the same spatial resolution and same time period.i.e., remote-sensing monitoring data of land use in China in 2015 (hereinafter referred to as LUCC 2015) and the tidal flat map derived from Murray et al. [21] (hereinafter referred to as Murray's).LUCC 2015 was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn),and derived

Comparison with Other Available Tidal Flat Maps
The resultant tidal flat map was compared with other available maps with the same spatial resolution and same time period.i.e., remote-sensing monitoring data of land use in China in 2015 (hereinafter referred to as LUCC 2015) and the tidal flat map derived from Murray et al. [21] (hereinafter referred to as Murray's).LUCC 2015 was provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn),and derived from visual interpretation based on Landsat 8 imagery.This dataset was updated from land-use data in 2010, covering cultivated land, woodland, grassland, water bodies, built-up land and unused land.Tidal flats were included in the water bodies, defined as offshore land between high-tide and low-tide level [18].Murray's was generated from GEE in collaboration with the RF machine learning algorithm [21].Both of tidal flat maps delineated the intertidal zone.
Ours, Murray's and LUCC 2015 maps were imported into ArcGIS for comparison and spatial analysis.In order to implement equal and objective analysis, we only extracted intertidal flats for comparison.The tidal flats area of three maps in terms of administrative division and the proportion of identical area to our map were calculated (Table 4).Murray's and our map represent a high correlation with R 2 value of 0.962 and a slope value of 1.532 as shown in Figure 10b, while the correlation between LUCC 2015 and our map is relatively low with R 2 value of 0.530 and a slope value of 0.410 (Figure 10a).Explicitly, LUCC 2015 delineated a smaller extent of tidal flats, where there is only 24.8% agreement to our map.By checking reference images, the ratio of identical area below 10% in coastal cities such as Hangzhou, Ningbo and Shaoxing, resulted from inconsistent detected area.In some cities like Jinzhou and Panjin situated in Liaohe Estuary, both LUCC 2015 and our map exhibited nearly equal areas, but showed low consistency.The dominant sources of disagreement are due to the different definition of tidal flats and potentially missing area.As shown in Figure 11, reclaimed tidal flats were included as tidal flats in LUCC 2015 (Figure 11b) but excluded in ours (Figure 11a), which are explicitly bounded by artificial marine facilities.In addition, a fairly large extent of intertidal flats in the east was not mapped in LUCC 2015.Supratidal barren flats and vegetated flats were also our target map classes.Although we note that there were some kinds of unused land in LUCC 2015, like saline-alkali land, swampland and bare soil delineating the same extent of supratidal barren flats and vegetated flats potentially, that was beyond the scope of this study.Another potentially missing area was in Jiangsu Province.The tidal flat area in Yancheng derived using the proposed method (1249.25 km 2 ) is much larger than LUCC 2015, while LUCC 2015 captured 441.02km 2 tidal flats.That mostly resulted from the intertidal flats called the radial sand ridges (RSR)(Figure 12a), which are located in region marked as A in Figure 8; the kind of tidal flats Explicitly, LUCC 2015 delineated a smaller extent of tidal flats, where there is only 24.8% agreement to our map.By checking reference images, the ratio of identical area below 10% in coastal cities such as Hangzhou, Ningbo and Shaoxing, resulted from inconsistent detected area.In some cities like Jinzhou and Panjin situated in Liaohe Estuary, both LUCC 2015 and our map exhibited nearly equal areas, but showed low consistency.The dominant sources of disagreement are due to the different definition of tidal flats and potentially missing area.As shown in Figure 11, reclaimed tidal flats were included as tidal flats in LUCC 2015 (Figure 11b) but excluded in ours (Figure 11a), which are explicitly bounded by artificial marine facilities.In addition, a fairly large extent of intertidal flats in the east was not mapped in LUCC 2015.Supratidal barren flats and vegetated flats were also our target map classes.Although we note that there were some kinds of unused land in LUCC 2015, like saline-alkali land, swampland and bare soil delineating the same extent of supratidal barren flats and vegetated flats potentially, that was beyond the scope of this study.Explicitly, LUCC 2015 delineated a smaller extent of tidal flats, where there is only 24.8% agreement to our map.By checking reference images, the ratio of identical area below 10% in coastal cities such as Hangzhou, Ningbo and Shaoxing, resulted from inconsistent detected area.In some cities like Jinzhou and Panjin situated in Liaohe Estuary, both LUCC 2015 and our map exhibited nearly equal areas, but showed low consistency.The dominant sources of disagreement are due to the different definition of tidal flats and potentially missing area.As shown in Figure 11, reclaimed tidal flats were included as tidal flats in LUCC 2015 (Figure 11b) but excluded in ours (Figure 11a), which are explicitly bounded by artificial marine facilities.In addition, a fairly large extent of intertidal flats in the east was not mapped in LUCC 2015.Supratidal barren flats and vegetated flats were also our target map classes.Although we note that there were some kinds of unused land in LUCC 2015, like saline-alkali land, swampland and bare soil delineating the same extent of supratidal barren flats and vegetated flats potentially, that was beyond the scope of this study.Another potentially missing area was in Jiangsu Province.The tidal flat area in Yancheng derived using the proposed method (1249.25 km 2 ) is much larger than LUCC 2015, while LUCC 2015 captured 441.02km 2 tidal flats.That mostly resulted from the intertidal flats called the radial sand ridges (RSR)(Figure 12a), which are located in region marked as A in Figure 8; the kind of tidal flats Another potentially missing area was in Jiangsu Province.The tidal flat area in Yancheng derived using the proposed method (1249.25 km 2 ) is much larger than LUCC 2015, while LUCC 2015 captured 441.02km 2 tidal flats.That mostly resulted from the intertidal flats called the radial sand ridges (RSR)(Figure 12a), which are located in region marked as A in Figure 8; the kind of tidal flats could be found in some studies [27,58,59].But LUCC 2015 (Figure 12b) did not delineate the Radial Sand Ridges within the map extent.
The aforementioned cases (Figure 11, Figure 12) did not occur when our map was compared with Murray's, where a significantly larger extent of tidal flats was delineated in Murray's within all the administrative regions.And from Table 4, the identical area between our map and Murray's is relatively higher; that is, both of our and Murray's methods showed high consistency in mapping intertidal flats.But from Figures 11c and 12c, some extent of supratidal flats was included in Murray's results, which we took as separate type different from intertidal flats.Additionally, tidal flats that had been reclaimed were also mapped in Murray's extent, but we removed those through morphological post-processing.In conclusion, the mapping extent of our map is more specific and detailed.could be found in some studies [27,58,59].But LUCC 2015 (Figure 12b) did not delineate the Radial Sand Ridges within the map extent.
The aforementioned cases (Figure 11,12) did not occur when our map was compared with Murray's, where a significantly larger extent of tidal flats was delineated in Murray's within all the administrative regions.And from Table 4, the identical area between our map and Murray's is relatively higher; that is, both of our and Murray's methods showed high consistency in mapping intertidal flats.But from Figure 11c and Figure 12c, some extent of supratidal flats was included in Murray's results, which we took as separate type different from intertidal flats.Additionally, tidal flats that had been reclaimed were also mapped in Murray's extent, but we removed those through morphological post-processing.In conclusion, the mapping extent of our map is more specific and detailed.

Limitation and Potential of GEE Cloud Platform for Tidal Flats Mapping
Ideally, all the time-series images capture every tidal level, including the nominal lowest and highest tide during observation period.However, Landsat satellites, just like other sun-synchronous satellites, can only partially obtain the information of full tidal range, and extreme low or high tide could be observed less frequently [28,60].So, actually, the extent of tidal flats derived from the proposed approach here demonstrated the maximum exposed extent of tidal flats during the observation period.
In this way, the resultant map here still has some limitations in accuracy.However, it is still certain that the proposed approach through GEE has the capability to delineate a reasonable and objective extent of tidal flats.Although we chose the three-year temporal window as investigating period, the extent of tidal flats can be generated year by year.For instance, the annual extent of tidal flats in 2013 and 2014 could be derived from images acquired from 1 January 2012 to 31 December 2014 and 1 January 2013 to 31 December 2015, respectively.That is, the proposed method in this study can be utilized for change analysis of the annual extent of tidal flats.
GEE witnessed an upward trend in mapping land use and cover currently [16,[29][30][31]33,61].Technically, GEE is an integrated toolbox, providing a flood of satellite images and data products as well as professional spatial analysis in collaboration with advanced machine learning algorithm.Users can freely implement spatial analysis or image classification through the online code editor (https://code.earthengine.google.com/); that is, the proposed approach can be easily tuned to time, place and data inputs.

Limitation and Potential of GEE Cloud Platform for Tidal Flats Mapping
Ideally, all the time-series images capture every tidal level, including the nominal lowest and highest tide during observation period.However, Landsat satellites, just like other sun-synchronous satellites, can only partially obtain the information of full tidal range, and extreme low or high tide could be observed less frequently [28,60].So, actually, the extent of tidal flats derived from the proposed approach here demonstrated the maximum exposed extent of tidal flats during the observation period.
In this way, the resultant map here still has some limitations in accuracy.However, it is still certain that the proposed approach through GEE has the capability to delineate a reasonable and objective extent of tidal flats.Although we chose the three-year temporal window as investigating period, the extent of tidal flats can be generated year by year.For instance, the annual extent of tidal flats in 2013 and 2014 could be derived from images acquired from 1 January 2012 to 31 December 2014 and 1 January 2013 to 31 December 2015, respectively.That is, the proposed method in this study can be utilized for change analysis of the annual extent of tidal flats.
GEE witnessed an upward trend in mapping land use and cover currently [16,[29][30][31]33,61].Technically, GEE is an integrated toolbox, providing a flood of satellite images and data products as well as professional spatial analysis in collaboration with advanced machine learning algorithm.Users can freely implement spatial analysis or image classification through the online code editor (https://code.earthengine.google.com/); that is, the proposed approach can be easily tuned to time, place and data inputs.
There is one thing that should be noted: the amazing computational efficiency of GEE.In our approach, it took less than 20 minutes to process 1802 Landsat images, including preprocessing, multi-temporal feature extraction and image classification.This is obviously superior to the traditional remote-sensing image analysis process using desktop softwares.This provides a fascinating tool for large-scale and long-term analysis of tidal flats, as well as other applications.
In this study, we demonstrated the feasibility of mapping tidal flats by combining Landsat OLI imagery archive together with a random forest algorithm through GEE.Along with a growing number of satellite images with high spatial resolution or from advanced sensors being ingested, we envisage that more elaborate and up-to-date map for tidal flats within a larger extent would be generated through the method developed here.

Conclusions
In this study, we have proposed a novel approach for tidal flats mapping using time-series images, including time-series analysis, feature selection, machine learning classification and morphological post-processing.The contribution of this work could be summarized as follows: Firstly, we have analyzed the characteristics of different tidal flats in remote-sensing images, and the time-series spectral behaviors.And then, we proposed a method based on statistical attributes for feature selection, and extract the tidal flat information with the selected features.In terms of time-series analysis in statistical features, we proceeded to select image features which were successfully applied for accurate tidal mapping.
Secondly, we have proposed a novel morphological processing for tidal flat mapping.Since intertidal flats have significantly different features in time-series images, the results are more reliable.The morphological process is based on the fact that intertidal flats, supratidal barren flats and vegetated flats are spatially adjacent.Thus, the process started with intertidal flats, and iteratively merged other tidal patches.Finally, the post-process could significantly eliminate the misclassified patches, and improve the visual effect.
By comparing with the credible LUCC 2015 and Murray's methods, it is found that the proposed method can obtain a more accurate and visually better tidal flat map.
The GEE platform was used for data management and processing.GEE is amazing in its capability, convenience and efficiency in time-series image analysis, especially for national and global applications.It showed significant advantages in tidal flats mapping as well as other geographical and environmental applications.We envisage that it will be an important tool in earth science and environmental science.

20 Figure 1 .
Figure 1.Sketch-map of tidal flat classification used in this study: intertidal flats, supratidal barren flats and vegetated flats.

Figure 1 .
Figure 1.Sketch-map of tidal flat classification used in this study: intertidal flats, supratidal barren flats and vegetated flats.

Figure 2 .
Figure 2. Location of study area, a 50 km coastal buffer for mapping tidal flats and tiles of Worldwide Reference System 2 (WRS-2) path/row covering the buffer area from Hangzhou Bay to Yalu River in China's eastern coastal zone.Points of interest (POIs, including training and validation sample points) were marked by colored points within this figure.

Figure 2 .
Figure 2. Location of study area, a 50 km coastal buffer for mapping tidal flats and tiles of Worldwide Reference System 2 (WRS-2) path/row covering the buffer area from Hangzhou Bay to Yalu River in China's eastern coastal zone.Points of interest (POIs, including training and validation sample points) were marked by colored points within this figure.

Figure 3 .
Figure 3. Workflow used in this study for tidal flats mapping.

Figure 3 .
Figure 3. Workflow used in this study for tidal flats mapping.

Figure 4 .
Figure 4. Satellite images interpretation of (a) intertidal flat, (b) supratidal barren flat and (c) vegetated flat as well as (d) offshore ponds shown in true color Landsat Operational Land Imager (OLI) images.Field photos taken in (e) Shanghai City and (f) Yancheng City, Jiangsu Province, which are corresponding to the sites of (a) intertidal flat and (c) vegetated flat, respectively.

Figure 4 .
Figure 4. Satellite images interpretation of (a) intertidal flat, (b) supratidal barren flat and (c) vegetated flat as well as (d) offshore ponds shown in true color Landsat Operational Land Imager (OLI) images.Field photos taken in (e) Shanghai City and (f) Yancheng City, Jiangsu Province, which are corresponding to the sites of (a) intertidal flat and (c) vegetated flat, respectively.
marks the intertidal (I), supratidal barren (II) and vegetated (III) flat within the false color Landsat 8 OLI image on 4 May 2015, and Figure5c-e present their time-series variation of spectral indices, respectively.For those areas, the variability of spectral behaviors has been recorded during a 3-year period.

Figure 5 .
Figure 5. Temporal profile of Enhanced Vegetation Index (EVI), Land Surface Water Index (LSWI), Modified Soil-Adjusted Vegetation Index (mSAVI), Normalized Difference Vegetation Index (NDVI), Modified Normalized Difference Water Index (mNDWI) and soil brightness derived from Landsat OLI imagery archive at the selected tidal flat sites (b) of Yellow River Delta (a); (c-e) are corresponding to intertidal flat (I), supratidal barren flat (II) and vegetated flat (III), respectively.

Figure 6 .
Figure 6.Boxplot of six land use classification knowledge captured to train random forest machine learning algorithm in (a) NDVI minimum and (b) LSWI maximum for POIs.

Figure 6 .
Figure 6.Boxplot of six land use classification knowledge captured to train random forest machine learning algorithm in (a) NDVI minimum and (b) LSWI maximum for POIs.2.3.3.Random Forest (RF) Machine Learning Algorithm via Google Earth Engine (GEE) About 80% of reference sample points were randomly selected as training dataset for RF machine learning algorithm through GEE, which consisted of 488 tidal flats samples (288 in intertidal flats, 131 in supratidal barren flats and 69 in vegetated flats) and 1145 non-tidal flats samples (Section 2.3.1).RF classifier consists of a number of decision trees, involving voting mechanism for improving accuracy[54], with high cost in computation when the number of trees is set very large.As the composite image with 48 bands was used as target image, the number of decision trees was limited to 100, keeping a trade-off between accuracy and efficiency.

Figure 7 .
Figure 7. Example of morphological process for tidal flats mapping: (a) is the preliminary result from GEE in collaboration with random forest algorithm, (b) manifests intertidal flat patches after using the kernel of eight-connected neighbors, (c) is the result filtered by minimum pixels count, and (d) is final map generated from the iterative eight-connected kernel and manually revised.

Figure 7 .
Figure 7. Example of morphological process for tidal flats mapping: (a) is the preliminary result from GEE in collaboration with random forest algorithm, (b) manifests intertidal flat patches after using the kernel of eight-connected neighbors, (c) is the result filtered by minimum pixels count, and (d) is final map generated from the iterative eight-connected kernel and manually revised.

Figure 8 .
Figure 8. Tidal flat map generated by proposed approach in this study.

Figure 8 .
Figure 8. Tidal flat map generated by proposed approach in this study.

Figure 9 .
Figure 9.The relationship between the overall accuracy and (a) different temporal windows as well as (b) the proportion of reference sample points for training, respectively.

Figure 9 .
Figure 9.The relationship between the overall accuracy and (a) different temporal windows as well as (b) the proportion of reference sample points for training, respectively.

Figure 10 .
Figure 10.Scatter plot and regression line between our tidal flat map and (a) LUCC 2015 as well as (b) Murray's in city level, respectively.

Figure 11 .
Figure 11.Detection of tidal flats in Liaohe Estuary within (a) our map, (b) LUCC 2015 and (c) Murray's, respectively.The base map from Landsat OLI acquired on 9 August 2018 is displayed in true color.

Figure 10 .
Figure 10.Scatter plot and regression line between our tidal flat map and (a) LUCC 2015 as well as (b) Murray's in city level, respectively.

20 Figure 10 .
Figure 10.Scatter plot and regression line between our tidal flat map and (a) LUCC 2015 as well as (b) Murray's in city level, respectively.

Figure 11 .
Figure 11.Detection of tidal flats in Liaohe Estuary within (a) our map, (b) LUCC 2015 and (c) Murray's, respectively.The base map from Landsat OLI acquired on 9 August 2018 is displayed in true color.

Figure 11 .
Figure 11.Detection of tidal flats in Liaohe Estuary within (a) our map, (b) LUCC 2015 and (c) Murray's, respectively.The base map from Landsat OLI acquired on 9 August 2018 is displayed in true color.

Figure 12 .
Figure 12.Detection of radial sand ridges (RSR) in Jiangsu Province within (a) our map, (b) LUCC 2015 and (c) Murray's, respectively.The base map from Landsat OLI acquired on 1 January 2016 is displayed in true color.

Figure 12 .
Figure 12.Detection of radial sand ridges (RSR) in Jiangsu Province within (a) our map, (b) LUCC 2015 and (c) Murray's, respectively.The base map from Landsat OLI acquired on 1 January 2016 is displayed in true color.

Table 1 .
Classification system used in this study and their corresponding sample numbers.

Table 1 .
Classification system used in this study and their corresponding sample numbers.

Table 2 .
Confusion matrix for assessing the performance of proposed approach here for mapping tidal flats.

Table 3 .
Statistics table of three categories among tidal flats in terms of administrative divisions.

Table 4 .
Statistical table for tidal flats area of administrative regions in the study area, derived from remote-sensing monitoring data of land use in China (LUCC 2015), Murray's and ours, respectively.Numbers in round brackets are the proportion of identical area to the intertidal flats area derived from our map."Proposed" means the resultant tidal flat map in this study.