Remote Sensing of Environment

The Australian Geoscience Data Cube — Foundations and lessons learned Adam Lewis ⁎, Simon Oliver , Leo Lymburner , Ben Evans , Lesley Wyborn , Norman Mueller , Gregory Raevksi , Jeremy Hooke , Rob Woodcock , Joshua Sixsmith , Wenjun Wu , Peter Tan , Fuqin Li , Brian Killough , Stuart Minchin , Dale Roberts , Damien Ayers , Biswajit Bala , John Dwyer , Arnold Dekker , Trevor Dhu , Andrew Hicks , Alex Ip , Matt Purss , Clare Richards , Stephen Sagar , Claire Trenham , Peter Wang , Lan-Wei Wang a


Introduction
Remote sensing can provide landscape-level analysis of vegetation dynamics (Goetz et al., 2009), aiding in broad-scale modeling of phenology, primary production and climate change effects Turner et al., 2004). In coastal marshes particularly, vegetation helps contribute to "blue carbon" in soils through both autochthonous production and sediment trapping (Mcleod et al., 2011), and remote monitoring can help quantify the carbon sequestration potential of these systems (Byrd et al., 2014;Stralberg et al., 2011). However, remote sensing in coastal marshes is complicated by periodic tidal flooding, which dampens spectral reflectance (Cho et al., 2008) and introduces noise that obscures the shape and magnitude of vegetation patterns. Such noise results from water and moist soil, which differentially attenuate the intensity of visible (VIS), near-infrared (NIR) and short-wave infrared (SWIR) spectra (Kearney et al., 2009;Mishra et al., 2012Turpie, 2013). Consequently, vegetation time-series based on traditional vegetation indices, such as the Normalized Difference Vegetation Index (NDVI) or the Wide Dynamic Range Vegetation Index (WDRVI) have a reduced red-NIR contrast and a lower overall range compared to terrestrial systems (Fig. 1). As a result, tidal effects can decouple vegetation indices from vegetation patterns (Cho et al., 2008). An inundation index is therefore needed to identify flooded observations so we can use remote sensing to effectively model true tidal marsh vegetation dynamics.
In order to separate water from other land-uses, researchers generally apply one of a suite of normalized difference water indices (NDWI); for a review see Ji et al. (2009). These NDWI-type indices are successful at mapping water bodies (Ji et al., 2009) and the best indices for mapping are those with similar index values for vegetation and soil (Ji et al., 2009). Optimization of thresholds for NDWI-type indices have generally been established with data from a single, short time period and were most successful at discriminating pixels with 100% water cover (Ji et al., 2009). We found these indices were less accurate for mapping periodic tidal inundation across multi-temporal vegetated marsh datasets. As the annual cycle progresses, the relative difference between water and vegetation also changes, requiring an index that accounts for seasonal phenology.
We built on the existing NDWI indices to create an optimized index for pre-processing MODIS (Moderate Resolution Imaging Spectroradiometer) satellite data in tidal marshes. Our goal was to develop an index that makes filtering flooded pixels from multispectral remote sensing data straightforward. We focused our efforts on the MODIS sensor, curated by NASA (the National Aeronautics and Space Administration). We used MODIS because it is freely available, has global coverage, a 1 day return frequency, and provides many useful products (modis.gsfc.nasa.gov), making this platform popular for broad-scale vegetation analysis. The daily 500-m MODIS surface reflectance product (MOD09GA) is especially useful along coasts, where frequent cloud cover interferes with space-to-ground spectral detection. Although MODIS surface reflectance-derived products are widely used in terrestrial systems, they are often inaccurate within tidal marshes because current compositing algorithms do not account for tidal flooding (Huete et al., 1999). These existing products include MOD15 (the leaf area index (LAI) and fractional photosynthetically active radiation (fPAR) product), MOD13 (NDVI and enhanced vegetation index (EVI) product), and MOD17 (the net and gross primary productivity product (NPP/GPP). This also means that independent wetland productivity models relying on these products as input data suffer from tidal fluctuation uncertainty (Barr et al., 2013). Additionally, studies that use raw surface reflectance for habitat mapping and vegetation dynamics are often forced to resort to complex and tedious procedures to filter inundated pixels (Knight et al., 2006;O'Donnell and Schalles, 2016;Wang et al., 2012).
In this paper, we describe the development and application of a new index for pre-processing MODIS data in coastal marshes, which we call the Tidal Marsh Inundation Index (TMII). We calibrated and validated the TMII using ground-truth data from a tidal Spartina alterniflora salt marsh located in Sapelo Island, GA, USA. We then applied the optimized TMII to two additional S. alterniflora MODIS pixels on Sapelo Remote Sensing of Environment 201 (2017) 34-46 Island as well as a S. alterniflora pixel on Plum Island, MA, USA. We demonstrated how to apply and optimize the TMII to produce 16-d composited vegetation time-series at the S. alterniflora sites in GA and MA, as well as two additional sites with new marsh species on the Gulf coast: a Spartina patens marsh in Louisiana and a Juncus roemerianus marsh in Mississippi. In all cases, we compared TMII-filtered composites to 16-d NDVI MOD13 composites. The TMII represents a step forward for coastal remote sensing that will be useful for improved phenology, biomass and carbon estimation in coastal marshes.

Study pixel footprints
The footprint of the MODIS pixel used for TMII development is located on Sapelo Island, GA, within the Georgia Coastal Ecosystems (GCE) Long Term Ecological Research (LTER) site and the Sapelo Island National Estuarine Research Reserve (NERR), hereafter referred to as "Sapelo Island" (lat: 31.441°, long: − 81.284°; Fig. 2, MODIS pixel A). Tides are semi-diurnal and high tides > 2 m are common (tidesandcurrents.noaa.gov). The marsh platform in this area is a S. alterniflora monoculture traversed by tidal channels. Within pixel A, the GCE-LTER maintains an eddy covariance carbon flux tower that measures exchange of carbon dioxide, water and energy between the marsh surface and the atmosphere.
Once the TMII was optimized, we applied it to two new MODIS pixels on Sapelo Island (Fig. 2,

Analytical approach
To create and apply the TMII, we followed the steps outlined in Fig. 3. First we obtained and pre-processed daily MODIS 500 m spectral reflectance data (steps 1a-1d). Next, we needed an accurate record of tidal inundation within the model development pixel footprint (step 2a). We then identified predictors (band combinations) that could explain the response variable (marsh flood status) to create a marsh flood classification model (steps 2b-2f). Next, we assessed the accuracy of TMII across testing and validation in comparison to existing indices (step 2g) and then applied the TMII to new MODIS pixels (steps 2h). Finally, we demonstrated the generality and steps to apply TMII by creating time-series vegetation composites in additional coastal marshes (steps 3a-3f). We used R version 3.3.3 to conduct all analyses and data processing steps.

Acquisition and pre-processing of daily MODIS imagery
Our first step was to develop a record of MODIS 500-m spectral reflectance that corresponded to the inundation data. We acquired daily MODIS TERRA version 5 products (Table 1) and followed steps 1a-1d in Fig. 3. Specifically, we extracted product layer values for each study pixel and then used the "reflectance data state" and "view zenith" layers from MODIS dataset A to set as missing data observations contaminated by clouds or cloud shadows, observations with poor data quality flags, or observations with high view zenith values (> 50°). We removed data with high view zenith values because Huete et al. (1999) suggested atmospherically corrected off-nadir pixels have vegetation indices that are artificially increased due to the coarse spatial resolution and spatial distortion. Thus, off-nadir pixels have more vegetation in the sensor instantaneous field of view. The MODIS science team flags high view zenith observations as poor data quality and minimizes their inclusion within vegetation composites. We used the satellite acquisition time ("Day view time" layer) from MODIS dataset B to temporally merge MODIS spectral reflectance datasets with our inundation groundtruth data. View time stored in this layer is the apparent solar time, based on the sun's position in the sky. We converted apparent solar time to local standard time (see Holbert, 2011 and Appendix A) and rounded to the nearest half-hour, matching the format of our inundation groundtruth data from the PhenoCam (described below). These satellite acquisition time data come from MOD11A1, whose 1-km spatial grid differs from the 500-m spatial grid of the MOD09GA spectral reflectance data we used to create the TMII. However, the acquisition time of MODIS pixels changes slowly at the landscape scale because MODIS collects measurements close to solar noon. One can therefore assume that acquisition time from the 1-km MOD11A1 pixel is representative of the acquisition time of all the 500-m MOD09GA pixels nested within it.

Accurate tidal marsh inundation ground-truth data
Next, we identified accurate inundation data to use for model development and ground-truthing (step 2a). We had access to data from ancillary sensors associated with the GCE-LTER flux tower, including a pressure transducer (Campbell Scientific Instruments Model CS455, Logan UT, USA), which measures water height in an adjacent tidal creek, and a 3D sonic anemometer (Campbell Scientific Instruments Model CSAT3, Logan UT, USA), which records wind direction and speed at 5 min intervals. The tower also has a PhenoCam (StarDot NetCam SC 5MP IR, StarDot Technologies, Buena Park, CA USA), which is part of the National PhenoCam network (phenocam.unh.edu). The PhenoCam is a digital camera that faces northeast and auto-collects fixed-scene oblique angle images of the marsh platform every half-hour during daylight (Fig. 4). Images from Sapelo Island are archived at phenocam. unh.edu under the site label "GCESapelo". Most of the flux tower datasets were available from 17 September 2013 through 31 December 2016, when this analysis was complete; PhenoCam data were missing from 16 July 2014 through 02 September 2014 because of a mechanical failure.
To estimate marsh flooding within the footprint of Sapelo Island pixel A, we first combined information from the PhenoCam with the creek water height and weather information. The PhenoCam provides a photographic record of the aerial percentage of marsh flooding during daylight hours. To extract flood information from PhenoCam images, we used our previously developed algorithm, called the "smart classifier", which quantifies percent cover and presence of water with 82-92% and 86-92% accuracy, respectively (O'Connell and Alber, 2016). We used generalized linear models (GLMs) with binomial error (Zuur et al., 2009) to relate smart-classifier estimated flooding extent to creek water height. We then refined this analysis by including creek water height, wind speed, wind direction, month and year as additional marsh flooding predictors. We used explained deviance, a log-likelihood proxy for explained variance (Zuur et al., 2009) to estimate the strength of the relationship. Explained deviance can be calculated as: ((Null deviance − Residual deviance) / Null deviance * 100) (Zuur et al., 2009). In the end, we found that degree of marsh flooding was not well-related to creek water height and that the PhenoCam record alone was the best variable for calibrating the MODIS TMII.

Identifying predictors and response variables for the TMII model
Before model fitting, we divided our model development data for 2013-2015 from Sapelo Island pixel A into training (70% of data) and testing (30% of data) sets to create an internally consistent model that  minimized over-fitting. We reserved 2016 as temporally independent validation data that could estimate true out-of-sample error. We then began building our model using the training data. To build the model from training data, first we identified which bands should be predictors (step 2b). We identified bands related to marsh inundation by plotting the mean and standard deviation of MODIS spectra from summer (Jun-Aug) and winter (Dec-Feb) against PhenoCam percent flooding data. These plots revealed which bands created the greatest separation between flooded and dry observations and whether there were seasonal differences that we needed to consider.
Next, we created a binary variable for flood status (0: dry; 1: flooded) that could be the response in our model (step 2c). We used a binomial GLM to relate percent cover of water to MODIS band 5. We chose band 5 because the spectral reflectance plots from step 2b showed that reflectance within this band was greatly reduced during flooding (see Results). The relationship between flooding and band 5 reflectance disappeared if percent cover of water was < 7%, so we used 8% as our detection limit for flood status. Based on the above analyses, we used the following rules to create the flood status variable: a) if percent of water in the PhenoCam image was > 8% during the MODIS acquisition window, we classified the observation as flooded; b) if percent of water was < 1%, we classified it as dry; c) if MODIS acquisition time was unknown but the minimum water cover during the possible acquisition window (10:00 to 12:000 for MODIS Terra) was > 8%, we classified the observation as flooded; d) if the MODIS acquisition time was unknown but water cover was always < 1% during the possible acquisition window, we classified it as dry; c) all other observations had flood status set as missing data.

