The Google Earth Engine Mangrove Mapping Methodology (GEEMMM)

: Mangroves are found globally throughout tropical and sub-tropical inter-tidal coastlines. These highly biodiverse and carbon-dense ecosystems have multi-faceted value, providing critical goods and services to millions living in coastal communities and making signiﬁcant contributions to global climate change mitigation through carbon sequestration and storage. Despite their many values, mangrove loss continues to be widespread in many regions due primarily to anthropogenic activities. Accessible, intuitive tools that enable coastal managers to map and monitor mangrove cover are needed to stem this loss. Remotely sensed data have a proven record for successfully mapping and monitoring mangroves, but conventional methods are limited by imagery availability, computing resources and accessibility. In addition, the variable tidal levels in mangroves presents a unique mapping challenge, particularly over geographically large extents. Here we present a new tool—the Google Earth Engine Mangrove Mapping Methodology (GEEMMM)—an intuitive, accessible and replicable approach which caters to a wide audience of non-specialist coastal managers and decision makers. The GEEMMM was developed based on a thorough review and incorporation of relevant mangrove remote sensing literature and harnesses the power of cloud computing including a simpliﬁed image-based tidal calibration approach. We demonstrate the tool for all of coastal Myanmar (Burma)—a global mangrove loss hotspot—including an assessment of multi-date mapping and dynamics outputs and a comparison of GEEMMM results to existing studies. Results—including both quantitative and qualitative accuracy assessments and comparisons to existing studies—indicate that the GEEMMM provides an accessible approach to map and monitor mangrove ecosystems anywhere within their global distribution.


Introduction
Mangroves are a species of woody plants which comprise unique, halophytic communities in the tropical and sub-tropical inter-tidal coastlines of the world [1]. When meeting accepted definitions based on attributes including height, diameter and canopy closure, mangroves can qualify as forest [2]. Areas not qualifying as forest are peripheral parts of wider mangrove ecosystems, including expanses audiences [34,37]. In addition, tools built for GEE and distributed over the Internet can facilitate methodological repeatability while providing opportunities for adaptability and customization [38].
To date, several studies have explored and demonstrated the utility of GEE for mapping mangroves yielding encouraging results and improvements over conventional methods [39][40][41][42][43]. While there is clear utility for mapping and monitoring mangrove ecosystems using GEE, published methodologies remain inaccessible to many would-be users. To replicate published methods requires an advanced level of specialized expertise with remote sensing, geospatial processing techniques, and/or coding. To date, no intuitive and accessible version of a mangrove mapping methodology within GEE has been proposed which caters to a wider audience of non-specialist conservation managers and decision makers. In addition, existing tools fail to fully capitalize on the wealth of local knowledge and understanding often held by coastal managers. Lastly, no single methodology comprehensively incorporates all of the best available options for mapping and monitoring mangrove ecosystems from across existing published studies and includes a widely applicable approach toward tidal calibration.
Herein we present a comprehensive, intuitive, accessible, and replicable methodology encapsulated in a new tool-the Google Earth Engine Mangrove Mapping Methodology (i.e., the GEEMMM). The GEEMMM was designed to provide a ready-to-go methodology for non-expert practitioners to map and monitor mangrove ecosystems, enabling them to combine their local knowledge with GEE's cloud computing capabilities. We developed the GEEMMM following a thorough review of mangrove remote sensing literature and incorporating the best available practices. In addition, our approach to tidal calibration operates completely within GEE based entirely on shoreline reflectance (i.e., image-based). To demonstrate the tool, we present an example of multi-date, desk-based (i.e., involving no field work) mapping and change assessment for Myanmar (Burma)-a global loss hotspot [19]. The GEEMMM-freely accessible to non-profit users-runs on detailed and well commented code within the GEE environment and is adaptable to any mangrove area of interest. GEEMMM outputs include multi-date classified maps, accuracies, and dynamic assessments. To set the stage for trailing the GEEMMM for Myanmar and contextualizing the outputs, and similar to methods detailed in Gandhi and Jones [19], all existing single-and multi-date mangrove maps for Myanmar were inventoried, described, and compared, with an emphasis on existing information on distribution and dynamics. We introduce the pilot area of interest (i.e., AOI), describe existing datasets, overview the GEEMMM tool, and compare the results to existing datasets.

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).

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]. 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.

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 offers 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.

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 buffer 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]. 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][84][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 buffered 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 difference water index (i.e., MNDWI) threshold established by Xu [87]. The modified normalized difference water index (MNDWI) was developed to detect water pixels by calculating the normalized difference 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 buffered at 5, 10, 15, 20, 25, 30, and 35 km intervals. 5 km intervals were used to ensure observable differences in buffer distances. 35 km was used as a maximum extent based on observations in several countries, including Myanmar. These buffer distances are used to calculate the area of known mangroves that falls within their respective bounds. The user either selects their buffer distance preference from a drop-down menu containing values in between, greater than, or less than the listed intervals.
In either case, the buffer 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 buffered 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 difference water index (NDWI), and was developed explicitly for detecting water and non-water pixels [87]. The shoreline is buffered 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 effort 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.  Enhanced Vegetation Index EVI 2.5*((NIR − red)/NIR + 6*Red − 7.5*Blue + 1)) Huete et al. [101] Mangrove Recognition Index MRI * |GVI(l) − GVI(h)|*GVI(l)* (WI(l) + WI(h)) Zhang and Tian [93] Submerged Mangrove Xia et al. [29] Land Surface Water Index LSWI Van Deventer et al. [103] Enhanced Built-up and Bareness Index EBBI (SWIR1 − NIR)/(10* √ (SWIR1 + LWIR)) As-syakur et al. [104] * denotes mangrove-specific spectral index.

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 coefficients 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 coefficients 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 GEEMMM in Myanmar, six classes were initially 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).

