2.1.1. Regional Context

The region encompassing south (S) Asia, southeast (SE) Asia, and Asia-Pacific is home to approximately 46% of the world's mangroves [5,44]. This region includes some of the world's most productive, oldest, and biodiverse mangrove forests [45]. Regional loss—the highest in the world—is driven by conversion to aquaculture ponds (i.e., shrimp and fish farms), oil palm plantations and rice paddies, coastal development, and over-extraction for wood [11,46–57]. Natural processes and phenomena (e.g., rising ocean temperatures and sea-levels, severe tropical storms, and natural disasters) also contribute to regional dynamics [48,53,56,58–66]. Notably, SE Asia is exceptionally biodiverse containing 51 of the world's 73 documented mangrove species, compared to 10 in the Americas and Africa [5,67]. SE Asia alone contains an estimated 34% of the world's mangroves [5,68]. Recent studies show that mangrove areas in SE Asia are experiencing the highest prevalence of anthropogenic activity in the world [68,69].

#### 2.1.2. Myanmar—A Regional and Global Loss Hotspot

Located within SE Asia, the preliminary AOI for this pilot study is all of coastal Myanmar (Figure 1). As confirmed by Gandhi and Jones [19], within SE Asia, mangrove loss is most notable in Myanmar, making the country both a regional and global loss hotspot. Giri et al. [70] reported a 35% decrease in mangrove extent from 1975–2005 whereas De Alban et al. [57] reported a 52% decrease from 1996–2016 [57,70]. According to De Alban et al. [57] and Estoque et al. [56], the primary anthropogenic drivers of this loss include conversion to rice paddies, oil palm and rubber plantations, and increasingly for aquaculture (e.g., shrimp, fish) [56,57]. Natural drivers include tsunamis triggered by seismic activity, and tropical storms [68,71,72]. Within Myanmar, according to Giri et al. [70], Saah et al. [73], Bunting et al. [74], De Alban et al. [57], and Clark Labs [75] sub-national loss hotspots include the northwestern (NW) coastline, much of the Ayeyarwady peninsula, and a smaller area slightly east of the Ayeyarwady peninsula (Figure 1).

#### 2.1.3. Myanmar—Inventory, Summary and Acquisition of Existing Datasets

All national-level mangrove datasets providing single- or multi-date coverage for Myanmar up to July 2020 were inventoried through an exhaustive online search and literature review. When available, datasets were obtained from online repositories or through contacting authors. When not available, datasets were described based on associated literature. All datasets were summarized based on producer/organization/reference, single- vs. multi-date, temporal and spatial extent, availability, imagery source(s), mapping methods, and whether discrete or continuous (Table 1).