Creating and optimizing the Tidal Marsh Inundation Index (TMII) model
The next step ( Fig. 3 step 2d) was to create a binomial GLM model that could serve as the basis for TMII, with flood status (0, 1) as the response variable and MODIS band combinations as predictors. We evaluated as predictors band combinations suggested by the spectral plot from step 2b, with a particular emphasis on bands that were maximally different among flooded and dry observations (for a full list of band combinations tried, see Appendix A). In addition, we also evaluated existing NDWI-type indices from the literature. Existing indices included mNDWI (Xu, 2006), NDPI (Lacaux et al., 2007) and several variations of NDWI (McFeeters, 1996;Gao, 1996;Rogers and Kearney, 2004;Ji et al., 2009; Appendix A).
The spectral plot from step 2b suggested seasonal differences in the spectral reflectance of marsh flooding. Similarly, many of the NDWItype indices and other band combinations exhibited seasonal differences. Therefore, we also explored the addition of a phenology index that could seasonally-adjust the classification model. Phenology parameters explored included vegetation indices, any NDWI indices that mirrored phenological signals in our training pixel (had the same relative curvature and temporal alignment as the vegetation indices), as well as a cosine fit of the phenology index point cloud (see Appendix 1 for a complete list of phenology predictors).
To select the best model for TMII, we used binomial GLMs to relate flood status to each possible predictor and selected the model that maximized explained deviance. We also used a contingency table approach to assess model output (Bennett et al., 2013;Kuhn, 2008), including accuracy (correctly classified flood presence and absence / total number of samples * 100), as well as sensitivity (true positives: correctly classified flood presence / total observed flood presence * 100) and specificity (true negatives: correctly classified flood absence / total observed flood absence * 100) against training, testing, and validation data. As an additional check on model selection, we created Receiver Operating Characteristic (ROC) curves and selected the best model as that which maximized the Area Under the Curve (AUC) (from R package "pROC") (DeLong et al., 1988; Hanley and McNeil, 1982;Perkins and Schisterman, 2006). ROC curves compare classifiers by plotting model specificity against 1 -sensitivity (i.e. false positives) for all possible decision boundaries (Hanley and McNeil, 1982). ROC curves distinguish among poor classifiers whose decision boundaries are not optimized vs those that are poor because the predictors do not explain the data (Hanley and McNeil, 1982). This approach recognizes that classification sensitivity and specificity are trade-offs, where often you increase one at the expense of the other, while accuracy is maximized at some intermediate value of each (Hanley and McNeil, 1982). Models that maximize the AUC have the greatest potential overall classification accuracy. The greatest possible accuracy for each model will be obtained with the decision boundary that creates the specificity and sensitivity of the point in the upper most left corner of the ROC curve. In this case, explained deviance and AUC always selected the same model as best.
For the best fit model, we optimized the decision boundary for separating flooded and dry pixels (Fig. 3, step 2e). Our optimization criteria were not to maximize classification accuracy of both wet and dry observations, but rather to maximize sensitivity (ability to detect flooding) while retaining > 70% specificity (ability to correctly identify dry observations). These criteria are purposely conservative because daily data have enough annual observations to allow an index with low specificity (dry images incorrectly identified as flooded). These criteria identify more of the truly flooded observations than a model that focuses on overall accuracy, resulting in a tide-free dataset that can be trusted for downstream analysis.
Finally, we used the best fit binomial GLM equation to create a simplified TMII index that could be applied to any MODIS pixel (Fig. 3,  step 2f). For this, we rounded the model parameter estimates to the nearest integer, applied the logit transformation to convert the best fit GLM equation from linearized model space back to the units of the predictors (see Zuur et al., 2009), and then reduced the resulting equation as much as possible. We then optimized the decision boundary again and compared the results to the original best fit GLM equation, ensuring that the simplified equation behaved identically.
We did several types of comparisons to show that the optimized TMII was a step forward for remote sensing of inundation status in tidal marshes. This involved comparing final model performance in the model development pixel (Sapelo Island pixel A) to the best NDWI index from the literature in terms of explained deviance, as well as calculating the various accuracy metrics and AUC (step 2g). We also used the TMII to identify flooded and dry conditions for MODIS pixels B and C on Sapelo Island and MODIS pixel D on Plum Island, and calculated sensitivity, specificity and accuracy in comparison to groundtruth data (step 2h).

