Study Area
Macajalar Bay is an embayment of the Mindanao Sea that encompassed the municipality of Laguindingan to Kinoguitan, Misamis Oriental (Fig. 1). A total of 12 coastal municipalities and 2 cities situated in Misamis Oriental, Philippines comprise the 470 square kilometer coastline of the bay (Canoy and Roa-Quiaoit, 2011). Its intertidal area is mainly occupied by mangroves forest stands, seagrass meadows, and coral reefs (Consolacion et al., 2021). Existing studies on mangrove forests found in the bay have put emphasis on diversity assessments and its importance to sustain corresponding faunal species (Olila and Relox, 2021; Taneo, 2021). However, important information that describes its current extent and historical change dynamics is not yet available. Moreover, recent reports on the Blue Carbon Ecosystems (BCE) studies in the Philippines elucidate that the mangroves forest located in Northern Mindanao is one of the least studied BCE sites in the Philippines (Corcino et al., 2023), which makes it imperative to conduct remote sensing-based studies to improve our understanding about the current state of mangroves present in the bay. Moreover, remote sensing-based studies that may depict the contemporary extent as well as disclose the historical changes are recommended to monitor development changes.
Dataset
Mangrove forest was estimated using National Mapping and Resource Information Authority (NAMRIA) topographic maps developed on aerial photographs from 1950 coupled with Landsat Surface Reflectance (SR) images from 1990, 2000, 2010, and 2020 to tract mangrove forest changes over a 70-year period. For the 1950 mangrove forest extent, topographic maps showing the underlying vegetation were georeferenced and projected to Universal Transverse Mercator (UTM) Zone 51N. The mangrove extent present in the georeferenced maps was digitized using an offline GIS software package. On the other hand, Landsat 5 Level 2, Collection 2, Tier 1 was utilized to map mangrove forests for 1990, 2000, and 2010 while Landsat 8 Level 2, Collection 2, Tier 1 was utilized to map the extent of mangroves for the year 2020, respectively. These image collections available at GEE data repository were radiometrically corrected using the methodologies described in Chander et al., (2009) and atmospherically corrected using Landsat Disturbance Adaptive Processing System (LEDAPS) algorithm version 3.4.0. Four bands (Band 3–5, and 7) were rescaled and included in the final image mosaics used to determine the existing mangrove extent (Potapov et al., 2012).
Our method employed a Landsat multi-date scene compositing to primarily overcome extreme cloud cover issues and to generate cloud-free satellite image mosaics representative of the target time epochs. For each epoch, image compositing was prepared independently where each of which was allocated with 2-year time intervals to ensure the production of cloud-free image mosaics. The selection of satellite images to map the existing mangrove forest in Macajalar Bay employed a temporal sampling design sensitive to the changing seasonal conditions. Green et al., (2005) mentioned that seasonal changes affect the phenological conditions of trees, and thus may subsequently result in different reflectance from month to month. In this study, image reflectance variability arising from external seasonal differences was minimized by selecting satellite images captured in the same season. We used the approach of Potapov et al., (2012) where images captured during the growing season of mangroves were used to create a cloud-free image mosaic of Macajalar Bay as they are more suitable conditions to map forest than that of the senescent condition. Growing seasons were defined geographically using the temperature and rainfall data observation of the Philippine Atmospheric, Geophysical, and Astronomical Services Administration (PAGASA) as a primary reference to delineate varying seasonal differences in the country. In this study, the growing season was from Julian Day 152 to Julian Day 334, respectively.
Creation of Finalized Region of Interest (FROI) to Delineate Potential Mangroves Sites
The mangrove forest ecosystem has a unique spatial pattern as it can only be observed at the coastal fringes of countries with tropical climate conditions. In this regard, several studies have shown that reducing the area extent of analysis may prevent spectral confusion with other existing image pixels (Jones et al., 2016; Yancho et al., 2020). Hence, the accuracy of mangrove forest maps critically depends on the demarcation of coastal fringes and the elimination of landscapes that may not utterly support mangrove survival.
Creation of the initial Region of Interest (ROI) was performed using the ‘drawing function’ tool of the GEE, with the high-resolution Google Earth base map used as the main reference to intuitively draw the initial ROI. A clip geoprocessing process between the initial ROI and Large-Scale International Boundary (LSIB) Polygon was employed to create an intersection polygon of the initial ROI and actual boundary of Misamis Oriental shoreline fronting the Macajalar Bay. Correspondingly, a buffer geoprocessing tool was executed at 5,10,20,25,30 kilometers (Yancho et al., 2020). We used the 25-kilometer buffer since it depicts the most conservative polygon that captures all mangrove forests present in 1950 within its boundary. The last procedure in creating the final ROI is to determine the elevation threshold limit to finalize the topographic data mask. In this study, we used 50 meters as the elevation threshold limit (Buitre et al., 2019). Areas that qualify for the aforementioned conditions were reproduced to a separate vector file format using a Boolean programming approach.
Collection of Classification Reference Points and Image Classification Schemes
The land cover classification was based on Santillan et al., (2019) but implemented in this study with modification which generally aggregated land covers that may fall within the same classes (e.g., different forest types). For this study, a total of six different land cover classes were included for the final land cover map analysis, including Forest, Sparse Vegetation, Water Bodies, Mangrove Forest, Built-up Areas, and Barren Land. A detailed class description as well as their corresponding numbers of Classification Reference Point (CRP) is shown in Table 1.
Table 1
Description of land cover classes and the number of Classification Reference Point (CRP) inputs used to generate the land cover map
Land Cover Classes | Description | Number of Classification Reference Point (CRP) |
1990 | 2000 | 2010 | 2020 |
Forest | Areas where terrestrial forest stands are observed. Pixels where FCD value is greater than 30%, and an NDVI value of > 0.6 | 86 | 74 | 55 | 56 |
Sparse Vegetation | Areas where terrestrial grasslands, and crops are observed. Pixels where FCD value is between 0–30%, and an NDVI value of < 0.6 | 136 | 74 | 124 | 99 |
Water Bodies | Areas occupied by water | 154 | 139 | 66 | 51 |
Mangrove Forest | Tall and short stand mangroves forest that is found within the elevation threshold mask. Pixels where FCD is > 30%, and MVI value is within the set threshold | 251 | 71 | 117 | 100 |
Built-up Areas | Areas where houses, roads, facilities are found | 50 | 50 | 50 | 50 |
Barren Land | Areas where exposed soil, sediments, and sands are found. This may also include agricultural areas under fallow period. | 50 | 50 | 50 | 50 |
To create an accurate land cover mapping model using GEE, an extensive set of CRP collections was prepared based on visual interpretations of mosaic satellite images and accompanied by information produced from band ratio calculation. Several external geographic datasets including Quickbird images from Google Earth, Planet Imagery footprints, and expert information were also considered for building the CRP collection. Identification of terrestrial forests was guided by the Normalized Difference Vegetation Index (NDVI) (Tucker, 1979) and Forest Canopy Density (FCD) models. A data mask that captures pixels that have an NDVI value of > 0.6 as well as an FCD value of > 30% was created to isolate forest pixels (Perez and Comiso, 2014). Pixels were carefully evaluated before digitization using High-Resolution maps of Google Earth as well as field-derived images. Pixels that qualify for the aforementioned conditions were manually digitized to represent forest areas for the image classification procedure. Pixels that record an NDVI value of < 0.6 and a canopy density of < 30% were digitized as sparse vegetation. Identification of bare ground was identified using the Modified Bare Soil Index (Nguyen et al., 2021).
Mangroves were identified based on two separate vegetation indices which aim to dissociate the mangroves from terrestrial vegetation and second, to separate mangroves from other sparse vegetation that may potentially co-exist in the mangrove forest ecosystem. First, a band ratio calculation called Mangrove Vegetation Index (MVI) was calculated to detect the potential location of mangroves within Macajalar Bay (Baloloy et al., 2020). Mangroves are quite distinct due to their inherent moisture compared to terrestrial vegetation and thus, it can be easily identified using MVI. The second approach concentrates on identifying thickly covered mangrove vegetation using FCD, which has been applied to classify mangrove forests based on their canopy cover (Sahana et al., 2015). In this study, we considered 30% as mangrove forest as it encompassed mangrove extent as viewed on high-resolution maps. Pixels that qualify for the aforementioned conditions were then manually interpreted and added as part of the CRP to represent mangrove forests.
Image Classification, Map Validation and Mangrove Extent Determination
Random Forest (RF) algorithm was utilized to perform a pixel-based image classification due to the proven success in producing highly-accurate mapping results for mangrove forests (Jhonnerie et al., 2015; Rogers et al., 2017). RF is non-parametric, and an ensemble learning that is based on multiple decision trees classifier. The number of decision trees was set to 100. CRP was segregated based on a 70:30 ratio to generate points for training and validating the model. In this study, 70% of the collected CRP were randomly selected as training data points, while the remaining 30% was used as a validation point dataset. The “.smileRandomForest” package implemented using the JavaScript environment of GEE was used to undertake image classification. The overall accuracy of the land cover map was determined using the error matrix “accuracy tool” native within the GEE. Moreover, the User’s Accuracy, Producer’s Accuracy, and Kappa Coefficient were calculated externally based on the error matrix produced by GEE.
Mangrove change dynamics were carried out in GEE as well. Detection of changes was applied in four time periods covered in this study: 1950–1990, 1990–2000, 2000–2010, and 2010–2020. Raster data showing the mangrove extent for each period were automatically saved as operational data that can be imported to a new code editor environment of GEE. Actual change dynamics analysis was performed by employing a sequential Boolean method that calculates total mangrove forest loss, gain, and persisted between periods. Results on the mangrove change dynamics mapping were superimposed with the land cover map over the finalized ROI to identify actual land cover conversion. The historical mangrove forest conversion was then summarized using the “davidsjoberg/ggSankey” package library generated using RStudio.
Discussions on the different drivers of mangrove forest loss have used historical inventories of mangrove stands and applied statistical techniques to elucidate the rate of clearing (Primavera, 1995; Primavera, 2000). In this study, a new way to systematically describe the underlying mangrove forest clearing driver was employed by detecting a minimum pixel neighboring approach. We implemented the “.connectedPixelCount” algorithm available at GEE to detect segments of the mangrove forest loss data with 34 connected pixels, or equivalent to 4 contiguous hectares of mangrove clearing. Moreover, conditional code formatting was generated to modify the code that detects large-scale clearing. Particularly, we added a line that aims to detect less than 34 connected pixels to quantitatively describe small-scale clearing loss in the entire Macajalar Bay, Misamis Oriental coastal area.