**Figure 1.** The preliminary region of interest (ROI) for the GEEMMM pilot representing coastal Myanmar; sub-national AOIs wherein qualitative accuracy assessments (QAAs) were untaken for existing maps (i.e., Baseline QAA AOIs); sub-national AOIs wherein GEEMMM QAAs were undertaken (i.e., GEEMMM QAA AOIs). Also shown are sub-national AOIs wherein classification reference areas (CRAs) were derived (i.e., CRA AOIs) and the location of known mangrove loss hotspots based on existing studies (i.e., Giri et al. [70], Saah et al. [73], GMW (Bunting et al. [74]), De Alban et al. [57], Clark Labs [75].

**Table 1.** Inventory and summary of existing national-level mangrove datasets for Myanmar—July 2020.

‡ Iterative Self-Organizing Clustering.

#### 2.1.4. Myanmar—Comparison of Existing Datasets and Baseline QAA

Once inventoried, all known datasets were compared based on mapped classes, mangrove distribution, accuracy, dynamics (when available), and known limitations. Where provided, mangrove distributions and dynamics were extracted from publications and supporting materials. If not readily apparent—and if the datasets were available—dynamics were calculated. Adding to the standard reported metrics, the accuracy was further qualitatively assessed for all available datasets through cross-checking in reference to high spatial resolution satellite imagery viewable through Google Earth Pro (GEP) [79]. This secondary qualitative accuracy assessment—or QAA—first reported in Gandhi and Jones [19], provides a more thorough understanding of existing mangrove datasets.

The QAA of existing maps (i.e., baseline QAA) was undertaken for the most recent entry in each discrete dataset, when available. Datasets were acquired in both raster and vector format, and in a range of coordinate systems, necessitating several pre-processing steps. For each baseline QAA, three 100 × 100 km sub-national AOIs were selected across Myanmar: in the north (Rakhine), in the center (Ayeyarwady Delta), and in the south (Tanintharyi) (Figure 1). Each baseline QAA AOI was divided into 10 × 10 km boxes, and working from NW to SE, every sixth box was selected for spot-checking, such that approximately 17% was systematically assessed. QAA of the Giri et al. [70] dataset was already conducted [19]. For the remaining datasets, each spot-check entailed comparing mangrove coverage to GEP imagery as close to the dataset's capture date as possible. In some instances, particularly in the southernmost Tanintharyi AOI, GEP imagery was partly/fully cloud-covered, limiting the ability to conduct QAA (limitations also noted by Estoque et al. [56]). A single mangrove class, representing the variability of canopy cover, height and stand structure in mangrove forests (as used in GEEMMM pilot classifications and defined below) was qualitatively assessed within each spot-check as either well-, under-, or over-represented. For each dataset, results help contextualize the representation of mangrove distribution and dynamics.

#### *2.2. The Google Earth Engine Mangrove Mapping Methodology (GEEMMM)*

The GEEMMM is intended to facilitate the mapping and monitoring of mangrove ecosystems anywhere in the world, without requiring a dedicated in-house geospatial expert. Intended users need basic computer skills and an understanding of the key steps required for mapping mangroves, but are not expected to hold advanced expertise in remote sensing, geospatial analysis, and/or coding. The interactive tool is broken into three modules—Module 1: defines customized region of interest (ROI) boundaries and generates multi-date imagery composites; Module 2: examines spectral separation between target map classes and undertakes multi-date classifications and accuracy assessments; Module 3: explores dynamics and o ffers an optional QAA. Each module is broken into thoroughly commented and referenced sections, bringing the user through all steps while making reference to this manuscript for full methodological details and context. Each module and the parameters used in this pilot study are described below. Table 2 provides a summary of all GEEMMM user inputs and variable selections for the Myanmar pilot study.

#### 2.2.1. Module 1—Defining the ROI and Compositing Imagery

In the first step of Module 1: Section 1, the user must identify key datasets to be used in the GEEMMM. The first user-defined dataset is a preliminary ROI. This is generated using the 'drawing tools' function built into GEE and clips all subsequent user-defined datasets. The second user-defined dataset is the known extent of mangroves which is used to calculate elevation and slope thresholds and shoreline bu ffer distance. The user can select the baseline GEE data set representing global mangrove distribution circa 2000 (i.e., Giri et al. [44]) or upload their own. The third user-defined dataset is coastline, for which the user can select the baseline Large Scale International Boundary Polygons [80] or upload their own. The fourth user-defined dataset concerns elevation which is required to generate topographic masks (i.e., elevation and slope). The user can select the GEE JAXA-ALOS satellite radar DSM (30 m) [81] or upload their own. For the Myanmar pilot, the preliminary ROI is shown in Figure 1, the GEE global mangrove distribution circa 2000 was used for known mangrove extent, the Global Administrative Boundaries database (GADM) Myanmar dataset (v3.6, www.gadm.org, an external source) for coastline, and the GEE JAXA-ALOS Global PALSAR-2/PALSAR Yearly Mosaic 25 m land-cover data for elevation [81,82].


**Table 2.** Summary of GEEMMM user inputs and selected variables used in the Myanmar pilot.

In the second step of Module 1: Section 1, the user defines input variables and sets how workflow thresholds are calculated. Table 2 lists all of the user variables and user inputs for the GEEMMM including those used in this pilot study. GEE provides unprecedented access to the Landsat catalog, offering approximately 1.3 M scenes from 1984 to present [34]. While it is certainly advantageous to have access to so many images, the choice of imagery based on parameters such as year(s) and time of year(s) must be considered carefully. Two variables define contemporary and historic year(s) of interest. There are two four-digit date inputs to bookend the historic and contemporary year windows. If the user wishes to isolate a single historic or contemporary year the same is selected for each book-end. Following the year(s) of interest, the month(s) of interest are selected. Seasonal variations can affect terrestrial vegetation adjacent to mangroves, and atmospheric conditions can change throughout the year, so the ability to target specific months is essential to generating optimal image composites [83–85]. The user identifies the month(s)-of-interest using two book-end numbers corresponding to the 12 months of the year; they may overlap the new year; e.g., "11" (Nov.) to "2" (Feb.). Next, the allowable cloud cover limit, an integer between 0 and 100, is used to filter the Landsat metadata [86]. Also related to cloud cover, the user decides whether to mask the imagery, and to what extent, i.e., setting a mild cloud mask using the USGS-provided (United States Geological Survey) quality band, or an aggressive cloud mask where pixels are excluded based on their 'whiteness' and a temperature band threshold [35]. For the sixth input, the approximate tidal zone—a numeric

input (in m) that represents the tidally active zone bu ffered inland from the coastline—is entered. Approximate tidal zone helps isolate the portion of images subject to reflectance changes from tidal variation, while reducing influence from other non-tidal variability. The default value is 1000 m. Next, the user chooses how water is masked out of the imagery, either using a mask developed from the water present in the contemporary imagery alone, or a combination mask based on pixels determined to be water in both historic and contemporary imagery. A pixel is determined to be water if its value was greater than the 0.09 modified normalized di fference water index (i.e., MNDWI) threshold established by Xu [87]. The modified normalized di fference water index (MNDWI) was developed to detect water pixels by calculating the normalized di fference between the green and short-wave infrared (e.g., Landsat 8 Operational Land Imager, 1.57–1.65 μm) bands, making it suitable for measuring the amount of water present in an acquisition. Topographic thresholds are set to generate masks based on elevation and slope. The user can either manually enter the elevation (m) and slope (%) thresholds, or have them automatically calculated based on the 99th percentile values extracted from within the known mangrove extent dataset. The user can further opt to search for inland-fringe mangroves, which have been documented as far as 85 km inland [75,88]. If inland-fringe mangroves are targets for the classification(s), the preliminary ROI is doubled for elevations lower than 5 m based on [89]. The last step in Module 1: Section 1 is the selection of spectral indices which the user would like to calculate for each image composite. After the workflow begins, the user chooses which indices they would like to calculate from a list of fourteen indices, including some which are mangrove-specific. The complete list of indices included in the GEEMMM can be found in Table 3. The contemporary and historic windows from which imagery was selected for the Myanmar pilot study were 2014–2018 and 2004–2008, respectively. The months of acquisition were limited to June through December, corresponding with the wet season and the months directly following that time [90]. The imagery was filtered using cloud cover information for each acquisition at a 30% threshold. All 14 spectral indices were selected for calculation.

Module 1: Section 2 determines the finalized ROI for processing. Numerous studies have demonstrated the utility of reducing the classification extent to the minimum required area—this approach helps reduce spectral confusion with unnecessary scene components [44,91]. The preliminary ROI is used to isolate a section of shoreline which is bu ffered at 5, 10, 15, 20, 25, 30, and 35 km intervals. 5 km intervals were used to ensure observable di fferences in bu ffer distances. 35 km was used as a maximum extent based on observations in several countries, including Myanmar. These bu ffer distances are used to calculate the area of known mangroves that falls within their respective bounds. The user either selects their bu ffer distance preference from a drop-down menu containing values in between, greater than, or less than the listed intervals.

In either case, the bu ffer distance is used to create the finalized ROI. This ROI is used to select Landsat path/row tiles and generate image composites, clip composite imagery and masks (i.e., elevation, slope, and water), define the classification and dynamics extent, and provide a visual aid for optional QAAs. The finalized ROI used in the pilot study was based on a 23 km bu ffered shoreline which represents the maximum observed distance between known mangrove extent (i.e., Giri et al. [44]) and Myanmar's coastline.

Module 1: Section 3 generates the imagery composites required for multi-date classifications. Given the daily dynamic nature of mangrove ecosystems—wherein tides inundate 2–3 times per day on average—tidal conditions and the associated presence (or lack thereof) of water must be considered—there are a growing number of mangrove detection indices which rely on the isolation of high and low tide imagery [29,92,93]. The GEEMMM uses an image-based approach to calibrate imagery based on high and low tide. For each available image, an MNDWI is generated and the land is masked out using JAXA-ALOS Global PALSAR-2/PALSAR Yearly Mosaic 25 m land-cover data. The MNDWI was selected as the key spectral index because it has been proven to be an improvement over the normalized di fference water index (NDWI), and was developed explicitly for detecting water and non-water pixels [87]. The shoreline is bu ffered to the user-defined tidal zone value and mean MNDWI is used to create a constant value band wherein the greater the MNDWI mean value, the more water present within the tidal zone, corresponding to higher-tidal conditions. A second value band is added to each available image by multiplying mean MNDWI by −1, isolating lower-tide conditions. Clouds, if present and opted to be, are masked prior to the calculation of mean MNDWI using only the pixel quality band or an aggressive approach where the three visible (red, greed, and blue) and thermal bands mask based on digital number reflectance thresholds. Under the aggressive filter, a pixel is considered to be a cloud if its visible spectrum bands digital number reflectance values are greater than or equal to 1850, and the thermal band (brightness temperature, Kelvin) digital number is less than or equal to 2955. For the Myanmar pilot, the aggressive cloud filter option was selected to filter the imagery in an e ffort to remove low-altitude clouds which were not correctly classified by the Landsat cloud detection algorithm. If/once clouds have been masked, all available images and their corresponding tidal value bands are used to create best available pixel-based highest observable tide (i.e., HOT) and lowest observable tide (i.e., LOT) composites. Composite generation works as if all available images were stacked and organized by desired tidal condition. For example, as the LOT composite is being generated, the imagery with the lowest observed tidal condition is placed on top, and any missing pixels in that image, e.g., clouds masked, would be filled by the next best tidal observation and so on until all the gaps are filled. This process takes place for both the contemporary and historic data sets, resulting in a maximum of four composites (i.e., HOT and LOT contemporary, HOT and LOT historic). Because tides are determined using value bands, it is possible that all of the pixels for HOT and/or LOT composites within a particular area may be from one image (e.g., if no clouds were present and that image represented best available tidal conditions). The GEEMMM employs USGS surface reflectance Landsat products, which are readily available within GEE [35,94].

Module 1: Section 4 calculates the user selected indices from Section 1 of Module 1 (Table 3). There are a growing number of Landsat-related spectral indices available, many of which relate directly to mangroves such as the submerged mangrove recognition index (SMRI) and the modular mangrove recognition index (MMRI) [29,43]. The GEEMMM provides the user with the option to select from 14 spectral indices, of which four are mangrove-specific. The selected indices are calculated for both contemporary and historic HOT and LOT composites and added as potential classification inputs. Figure 2 compares the appearance of a typical mangrove-dominated area in Myanmar across all of the available mangrove-specific spectral indices (i.e., combined mangrove recognition index (CMRI), MMRI, SMRI, MRI) in the GEEMMM [29,43,93,98].

In Module 1: Section 5 the classification extent is further reduced through masking. In accordance with numerous mangrove mapping studies (e.g., Jones et al. [91], Thomas et al. [68], and Weber et al. [105]), the GEEMMM incorporates cloud, water, slope, and elevation masks to produce a finalized AOI. The cloud mask is generated and applied before composites are produced. The water mask is calculated for each composite using the methodology established in Xu [87], where the MNDWI layer for historic and contemporary LOT composites are generated and then a threshold is applied. Pixels with a value greater than 0.09 are considered to be water and a binary mask is produced. Depending on user selection, the water mask is finalized by either using just contemporary or combining the historic and contemporary and selecting only pixels determined to be water in both composites. This pilot study used the combined water mask. The two topographic masks are generated through user-defined thresholds or automatically determined using the 99th percentile of elevation and slope for known mangroves. The Myanmar pilot study used the known mangrove extent to generate topographic masks based on elevation values > 39 m and slope values > 16%. Noting how minor elevation is within mangrove ecosystems, the elevation threshold actually represents an approximate combined elevation + canopy height past which mangroves are not found. The generated masks are combined to create a binary, single unified final mask which is applied to all composites within the finalized ROI.

Module outputs include: (1) HOT contemporary composite, (2) LOT contemporary composite, (3) HOT historic composite, (4) LOT historic composite, (5) finalized ROI, and (6) Finalized Mask.


**Table 3.** List of all spectral indices available in the GEEMMM including mangrove-specific.

\* denotes mangrove-specific spectral index.

**Figure 2.** The appearance of a typical mangrove-dominated area in Myanmar across all of the available mangrove-specific spectral indices (i.e., CMRI, MMRI, SMRI, MRI) in the GEEMMM [91].

#### 2.2.2. Module 2—Spectral Separability, Classifications and Accuracy Assessment

For Module 2: Section 1, user inputs address classification variables and settings. The user enters the asset path for historic and contemporary classification reference areas (CRAs) (i.e., the user-defined examples of target map classes) and identifies the unique column labels for class names and numeric codes. Next, the user identifies whether CRAs are spatio-temporally invariant (i.e., each CRA represents a class example in both contemporary and historic imagery). If the CRAs are not spatio-temporally invariant, the spectral properties of the contemporary CRAs are extracted and used to define class boundaries in the historic classification(s). For classification algorithm the single option is currently random forest [106]. The user determines how many trees are employed. The final input determines classification outputs. Users have the option to select outputs from either HOT or LOT composites for contemporary and historic inputs (i.e., four possible outputs), and/or a combined classification where HOT and LOT composites are merged to create single outputs (i.e., two more possible outputs), totaling six possible classification outputs. Zhang and Tian [93] demonstrated the utility of using combined HOT and LOT image composites as classification inputs. For the Myanmar pilot, 200 trees were selected with outputs based on combined (i.e., HOT and LOT) historic and contemporary classifications (i.e., two classifications).

In Module 2: Section 2, the user can examine correlation between potential spectral indices and the spectral separability of CRAs across all potential classification inputs. The Pearson's correlation is calculated for each selected index to all others and these values are used to generate a correlation matrix with values ranging from −1 to 1 [107]. A value of 1 means that the potential inputs have a perfect, positive, linear correlation, and a value of −1 indicates that the indices have a perfect, negative, linear correlation. Users are encouraged to select indices that are not highly correlated indicating that they provide unique information. As a general rule, correlation coe fficients with absolute values greater than 0.7 are considered moderately to strongly correlated and thus present similar information [107]. Users are advised to consider that correlation coe fficients are also impacted by the amount of variability in the data, the shapes of distributions, and the presence of outliers among other factors [108].

The spectral separability between target map classes as represented by CRAs is explored through the generation of three types of graphs. First, the user can view the spectral separability between each target class and each Landsat band—the user has the option to view this output for each of the four imagery composites. Box-and-whisker plots show the min, max, and inter-quartile range for each band and each map class. The second set of graphs is similar to the first, except that spectral separability is shown for individual indices across all of the target classes, showing only one index at a time. The final graph shows spectral feature space, where the x and y axes are user selected bands or indices. For the pilot study, and based on previously established precedents in Jones et al. [91], we included the visible, NIR and SWIR Landsat bands. Based on the correlation matrices and further the spectral separability they provided, the MNDWI, CMRI, MMRI, enhanced vegetation index (EVI), and Land Surface Water Index (LSWI) indices were selected as additional classification inputs.

For piloting the GEEMMMinMyanmar, six classes wereinitially targeted, including, (1) closed-canopy mangrove, (2) open-canopy mangrove, (3) terrestrial forest, (4) non-forest vegetation, (5) exposed/barren, and (6) residual water. Table 4 provides class descriptions and an overview of how many CRAs were digitized per class. CRAs can be derived within the GEE environment or externally. For this pilot, 90 × 90 m (i.e., 3 × 3 Landsat pixels) CRAs were derived externally. To ensure that internal class variability was captured for each class and across the AOI, three sub-national AOIs were used to define CRAs (Figure 1). CRAs were derived referring to finer spatial resolution satellite imagery viewable in Google Earth Pro (Google, Mountain View, CA, USA), existing contemporary land-cover maps for Myanmar (i.e., Giri et al. [44], Saah et al. [73], and De Alban et al. [57]), and expert interpretive knowledge gained with mapping mangroves in other regions of the world. Two mangrove classes were defined to ensure that the internal variability of mangrove forests based on stature, canopy cover and density was captured. Figure 3 shows examples of all targeted classes in HOT, LOT, a key spectral index, and finer spatial resolution imagery viewable in Google Earth Pro (Google, Mountain View, CA, USA) [79].


**Table 4.** Names and description of classes and numbers of classification reference areas (CRAs). Also shown is how many CRAs were derived within each sub-national CRA AOI (Figure 1).

In Module 2: Section 3: once the user confirms their final choice of classification inputs and target classes, classification—the process by which remotely sensed data is assigned land-cover classes—can occur [109,110]. There are many established algorithms for classifying Landsat data to produce maps of mangrove distribution, including classification and regression trees (CART), support vector machines (SVM), unsupervised k-means, decision trees, and maximum likelihood (ML) [28,29,42,111,112]. Many of these algorithms are available to use within the GEE environment; however, random forests—also available in GEE—is well established and used to map mangroves across the world, with distinct success within the GEE environment [27,41,43,106,113]. The inputs for random forest include an imagery data set (i.e., selected Landsat bands and spectral indices), training data (i.e., randomly selected 70% of CRAs), and a numeric parameter determining the number of 'trees' to be employed. For each classification the output is a single band raster with the same spatial resolution as the input data (30 m), with each pixel assigned a map value based on target classes. Following classification, the user can choose to merge map classes—this is particularly advantageous in scenarios where initial map classes were used to capture variability, but for which confidence in class boundaries may be lacking. For example, in the Myanmar pilot, we merged the two mangrove classes (i.e., closed- and open-canopy) post-classification. This ensured capturing mangrove variability while not having to draw a distinct boundary between these potentially overlapping classes in the final map.

Classification accuracy—defined as "a comparison of the derived product to ground condition"—is not reported in numerous studies involving mangrove mapping [114,115]. Following classification and optional class merging, in Module 2: Section 3, the GEEMMM automatically produces resubstitution and error matrices for all output classifications [116]. The resubstitution matrices determine end land-cover class for the CRAs used for training the classifier. The error matrices use 30% of CRAs held back from classification to independently evaluate map accuracies. The overall accuracy is reported using the error matrix 'accuracy' tool, found within the GEE library. Overall accuracy is printed below both the error and resubstitution matrices. By reviewing the error matrices and visually inspecting the output maps the user may wish to collapse/further collapse classes (e.g., if two classes are very confused). If the user combines classes, they can opt to re-calculate accuracy, re-generating resubstitution and error matrices. The final step for all users to exporting the classification maps to their assets. Module 2 outputs include, (1) correlation and spectral separability graphs, (2) classified maps, and (3) accuracy assessments.

**Figure 3.** The appearance of all targeted classes in highest observable tide (HOT), lowest observable tide (LOT), key spectral indices, and fine spatial resolution satellite imagery viewable in Google Earth Pro (Google, Mountain View, CA, USA) [79]. The HOT and LOT composites represent 432 (R: NIR, G: red, B: green) or 453 (R: NIR, G: SWIR, B: red) false color. The spectral indices include enhanced vegetation index (EVI—[101]), combined mangrove recognition index (CMRI—[98]) and modified normalized difference water index (MNDWI—[87]).

#### 2.2.3. Module 3—Dynamics and QAA

In Module 3: Section 1, the user indicates which classification(s) will be used to calculate dynamics and/or assess optional QAA. If desired, the user can further clip classifications to a country's boundary— if pertinent—using the GEE Large Country Boundary Polygons, or by uploading an external dataset. For the

Myanmar pilot we further clipped using a uniquely uploaded boundary (GADM) and exclusive economic zone (EEZ) from Marine Regions (v10 World EEZ,) [117]. For the QAA, the user enters CRA information (e.g., asset path, class names, and unique class numbers).

In Module 3: Section 2, multi-date outputs are used to quantify dynamics. This is foundational to understanding long-term trends and the effectiveness of conservation efforts. The user selects which map class they would like to view, and loss, persistence, and gain (i.e., LPG) are calculated. The automatically produced, self-masked layers are added to the GEE-GUI map interface. The resulting area for each dynamic assessment is printed to the console, expressed in hectares. Building on the inventory, description, acquisition and comparison of existing datasets, the dynamics resulting from this GEEMMM pilot were also compared to published values.

Module 3: Section 3—building on the previously referenced methods detailed in Gandhi and Jones [19]—facilitates an optional QAA. For this GEEMMM QAA, an interactive map is divided into three linked maps (Figure 1). In each map, two sets of grids are automatically generated, (1) 100 km by 100 km grids, and (2) within each of those cells, sub-divided 10 km by 10 km grids. The 100 km × 100 km grid cells are randomly selected, retaining 50% of the grid cells that intersect the ROI. In slight contrast to the baseline QAA described in Section 2.1.4, for the QAA tool in the GEEMMM, within each selected grid cell, 20% of the sub-grid cells are selected. The tool works by cycling through the sub grid cells, and giving the user the option to view simultaneously on linked maps showing Landsat composites where the date can be changed at the user's preference, the classifications produced in Module 2, and the imagery used for the classifications generated in Module 1. The user then has the ability to record in the GUI whether each map class is under-, well-, or over-represented, and record 'free comments' for each sub-cell.

Module 3 outputs include: automatically generated LPG as raster and—if performed—QAA grid (for viewing outside of GEE). The user also has the option to export the QAA table (containing the under, over, and well representation statistics, and the free comments) as a CSV (i.e., comma separated values) file at any point during the QAA.

#### **3. Results and Discussion**

#### *3.1. Myanmar—Comparison of Existing Datasets*

Table 5 provides a comparison of all single- and multi-date datasets based on dataset/authors, year, extent (ha), dynamics (ha and %), whether discrete or continuous, mapped classes, accuracy, and known limitations. Figure 4 provides a comparison of all distributions across time across all datasets. Results show that Myanmar's mangrove distribution ranged from 851,452–1,323,300 ha circa 1975–1987 (i.e., historic) to 475,637–1,002,098 ha circa 2014–2018 (i.e., contemporary). Of the 11 existing studies, only five provided quantitative accuracy assessments, with overall accuracies ranging from 76% (i.e., Saah et al. [73]) to 97% (i.e., Estoque et al. [56]), mangrove producer's accuracies ranging from 75% (i.e., De Alban et al. [57]) to 93.1% (i.e., also De Alban et al. [57]), and mangrove user's accuracies ranging from 92.3% (i.e., De Alban et al. [57]) to 98.1% (i.e., Clark Labs [75]). Of the existing studies, eight provided dynamics, including a loss of 300,091 ha/35.2% from 1975–2005 (Giri et al. [70]), 195,227 ha/16.3% from 1987–2018 (Saah et al. [73]), 43,208 ha/8.0% from 1996–2016 (Bunting et al. [74]), 694,600 ha/52.5% from 1996–2016 (De Alban et al. [57]), 76,465 ha/10.9% from 1999–2018 (Clark Labs [75]), 27,064 ha/9.7% from 2000–2014 (Hamilton and Casey [78]), 27,770 ha/5.5% from 2000–2012 (Richards and Friess [47]), and 191,122 ha/28.7% from 2000–2014 (Estoque et al. [56]). Two reported specifically on sub-national loss hotspots (i.e., De Alban et al. [57] and Estoque et al. [56]). According to De Alban et al. [57], Bago, Mon, Yangon—the three states immediately to the east of the Ayeyarwady delta—suffered greatest proportionate loss from 1996–2016 totaling more than 80% of their extents. In terms of absolute loss, from 2000–2014, Estoque et al. [56] reported Rakhine as the state with the greatest loss (75,494 ha/39.5% of Myanmar's total loss), followed by Ayeyarwady experiencing 69,431 ha/36.3% of Myanmar's total loss.