Demonstration of how new users can apply TMII to create TMIIfiltered composited vegetation time-series
The primary purpose of the TMII is to assist with creating filtered time-series that track vegetation dynamics in coastal marshes. We developed a workflow that follows established compositing procedure as much as possible (Huete et al., 1999), except that we incorporated the TMII as a tide-filtering step before compositing (steps 3a-f). Thus, from daily MODIS data, we subset to cloud-free data, data with good quality flags and view zenith angles < 50°(step 3a), calculated TMII for the remaining pixels (step 3b) and further removed pixels with TMII index values greater than our optimized decision boundary (step 3c). Then, following the general guidance outlined in Huete et al. (1999), we created a custom R function to generate an NDVI time-series composite (step 3d). The custom function composited NDVI with the following rules: for each 16-d window, order NDVI observations by view zenith angle. Up to five observations, in order of lowest view zenith, were used to create the average NDVI for the window, provided that these observations all had view zenith angles < 35°. If low view zenith observations were not available, the function selected the single NDVI value corresponding to the lowest view zenith angle < 50°. If no low view zenith values were available < 50°within the window, then the window NDVI value was set as missing data.
We demonstrated the transferability of TMII by applying our compositing procedure to generate tide-filtered 16-d NDVI time-series composites for 2013-2016 for Sapelo Island (pixels A-C), Plum Island (pixel D), Bayou Sauvage (pixel E) and Grand Bay (pixel F). In each case we compared the resulting vegetation time-series with those obtained by using MOD13Q1, a 16-d MODIS NDVI composite product.
Next, we provided a demonstration of how a new user can validate TMII-filtered vegetation composite time-series (step 3e). Validation steps may require creativity because accurate marsh inundation data often are rare (see Results and Discussion). An alternative approach is to validate the TMII-filtered vegetation-time-series, rather than validate TMII itself. This is effective because appropriate filtering steps result in accurate downstream datasets. Therefore, end-user validations options include 1) compare the TMII to accurate marsh inundation data, though data may be rare, 2) compare the resulting TMII-filtered vegetation time-series with ground-truth data such as in situ spectral reflectance or biophysical vegetation data. For this latter option, data availability varies but should be easy to acquire with planning. Lastly, 3) we can compare the time-series with known patterns in the biophysical data, such as might be derived from the literature. We expect this last option to be the most widely available, but care needs to be taken because previously published literature patterns may not match local conditions. However literature patterns can still provide a coarse litmus test to decide if the data processing steps give broadly reasonable results. For example, Spartina alterniflora typically has maximum standing biomass in mid to late summer (Giurgevich and Dunn, 1982). If our reflectance vegetation time-series has a peak outside this range without a corresponding environmental explanation, than we should be skeptical.
At new sites, one can use the decision boundary derived from this paper as a first cut. Following this, the TMII decision boundary can be adjusted until it 1) classifies flooded observations with the specificity and sensitivity that matches user criteria against ground-truth data or, 2) maximizes the correlation between known patterns in field data (step 3f) and repeat steps 3e-d as necessary.
To validate the time-series data for pixels other than pixel A, we used all of the above methods. First, we sought inundation data for the other MODIS pixels analyzed. The PhenoCam field of view spans the Sapelo Island MODIS pixel footprints we analyzed, with pixels A and B in the foreground and pixel C in the distance (Fig. 4). We therefore assumed that the flooding regime for each of these pixels was identical. For Plum Island, the PIE-LTER maintains a tide gauge near to study pixel D (Fig. 2). We assumed that the lag between the tide gauge and marsh flooding was 30 min based on observations by PIE-LTER researchers (I. Forbrich, pers. comm.). Therefore, we only retained observations corresponding to extremely high or low tides (i.e., when we could be fairly confident that the marsh was dry or flooded).
Independent marsh inundation data were not available for the Gulf Coast (pixels E-F) and our initial exploration of nearby creek gauges yielded improbable results (data not shown). We therefore could not directly assess TMII wet/dry classifications of these pixels. However, we had green biomass data that we compared to the TMII-filtered vegetation series. Green biomass typically has a moderate relationship with NDVI (Gamon et al., 1995;Ghosh et al., 2016;Mutanga and Skidmore, 2004). Aboveground green biomass data were available for both of the Gulf Coast sites every other month from January-October 2015, and were collected from eight 45.7 cm 2 area clip plots within MODIS pixel footprints E and F, oven dried at 40°C and scaled to a meter-squared basis. Biomass data from a given sample date were averaged to a single value within each pixel footprint. We then compared the correlation coefficient between biomass and each NDVI time-series. The best MODIS time-series should have a correlation coefficient with biomass that resembles the relationship among in situ NDVI and biomass.
In the Gulf Coast, we also had in situ spectral reflectance data for each biomass plot that we compared to the composited vegetation reflectance data. In situ reflectance was collected with a SVC HR-1024i spectroradiometer (Spectral Vista Corporation, Poughkeepsie NY). The spectroradiometer collected data across 350-2500 nm (1.5-3.8 nm band widths) as an average of 10 measurements. We used 4°field of view foreoptics at nadir 1 m above the vegetation canopy. We calculated canopy reflectance as the ratio of canopy radiance to radiance measured from a calibrated white reference (Spectralon® Labsphere, Inc.). White reference measurements were collected every 10 min. We used R to preprocess spectroradiometer reflectance by removing wavelengths with high noise from water absorption features and then resampled the data to average reflectance with the same band widths as MODIS (e.g. as in Table 1). Finally, as with the biomass data, reflectance samples from a given date that fell within the same MODIS footprint were averaged together (eight samples each per pixel per sample date). We then calculated NDVI using the same formula described in Appendix 1. To compare in situ NDVI to composite MODIS NDVI, we needed to estimate satellite-derived NDVI on the day the in situ NDVI was collected. Thus, we fit a smoothing spline through each composite MODIS MOD13 and TMII-filtered composited NDVI timeseries and predicted NDVI on the date of field sampling. As a measure of accuracy, we calculated root mean square error (RMSE) between field and composited MODIS NDVI as (Σ (in situ NDVI t − composited NDVI t ) 2 /N) −(1/2) for each time-series, where NDVI t is the NDVI value at time t and N is the sample size (Hyndman and Koehler, 2006). Predicted MODIS NDVI should be much lower than field-measured NDVI if dry and inundated observations were composited together.

