Gap-Free Monitoring of Annual Mangrove Forest Dynamics in Ca Mau Province, Vietnamese Mekong Delta, Using the Landsat-7-8 Archives and Post-Classiﬁcation Temporal Optimization

: Ecosystem services o ﬀ ered by mangrove forests are facing severe risks, particularly through land use change driven by human development. Remote sensing has become a primary instrument to monitor the land use dynamics surrounding mangrove ecosystems. Where studies formerly relied on bi-temporal assessments of change, the practical limitations concerning data-availability and processing power are slowly disappearing with the onset of high-performance computing (HPC) and cloud-computing services, such as in the Google Earth Engine (GEE). This paper combines the capabilities of GEE, including its entire Landsat-7 and Landsat-8 archives and state-of-the-art classiﬁcation approaches, with a post-classiﬁcation temporal analysis to optimize land use classiﬁcation results into gap-free and consistent information. The results demonstrate its application and value to uncover the spatio-temporal dynamics of mangrove forests and land use changes in Ngoc Hien District, Ca Mau province, Vietnamese Mekong delta. The combination of repeated GEE classiﬁcation output and post-classiﬁcation optimization provides valid spatial classiﬁcation (94–96% accuracy) and temporal interpolation (87–92% accuracy). The ﬁndings reveal that the net change of mangroves forests over the 2001–2019 period equals − 0.01% annually. The annual gap-free maps enable spatial identiﬁcation of hotspots of mangrove forest changes, including deforestation and degradation. Post-classiﬁcation temporal optimization allows for an exploitation of temporal patterns to synthesize and enhance independent classiﬁcations towards more robust gap-free spatial maps that are temporally consistent with logical land use transitions. The study contributes to a growing body of work advocating full exploitation of temporal information in optimizing land cover classiﬁcation and demonstrates its use for mangrove forest monitoring.