**Table 5.** Comparison of single- and multi-date datasets based on mapped classes, accuracy, mangrove distribution (ha), dynamics, and known limitations. Accuracy:


\* WAM data over Myanmar from Ministry of Forestry's Remote Sensing and GIS Section, derived from Landsat imagery 2000–2007.

#### *Remote Sens.* **2020**, *12*, 3758

**Table 5.** *Cont.*

**Figure 4.** Comparison of distribution for all existing single- and multi-date mangrove distribution maps for Myanmar, including results of GEEEMMM pilot.

Direct comparisons of existing datasets are challenging due to differences in temporal coverage, methodologies, and imagery sources. Although most studies use optical imagery (typically medium-resolution Landsat), some of the more recent studies combine optical with radar imagery (e.g., Bunting et al. [74]; De Alban et al. [57]). Several different mapping techniques are employed, while two of the datasets (Hamilton and Casey [78]; Richards and Friess [47]) calculate and present continuous measures of mangrove canopy cover, rather than discrete (i.e., presence vs no presence). Interpreting continuous datasets for areal mangrove extent is problematic as pixels containing just 0.01% canopy cover are included as mangrove falling well below commonly used minimum definitions mangrove forest (e.g., 30%) [18,78,91].

Of the five datasets reporting, all achieve overall accuracies of >75%, with four >85% [56,57,74,75]. QAAs further identified Clark Labs [75] as mapping mangroves in Myanmar most consistently. Mangroves were under-represented in the remaining five datasets assessed, particularly in Giri et al. [70] and Saah et al. [73], but also in Bunting et al. [74] (Table 6).