Acquisition and pre-processing daily MODIS of imagery
We acquired imagery from 1 Jan 2013 through 31 December April 2016 for all six pixels included in this study (mean: 1445 observations per pixel). We subset these observations to those with high band quality, zero cloud cover, and view zenith < 50°, which removed 77% of observations, leaving an average of 337 observations per pixel across all years.

Accurate tidal marsh inundation ground-truth data
One finding of this investigation was that creek water height was not well-related to the PhenoCam-derived marsh percent flooding (explained deviance: 53%). At intermediate water heights, for example (0.9-1.0 m), the PhenoCam record suggested that percent marsh flooding ranged from 0-99%, with a mean and standard deviation of 46% ± 35% (Fig. 5). Adding information on wind speed, wind direction, month and year improved the relationship only slightly (explained deviance: 58%). We concluded that measurements of water height do not necessarily provide an accurate indication of marsh flooding, and ultimately used the PhenoCam record to develop and assess the TMII.

Identifying predictors and response variables for the TMII model
MODIS spectral plots showed differences in terms of both seasons and flood status (Fig. 6). Mean reflectance was greater during summer and lower during winter, especially in the visible through NIR spectra. Flooding depressed mean reflectance throughout the NIR and SWIR spectra (bands 2-859 nm, 5-1240 nm, 6-1640 nm, and 7-2130 nm), particularly when compared to the non-flooded pixels from the same season. Among these longer wavelength bands, band 7 generated the least separation among flooded and non-flooded observations. However, flooding hardly influenced visible wavelengths (Fig. 6). These results suggested broad groupings, such that bands 1, 3 and 4 (VIS) and bands 2, 5 and 6 (NIR and SWIR) each were similar and probably could be used somewhat interchangeably in a flood classification model (Fig. 6).
There also were seasonal differences. Flooded conditions during winter produced the lowest reflectance in NIR and SWIR bands. Winter dry conditions also overlapped with summer flooded conditions. Flooding reflectance may have changed seasonally because flooded marshes consisted of mixed reflectance from both vegetation and water. The overlap between winter vegetation and summer flooding suggested winter vegetation without pigment scattered light in a manner that mirrored the reduced spectral reflectance (increased absorption) of partially inundated summer vegetation. When vegetation was both flooded and lacked pigment, reflectance in NIR and SWIR was lowest. These observations indicated that flood status indices need to account for seasonal variation across some wavelengths. However the overlap in spectral reflectance, both seasonally and with inundation status suggests a perfect classifier might be difficult to create.
Binomial GLMs suggested percent cover of water had no relationship (P > 0.1) to MODIS bands when water was < 7% of the PhenoCam image. As described in the Materials and methods, we used this information to set detection limits for the binary flood status response variable in the TMII models, such that a flood status of 1 indicates > 8% cover of water within the pixel.

