Using Continuous Change Detection and Classification of Landsat Data to Investigate Long-Term Mangrove Dynamics in the Sundarbans Region

Abstract

Mangrove forests play a global role in providing ecosystem goods and services in addition to acting as carbon sinks, and are particularly vulnerable to climate change effects such as rising sea levels and increased salinity. For this reason, accurate long-term monitoring of mangrove ecosystems is vital. However, these ecosystems are extremely dynamic and data frequency is often reduced by cloud cover. The Continuous Change Detection and Classification (CCDC) method has the potential to overcome this by utilising every available observation on a per-pixel basis to build stable season-trend models of the underlying phenology. These models can then be used for land cover classification and to determine greening and browning trends. To demonstrate the utility of this approach, CCDC was applied to a 30-year time series of Landsat data covering an area of mangrove forest known as the Sundarbans. Spanning the delta formed by the confluence of the Ganges, Brahmaputra and Meghna river systems, the Sundarbans is the largest contiguous mangrove forest in the world. CCDC achieved an overall classification accuracy of 94.5% with a 99% confidence of being between 94.2% and 94.8%. Results showed that while mangrove extent in the Sundarbans has remained stable, around 25% of the area experienced an overall negative trend, probably due to the effect of die-back on Heritiera fomes. In addition, dates and magnitudes of change derived from CCDC were used to investigate damage and recovery from a major cyclone; 11% of the Sundarbans was found to have been affected by Cyclone Sidr in 2007, 47.6% of which had not recovered by mid-2018. The results indicate that while the Sundarbans forest is resilient to cyclone events, the long-term degrading effects of climate change could reduce this resilience to critical levels. The proposed methodology, while computationally expensive, also offers means by which the full Landsat archive can be analyzed and interpreted and should be considered for global application to mangrove monitoring.

1. Introduction

Mangroves are salt-tolerant trees which occupy the intertidal zone. Globally mangrove forests cover an area of nearly 140,000 km² [1,2] in over 120 countries and territories [3]. These forests play an important role as carbon sinks [4] in addition to providing many ecosystem goods and services [5]. Occupying a narrow ecological niche, mangroves are particularly vulnerable to climate change effects such as sea level rise, changing ocean currents, and increasing temperatures, which lead to greater erosion, increased salinity, and reduced sediment deposition [6]. Globally, many mangrove ecosystems remain threatened by anthropogenic activity [7]. Understanding the dynamic nature of mangrove ecosystems is vital for both preservation and for utilization of theses ecosystems as climate change markers.