**Table 6.** Results of QAA for available/acquired datasets (1 = best).

Existing studies clearly establish that Myanmar has experienced consequential mangrove loss; however, baseline distributions and dynamics (when available) are highly variable. These discrepancies are likely attributed to the di fferences highlighted in Table 5. In addition, the definitions for mangroves and surrounding land-cover classes and the actual examples used for classification (i.e., CRAs) likely further account for di fferences. Only with agreed upon conventions for defining mangroves and providing examples as CRAs can cross-study comparisons become standardized and optimized. Falling short of this, discrepancies will remain common.

#### *3.2. Results of the Google Earth Engine Mangrove Mapping Methodology (GEEMMM)*

#### 3.2.1. Module 1—Defining AOI and Compositing Imagery

As confirmed through qualitative ye<sup>t</sup> systematic spot-checks, the imagery generated from the Myanmar pilot reflects the selected inputs well—both historic and contemporary composites are mostly cloud- and artifact-free, and clearly represent distinct HOT and LOT conditions (Figure 5). Figure 6 shows a national overview of the AOI including contemporary and historic HOT and LOT composites. The challenges that have been identified can be attributed to the extent of the study area and trying to capture a long, complex coastline in a series of contiguous composites. The most notable challenge relates to seasonal variability observed primarily in areas where large clouds were masked out of one image and the pixels selected to fill captured seasonally different land-cover conditions. Notably, this issue was almost entirely associated with areas which undergo significant changes throughout the year, i.e., agricultural mosaics and non-forest vegetation. Even within the defined seasonal window, variability was observed. Users are advised to select meaningful seasonal windows that restrict such variability while still offering enough imagery to make optimal composites—this is constrained by the size of the AOI.