Creating and optimizing the Tidal Marsh Inundation Index (TMII) model
To predict flood status, we tried many band combinations suggested by Fig. 6 with varied success (see Appendix 1). NDWI-type indices also were not perfect predictors of marsh flood status, such that explained deviance ranged from 9-44% when applied to our testing data from Sapelo pixel A. In general, band combinations composed of both visible and non-visible spectra performed best because the difference of VIS wavelengths and NIR/SWIR wavelengths was greatest between flooded and dry plots. The best result (explained deviance 44%, Table 2, Appendix 1) was obtained through an NDWI built from band 6 (the SWIR band centered at 1640 nm) and band 4 (the green band centered at 555 nm) (band4 − band6) / (band4 + band6), similar to the index proposed by Xu (2006).
Combining NDWI 4,6 with a phenology variable improved explained deviance from 44% to 47% (Table 2, Fig. 7). The best phenology variable to include was the one least influenced by flooding. This allows the phenology parameter to track vegetation development and scale TMII seasonally. Unlike the flooding variable, the most appropriate  Table 2 Goodness of fit measures for the Tidal Marsh Inundation Index (TMII) and the best Normalized Difference Water Index from the literature (Xu, 2006), based on MODIS band 4 and band 6 (NDWI 4,6 ). ED is explained deviance and AUC is Area Under the Curve. Sensitivity and specificity are the accuracy of classifying flooded and dry observations, respectively; accuracy is overall classification accuracy; N is sample size.  phenology variable to incorporate into TMII was one composed of only longer wavelengths, because the relative difference between long wavelengths was similar in both flooded and non-flooded pixels (Fig. 6). Conversely, potential phenology variables composed of combinations of visible and longer wavelengths had higher noise and less distinctions among winter and summer values (Fig. 7c). The phenology variable that provided the best signal (highest peak in summer, lowest peak in winter, and low local variance) was NDWI 2,5 , composed of the NIR band at 859 nm and SWIR band at 1240 nm (band2 − band5)/ (band5 + band2). We further reduced local variance in this phenology variable by using its rolling mean from a time-series window 40 observations wide, centered on the day. The soil adjusted vegetation index (SAVI), an index which attempts to reduce the influence of soil moisture, was also a potentially valid phenology candidate, but had lower explained deviance (see Appendix 1).
This best model explained 47% of model deviance (Table 2). We transformed Eq. (1) The TMII had a higher explained deviance and AUC than the best NDWI index, and performed slightly better in almost all cases, particularly in terms of sensitivity to flooding and overall accuracy (Table 2; Figs. 7 and 8).
The optimized decision boundary for TMII in Eq.
(2) was 0.2 (index value > 0.2 = flooded, < 0.2 = dry). This decision boundary generated a classification accuracy, sensitivity and specificity equivalent to the underlying GLM used to estimate flood status (e.g., Eq. (1)). AUC for the simplified index equation was 0.91 (Eq. (2); Table 2; Fig. 8). Applying Eq. (2) and the optimized decision boundary to the training set yielded a sensitivity of 81%, specificity of 79% and overall accuracy of 79%. After all filtering steps, Sapelo Island pixel A from 2013-2016 was reduced from 1447 observations (full dataset), to 324 (cloud, view zenith and quality control filtering), to 223 (tide-filtered), representing an 85% data reduction. The fully-filtered dataset averaged 56 observations per year, and across all years, averaged 5 observations per month (SD = 2, range 1-9). These numbers are the numbers a typical end-user will have to work with, but differ from sample sizes in Table 1 because Table 1 represents data where we also have ground-truth observations from the PhenoCam.