Class
Class Description Contemporary Historic

Non-Forest Vegetation
Grass and/or shrubs dominate; some exposed soil + scattered trees; canopy < 30% closed; active cropland, vegetation appears green 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.

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 boundaryif 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. 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] [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.   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].

Myanmar-Comparison of Existing Datasets
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). 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 differences 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 differences. 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.

Module 1-Defining AOI and Compositing Imagery
As confirmed through qualitative yet 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.

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.    (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). 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.

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, effectively 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-offers 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 offers >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 offer the choice of Sentinel imagery to users as an option.

Conclusions
We present a new tool-the GEEMMM-for mapping and monitoring mangrove ecosystems. By leveraging GEE, this new tool circumvents many traditional barriers to conventional methods. In addition, it presents an internal, image-based approach for tidal calibration. The GEEMMM-including the well commented source code-is available online and is ready to be used by practitioners anywhere mangrove ecosystems exist; please see information in Supplementary Material Section on how to access the GEEMMM.
While operational, the GEEMMM is not without its limitations: the larger the area the more complex the mapping task, particularly when it comes to creating optimal imagery composites within defined seasonal windows. In addition, the upper limits of GEE and internet connectivity present a challenge in terms of the time associated with and reliability of running the GEEMMM; however, when compared to the conventional processing times associated with standalone workstations it remains much faster, and once a part of the GEEMMM starts running it will continue to run even if the internet connection is lost. In any application, the resulting maps and dynamics assessments will only ever be as good as the examples of target map classes provided. Coastal managers will normally have such information available to them and GEEMMM provides them with a framework through which to capitalizes on this local knowledge, rather than relying on external datasets, which allow little to no customization, to map and monitor their mangroves.
The GEEMMM makes a significant and ready-to-go contribution toward accessible mangrove mapping and monitoring. It also remains a living tool wherein non-profit users are encouraged by the authors to make useful suggestions for modifications or additions, or modify the tool directly themselves to meet their own customized needs. While piloting the GEEMMM for Myanmar is an important first step, additional applications and tests are required, particularly for smaller areas of interest, wherein the GEEMMM can help fill a critical sub-national mapping gap. The authors welcome the opportunity to receive feedback from and work with users to more comprehensively assess the tool and gauge areas for improvement. A series of in-person and online instructional materials will go a long way toward ensuring the maximum and optimal utility of the GEEMMM. This first iteration of the GEEMMM further sets the stage for a comparatively more automated and even more accessible version to be deployable completely on mobile devices.