Introduction
Ecosystem services offered by mangrove forests are facing severe risks. Within the transition zone of land and sea of (sub)tropical coastal regions, mangroves have carved out a distinct niche to flourish and thereby provide vital services to mankind. Specifically, mangrove ecosystems have shown to be one of the world's most productive in terms of carbon sequestration, shelter and breeding grounds for aquatic species, and as important physical barriers against tides and ocean surges [1][2][3][4]. Despite the multitude of crucial ecosystem services these coastal forests offer to communities in coastal regions of more than 124 countries, the status of mangrove forests in many regions is under pressure due to forest loss and land degradation, caused by overexploitation and land use change driven by human development [5][6][7]. Due to the inaccessible, ever-changing, and extensive nature of these mangroves, remote sensing has become a primary instrument to monitor the health and dynamics of these ecosystems [8][9][10].
The field of Satellite Remote Sensing has moved into an era in which a tremendous wealth of earth observation (EO) data are gathered at increasing spectral, spatial, and temporal resolutions-supporting the wide-spread application of satellite data for studying global changes [11]. Orbiting EO satellites allow us to repeatedly revisit areas of interest to study temporal changes and facilitate time series analysis. The iconic Landsat-7 and Landsat-8 missions both offer average revisit intervals of 16 days and observations that go back as early as the year 2000. The later Landsat-8 mission collected over 3. 35 Petabyte of scenes over the course of a single year in 2017 [12]. These data collections hold great potential to improve our monitoring efforts of mangrove ecosystems and study changes over time.
A critical review by Younes Cardenas et al. (2017) on using satellite remote sensing to monitor mangrove ecosystems points out that the majority of studies conducted-reviewing 55 recent peer-reviewed articles using Landsat/Aster imagery-are not making full use of the wealth of EO data available [13]. The authors specify that most studies between 2001-2016 used fewer than 10 images and longitudinal studies often analyze temporal changes with 7-11 years between scenes which leaves much of the potential of current satellite archives unlocked [13]. Yet, mangrove forests are frequently part of fast-changing landscapes driven by land use change at the interplay of volatile aquaculture markets, policy-making, and the biophysical dynamics of erosion, sedimentation, and changing tides [14][15][16]. This raises the question of how we can better unlock the potential of available satellite imagery archives to facilitate high temporal resolution monitoring of the fast-paced land use processes surrounding mangrove forest ecosystems.
The advances in high-performance computing (HPC) in combination with cloud-computing services, such as provided by the Google Earth Engine platform (GEE), allow us to address the major challenges of processing and handling enormous EO datasets and turning these into comprehensible information [13,[17][18][19][20][21]. The GEE platform provides straightforward HPC cloud access to many of the major satellite archives as well as numerous image classifiers for mapping applications, including Classification and Regression Trees (CART) and Random Forests (RF) approaches. Illustrative of its capabilities, Hansen et al. (2013) mapped global forest cover change products from over 650 thousand Landsat-7 scenes [22]. Following this, a large body of regional studies has demonstrated high mapping accuracies using GEE's land use classifiers (CART) with Landsat images [19,23,24]. More specifically, we observe an increasing use and successful implementation (accuracies between 92% and 97%) of GEE-based land use classification for mangrove mapping [25][26][27].
Through GEE, we can efficiently organize longitudinal time series from satellite observations and independent classification efforts can be repeated over time with ease. These time series can be valuable to study and monitor temporal changes in land use. Conversely, the temporal dependencies of each time point in the series can also be used to further optimize the time series in terms of missing information and consistency. In other words, the temporal domain can facilitate post-classification optimization of GEE output towards maps that are gap-free and temporally consistent with logical land use transitions as well as provide a means of cross-validation. This is particularly crucial in cases Remote Sens. 2020, 12, 3729 3 of 16 with hampering climatic conditions (clouds, snow, dust, and aerosols), instrumentation errors, losses of image data during data transmission, or high uncertainties in information processing [28].
Temporal gap-filling and smoothing approaches are common practice in remote sensing of phenology and cropping cycles through continuous parameters, such as vegetation indices (e.g., NDVI, EVI) and surface parameters (Land Surface Temperature) [28][29][30][31][32]. However, in discrete land cover classification exercises, this practice remains less common, including in combination with the GEE platform [33,34]. Current studies tend to focus on gap-filling based on spatially neighboring pixels [35,36], spectral similarity, and/or multi-sensor (source) data fusion [34,37,38], rather than temporal integration. As such, few land use studies have taken full advantage of temporal dependencies to reduce both information gaps and inconsistent land use transitions [13,[39][40][41]. This is a particularly rare undertaking for the monitoring of mangrove forests land use changes, whereas consistent and gap-free time series are crucial to closely monitor mangrove deforestation, degradation, and disturbance [13,15]. Land use changes tend to follow logical temporal land use transitions which can guide the optimization of classification maps [13,40].
The main objective of our study is to deploy high-performance computing techniques to monitor mangrove forest cover changes in our case study area; the mangrove-rich Ngoc Hien District, Ca Mau province in the Vietnamese Mekong delta. Rather than a single land use classification approach, we demonstrate how independent land use classifications conducted in GEE can be combined to optimize classification results in terms of completeness and consistency. As such, the study exploits both; (1) the computational capacity of GEE to deal with the entire Landsat-7 and -8 archives and (2) the temporal element of a longitudinal time series to optimize land use classification results into "gap-free" and temporally consistent information. This can help us better understand the spatio-temporal dynamics of mangrove forests, in terms of extent, distribution, and land use change and disturbances that threaten their conservation.

Study Area
The study area focuses on Vietnam's southernmost district, Ngoc Hien, Ca Mau province, located in the Southern Mekong Delta between latitude 8 • 33 -8 • 45 N and longitude 104 • 42 45"-105 • 3 54"E, spanning an area of 743 km 2 (See Figure 1). The district has been well-studied for its importance as a major aquaculture hub and its significant reserves of Vietnam's largest and last remaining old-growth mangrove forests, including the internationally acknowledged RAMSAR site of Mui Ca Mau (2012) and UNESCO Biosphere Reserve (2009) [42][43][44]. The landscape supports both ecologically important mangrove ecosystems and the economic livelihoods based on aquaculture.

Remote Sensing Data Pre-Processing
This study makes use of the archives of Landsat's later missions embodied by the Landsat-7 and Landsat-8 multispectral imagery available through GEE's public data catalogue of atmospherically corrected surface reflectance data. We have made use of all available 30 m spatial resolution bands of both missions, this implies: two short-wave infrared (SWIR) bands and four/five visible and near-infrared (VNIR) bands for Landsat-7 and Landsat-8, respectively. The study area is centered within the path-rows of 125-054 and 126-064 with an average 16-day revisit time. The Landsat-7 and -8 Quality Assessment bands and calculated F-mask were used to filter out pixels containing clouds, cirrus, cloud shadow, and atmospheric contamination of the reflectance signal. Table 1 gives an overview of the available images for each year. Based on the spectral reflectance of all available images annual median composites are compiled that provide consistent cloud-free median images of spectral reflectance [45]. The malfunction of the Scan Line Corrector (SLC) of the Landsat-7 imager has resulted in that on average around a quarter of the data in a Landsat-7 scene is missing from the 31st of May 2003 onwards [46]. These products hence have considerable data gaps, Remote Sens. 2020, 12, 3729 4 of 16 but still maintain the same radiometric and geometric corrections as data collected prior to the SLC failure. The combination of high cloudiness of the region with SLC failure results in limited usability of the years 2003-2012. Therefore, our analysis deploys a cautionary interpretation of these years. Furthermore, years characterized by SLC failure are combined bi-annually to increase data availability, coverage of the region, and to lower uncertainties.

Study Area
The study area focuses on Vietnam's southernmost district, Ngoc Hien, Ca Mau province, located in the Southern Mekong Delta between latitude 8°33′-8°45′N and longitude 104°42′45″-105°3′54″E, spanning an area of 743 km 2 (See Figure 1). The district has been well-studied for its importance as a major aquaculture hub and its significant reserves of Vietnam's largest and last remaining old-growth mangrove forests, including the internationally acknowledged RAMSAR site of Mui Ca Mau (2012) and UNESCO Biosphere Reserve (2009) [42][43][44]. The landscape supports both ecologically important mangrove ecosystems and the economic livelihoods based on aquaculture.

Land Cover Classification
After pre-processing, the resulting cloud-free median multispectral annual composites are used to characterize land use, and the land use changes over time. The land use classification scheme of our study takes into consideration four dominant land uses within the Ngoc Hien district, namely (1) Dense Mangrove Forest, (2) Sparse Mangroves, (3) Aquaculture/Waterbodies, and (4) Built-up and Barren lands. Dense mangrove forest is defined by a minimum of 30% canopy cover. Vegetated mangrove areas that are 10-30% crown cover are classified as sparse mangroves.
We conducted a supervised classification to develop land use maps. There are several classification algorithms available within GEE, including; Classification and Regression Trees (CART), Random Forest (RF), Naïve Bayes, and Support-Vector Machine (SVM). Our study opted for the commonly used CART classifier which has produced relatively high accuracies when applied to Landsat data in numerous settings [19,23,24,26]. More specifically, several studies have reported the highest accuracy for CART land use classification of coastal wetlands and mangroves using GEE compared to other classifiers [25,27]. Most importantly, we ran trails in the study area for both CART and RF in which the first yielded the highest classification accuracy (94-96% for CART, against 89-94% for RF, respectively). GEE code implementations of both approaches and its validation against test data can be found in the Supplementary Materials Table S1.
Within CART, a decision tree (DT) classifier was instantiated and trained on field data using GEE's default parameters. The CART algorithm runs through a series of nodes that recursively split the input data in such a way that there is a decrease in entropy and an increase in information gain after the split. GEE's CART uses the Gini Impurity Index to decide the input features which will provide the best split at each node. A tabular overview of the exact decision rules for building the model can be found in Supplementary Materials Table S2. One disadvantage of the DT classifier is the considerable sensitivity to the training dataset. A small change to the training data can result in a very different set of subsets and can result in overfitting [19,47]. Nevertheless, our training and validation data relies on extensive fieldwork, including 514 georeferenced points gathered in-situ across the Ngoc Hien district in 2015, subdivided into the four classes; dense mangroves (n = 247), sparse mangroves (n = 72), waterbodies/aquaculture ponds (n = 120), and built-up and barren lands (n = 75). We used 70% of the field data for training and the remaining 30% for validation, thereby estimating the classification errors independently.
Following the initial training of the classifier, it is then deployed backward (LS-7) or forward (LS-8) through the time series based on spectral/change information of the surface reflectance data available in the composite datasets. Based on this method, land cover maps are generated from the surface reflectance of pre-processed yearly median composites between 2001 and 2019. The workflow of GEE pre-processing and land use classification is presented schematically in Figure 2. GEE code can be accessed through the URLs published in Supplementary Materials Table S1.

Post-Classification Optimization through Time Series Temporal Data Fusion
The longitudinal temporal data of the Landsat archives enabled the use of neighboring time points to cross-validate findings, fill in missing data through temporal data fusion, detect and revise illogical land use changes in the post-classification analysis [33,[39][40][41]46,48].
Gap-filling through consideration of the temporal dimension in years preceding and following a missing value has allowed us to interpolate missing pixels to assure spatial and temporal continuity over the study area. The approach for temporal interpolation follows an adaptation of inverse distance weighting methods applied to discrete land use classification data, the basic assumption is that a temporally distributed variable at short distance is generally more similar than at larger Remote Sens. 2020, 12, 3729 6 of 16 distance [49,50]. The applied approach integrates land use classes of adjacent years weighted by a power (p = 1.5) of the distance (d) to the year of interpolation, formulated in an equation as: in which i = 0 indicates the time point to be interpolated with predicted land use (ẑ), and which i = n indicate the n years adjacent where land use (z) has been observed. As such, pixels missing valid observations are estimated by taking into account a seven-year pattern and scoring neighboring time points, in which observations nearest in time weight the heaviest. The seven-year time window corresponds with the consecutive years with full data availability (2013-2019). The land use class scored the highest will be used for gap-filling.
Within CART, a decision tree (DT) classifier was instantiated and trained on field data using GEE's default parameters. The CART algorithm runs through a series of nodes that recursively split the input data in such a way that there is a decrease in entropy and an increase in information gain after the split. GEE's CART uses the Gini Impurity Index to decide the input features which will provide the best split at each node. A tabular overview of the exact decision rules for building the model can be found in Suppl. Mater. S2. One disadvantage of the DT classifier is the considerable sensitivity to the training dataset. A small change to the training data can result in a very different set of subsets and can result in overfitting [19,47]. Nevertheless, our training and validation data relies on extensive fieldwork, including 514 georeferenced points gathered in-situ across the Ngoc Hien district in 2015, subdivided into the four classes; dense mangroves (n = 247), sparse mangroves (n = 72), waterbodies/aquaculture ponds (n = 120), and built-up and barren lands (n = 75). We used 70% of the field data for training and the remaining 30% for validation, thereby estimating the classification errors independently.
Following the initial training of the classifier, it is then deployed backward (LS-7) or forward (LS-8) through the time series based on spectral/change information of the surface reflectance data available in the composite datasets. Based on this method, land cover maps are generated from the surface reflectance of pre-processed yearly median composites between 2001 and 2019. The workflow of GEE pre-processing and land use classification is presented schematically in Figure 2. GEE code can be accessed through the URLs published in Suppl. Mater. S1.  Similarly, further optimization of classification results can be achieved by taking into account that land use changes usually occur characterized by a logical transition [40,41]. Land use changes and transitions follow ecological rules [40,41]. For instance, the growth of dense mangrove forests takes at least multiple years. Understanding these land use transitions can help setting rules determined by ecology and feasibility to detect illogical land use transition from remotely sensed time series.
In our study, we have opted for a post-classification approach to detect and revise illogical land use (See Table 2). The assessment of uncertainties and reclassification by the CART classifier based on the temporal exclusion of certain illogical transition rules (e.g., maximum a posteriori (MAP) classification rule in a Bayesian framework [41]) has not been embedded within the GEE platform. Therefore, a post-classification approach that scores the likelihood of land use classes to occur based on temporal context is deployed. Table 2 provides an overview of the anomalies detected. Similar to the gap-filling of missing pixels, illogical land use transitions are revised by distance weighting of neighboring time points (analyzing a seven-year pattern) of the pixel under scrutiny to determine what land cover time is most probable to be in place. The entire workflow from pre-processing, land cover classification, post-classification optimization is presented schematically in Figure 2's diagram.

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accurac confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-fillin algorithm's performance to predict land use type based on the seven-year pattern surroundin missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran withi separately for both Landsat-7, and -8 on collected field data, produces 95.6% confiners for LS-7 and LS-8, respectively, for the reference year 2015 (Table algorithm's performance to predict land use type based on the seven-year p missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accurac confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-fillin algorithm's performance to predict land use type based on the seven-year pattern surroundin missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran withi separately for both Landsat-7, and -8 on collected field data, produces 95.6% confiners for LS-7 and LS-8, respectively, for the reference year 2015 (Table algorithm's performance to predict land use type based on the seven-year p missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accurac confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-fillin algorithm's performance to predict land use type based on the seven-year pattern surroundin missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran withi separately for both Landsat-7, and -8 on collected field data, produces 95.6% confiners for LS-7 and LS-8, respectively, for the reference year 2015 (Table algorithm's performance to predict land use type based on the seven-year p missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).  We tested the gap-filling algorithm's predictive performance using a k-fold leave-one-out strategy. In these tests, we purposely removed a datapoint to run the gap-filling algorithm to predict the land use based on temporal neighbors. We assessed the percentage of correctly predicted land use classifications. We have conducted this for all available pixels and years in the analysis.
Taken together, temporal integration of GEE's individual land use classification maps into a gap-free sequence of maps that follow logical land use transitions enabled us to analyze the land use trends occurring throughout the available data. Land use changes were grouped into four categories illustrated through the color scheme in Table 2. Only mangrove forest changes were considered in the trend maps (other land use changes remain white). Moreover, we used 150 × 150 m box averages to highlight general trends and reduce uncertainties by spatial regularization.

Temporal and Spatial Accuracy Assessment
Classification of surface reflectance, using a CART classifier ran within GEE and trained separately for both Landsat-7, and -8 on collected field data, produces 95.6% and 94.1% accuracy confiners for LS-7 and LS-8, respectively, for the reference year 2015 ( Table 3). The gap-filling algorithm's performance to predict land use type based on the seven-year pattern surrounding missing pixels hovers between 87 to 92% depending on the year (Table 4).  An overview of all land use classification maps after gap-filling and post-processing for logical land use transitions is presented in Figure 3. The maps depict all years with full data-availability and therefore the highest reliability ranging from 2001-2002, and 2013-2019. A visual inspection of the temporal dynamics indicates that the central regions in Ngoc Hien have been subject to a high frequency of land use changes whereas towards the coast, specifically the Western Cape, stable havens of dense mangrove forest have to large extent been minimally subjected to forest clearance and land use change. Over the years, three core mangrove forests areas can be identified (See Figure 1, for commune locations); (1) the cluster in the south-east of Ngoc Hien, including the eastern part of Tan An commune and adjacent southern Tam Giang Tay commune, (2) the mangrove forest strip along the southern coastline, and (3) the largest core mangrove area is located in the western tips of Ngoc Hien district; Dat Mui commune which encloses the Ca Mau Cape National Park and the north-western part of Vien An commune.
The total land use changes over time are highlighted in Figure 4. These results show that the total shares of land use cover between four classes have remained relatively stable with yearly fluctuations. In the span of two decades, we find sparse mangrove cover fluctuating between 33,000 ha and 40,000 ha and dense mangrove forests hovering between 18,000 ha and 24,000 ha. The average net change of dense mangroves over the 2001-2019 equals −0.01% annually. In recent years since 2013, we observe upward trends in dense mangrove forest cover as well as built environment and mudflats. Sparse mangrove cover has been decreasing in that same period. Waterbodies/ponds have remained more or less equal.
Despite these overall trends, the figure conceals the spatial dynamics behind the distribution of mangrove forests and land use changes that characterize the region. To gain insight into the spatial distribution of land use changes, an aggregated trend map was created ( Figure 5) to identify hotspots of mangrove forest change while leaving changes of non-forest classes out of consideration. Generally, the dense mangrove patches in the Eastern regions have shrunk in size and extent, whereas the core mangrove forests on the Western Cape have seen an expansion through sedimentation and seaward colonization. Ngoc Hien's central communes have seen an increased integration of dense mangroves patches in the mosaic of waterbodies/ponds and sparse mangrove cover. Areas most prone to mangrove forest loss and degradation are the eastern central regions, possibly due to its vicinity to districts with higher urban development, infrastructure, and aquaculture production [51].
The trend map presented in Figure 5 gives insight into the spatial distribution of temporal dynamics of mangrove forest loss and gains. The map highlights how coastal erosion along the southward coastline and sedimentation along the westward frontier have resulted in losses and gains of mangrove forests. The strip of dense mangrove forests serving as a protective buffer along Ngoc Hien's southern coastline has seen shifting land use changes; on the one hand, inland regrowth and increased connectivity are observed, whereas on the other hand coastal erosion is becoming increasingly severe leading to seaward losses of dense mangrove forest. In Ngoc Hien's Western regions, Dat Mui and Vien An communes (see Figure 1, for commune locations), the remaining core mangrove forests have seen relatively few land use changes except for seaward mangrove expansion and colonization caused by coastal sedimentation (Figure 5).
land use transitions is presented in Figure 3. The maps depict all years with full data-availability and therefore the highest reliability ranging from 2001-2002, and 2013-2019. A visual inspection of the temporal dynamics indicates that the central regions in Ngoc Hien have been subject to a high frequency of land use changes whereas towards the coast, specifically the Western Cape, stable havens of dense mangrove forest have to large extent been minimally subjected to forest clearance and land use change. Over the years, three core mangrove forests areas can be identified (See Figure  1, for commune locations); (1) the cluster in the south-east of Ngoc Hien, including the eastern part of Tan An commune and adjacent southern Tam Giang Tay commune, (2) the mangrove forest strip along the southern coastline, and (3) the largest core mangrove area is located in the western tips of Ngoc Hien district; Dat Mui commune which encloses the Ca Mau Cape National Park and the northwestern part of Vien An commune.  The total land use changes over time are highlighted in Figure 4. These results show that the total shares of land use cover between four classes have remained relatively stable with yearly fluctuations. In the span of two decades, we find sparse mangrove cover fluctuating between 33,000 ha and 40,000 ha and dense mangrove forests hovering between 18,000 ha and 24,000 ha. The average net change of dense mangroves over the 2001-2019 equals −0.01% annually. In recent years since 2013, we observe upward trends in dense mangrove forest cover as well as built environment and mudflats. Sparse mangrove cover has been decreasing in that same period. Waterbodies/ponds have remained more or less equal. Despite these overall trends, the figure conceals the spatial dynamics behind the distribution of mangrove forests and land use changes that characterize the region. To gain insight into the spatial distribution of land use changes, an aggregated trend map was created ( Figure 5) to identify hotspots of mangrove forest change while leaving changes of non-forest classes out of consideration. Generally, the dense mangrove patches in the Eastern regions have shrunk in size and extent, whereas the core mangrove forests on the Western Cape have seen an expansion through sedimentation and seaward colonization. Ngoc Hien's central communes have seen an increased integration of dense mangroves patches in the mosaic of waterbodies/ponds and sparse mangrove cover. Areas most prone to mangrove forest loss and degradation are the eastern central regions, possibly due to its vicinity to districts with higher urban development, infrastructure, and aquaculture production [51].
The trend map presented in Figure 5 gives insight into the spatial distribution of temporal dynamics of mangrove forest loss and gains. The map highlights how coastal erosion along the southward coastline and sedimentation along the westward frontier have resulted in losses and gains of mangrove forests. The strip of dense mangrove forests serving as a protective buffer along Ngoc Hien's southern coastline has seen shifting land use changes; on the one hand, inland regrowth and increased connectivity are observed, whereas on the other hand coastal erosion is becoming increasingly severe leading to seaward losses of dense mangrove forest. In Ngoc Hien's Western regions, Dat Mui and Vien An communes (see Figure 1, for commune locations), the remaining core mangrove forests have seen relatively few land use changes except for seaward mangrove expansion and colonization caused by coastal sedimentation (Figure 5).

Overcoming Observation Gaps in Mangrove Monitoring
In the critical review by Younes Cárdenas et al. (2017) based on 55 recent peer-reviewed journal articles using Landsat or ASTER images to monitor mangrove forests, the authors conclude that the majority of multi-temporal studies focus on only a fraction of available satellite imagery with on average 7-11 years between scenes of multitemporal analysis [13]. Yet, high temporal change detection of mangrove forests is vital in mapping threats from aquaculture expansion and coastal development as well as to understand cyclic processes such as logging in production areas and seasonal biomass fluctuations [15]. In our study, 446 unique scenes were processed into annual median composites to study the temporal dynamics of mangrove forests in one of Vietnam's most prominent regions for mangrove conservation, Ngoc Hien district in the Ca Mau province. Despite a large number of scenes per composite (15.4 scenes on average), it still resulted in a high number of pixels with no observations and/or illogical land use transitions for different years between 2001-2019 (See Table 1). Following temporal gap-filling and post-classification optimization, the resulting optimized time series allowed us to better monitor the state of mangrove forests in Ngoc Hien with spatial and temporal continuity and logical consistency in transitions.

Overcoming Observation Gaps in Mangrove Monitoring
In the critical review by Younes Cárdenas et al. (2017) based on 55 recent peer-reviewed journal articles using Landsat or ASTER images to monitor mangrove forests, the authors conclude that the majority of multi-temporal studies focus on only a fraction of available satellite imagery with on average 7-11 years between scenes of multitemporal analysis [13]. Yet, high temporal change detection of mangrove forests is vital in mapping threats from aquaculture expansion and coastal development as well as to understand cyclic processes such as logging in production areas and seasonal biomass fluctuations [15]. In our study, 446 unique scenes were processed into annual median composites to study the temporal dynamics of mangrove forests in one of Vietnam's most prominent regions for mangrove conservation, Ngoc Hien district in the Ca Mau province. Despite a large number of scenes per composite (15.4 scenes on average), it still resulted in a high number of pixels with no observations and/or illogical land use transitions for different years between 2001-2019 (See Table 1). Following temporal gap-filling and post-classification optimization, the resulting optimized time series allowed us to better monitor the state of mangrove forests in Ngoc Hien with spatial and temporal continuity and logical consistency in transitions.

Post-Classification Temporal Optimization
Accurate land use classification remains a challenging undertaking in landscapes dominated by aquaculture and mangrove forest land use [13]. In our study, the accuracy assessment for the reference year 2015 yielded assuring accuracy confiners (Table 3). We used 514 single-year reference data points observed in-situ for training the classifier and validation of the results. Ideally, we would have such data available for multiple years; however, this is difficult to acquire and organize. With limited availability of ground-truth data, reducing classification uncertainties and increasing temporal consistency is key to provide high-quality land use maps.
The use of yearly median composites allowed for fast computational processing and a comprehensible annual time series output for policy and decision-makers. Nevertheless, the application of median composites also poses disadvantages. Adequate composites still require the presence of sufficient high-quality observations. Moreover, a year-round even temporal distribution of scenes is required to facilitate an adequate median composite that is representative for the entire year. Knowledge regarding yearly phenology, the impact of tides of reflectance signals, trends in biophysical parameters (functional traits) of mangrove ecosystems is for a large part still lacking to make appropriate judgments on possible biases in median composites [13]. In other words, gaining more understanding regarding these temporal processes, also within yearly cycles, will help gain insight into the robustness of median composites for mangrove forest ecosystems.
Further challenges in accurate classifications of mangrove forests are raised by; (1) the fine-grained landscape mosaic with mangrove plots and aquaculture ponds often sized at sub-pixel (30 m × 30 m) measures, (2) the unknown implications of tidal effects on spectral signals and the high level of water vapor observed in these coastal regions, and (3) recent trends towards integrated mangrove-shrimp farming production systems which have made discrimination between mangrove, aquaculture paddy land use classes more ambivalent [9,13,15,52,53]. These challenges highlight the importance of making the most of the temporal information available to lower uncertainties in the final classification product. This is particularly important when noise in remote sensing signals is high-which is commonly the case in cloud-covered mangrove regions-and when multi-annual validation/training data are scarce. Several situations can cause illogical changes in a land use time series, such as classification errors, reflectance signal noise, and imperfect image co-registration.
Here, we demonstrated how the exploitation of information in the temporal domain allowed for additional optimization and a means of cross-validation of the GEE classification outputs. Studying temporal patterns and cross-validating land use changes in relation to the previous and following years help increase the robustness to noise and credibility of classification efforts. Specifically, the temporal information of land cover maps in a time series has been used to detect illogical land use changes and improve classification results [48]. The approach for temporal mangrove monitoring outlined is relatively easy to implement using GEE output and post-classification optimization. At the same time, it provides valid spatial classification (Table 3; 94-96% accuracy) and temporal interpolation (Table 4; 87 to 92% accuracy). Temporal interpolation follows a simple discrete inverse distance weighting algorithm, however other and more advanced statistical learning approaches can potentially be interesting alternatives [41,49,50,[54][55][56]. Figure 5 illustrates the further implementation of our gap-free time series in studying multi-temporal mangrove land use trends in Ngoc Hien. Bi-temporal approaches risk highlighting observations that result from isolated instances that can introduce inconsistencies and uncertainties in classification, especially considering the, at times, unfavorable signal-to-noise ratios found in satellite remote sensing [13]. Instead of comparing two single timestamps, the integration of yearly gap-free land use classification enables temporal cross-validation along logical land use transitions and gap-filling based on temporally neighboring information. This temporally dense analysis allows us to fully assess the direction and frequency of land use changes affecting mangrove forests [57,58]. This is also important to ensure that the unchanged forest between two time points has remained undisturbed in the years in between. Moreover, the multi-annual approach allows us to assess and quantify observed land use changes, e.g., forest disturbances, at a high temporal frequency, thereby opening venues to better monitor and study mangrove forest disturbance regimes and mangrove degradation processes as compared to bi-temporal land use change.
The land use trend map ( Figure 5) facilitates the detection of hotspots for mangrove forest change; deforestation, degradation, and regrowth. While the map is informative for deforestation and the land use change drivers behind it, the assessment of forest degradation remains arbitrary inherent to the operationalization of the land use classification scheme. This ties in with the challenges and complexity of defining forest degradation [59]. A large variety of existing definitions of forest degradation require different methods for assessment based on the objectives of the intended study [60]. The gap-free land use change maps presented here may help flag changes in mangrove forest extent. However, to overcome the arbitrariness of classes, complimentary maps on quantitative canopy/leaf traits and biomass can further enhance our ability to assess forest degradation and forest regrowth along a spectrum of ecologically relevant indicators [61][62][63].

Future Implementation
The case study presented here builds on and contributes to a growing body of work advocating that approaches that account for temporal information in optimizing land cover classification are superior to temporally independent time series classifications and other single-date methods [39,64]. The application of the post-classification optimization to stabilize land cover trajectories, mitigate unrealistic land cover transitions and overcome the limitations and costs of obtaining ancillary and field data are still rarely applied in land use time series studies [39,64]. GEE implementation of these methods within its standard functionality offers an opportunity to improve temporal consistency and promote temporally dense analyses as a common practice. In addition to the post-classification method presented here, other opportunities to fully benefit from the temporal information available would be a GEE implementation of methods combining the probability functions of land use trajectories and the propagation of classification uncertainties temporally to decide on the final classification outcome [33,[39][40][41]46,54,65]. Furthermore, data fusion of multiple sensors have shown large potential to further increase the temporal resolution of analysis [34,37,38]. In particular, the integration of optical and LiDAR or synthetic aperture radar imaging offers opportunities for high-frequency temporal analyses that are independent of the weather and cloud cover [66].

Conclusions
Our study aimed to advance temporal mangrove forest monitoring efforts to benefit from the potential that currently available satellite earth observation data and cloud-based high-performance computing can offer. The temporal domain of these information-dense datasets opens opportunities to apply data fusion principles to optimize classification outputs to be gap-free and temporally consistent with logical land use transitions as well as provide a means of cross-validation. The results of our case study on mangrove forests demonstrate how this information can be valuable in understanding the spatio-temporal dynamics, processes, and trends of land use changes and improve decision-making with detailed information. Thereby, our study builds on and contributes to a growing body of work advocating that accounting for temporal information in optimizing land cover classification is superior to temporally independent classifications in time series and other single-date methods.
Despite growing awareness, most mangrove forest cover classification studies are yet to take full advantage of Earth observation's potential and the rich temporal information available from time series data. Implementation of temporal optimization, either post-classification or during the classification process, within future implementations or that can be automated on top of GEE's output as presented here, can hopefully contribute to advance mangrove monitoring studies towards fully unlocking the potential of data available as the field of earth observation keeps evolving. The integration of synthetic aperture radar remote sensing in addition to optical observations will be key in advancing the approach presented here and increase the temporal resolution by overcoming data gaps due to weather or lighting conditions. Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/22/3729/s1, Table S1: Links to the GEE code implementation used for this study, Table S2: Decision rules of the regression tree (CART) build for land use classification (2015) by Landsat-7 and Landsat-8 respectively.