Demonstration of how new users can apply TMII to create TMII-filtered composited vegetation time-series
Incorporating the optimized TMII into MODIS data processing steps allowed us to create composited vegetation time-series for 2013-2016 (Fig. 9). We validated these using available ground-truth data and by comparing the resulting phenologies to expected patterns. In Sapelo Island pixel A, the model development pixel, scenes classified as flooded based on the TMII generally corresponded to flooded scenes recorded by the PhenoCam, which was available starting in September 2013 (Fig. 9b). The optimal TMII-filtered composite had consistent growing season maxima in September (2 Sept 2013; 5 Sept 2014; 24 Sept 2015; 10 Sept 2016) and late winter/early spring minima (14 Apr 2014;16 Mar 2015, 19 Apr 2016) (Fig. 9c), similar to known phenology patterns. This contrasts with the existing MOD13Q1 time-series product, which for Sapelo pixel A had a weaker seasonal cycle and phenological timing that did not correspond with ground observations (i.e. an inaccurate phenology peak during November 2013). The TMII-filtered composite was also better-aligned with known patterns for Sapelo Island pixels B-C, where growing season maxima ranged from July through September and minima ranged from January through early March (Fig. 10). For these pixels, we could also validate TMII against flood data. Specificity and accuracy measures for the two additional pixels on Sapelo Island ranged from 73 to 89%, whereas sensitivity was lower (Table 3).
Applying the TMII to additional marsh sites also resulted in phenology peaks and minima more closely aligned with expected time periods than MOD13. At the S. alterniflora salt marsh on Plum Island (pixel D), the MOD13 composite and the TMII composite were superficially similar, with both series showing phenology peaks in July and  August (although the peaks occurred slightly later in the TMII-filtered data) (Fig. 11). However, the MOD13 time-series had minima in March-April, whereas in the TMII time-series, the minima occurred in January through February (Fig. 11a-c). The Plum Island data set only included very high and low tides (see Materials and methods), and under these conditions the TMII itself was approximately 80% accurate with a sensitivity of 60% (Table 3). Adjusting the TMII decision boundary did not improve accuracy at this site, so we used the optimized cut-off of 0.2 identified in the previous section.
In the S. patens marsh at Bayou Sauvage, where marsh flood groundtruth data were not available, the TMII composite NDVI had a lower RMSE when compared with in situ NDVI than did the MOD13 composite when using the decision boundary from Sapelo marshes (Fig. 12a) (RMSE TMII : 0.094; RMSE MOD13 : 0.183). At this site, as a demonstration of site specific optimization, we improved RMSE somewhat by adjusting the TMII cut-off to 0.01 (after adjustment, RMSE TMII : 0.084). The resulting rebuilt TMII composite was the only one where the MOD13 and TMII composites were nearly identical, with maxima in mid-August through early September in all cases and minima in January and February (Fig. 11d-f). However, the MOD13 did not have as high a peak in summer as was observed in both the TMII composite and the in situ NDVI (Fig. 12a, b). The biomass data collected in this pixel did not perfectly resemble any NDVI time-series (Fig. 12b). However the correlation coefficient between S. patens biomass and NDVI was similar for in situ NDVI and the TMII composite (r: 0.14; and 0.21, respectively), while the MOD13 NDVI composite had a negative correlation with biomass (r: −0.15).
In the J. roemerianus marsh pixel at Grand Bay, marsh flood ground-truth data again were not available. As observed at Bayou Sauvage, the in situ to TMII composite NDVI had a lower RMSE than did the in situ to MOD13 NDVI composite, when using the decision boundary from Sapelo marshes (Fig. 12a) (RMSE TMII : 0.046; RMSE MOD13 : 0.052). Again, as a demonstration of site optimization, we improved RMSE by adjusting the TMII cut-off to 0.01 (after adjustment, RMSE TMII : 0.033), and then rebuilt the composite. The resulting vegetation time-series had a seasonal cycle with a lower amplitude than observed in the other pixels evaluated (Fig. 11g-i). The TMII composite also had a more reasonable phenology than MOD13. For example, the MOD13 composite indicated a phenology peak during September in 2013, both an April and November peak in 2014, and a February and December peak in 2015 (Fig. 11g). Minima in the MOD13 composite were similarly scattered throughout the calendar year. The phenology peaks in the TMII composite were more consistent, with September peaks in both 2013 and 2014, a December peak in 2015, and a July peak in 2016 and minima in February through April (Fig. 11i). The biomass data collected in this pixel also mirrors both in situ and TMII composite NDVI more closely than the MOD13 composite (Fig. 12c). The correlation coefficient between J. roemerianus biomass and NDVI was similar for in situ NDVI and the TMII composite (r: 0.77 and 0.74, respectively), and was far less for the MOD13 NDVI composite (r: 0.21).