While many studies exist which aim to monitor long-term change in mangrove ecosystems (e.g., [8,9,10,11,12], the focus has remained on single-image classification techniques, typically at yearly resolution or less. In recent years land cover monitoring has seen a move away from such approaches and towards methods which make use of all available observations [13]. Method such as Continuous Change Detection and Classification (CCDC) [14,15] use a per-pixel model fitting approach to capture the seasonal dynamics of land cover types in addition to their inter-year greening and browning trends. One advantage of model fitting is that it reduces reliance on individual observations, instead attempting to capture the broad phenological cycle of the underlying vegetation. It also reduces the need for whole cloud-free images, instead relying on each individual pixel. This is particularly useful for monitoring mangrove ecosystems, which typically exist in tropical and sub-tropical regions where high cloud cover reduces observation frequency [3]. Changes in land cover can then be deduced by finding where in time the fitted season-trend model breaks down. In addition, a classification of the different models can be made based on their various parameters. Unlike a single-image classification approach, this type of classification takes into account the temporal signature of the vegetation rather than spectral values at a single point in time, potentially leading to higher classification accuracy.

Spanning parts of both Bangladesh and India, the Sundarbans is the largest contiguous mangrove forest in the world [16,17,18]. It is an area of high biodiversity [19] which provides many ecosystem goods and services, such as timber for building and fuel [20]; materials for medicines [19]; stabilization of the coastline and protection against extreme events such as cyclones and tsunamis [21]; maintaining water quality [22]; and habitat for a variety of endangered species. More than 5 million people depend on the Sundarbans mangroves for resources [19]. While much of the Sundarbans has been designated a United Nations Educational, Scientific and Cultural Organization (UNESCO) World Heritage Site [23], it remains threatened by anthropogenic activities such as over-exploitation, illegal logging, pollution, and expanding industries such as shrimp farming. In addition, rising sea levels are increasing salinity in the delta, with adverse effects on the dominant mangrove species, Heritiera fomes [24,25]. There is also evidence that the frequency of tropical cyclones making landfall in Bangladesh is increasing [26]. Damage from such cyclones, in combination with rising sea levels, means that the protective capacity of the Sundarbans could be reduced in future years. High importance must therefore be placed on understanding the potential long-term impact of cyclones on the Sundarbans mangroves in addition to tracking mangrove extent.

Given the importance of the Sundarbans, both in terms of providing resources and protecting coastal communities, adequate monitoring of mangrove abundance and health is vital. However, as for many mangrove forests, field studies are expensive and time consuming due to the size of the area covered and difficulty of access [8,24]. The Sundarbans is also an extremely dynamic region which experiences rapid changes due to tides and flooding from tidal surges and seasonal rains, in addition to long-term erosion and accretion of land [8]. These dynamics are difficult to capture on the ground, where a simultaneous snapshot of the entire area is not possible. Remote sensing offers a solution to this. Datasets such as the Landsat archive can provide data going back over multiple decades at a high enough spatiotemporal resolution to identify and track localized changes, both current and historical. In addition, as the price of computing infrastructure continues to fall, the processing of large areas such as the Sundarbans is becoming faster and more viable.

Many studies have used satellite imagery to monitor dynamics in the Sundarbans in recent years. These have covered a wide range of topics including overall mangrove extent [8,9], land cover and species level change [22,24,27], mangrove phenology [28], and coastline change [29]. However, these studies are limited in that they only use one or two images per year [9,22,24,27], do not cover the whole extent of mangroves in the Sundarbans region [22,29], or use data at lower spatial resolution than Landsat [28]. While several studies attempt to quantify damage and recovery from cyclones (e.g., [17,30,31]), the focus on a single image or comparison between very few before/after images is limiting in that it does not take into account the wider long-term variability of the Sundarbans area.

Mangroves typically exists in tropical and sub-tropical regions which are often covered by clouds. This creates difficulty in obtaining cloud-free images from the same time of year which can be compared. For this reason, Cornforth et al. [17] used Synthetic Aperture Radar (SAR) data to study damage and recovery from a major cyclone in the Sundarbans. However, there are disadvantages to the use of SAR imagery. Firstly, no single SAR dataset offers the continuous decades-long coverage of Landsat, making long-term monitoring difficult. Secondly, while SAR data can provide information about tree canopy structure which is absent from optical data, there is evidence that optical data can provide better discrimination between mangroves and other similar tropical vegetation than SAR data [2]. In central areas of mangrove forest where vegetation is overwhelmingly likely to be mangroves, discrimination between land, water, and vegetation is adequate. However, for detecting smaller clusters of mangroves over widers area a more comprehensive classification is needed.

There is clearly scope for a more in-depth study which makes full use of the temporal record for mangrove ecosystem monitoring. In this study, mangrove extent and greening and browning trends in the entire Sundarbans region are examined over thirty years using the Landsat data archive and CCDC. Damage extent and recovery from a major cyclone, Sidr, is also investigated. By taking advantage of all available observations at a per-pixel level, this study aims to provide a more comprehensive and accurate view of mangrove dynamics in the Sundarbans. In doing so, the utility of using a dense time series method for both mangrove classification and monitoring over a large area is demonstrated.

2. Materials and Methods

2.1. Study Area

The study area covers the region from 20°8′ to 24°0′N and 87°4′ to 93°0′E (Figure 1). This covers the entirety of the Sundarbans, approximately two-thirds of which is located in Bangladesh and the remaining third in India [20]. Located in the Bay of Bengal at the confluence of several major river systems (the Ganges, Brahmaputra and Meghna), the Sundarbans is a low-lying area of dense mangrove forest intersected by a complex network of river channels, islands, and mudflats. Elevation is between 0.9 and 2.1 m above Mean Sea Level (MSL) [32]. The climate is tropical, with the monsoon season extending from May to October [24]. The area is periodically hit by cyclones [26], many of which cause substantial losses of both lives and property [33].

The Sundarbans is an area of high biodiversity, home to many rare, threatened, or endangered species [34]. Three main mangrove species dominate: Heritiera fomes, which can only tolerate low water salinity; Excoecaria agallocha, which can grow in moderately saline water; and Ceriops decandra, which can tolerate high salinity [28]. These species have been estimated to make up 33.4%, 30.2%, and 32.4% of mangroves in the Sundarbans, respectively [24]. While there is evidence that numbers of Heritiera fomes and Excoecaria agallocha trees are declining, numbers of Ceriops decandra may be increasing [24].

2.2. Data and Pre-Processing

All United States Geological Survey (USGS) Collection 1 Tier 1 Landsat 4-5 TM, Landsat 7 Enhanced Thematic Mapper Plus (ETM+) and Landsat 8 Operational Land Imager (OLI) data covering the Sundarbans national park region for the period from January 1988 to June 2018 was downloaded from the Google public repository. Collection 1 Tier 1 provides high quality data which has been georegistered and inter-calibrated across the Landsat instruments, and is considered suitable for time series analysis [35]. The Sundarbans is covered by six Landsat scenes: WRS-2 (World Reference System 2) Row 45/Path 137, Row 44/Path 137, Row 44/Path 138, Row 45/Path 138, Row 45/Path 136, and Row 44/Path 136 (Figure 1). The downloaded images were atmospherically and radiometrically corrected converted to analysis ready surface reflectance data using the Atmospheric and Radiometric Correction of Satellite Imagery (ARCSI) Python package [36]. As part of this process cloud masks were also created using the Function of mask (Fmask) algorithm [37].

The data were then ingested into a data cube on Super Computing Wales (SCW). The Open Data Cube (ODC) is an open source program which uses a database along with a Python interface to simplify the processing and organization of geospatial data [38]. Data is tiled and spatially aligned to allow for easy analysis of per-pixel time series. Processing can also be parallelized to facilitate analysis of large scale and even continental scale datasets [39]. As part of the process of data cube ingestion, the data were divided into 105 × 105 km tiles (roughly 1° by 1°) and re-projected to the Bangladesh Transverse Mercator projection (EPSG 3106). It was then stored as Network Common Data Form (NetCDF) files on SCW ready for analysis. Based on work by Zhu et al. [15] only the red, green, near-infrared (NIR), and shortwave infrared (SWIR) bands were used for analysis.

2.3. Classification of Mangroves Using CCDC

2.3.1. The CCDC Algorithm

Once data had been pre-processed, we applied the Continuous Change Detection and Classification (CCDC) method to the Sundarbans area to provide data on the timing and magnitude of class changes (e.g., mangrove to non-mangrove) in addition to details on trend over time (Figure 2). CCDC is designed to work with multi-band Landsat data and focuses on changes in land cover class whereby breaks in the time series are identified and each segment (period between breaks) is classified independently [14]. Breaks are found by fitting a linear model to a stable history period of up to two years. New observations are then added and their residuals compared to the RMSE (Root Mean Square Error) of the history period. If the error of a new observation is too high, it is flagged as a potential change. Once six such observations are flagged in a row a change is identified and a new stable model is initialized. Six observations are used to provide a high confidence of change having occurred [15]. Tropical forests present a challenge for change detection analysis because the prevalence of clouds reduces data frequency. While CCDC has been applied to tropical forest monitoring [40], it has not previously been used for analysis of long-term trends in tropical vegetation or for classifying mangrove forest. However, given the length of the time series used some effects of missing data should be mitigated because observations from different years are unlikely to fall on the same day of year.

CCDC is a computationally expensive method [14]. The area processed is approximately 600 × 380 km in size constituting nearly 250,000,000 pixels. To process the area in a reasonable time frame, the SCW facility was used. This provided access to up to 600 cores simultaneously allowing for parallel processing of pixels. Running CCDC over the whole area took around two weeks, however, 600 cores were not continually available during this period.

The version of CCDC used was that proposed in [15] and was implemented in Python 3.6.8. This version fits a season-trend model of the form given in Equation (1), where $\widehat{\rho }(i,x)$ is the predicted value for the $ith$ Landsat band at Julian date x, $a_{0,i}$ is the coefficient for the mean of the $ith$ Landsat band, $a_{1,i}$ and $b_{1,i}$ are coefficients representing intra-annual change, and $c_{1,i}x$ is the coefficient representing the inter-annual change, or trend [14]. $T=365.25$ (the number of days in a year). CCDC is adaptive, including up to three harmonic (sine/cosine) terms in the model depending on available data quantity [15]. The model is fitted separately to the data for each Landsat band using Lasso regression. Lasso regression minimizes overfitting through regularization by limiting the magnitude of the model coefficients [15]. The degree of regularization depends on the value of the parameter $\lambda$ where $0<\lambda <\infty$. If the value of $\lambda$ is too high, the model will underfit and fail to capture the seasonal nature of the data. If $\lambda$ is too low, the model could overfit, resulting in a higher rate of false changes being detected. The value of $\lambda$ was decided through manual interpretation of models fitted to a sample of mangrove time series where $\lambda$ = 0.01, 0.1, 1, 5, and 20. A value of 1 was chosen in this case because it resulted in models which adequately captured the seasonal cycle of the data while minimizing overfitting.

$$\widehat{\rho }(i,x)=a_{0,i}+a_{1,i}cos\left(\frac{2\pi }{T}x\right)+b_{1,i}sin\left(\frac{2\pi }{T}x\right)+c_{1,i}x$$

(1)

2.3.2. Generation of Yearly Classification Maps

Classification maps identifying pixels as either mangrove, other terrestrial, or water were generated for each year from 1988 to 2017 (inclusive). A classification map was not created for 2018 because data had only been processed up to June 2018. A full year of data was, therefore, not available.

Model Classification

In addition to dates of change for each pixel through time, the algorithm outputs were per-band model coefficients ($a_1,b_1\dots a_3,b_3$), RMSE, and an overall value for each model calculated using the slope and intercept [14]. Given an input of five image bands and a third-order harmonic model, this resulted in a set of 45 variables for each model. These variables were used to classify each model as mangrove, water, or other terrestrial using the method provided by Zhu et al. [14]. The assumption of this method is that different land cover types have different seasonal cycles, resulting in different model coefficients and different model fits. A classifier can therefore be trained to identify different land cover types based on the fitted season-trend model.

First training data was obtained using the Global Mangrove Watch (GMW) Version 2.0 dataset. The GMW provides a highly accurate global mangrove baseline for 2010 [2]. A sample of 50,000 pixels was taken randomly from within the GMW 2010 mangrove area. To obtain a non-mangrove sample, ten polygons of approximately 3000 km² in area were then chosen from around the Sundarbans region to maximize the distribution of the training data (Figure 1). Data were not taken from the entire non-GMW region due to the risk of including mangrove pixels not identified by the GMW baseline. A sample of 50,000 pixels was then randomly selected from each polygon to obtain a total sample of 500,000 pixels. The sample size was large enough to ensure that most dominant land cover classes were represented. Once the training pixels had been selected, the results of the change detection step were used to identify whether those pixels had fitted models suitable for use as training data.

Each model generated by CCDC covers a specific period of time for a specific pixel. Therefore, in order to classify models based on the 2010 dataset, only models which covered the whole of 2010 were chosen as training data because a break detected in 2010 could signal a change in land cover type. There would then be no way of knowing which model aligned with the mangrove/non-mangrove classification provided by GMW. If a model covered a period starting before 2010 and ending after 2010, that model likely represented a mangrove/non-mangrove signal as defined by the 2010 data. The non-mangrove training models were further divided into land and water models using a land/water mask for 2010. This gave final final training datasets of 49,923 mangrove models, 87,736 water models, and 403,088 other terrestrial models.

Once a set of stable models was identified, the model parameters were used to train a Random Forest (RF) classifier implemented using the scikit-learn Python library [41]. The same model parameters suggested by Zhu et al. [14] were used, however in almost all cases CCDC fitted a third-order harmonic model resulting in a set of nine variables per band or 45 total variables used for classification per model. In the rare cases where a less complex model was fitted the extra coefficients were set to zero to maintain compatibility with the classifier. To find the best hyperparameters for the classifier, a random grid search with cross validation was used to narrow down the potential range of values. A second grid search was then performed to produce the final set of hyperparameters. Once the classifier was trained, all models produced as output from running the CCDC algorithm per-pixel were classified.

Yearly Maps

The output of the CCDC algorithm resulted in a set of classified models per-pixel, each with a start and end date. Some models covered only a few years whereas others covered the entire time period. To generate yearly classification maps it was necessary to summarize these models into a yearly class for each pixel. To achieve this the majority class within a given year was used. If there was no majority, the pixel was given a value of 0 (not enough data) for that year. Mangrove pixels were given a value of 1, other terrestrial pixels a value of 2, and water pixels a value of 3. The resulting yearly maps were then processed using the Shuttle Radar Topography Mission (SRTM) global 1 arc second product to remove any mangrove pixels above 30 m. Mangroves only grow at low elevations in tidal and intertidal zones [9], so any pixels above this height were highly likely to be miss-classifications. Finally, any clumps of land cover less than 10 pixels (0.9 ha) in size were removed and replaced with the value of the largest neighboring class using the Geospatial Data Abstraction Library (GDAL) [42], in line with recommendations by Bunting et al. for reducing error associated with small-scale features [2]. Yearly mangrove extent in km² was then calculated based on the number of pixels in each yearly map assigned to class 1.

Validation

Due to the spatiotemporal extent of the data, validation using field data or a higher resolution dataset was not possible. No other dataset is available which covers the same time period as Landsat at similar or higher spatial resolution. Following the same methodology as the GMW [2], a validation set was therefore generated by taking a randomly selected set of 47,000 pixels across both space and time and interpreting them manually using Landsat data. Reference to up to date high resolution Google Earth imagery was made and the Landsat imagery was displayed using the NIR, SWIR1 and Red colour composite, in which the mangroves are spectrally distinct (see Figures 4 and 5), to aid identification of mangroves from other terrestrial land covers.

The CCDC classification resulted in some pixels for which no overall land cover class could be determined (i.e., there was no majority class in the year) so the final number of validation points used was 45,440. User’s and producer’s accuracy were calculated for each of the three classes along with overall accuracy. Quantity disagreement and allocation disagreement were calculated as described by Pontius and Millones [43].

2.4. Long Term Vegetation Trends in the Sundarbans

As well as determining yearly classifications, the models created by the CCDC process were used to generate a mangrove trend map covering January 1988 to June 2018. The purpose of this was to summarize greening and browning trends for the Sundarbans mangroves both spatially and temporally. Using the yearly classification maps, all pixels which had been classified as mangrove at any point between 1988 and 2017 were identified. For each pixel, the overall Normalized Difference Vegetation Index (NDVI) trend was then calculated using the method described by Zhu et al. [44] over all models for that pixel. Briefly, this method uses the intercept and slope of the fitted seasonal models for the red and NIR bands to calculate an overall NDVI value for the start and end date of each model. Overall trend is then calculated by adding up the within-model and between-model differences.

2.5. Investigation of Dynamics around Cyclone Sidr

In addition to investigating the accuracy of CCDC for mangrove classification in the Sundarbans, a specific use case was also desirable to illustrate the potential for more in-depth investigation of vegetation dynamics. Given the computational complexity of CCDC, its usefulness for many applications could be increased if multiple analyses can be carried out using the same results. In this case disturbance and recovery dynamics around Cyclone Sidr were investigated as the final stage of analysis (Figure 2).

Cyclone Sidr was a category 5 storm which made landfall in Bangladesh on the 15th of November 2007. It caused a 3 to 4 m tidal surge along with wind speeds up to 220 km/h [30]. Sidr was chosen as a use case for several reasons. Firstly, the cyclone hit the Sundarbans area directly [33]. Secondly, Sidr hit in 2007, meaning that data were available for multiple years both before and after it made landfall. Thirdly, estimates of damage done by Sidr around the time of impact range widely from around 22% to 45% of mangroves [17,31] and no research was found into the long-term effects of cyclone Sidr beyond 2010.

In order to investigate the possible long-term effects of Sidr, all pixels classified by CCDC as mangrove in 2006 were selected. From this set all pixels were selected for which CCDC recorded a break in the two months following Sidr. Two months was chosen to provide a reasonable level of certainty that the change was attributable to Sidr rather than other causes, given the lack of available validation data. The magnitude of the reported break was recorded along with the classification of the first post-Sidr model and the number of days between the end of the pre-Sidr model and the start of the first post-Sidr model. The trend method outlined in Section 2.4 was then used to estimate an overall NDVI value for the end of the final pre-Sidr model, i.e., a greenness value for the pixel before the cyclone hit. The trend for each subsequent model was then tracked to find the first year in which the overall NDVI value reached or exceeded the final pre-Sidr value. This year was then recorded as the year of recovery from Sidr. If the pixel had not recovered by the end of the available data in mid-2018, it was recorded as Not Yet Recovered (NYR).

3. Results

3.1. Classification of Mangroves Using CCDC

The classification achieved an overall accuracy of 94.5% with a 99% confidence of being between 94.2% and 94.8% (

Table 1. Results of spatiotemporal accuracy assessment for CCDC.

		Reference
		Mangrove	Water	Other	Total	User’s (%)
Classifier	Mangrove	11,835	531	219	12,585	94.0
	Water	174	22,422	367	22,963	97.6
	Other	684	609	8599	9892	86.9
	Total	12,693	23,562	9185	45,440
	Producer’s (%)	93.2	95.2	93.6		94.5

). Quantity disagreement was 0.02 and allocation disagreement was 0.04. Mangroves were most frequently confused with other terrestrial pixels. Other terrestrial pixels were most often confused with water.

Based on the yearly classification maps, the year with the highest mangrove extent was 2004 with mangroves covering approximately 6672.6 km². The year with the lowest mangrove extent was 1988 when mangroves covered approximately 6368.3 km². The lower values for 1988 and 1989 (Figure 3) are probably because Bangladesh experienced both catastrophic flooding and a severe cyclone in 1988. These events caused rapid changes in land cover (i.e., water inundation) and destruction to vegetation making it difficult to fit a stable model. Mangrove extent also increased substantially between 1989 and 1990 (Figure 3), possibly because the forest was still recovering. Excluding 1988 and 1989, the year of lowest extent was 2007 (6576.0 km²). Mangrove extent increased by 15.1 km² between 1996 and 2006 (0.23%), the ten years before Sidr, and decreased by 40.9 km² between 2008 and 2017 (0.61%). However, the majority of this decrease was between 2013 and 2017, with extent only decreasing by 0.03% between 2008 and 2013. In the last 20 years, between 1997 and 2017, mangrove extent reduced by 38.4 km² (0.58%). All of these figures are lower than the classification uncertainty and should therefore be interpreted cautiously. A full table of yearly extent change can be found in

Table A1. Mangrove area lost and gained for each year in comparison to the previous year (km²).

Year	Change	Net Loss	Net Gain	Total
1988	-	-	-	6368.3
1989	39.9	86.7	126.6	6408.2
1990	171.8	19.8	191.6	6580.0
1991	32.6	11.9	44.5	6612.6
1992	15.2	5.6	20.8	6627.8
1993	7.7	3.7	11.4	6635.5
1994	2.9	4.5	7.4	6638.4
1995	10.7	4.7	15.4	6649.1
1996	2.0	5.7	7.7	6651.1
1997	5.0	7.0	12.0	6656.1
1998	4.4	6.7	11.1	6660.5
1999	2.9	7.0	9.9	6663.4
2000	2.5	6.5	9.0	6665.9
2001	0.4	10.4	10.8	6666.3
2002	−0.7	12.3	11.6	6665.6
2003	1.5	12.9	14.4	6667.1
2004	5.5	11.9	17.4	6672.6
2005	−3.5	12.4	8.9	6669.1
2006	−3.0	10.6	7.6	6666.1
2007	−90.1	105.0	14.9	6576.0
2008	82.5	21.7	104.2	6658.5
2009	−2.5	17.4	14.9	6656.0
2010	9.2	9.6	18.8	6665.2
2011	−2.0	6.4	4.4	6663.2
2012	−3.2	7.0	3.8	6660.0
2013	−3.8	9.6	5.8	6656.2
2014	−7.1	14.1	7.0	6649.1
2015	−3.4	10.0	6.6	6645.7
2016	−40.6	43.8	3.2	6605.1
2017	12.5	10.7	23.2	6617.6

.

For 2010 there was substantial overlap between CCDC and the GMW. Only 1% of mangrove area captured by the GMW was not captured by CCDC. However, mangrove extent derived using CCDC for 2010 (6665.2 km²) was much higher than estimated by the GMW (6016.4 km²). 10.7% of the mangrove area captured by CCDC was not classified as mangroves by the GMW, indicating that the CCDC approach was able to capture areas of mangrove which were missed by the GMW (Figure 4). The CCDC classification was also able to capture long-term changes in mangrove extent caused by land erosion and accretion (Figure 5). A specific example is shown in Figure 6D of a pixel which developed from water to mangrove and was robust enough to recover from a major disturbance event.

3.2. Long Term Vegetation Trends in the Sundarbans

Approximately 0.002% of pixels reported an unrealistic overall trend (trend < −2 or trend > 2) and were removed from the analysis. These trends were probably caused by outliers (e.g., very high or negative values) in the red and NIR bands. Of all remaining pixels classified as mangrove at any point between 1988 and 2017 (6865.6 km²), 73.5% experienced an overall positive trend in NDVI (5044.1 km²) and 26.5% experienced an overall negative trend (1821.3 km²). The median positive trend value was 0.06 with 24.6% of positive trend pixels (18.1% of all pixels) exhibiting a trend >= 0.1 and 4.5% of positive trend pixels (3.3% of all pixels) exhibiting a trend >= 0.5. The median negative trend value was −0.03 with 19.5% of negative trend pixels (5.2% of all pixels) exhibiting a trend <= −0.1 and 7.8% of negative trend pixels (2.1% of all pixels) exhibiting a trend <= −0.5. These figures indicate that more greening than browning is occurring across the Sundarbans mangroves. Spatially, stronger greening and browning trends occur around the coastline with a more generalized area of browning also apparent to the north (Figure 7).

3.3. Investigation of Dynamics around Cyclone Sidr

An estimated 10.7% of mangroves were damaged by Cyclone Sidr based on change from the 2006 classification. Damaged areas were clustered around the coastline and to the west of the cyclone’s path (Figure 8). The median NDVI break magnitude was −0.1 with an interquartile range of 0.05. Break magnitude cannot be directly related to the amount of damage caused on the ground. However, 51.9% of pixels registered a break magnitude $>=-0.1$, indicating that the majority of disturbances were small (Figure 6C). Median days disturbed (days before CCDC was able to fit a new stable model) was 80, indicating that land cover often took more than two months to recover to a a level stable enough for a new model to be fitted. Minimum days disturbed was one. This was possible because a small amount of overlap between Landsat 5 and Landsat 7 scenes in January 2008 meant that some pixels had observations on two consecutive dates. Maximum days disturbed was 1040. 88.6% of pixels had a new model fitted within 180 days.

95.7% of damaged pixels were classified as mangrove once a new model had been fitted (684.0 km²) (e.g., Figure 6B), 4% were classified as other terrestrial (28.3 km²), and 0.3% were classified as water (2.4 km²) (e.g., Figure 6A). This suggests that the majority of damaged areas remained dominated by mangroves despite the storm damage.

47.6% of damaged pixels had not reached pre-Sidr NDVI levels by mid-2018 (Figure 8) with the majority of recovery occurring between 2013 and 2018 (Figure 9). Data were not normally distributed so a Spearman’s rank test was used to test for correlation between break magnitude and year of recovery. This resulted in a significant trend ($p<0.01$) with a value of $Rho=-0.24$, indicating a small negative trend (Figure 9). The trend remained similar with and without inclusion of the NYR category, which contained nearly 50% of the data.

4. Discussion

4.1. Classification of Mangroves Using CCDC

The CCDC method proved to be highly accurate for mangrove classification. This demonstrates that CCDC can provide a tractable solution for quantifying change in highly dynamic regions, and can be applied to areas such as the tropics which have high cloud prevalence. Allocation disagreement was higher than quantity disagreement, indicating more error in the spatial distribution of land cover classes than in the quantity mapped as each class. The main source of error was confusion with other terrestrial vegetation. This is probably because along with other types of tropical vegetation, mangroves do not have a very pronounced seasonal cycle. However, other comparable methods using optical imagery and single-date classifications report overall accuracies in the range of 64–88% [9,31,45]. Including temporal information in the classification process in the form of time series models therefore represents a reasonable improvement over single-image approaches. Possible further improvements could be made by including auxiliary data such as elevation or distance from water in the classification step, rather than in post-processing.

The lack of an independent dataset for verification is a limitation of this study. However, CCDC has previously been shown to be a accurate method of determining land cover type [14,44]. Substantial overlap was also found between the 2010 classification map generated by this study and the GMW estimate, and we found evidence that CCDC was also able to distinguish some smaller areas of mangroves missed by the GMW. Our classification estimates the total area of mangroves in the Sundarbans in the present day to be around 6600 km². While the total area of the Sundarbans is often cited as around 10,000 km² [8,16,24], recent studies using Landsat data place the estimate of actual mangrove forest closer to 6000–7000 km² [8,9,24]. Our estimate is, therefore, in line with previous studies, and demonstrates that evaluation of baseline extents of important habitats can be more accurate with new techniques, particularly if the full data archive is exploited.

The lower estimates of mangrove extent for 1988 and 1989 indicate that not enough data were available to determine land cover for many pixels. There are two factors which could have contributed to this. Firstly, 1988 was the first year of data used and classifications for this year might have been affected by a spin-up effect of the model fitting process. If a stable history period could not be initialized by mid-1988 then a pixel would not have a majority class for that year. Secondly, Bangladesh experienced extreme flooding in both 1987 and 1988 [46] exacerbated by a category 3 tropical cyclone making landfall in November 1988 which directly impacted the Sundarbans [47]. This may be why mangrove extent was also lower for 1989. These extreme weather events would have caused rapid changes in land cover and made it more difficult to fit a stable model. The drop in mangrove extent between 2006 and 2007 is discussed in Section 4.3.

The classification also indicates that there has been a slight decrease in mangrove extent over the last 20 years. This result is broadly in line with a study by Giri et al. [9] which found a slight decrease in mangrove forest in the Sundarbans between the 1990s and the 2000s and a slight decrease in mangrove forest over all of South Asia between 2000 and 2012. Other studies also report little or no long-term decrease in mangrove extent [8,24]. In both this study and that by Giri et al. the detected decrease is within the margin of error of the classification and therefore must be interpreted cautiously. Mangrove dynamics in the Sundarbans are also complex, with both land erosion and accretion contributing to mangrove numbers [9,28]. While evidence was found that mangrove quantity was increasing before Sidr but decreasing afterwards, the majority of that reduction occurred in the last five years. In addition, the number of pixels classed as non-mangrove after Sidr which were previously classed as mangrove was 4.3%, much higher than the estimated reduction in mangrove extent between 2008 and 2017. This suggests that even where mangroves were damaged badly enough to be classified as non-mangrove for a time, most either eventually re-established as mangroves or the loss was offset by gains elsewhere. This demonstrates the resilience of these habitats to events such as cyclones. However, long-term degradation and subsequent fragmentation could reduce this resilience to critical levels.

4.2. Long Term Vegetation Trends in the Sundarbans

Vegetation dynamics in the Sundarbans are complex, particularly when examined over a long period of time. By examining overall greening and browning trends, it was hoped that areas of the Sundarbans at potential risk of degradation in the future would be highlighted. Results showed that three-quarters of pixels identified as mangroves at some point during the 30-year period from 1988 to 2017 exhibited a positive trend over their entire time series. Of these areas, some are clearly regions of land accretion where new stands of mangroves have become established (Figure 7). Given that mangrove extent changed little over the study period, many of these areas of accretion must be being counterbalanced by mangrove loss. However, approximately three times more strong positive trends were found than strong negative trends. This suggests that over the majority of the Sundarbans, mangrove canopy cover has increased over the past 30 years. Excluding establishing mangrove stands, mangroves exhibiting a net increase in greenness could fall into one of several categories. Firstly, they could have experienced damage (e.g., from a cyclone) but recovered past their pre-damage greenness level. Secondly, they could have experienced damage prior to 1988. In this case, even if the mangroves did not recover to their pre-damage level, the trend would still present as being net positive. Finally, an increase in NDVI could be the result of a change in mangrove species composition. C. decandra, E. agollocha, and H. fomes often coexist in the Sundarbans, with the shorter C. decandra often forming as an understory to E. agollocha [24]. Decline in one species could therefore be masked in terms of NDVI trend if another species rapidly takes over and establishes.

Approximately one quarter of pixels classified as mangroves at some point in the 30-year period were found to have an overall negative trend. This suggests that even though there is very little net loss of mangroves, a substantial proportion may be at risk. In particular, there is a clear area to the north where negative trends are prevalent. A study into long-term mangrove species dynamics by Ghosh et al. found that this area has predominately been occupied by H. fomes, which has a low tolerance for saltwater. Die-back disease has been a large-scale problem since 1980 and was estimated to affect 5–6% of H. fomes in 2010. It is, therefore, possible that a large proportion of the negative trends found were in H. fomes. This highlights the particular vulnerability of these important habitats to global climate change and subsequent sea level rise, as well as increasing number of extreme weather events. As such, monitoring mangroves in this way can provide a barometer for regional and global vegetation changes.

4.3. Investigation of Dynamics around Cyclone Sidr

4.3.1. Comparison to Previous Damage Estimates

Results showed that around 11% of the Sundarbans mangroves were damaged by cyclone Sidr. This is lower than other estimates. However, it is important to highlight that the approach used in this study is substantially different to other studies, which used relatively few images (e.g., [17,30,31]). By fitting temporal models rather than relying on single-image classifications [30] or differencing between images [31], mangrove dynamics both before and after the cyclone could be more accurately captured. This creates difficulty in comparing the damage estimates between CCDC and the other methods. In addition, this study covered the entire Sundarbans region, including many outlying areas not covered by other studies. This likely reduced the damage estimate as a percentage.

Akhter et al. [30] used Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery from November 2007 to asses damage immediately following the cyclone. Based on field visits, the study area was classified into regions which were highly, moderately, and severely affected by Cyclone Sidr with the total damaged area estimated to be 22.2%. However, this study only focused on a small area of the forest which was directly impacted. Work by Bhowmik and Cabral [31] classified four Landsat 7 scenes from February 2007–2010 to distinguish between different mangrove species before using NDVI values to track disturbance and recovery from Sidr. This placed the estimate of damage at around 45% of the Sundarbans, with recovery being complete by 2010. However, this assessment was based on more generalized reductions and increases in relative NDVI values across images, and still only covered the area directly impacted by Sidr. Cornforth et al. [17] used backscatter values from radar imagery as a proxy for above ground biomass, where a reduction in backscatter indicates a reduction in biomass and therefore a change in forest condition. Imagery from 2007 was compared to images from 2008 and 2009, though due to a lack of available imagery, less than half of the total area of the Sundarbans was assessed. This study reported lower damage estimates than other studies which used optical data. However, their approach was still based on differencing between single images from different years and using a threshold to identify damaged areas. This meant that determination of changed areas was entirely dependent on relative change from a single previous time point. By fitting a stable temporal model to optical data, some of the problems caused by noisy or missing data can be eliminated by capturing the broader seasonal dynamics of a land cover type. Change analysis carried out this way is far less dependent on the value of individual observations.

This study also took a conservative approach to damage estimation by only looking for changes in the two months following Sidr. This was based on the fact that it is difficult to validate what truly caused a change. In addition, there were limitations on allocating the true date of change due to data availability. The only data available between September 2007 to February 2008 for the area directly impacted by Sidr was from the Landsat 7 satellite, and no data at all was available for June, July, or August 2007. This meant that the SLC failure on board Landsat 7 contributed to the lower estimate given by CCDC. While season-trend modelling approaches such as CCDC can interpolate between observations to fill in gaps in the data record, change cannot be assigned until an actual observation is available. For pixels with missing data due to the SLC failure, no observation was available until March 2008, meaning that if those pixels were damaged by Sidr, the change would not be detected by CCDC until nearly four months after Sidr made landfall. As two months was used as the cut-off for changes to be reliably assigned to Sidr, these damaged pixels were missed leading to gaps in the estimate of damage extent (Figure A1). This problem was difficult to avoid without substantially extending the range within which a change was considered attributable to Sidr, and therefore decreasing the likelihood of Sidr being the true cause of change. 11% is, therefore, a robust minimum estimate of the area of mangroves damaged by Sidr. Estimates based around comparing yearly images (e.g., [31]) are more likely to attribute changes to Sidr which may have other causes, and to be more affected by noise and cloud contamination.

4.3.2. Impact of Sidr on Mangrove Extent

The impact of Sidr contributed to a drop of 90.1 km² in mangrove extent between 2006 and 2007 (Table A1), despite the fact that the vast majority of damaged pixels were classified as mangrove before and after Sidr. This is because classification is based on a majority class in a given year. For some pixels damaged by Sidr, valid data was only available for May 2007 and March 2008. In these cases, a change would be detected in March 2008 once Landsat 5 data was available, with the end date of the last stable model being May 2007. Since this model covers less than half the year, these pixels would not have a majority class for 2007.

The attributed class of most pixels damaged by Sidr did not change, but the damage done was substantial enough to be flagged by the CCDC method. This suggests that the damage done by Sidr was extreme enough to emulate a change in land cover class, but that pixel trajectories after Sidr were still similar enough to other mangrove pixels to still be classified as mangrove. Therefore the CCDC method is probably capable of further separating the mangrove class into stable, recovering, and degrading mangroves. It is possible that as well as experiencing an increasing or decreasing trend, non-stable mangroves also have more subtle differences in seasonality such as smaller seasonal peaks which would make them distinguishable from the other classes. However, this type of classification would require accurate training data on recovering and declining mangrove populations and was beyond the scope of this study.

Spatially the areas that experienced the highest magnitude of damage were along the coast (Figure 8). In particular, the edges of more fragmented areas and smaller islands were more likely to experience greater damage. These areas would have been vulnerable to the storm surge and more likely to be already experiencing erosion or to consist of land which was still in the process of establishing. Given that the majority of breaks were fairly small in magnitude, fairly rapid recovery would be expected. Bhowmik and Cabral [31] investigated the impact of Sidr on biodiversity in the Sundarbans and concluded that the Sundarbans had recovered to a satisfactory level by 2010. Mangroves are also known to recover quickly once primary regeneration has taken place [31] and to be very resilient to many ecological changes [3].

4.3.3. Estimation of Recovery from Sidr

Median days disturbed was 80 indicating that the CCDC algorithm was able to fit a new stable model to many pixels within three months. The majority were stable within six months. Given that CCDC is capable of identifying a recovering mangrove pixel just as well as a stable mangrove pixel, these periods of disturbance indicate that many mangroves were damaged enough by Sidr that several months passed before they were once again recognizable as mangroves. This also suggests that fairly small magnitude changes can represent substantial disturbances on the ground.

However, even once new models were fitted, nearly 50% of damaged pixels had not recovered to their pre-Sidr overall NDVI value by mid-2018. While mangrove extent did drop between 2008 and 2017, the decrease was small, indicating that the vast majority of these pixels remained mangrove throughout that period. The number of pixels recovered increased much more rapidly between 2013 and 2018 than between 2007 and 2013 (Figure 9). A possible reason for this could be a protective effect whereby once some mangroves have recovered others are less exposed and recover more quickly.

Another possible explanation is bias in the Landsat 8 NDVI trajectories. There is evidence that Landsat 8 NDVI is positively biased compared to Landsat 5 and 7 data [44,48]. Zhu et al. [44] concluded that this was due to the atmospheric correction method used for Landsat 8. Given that all data was pre-processed using the same method (ARCSI) this was not expected to be an issue. However, the increase in recovered pixels does present around the time of Landsat 8’s launch and differences between the sensors cannot be discounted as a possibility. Even if bias exists in the NDVI trends, it is unexpected that so few mangroves appear to have recovered fully by mid-2018. This suggests that although mangroves generally recover quickly, and may appear to have recovered fully based on yearly observations, full recovery following a disturbance event is much slower. If 11% is assumed to be a reasonable minimum area of mangroves damaged by Sidr, the recovery trends indicate that around 5–6% of the Sundarbans mangroves are still recovering from the effects of cyclone Sidr in some way. These mangroves may be more vulnerable to further cyclones and to other environmental changes such as rising sea levels.

5. Conclusions

The methods outlined in this paper provide a tractable solution to exploiting full time series of available Landsat data to map mangrove extent and condition at large scale. Specifically, this was demonstrated for the first time in the Sundarbans, a highly dynamic ecosystem with long term trend changes as well as step changes and which is in a region prone to persistent cloud cover and associated missing data issues. In this example, mangroves were shown to be resilient in their recovery from large cyclonic events, but the approach was also able to detect longer-term decreases in greenness potentially related to climate change effects such as increased salinity caused by rising sea levels. Monitoring mangroves in this manner can therefore present important evidence as a barometer for climate change. By applying this approach to the the entire Sundarbans region, this study demonstrates that CCDC can provide data to help drive regional or even countrywide scale policies as part of strategies such as REDD+ [49] and blue carbon initiatives, informing future decisions as to where conservation efforts in the Sundarbans and other mangrove ecosystems should be focused.