**Figure 5.** Examples of image composite outputs from the GEEMMM showing lowest observable tide (LOT), panel (**a**) and highest observable tide (HOT), panel (**b**). The north oriented, false colour (R: NIR G: SWIR B: Red) Landsat image is over Kaingthaung Island, Ayeyarwady Region, Myanmar.

**Figure 6.** National overview of image composite outputs from the GEEMMM showing highest observable tide (HOT) and lowest observable tide (LOT) (false color composites, R: NIR, G: SWIR, B: Red). The composites were further reduced in area using topographic and combined water masks. (**A**) Contemporary HOT; (**B**) Contemporary LOT; (**C**) Historic HOT; (**D**) Historic LOT.

#### 3.2.2. Module 2—Spectral Separability, Classifications and Accuracy Assessment

Based on correlation analysis of all available spectral indices, five (i.e., MNDWI, CMRI, MMRI, EVI, and LSWI) stood out as not correlated (Appendix A) and were selected as classification inputs. Using the spectral separability tools, all target classes as represented by CRAs were assessed across all non-thermal (red, green, blue, NIR, SWIR1, and SWIR2) Landsat bands (Figure 7), and each selected index was further evaluated to confirm that it provided additional separation for one or more classes (e.g., MMRI: Figure 8). Results indicate that bands SWIR1 and SWIR2 were particularly helpful in separating non-forest vegetation. Non-forest vegetation was most confused with other vegetation classes in the visible spectrum and MNDWI. For terrestrial forest, NIR, MNDWI, and MMRI provided separability. In particular, MNDWI provided good separation from mangroves; whereas within the visible spectrum and CMRI the most confusion was noted, particularly with other vegetation classes. Closed-canopy mangrove was best distinguished by LSWI, MMRI, and to a limited extent bands SWIR1 and SWIR2. In contrast, open-canopy mangrove was best distinguished by CMRI,MNDWI, NIR, SWIR1, and SWIR2. While there are meaningful and distinct differences between the two canopy-based mangrove classes, there is spectral overlap—this speaks to the advantage of capturing the variability within mangrove forests while subsequently merging into a single class post-classification. Field work is required to confidently define the boundaries between these sub-mangrove types to make them final map classes—following classification and prior to validation, mangroves were merged into a single class (i.e., mangrove). Taken together, the combined mangrove class exhibited some confusion with terrestrial forest and non-forest vegetation classes in EVI, the visible bands, and SWIR1 and SWIR2. The exposed/barren class had the most separability in indices CMRI, MMRI, and LSWI, and the most confusion with non-forest vegetation and terrestrial forest notably in MNDWI and residual water in the visible bands. Residual water was easily distinguished with MNDWI, and the non-visible bands, but was confused with exposed/barren in the visible bands, non-forest vegetation within CMRI, and all classes within EVI.