Discussion
This paper describes the development of a new Tidal Marsh Inundation Index (TMII) for pre-processing MODIS imagery for vegetation analysis in coastal wetlands. Once developed, TMII identified tidal flooding in S. alterniflora marshes with 72-81% accuracy. The TMII also generally aligned NDVI marsh vegetation time-series with expected patterns, even where inundation ground-truth data were not available. Identifying and filtering flooded scenes was important because flooding dampens spectral reflectance (Adam et al., 2010;Cho et al., 2008;Kearney et al., 2009;, usually resulting in lower NDVI values in flooded scenes and potentially misrepresenting biomass (Kearney et al., 2009). Additionally, the inclusion of flooded scenes adds scatter to the time-series, making it difficult to discern seasonal cycles and other patterns in the underlying vegetation data (Kearney et al., 2009;.

Building, optimizing and accessing accuracy of TMII
The flood status indicator used to develop the TMII was calibrated with digital imagery from the "GCESapelo" PhenoCam, mounted on the GCE-LTER flux tower within the model development pixel. This was fortuitous, as our original intent was to predict marsh inundation based on water level measured in an adjacent creek. However, we observed large variation in the extent of marsh flooding at intermediate creek water heights (Fig. 5). Marsh flooding was likely influenced by a combination of factors, including wind speed and direction. Changes in local topography and drainage patterns also are common on the marsh surface (Hughes et al., 2009;Pennings and Richards, 1998), and can influence flooding. Vegetation development also may have been important, because the amount of water required to cover vegetation varied with plant height. We found creek water height combined with climate variables and temporal predictors estimated marsh flooding better than water height alone. However, creek water height still failed to explain much of the observed variation. External datasets of water height, such as creek water gauges or NOAA tide height data, may not be ideal for estimating marsh flooding. To measure surface inundation, we suggest digital imagery and/or water level gauges placed directly on the marsh surface would be more accurate than gauges placed in adjacent water bodies. An important feature of the TMII, however, is that once calibrated it can classify MODIS pixels without additional external water level information. As demonstrated here, our calibration from Sapelo Island performed fairly well in a wide range of situations.  In our analysis, the best existing NDWI-type index used the normalized difference of bands 4 (green) and 6 (SWIR), which is similar to the water mapping index proposed by Xu (2006). Ji et al. (2009), who reviewed the many NDWI indices in the literature, also found that this index performed well. The TMII builds on this by adding a vegetation phenology variable. The addition of phenology was necessary because the spectral reflectance of marsh vegetation changed seasonally such that winter vegetation overlapped with flooded observations from summer. Therefore, decision boundaries could not be optimized across winter and summer conditions. The incorporation of a phenology variable helps to scale the TMII relative to vegetation growth. The best phenology parameter was the localized mean of the normalized difference of MODIS band 2 (NIR) and band 5 (SWIR), an index similar to that proposed by Gao (1996). Though called NDWI, Gao's NDWI index was created as a vegetation index for tracking water content of vegetation canopies. Similar to Gao, we found that the band combination in this index was less noisy for tracking wetland vegetation than more traditional vegetation indices, such as NDVI or EVI, perhaps because it did not combine visible and longer wavelengths. Cho et al. (2008) also observed that NIR reflectance of submerged aquatic vegetation was reduced more rapidly than VIS as water depths incrementally increased. Thus indices composed of visible and longer wavelengths may be more influenced by variation in inundation and soil moisture, pervasive features of salt marshes that can mask vegetation dynamics. NDWI 2,5 was a good pretide filter phenology indicator because it followed the same shape as the other potential phenology parameters (e.g. NDVI), but was the least influenced by inundation.
Overall accuracy for the TMII ranged from 72-81% across training, testing and validation data. We could not create a perfect classifier because spectral reflectance for flooded and dry scenes overlapped. We reduced sources of error for the index by including quality control steps as part of the MODIS data workflow, effectively removing cloudy pixels, observations resulting from poor detector/sensor readings, and observations of different radiometric quality such as caused by variation in view zenith (Klemas, 2011). However, some error remained, resulting from the remaining variation in view zenith, imperfect atmospheric correction (Huete et al., 1999), and environmental variation such as water, chlorophyll and suspended sediment, all of which alter spectral reflectance (Han et al., 1994;Mittenzwey et al., 1992). Additionally, marsh surface flooding was uncommon and the magnitude of flooding varied, resulting in an unequal variance structure among flooded and dry conditions, which made it difficult to create a highly flood sensitive classifier. When we tuned the index to increase flood sensitivity (> 90%) in one data set, it was less generalizable to new data. TMII represents a compromise in that it identifies the majority of observations that create noise in wetland vegetation time-series, does not misclassify a majority of dry observations as flooded, and can be applied to new data with minimized accuracy losses.
The final TMII was fairly successful at identifying marsh flooding in the model development pixel, and performed better than the NDWI index alone (Table 2). It also predicted flooding at the other Sapelo Island pixels (B and C) (73-77% accuracy) as well as on Plum Island (D) (80% accuracy), where we had independent ground-truth information (Table 3). Although applying TMII removed nearly 25% of observations from the pre-processed pixel A data, MODIS is an ideal sensor for implementing tide-filtering because of its high temporal resolution. Even with a significant data loss during all pre-processing and filtering steps, data density was still high in pixel A (N = 223 in total, with an average N of 5 ± 2 per month). For comparison, the maximum number of observations from another common multi-spectral satellite, Landsat, is 2 per month (landsat.usgs.gov).

Demonstration of how new users can apply TMII to create TMII-filtered composited vegetation time-series
When TMII was included in the MODIS workflow, vegetation composites developed for S. alterniflora pixels were generally consistent with expected patterns and better fit field observations than existing MODIS products (e.g. MOD13). The difference between TMII-filtered and MOD13 composites was greatest for Sapelo Island. This was likely because Sapelo pixels were almost entirely comprised of short and medium form S. alterniflora (Hladik et al., 2013). Inundation should be detected from nadir views more often when low-statured, sparse vegetation is combined with high tidal amplitudes, increasing the need for tidal filtering. Although the vegetation and tidal range were similar on Plum Island (Redfield, 1972), the salt marsh extent there was small, had greater edge habitat and thus probably included more tall form S. alterniflora, and it was difficult to select pixels without some portion of upland cover. Mixed pixel cover types may have increased NDVI and reduced the influence of tidal flooding at this site. However, TMII filtering still improved the phenology signal at Plum Island by suggesting earlier and more appropriate phenology minima than MOD13 (e.g., January and February rather than March-April).
The improvement over MOD13 was less dramatic when TMII was applied to the two Gulf Coast sites, which may partly be due to the lower tidal range (Orson et al., 1985) and higher stem densities of S. patens and J. roemerianus (unpublished data). Tides likely influence Gulf Coast marsh pixels less than Atlantic Coast pixels. However, the RMSE between satellite and in situ NDVI was still lower after TMII filtering than that for MOD13. Although the optimized cut-off from Sapelo Island improved the Gulf Coast vegetation time-series, greater improvements were achieved by adjusting the decision boundary until RMSE between in situ and satellite NDVI was minimized. At both Gulf Coast sites, we ultimately used a decision boundary of 0.012. This adjustment probably improved the outcome because of the higher stem densities and lower tidal amplitude. The resulting TMII-filtered vegetation timeseries had more peaks and minima during expected times than those observed in the MOD13 composite and the TMII-composites also more closely followed the patterns in the green biomass data and in situ NDVI. Within Bayou Sauvage, in situ NDVI was not perfectly matched by either satellite composite, though the TMII composite was a closer match than the MOD13. We expect the cover types to differ somewhat between the satellite scale and the plot scale, particularly in terms of increased bare soil and dead vegetation at the larger scale, likely causing this discrepancy.
In the J. roemerianus marsh at Grand Bay, neither the vegetation time-series, the in situ NDVI, nor the biomass data reflected the sinusoidal vegetation pattern we observed elsewhere. The different phenological pattern of Juncus may reflect its physiology: Juncus is a C3 species while Spartina exhibits the C4 photosynthetic/production strategy (Giurgevich and Dunn, 1982). Although species variation is high, C3 plants are often associated with cool season growth whereas C4 plants generally have more efficient photosynthesis, greater tolerance for high temperature, and warm season growth (von Fischer et al., 2008). Thus, J. roemerianus typically has a fairly constant standing crop, resulting from even daily production with less seasonal variation than S. alterniflora (Giurgevich and Dunn, 1982).
An advantage of TMII is that it requires little additional data to apply and is relatively easy to incorporate into the MODIS workflow. Previous methods for filtering inundated pixels are more complex and may be less accurate. Such methods include modeling basin topography to create a submersion time index to weight vegetation observations (Wang et al., 2012), using complex algorithms and band rules to identify seasonal inundation (Knight et al., 2006), or modeling the movement of the tidal stage across the landscape from tide gauges and marsh elevation estimates (O'Donnell and Schalles, 2016). In contrast, Fig. 12. A) Comparison of Gulf Coast in situ NDVI collected with a field spectroradiometer vs MODIS-500 m scale MOD13 and TMII-filtered 16-d composites. In this panel, the 1:1 line is drawn for reference. Points close to the 1:1 line are more desirable. B) Bayou Sauvage S. patens in situ NDVI collected with a field spectroradiometer, MODIS-500 m scale MOD13 and TMII-filtered 16-d composites and aboveground biomass plotted over time. C) Grand Bay J. roemerianus in situ NDVI collected with a field spectroradiometer, MODIS-500 m scale MOD13 and TMII-filtered 16-d composites and aboveground biomass plotted over time.
TMII relies on a single index value calculation (Eq. (2)), from which observations are classified as flooded or dry based on an agreed decision boundary. In this paper we applied a decision boundary of 0.2, derived from a Sapelo Island S. alterniflora marsh, to all S. alterniflora pixels. The 0.2 cut-off also worked for other vegetation pixels, but we could improve the result in the Gulf Coast by validating against field data and optimizing the decision boundary until the correspondence between field and satellite data was maximized. In all cases, the TMIIfiltered NDVI composites compared favorably to MOD13. Additionally, we likewise expect improvements if we use TMII to create other composites similar to MODIS biophysical products such as MOD15 (LAI and fPAR product) and MOD17 (GPP/NPP product). The TMII also may be useful in coastal habitat mapping studies although for these studies, the best solution might be to use imagery from the lowest tides. The next step is to scale-up from the pixel scale to broad-scale landscape analyses, where TMII-filtering during preprocessing should make broadscale trends more apparent.

Conclusions
TMII filtering can help improve MODIS models of wetland vegetation, whether used alone or as part of a time-series compositing procedure as we demonstrated here. This can allow a more accurate view of vegetation dynamics in coastal ecosystems, and may also provide the ability to detect atypical patterns. All of this is highly relevant for studies of phenology, primary production and carbon storage in coastal wetlands. Incorporating TMII as a pre-processing step before constructing long-term time-series wetland products will not only yield more accurate phenological trends but can also be used in GPP prediction models. We recommend that TMII be routinely applied to preprocess MODIS time-series in coastal tidal marshes. Ultimately, TMII filtering in combination with vegetation reflectance analysis can assist with upscaling to understand broad-scale patterns.