For both historic and contemporary classifications, resubstitution accuracies were 100%, indicating all training data was assigned to the correct land-cover class. Based on accuracy assessments using independent validation data, overall accuracies for historic and contemporary classifications were 97.0 and 98.5%, respectively (Table 7). For the contemporary classification, there was slight confusion between terrestrial forest and mangroves. Additionally, there was a small amount of two-way confusion between non-forest vegetation and terrestrial forest. The greatest source of error for the historic classification was the non-forest vegetation class, which was at-times confused with mangroves and the exposed/barren class.

bands. The set of bar and whisker plots shows the min, max, and interquartile range.



#### 3.2.3. Module 3—Dynamics and QAA

Classification results indicate that circa 2004–2008, Myanmar contained 995,412 ha of mangroves. In contrast, by 2014–2018, Myanmar contained 642,659 ha of mangroves. These results sugges<sup>t</sup> that from 2004–2008 to 2014–2018 there was 551,570.99 ha of loss and 198,818.42 ha of gain (i.e., net loss 352,752.57 ha or 35.4%) (Figures 4 and 9, Table 5). As compared to Estoque et al. [56] and Giri et al. [70], estimated rates of loss are within reported trends and ranges; however, other studies reported lower rates of loss often coinciding with lower total estimates of mangrove cover (i.e., Bunting et al. [74]; Hamilton and Casey [78]; Richards and Friess [47]; Estoque et al. [56]). Figure 9 shows LPG from GEEMMM results within the loss hotspots identified through existing literature (i.e., Figure 1).

**Figure 9.** (**Left**) panels: Known mangrove loss hotspots (Figure 1). (**Top left**) shows loss, persistence, and gain (LPG) from 2004–2008 to 2014–2018 in Rakhine State; (**middle left**) panel shows the Ayeyarwady Region; (**bottom left**) shows Tanintharyi Region. (**Right**) Panel: contemporary high tide (HOT) image composite, false colour (R: NIR G: SWIR B: Red) with boundaries of left panels highlighted in cyan.

While there was a substantial net loss based on GEEMMM results, the reported gain seems relatively high. Portions of this likely reflect actual natural processes and increases in mangrove extent; however, the overall gain estimate is likely an overestimation. Exaggerated gain likely reflects the desk-based process of deriving CRAs. Clearly *any* classification is only as good as the examples used to calibrate the algorithm, and a limitation of this pilot was no direct access to field observations or ground truth, and constrained access to historical high spatial resolution satellite imagery. Disproportionate mangrove gain therefore likely reflects an underrepresentation of lower stature, less dense mangroves in the historic classification, which in turn exaggerates the amount of supposed gain (i.e., many of these areas were likely actually mangroves in both dates). Extensive field work and ground verification is required to confirm.

The GEEMMM QAA was conducted for the contemporary map, then repeated for the historic map. As part of the contemporary QAA, spot-checks were conducted over 108 sub-grid cells across Myanmar (Figure 1). The mangrove class was generally well-represented; however, at-times under-represented in favor of classes depicting portions of areas in the variable agricultural mosaic, i.e., non-forest vegetation, and exposed/barren. In both the contemporary, and less so the historic map, the agricultural mosaic was depicted as a patchwork of these two classes, on a pixel-by-pixel basis, given the inherent variability within the seasonal window. This resulted in some confusion between the two classes, and to some extent an under-representation of mangroves. In the contemporary map, sparser mangroves at the ecosystem periphery were at-times misclassified as non-forest vegetation, thereby under-representing mangrove and over-representing non-forest vegetation. Terrestrial forest was also at-times over-represented, occasionally at the expense of actual mangrove areas. Overall, the contemporary classification appeared to best represent Myanmar's south (i.e., the Tanintharyi coastline). The historic QAA, while not quite as comprehensive as the contemporary QAA (mainly due to the absence of historic imagery in GEE), found the mangrove class to be generally well-represented, though at-times over-represented at the expense of classes depicting the agricultural mosaic—an inverse to the contemporary map. Some portions of the agricultural mosaic were also found to be misclassified as terrestrial forest. As with the contemporary map, Myanmar's southern Tanintharyi coastline seemed best represented. Notably, most existing studies did not provide standard quantitative accuracy assessments, and no existing studies went beyond these and further qualitatively assessed resulting maps. While quantitative accuracy assessments should be a standard part of reporting, QAAs also help further assess resulting maps and identify areas for improvement. As such, the GEEMMM goes beyond standard accuracy—for which GEEMMM results were very high—and allows users to more closely examine actual distributions and subsequently dynamics.

#### 3.2.4. Dissemination and Improvement

The code is available in a GitHub repository (see Supplementary Materials Section), with a GNU GPLv3 license permitting free use, modification, and sharing, provided that the source is disclosed and not used for commercial purposes. The code runs based on provided links, or is copied-and-pasted into GEE, which remains available for free non-profit and educational use. The tool itself continues to be adjusted and updated, as the GEE library evolves and as new mangrove remote sensing techniques become available.

While the tool performs well there are always potential improvements. Notably, CRAs are a key input for the workflow, and highly influence the outcome of the classifications. Future applications of the GEEMMM would benefit from direct access to field-based ground truth when deriving CRAs, particularly when it comes to confidently using mangrove sub-types as final map classes. While the need for and merit of isolating tidal conditions is proven, which tidal conditions are best requires further exploration—we used combined HOT and LOT in this pilot, whereas HOT or LOT on its own could also be employed. Furthermore, the choice of tidal condition depends on the intended application. For example, mangrove carbon projects may favor using HOT composites on their own for more conservative estimates of mangrove extent and change.

While going beyond standard accuracy metrics, the QAA is a somewhat complex component requiring significant user interaction; however, it too will evolve as the GEEMMM is further tested with other settings and applied to other AOIs. GEE itself also has notable limitations: the AOI can be as large or as small as the user requires but GEE has computational limits. Google shares its cloud processing among all GEE users, which means that if the task requested to process is too large (e.g., a long complex coastline, with collections containing hundreds of images) the user's allocated capacity may be exceeded and error(s) returned. Additionally, the functioning of this tool requires a relatively stable and reasonably strong internet connection, especially to view images and products within the GUI. If internet connectivity is limited, there may be latency issues loading data or even time-out errors. One of the benefits of working within the GEE environment however is that once a data product export has begun it will be completed on Google's server side. This means that internet access can be interrupted while using the tool, and it will continue to run. It was this feature of GEE, the server-side image/vector data exporting that drove the current configuration of three modules, where intermediate data products are exported to the user's assets, e ffectively saving their progress through the tool.

The GEEMMM is currently designed around the use of Landsat data—this was a conscious choice based around data availability. Sentinel imagery—which is also available through GEE—o ffers an increased revisit time (i.e., higher temporal resolution) and finer spatial resolution; however, it remains limited by a 2015-present temporal window. In contrast, the Landsat archive in GEE o ffers >35 years of imagery which facilitates more historically meaningful and robust dynamics assessments while also providing enough imagery to draw from multiple years to produce composites within preferred

seasonal windows. Given the added benefits—especially once the archive spans 10+ years—future versions of the GEEMMM should also o ffer the choice of Sentinel imagery to users as an option.
