Influence of Pixel Quality, Land Cover, and Hydroclimatic Cycle on Moderate Resolution Imaging Spectroradiometer Inundation Monitoring Performance in the Pantanal, Brazil ()
1. Introduction
Large wetlands provide superregional ecological services and are highly vulnerable (Tockner & Stanford, 2002; Mitsch & Gosselink, 2000) . Flooding dynamics are the key driver of their biodiversity, ecological processes, and matter and energy fluxes (Junk et al., 2013) . The spatial and temporal representative monitoring of inundation by field methods, however, is mostly not viable because of their extension and logistic reasons, so remote-sensing methods have been widely applied (Bergé-Nguyen & Crétaux, 2015) .
Passive and active microwaves (Ward et al., 2014; Schroeder et al., 2015) , as well as optical sensors (Sakamoto et al., 2007; Huang et al., 2014) , have been found suitable for inundation detection in floodplains.
However, each of these approaches has shortcomings. Passive microwave remote sensing (RS) imagery has spatial resolution lower than 10 km, whereas available space-borne, active synthetic aperture radar (SAR) systems provide only short time series (TS), have temporal resolutions lower than 20 days and shortcomings in the distinction between flooded and non-flooded grasslands (Di Baldassarre et al., 2011) .
Revisit frequency of spatial mid/high-resolution optical RS systems with long-term coverage such as Landsat do not provide cloud-free datasets to capture sub-seasonal dynamics (Aires et al., 2018) . More recent systems, such as the Sentinel-2 constellation can now provide imagery with a five-day repeat cycle (DeVries et al., 2020) . Nevertheless, such temporal resolution continues low to represent intra-annual flooding regimes over large tropical wetlands, which are heavily affected by cloud cover. Furthermore, available TS are short yet to adequately assess the interannual range in flooding extents and their trends.
With this background, middle- and low-spatial-resolution optical remote-sensing imagery, such as the moderate resolution imaging spectroradiometer (MODIS), have been found to provide a compromise between imaging frequency and spatial resolution to assess inundation and recession patterns in large floodplains (Sakamoto et al., 2007; Ordoyne & Friedl, 2008; Almeida et al., 2015) . Importantly, MODIS provides TS long enough to assess long-term variations of flooding dynamics, with the coverage to operationally monitor on a global scale and continues to have a unique value for studying short and long-term inundation changes (Alonso et al., 2020) .
Shortcomings in MODIS flood-monitoring performance has been attributed to interference from pixel quality as a function of varying atmospheric conditions, sun-sensor-surface viewing geometries (Carreiras et al., 2003; Huete et al., 2002; Kobayashi & Dye, 2005) , flooding stage (Sakamoto et al., 2007; Islam et al., 2010; Almeida et al., 2015) and vegetation cover (Ordoyne & Friedl, 2008) .
Large tropical wetlands are located in humid (Amazon and Ganges Delta) or seasonal (Llanos del Orinoco, Okovango, Mekong Delta, Everglades, and Pantanal) climates, where, during at least parts of the year, passive remote-sensing methods are affected by cloud cover. Therefore, it is supposed that classification performance is subject to climatic-driven seasonality, implying cyclical variations of mapping accuracies during different flooding and retraction phases. Because of their large extent and low slopes, the timing of flood pulses varies inside the floodplains. In the case of the Pantanal, in the northern and western parts and near the non-flooded surrounding contribution areas, the high water stands coincide with the rainy season and dense cloud cover, whereas flooding in the central and southern parts is delayed by as much as five months. In these parts, high water stands occur in periods with minimal cloud cover. In floodplain areas, where precipitation and flooding cycles coincide, the mapping of maximum flooding extents is supposed to be more strongly biased than that in the terrestrial phases. The detection of other important phases in floodplain ecology—for example, when important shifts in river-floodplain connectivity occur at the beginning and end of the rainy season—may demonstrate poorer-than-average classification accuracies because of atmospheric conditions.
The fact that the vegetation overstore obscuring standing water affects the flood-mapping performance with passive remote-sensing methods is well known (Gibbs et al., 2016) ; however, it was not quantified for MOD13Q1 TS. To minimize such effects, classification approaches apply spectral water indices, such as the land surface water index (LSWI) (Gao, 1996; Chandrasekar et al., 2010) in conjunction with vegetation indices (VIs) as covariates (Huang et al., 2014) . Independently, the flooding signal has been found to be increasingly mimicked with densifying vegetation cover (Ordoyne & Friedl, 2008; Gibbs et al., 2016; Garcia et al., 2015) . The influence of canopy structure on flood mapping has not been systematically studied. It has been shown, however, that vegetation formations with very different canopy structures can have similar VI values and that empirical models to estimate biochemical and biophysical vegetation properties from VIs, for example, may strongly vary with canopy structure (Croft et al., 2014) .
Different approaches have been proposed for flood monitoring with RS TS, including spectral unmixing (Padovani, 2010) , dynamic thresholding (Ji et al., 2009) , sensor fusion (DeVries et al., 2020) , and the use of ancillary datasets, such as surface elevation (Lant, 2013) . Inundation and flood monitoring is a dichotomous classification problem. Thus, frequently, logistic regression (LR) models have been applied to estimate the probability of a candidate pixel being flooded or not being flooded (Huang et al., 2014; Ordoyne & Friedl, 2008) . This approach was followed in this study, principally because the influence of data quality and geoenvironmental conditions on flood-mapping performance was principally examined. Under such a framework, it is advantageous for LR validation procedures to be straightforward and for model outcomes to be easier to reproduce and influenced less by model parameters than in more-complex data-driven classification approaches (e.g., neural networks) or nonlinear approaches, such as neural networks, vector machines, or techniques that involve boosting (boosted regression trees, etc.).
To improve information extracting from multitemporal imagery, TS analysis (TSA) techniques in the time and spectral domains have been applied for smoothing, filtering, and gap filling of degraded pixels (e.g., Beck et al., 2007; Bradley et al., 2007; Jönsson & Eklundh, 2002, 2004 ). Because of the diversity of tested techniques, slight performance differences in many case studies, and differences in the environmental conditions (type of precipitation regime, cloud cover, land use and cover, duration and depth of flooding, etc.), which interfere with the performance of TSA, this review could not identify a possible best TSA method for flood monitoring. Because comparisons of candidate techniques are not proposed, for gap filling and smoothing the Savitzky-Golay (SG) filter, a standard time domain TSA approach was used. SG results in local smoothing of a TS by an iterative weighted moving average filter, with weighting given as a polynomial of a particular degree (Chen et al., 2004) . SG has been applied successfully for preprocessing of water and VIs in heterogeneous landscapes. SG was found to perform well in overall noise reduction, under diverse characteristics of TS degradation. Furthermore, SG can adequately represent asymmetric TS profiles (Hird & McDermind, 2009; Atkinson et al., 2012) , as found in floodplains, which frequently have a steep increase of flooded areas after embankment overflow and a slow reduction of flooding extent during drying process.
Because of the ongoing importance of mid-resolution optical RS TS for inundation monitoring in large tropical floodplains and the performance shortcomings in approaches targeting the operational generation of worldwide flood products, in the present study, the following ongoing research issues are addressed. 1) What is the influence of different quality assurance (QA) levels on flood-mapping performance in quantitative terms? 2) What pixel filter strategies based on available QA values should be used to optimize TS classification performance? 3) What are the influences of vegetation cover and period of flooding cycle on the classification performance? Furthermore, the reliability of ground truth influences the construction of empirical classification models (Raudys & Jain, 1991) and determination of classification accuracy metrics (Foody, 2010; Olofsson et al., 2014) . Previous studies on land cover estimation showed, for example, that internal validation outcomes may overestimate in more than 30% of accuracy metrics (Oldeland et al., 2010) . Many reviewed studies on MODIS-based flood mapping in large wetlands did not count with sufficient field datasets to sustain external validation. Therefore, validation has been conducted using remote-sensing imagery with higher spatial resolution, such as Landsat TM (Klein et al., 2017; Dao & Liou, 2015; Klein et al., 2014; Sun et al., 2012; Chen et al., 2013) or SAR (Pope et al., 1994; Smith, 1997; Hess et al., 1995) . Because of the known and significant shortcomings of these systems used for ground truth in flood mapping (Gibbs et al., 2016; Frazier & Page, 2000) , it must be assumed that such validation efforts are biased because of ground reference data errors (Foody, 2010) . Therefore, particular care was taken in the elaboration of multiyear training and validation ground truth sets to enable the determination of true error rates and compare them with accuracies obtained from a k-fold cross validation and internal validation.
2. Materials and Methods
2.1. Study Area and Field Data Acquisition
With an area of approximately 140,000 km2 only in Brazil, the Pantanal is the world’s largest continuous wetland (Nunes da Cunha & Junk, 2010) . Its unimodal flooding cycle is driven by the semihumid tropical climate in the upper Paraguay river basin with annual precipitations varying between approximately 950 mm in the central western part of the floodplain and 2000 mm in the northeastern uplands. In the northern part of the Pantanal studied here and in regions near the surrounding highlands, flooding peaks coincide with the highest rainfalls between January and March or are delayed by approximately two months at maximum in the extreme southwest of the study area. The Pantanal is covered by a heterogeneous vegetation mosaic. Semideciduous and evergreen broadleaf forests dominate on holozoic deposits, principally along permanent watercourses. Cerrado savannah formations with varying densities of tree layers are common on higher and older alluvial fans. Grasslands predominate over lower parts of the lobes, depression zones between lobes, and stands with poorly drained soils (Zeilhofer & Schessl, 2000) .
In situ flood measurements from between 2005 and 2009 were available for two test sites, which well characterize the hydro-ecological amplitude of the Northern Pantanal. The northern test site of the Pirizal long-term sampling site (LTSS) (Fantin-Cruz et al., 2010; Signor et al., 2010) is built up by a 25-km2 grid of 40 limnimetric stages and covers an alluvial lobe with a longer-flooded depression zone in its southeastern part (Figure 1).
Vegetation formations include broadleaf evergreen and semideciduous forests, grasslands, and pastures with or without tree or scrub layers of different densities (Cerrado savannah), both intercalated by non-flooded Cerrado woodland on ancient levees or relict dunes. The second site is a 12-km transect inside the Private Reserve of Natural Patrimony (RPPN) of SESC Pantanal and is composed of
![]()
Figure 1. Study area in a SPOT 4/3/2 composite from the dry season of 2007, field monitoring sites at the Pirizal LTSS and the SESC Pantanal RPPN transect, monthly precipitation averaged for the study area from the Tropical Rainfall Measuring Mission (TRMM), and 3B43 product and daily Cuiaba river discharge at the Porto Cercado gauge (https://www.hidroweb.ana.gov).
11 limnimetric stages. The northern parts of the Pirizal site are flooded at the beginning of the wet season by stagnant rainwater, whereas the SESC transect is exclusively flooded by dam overflow of the Cuiabá river and by the more sediment rich São Lourenco river in the extreme south (Fantin-Cruz et al., 2011) . Here, grassland and scrublands intercalate with periodically flooded, monodominant Vochysia divergens (“Cambarazal”) forests and strips of riverine forests along the Cuiabá river and its secondary channels.
2.2. MODIS Products
For most analyses, a total of 299 images of the 250 m resolution, 16-day composite MOD13Q1 VI product from 2000 to 2012 were acquired and reduced to 100 images available for the period between January 2005 and August 2009 (Figure 1), for which in situ water-level measurements were available.
The MOD13Q1 product was built from two spectral VIs—normalized difference VI (NDVI) and enhanced VI (EVI)—four reflectance bands of the visible (blue, red) and infrared (near infrared, midinfrared), the composite day of the year band, which provides each pixel acquisition date, and comprehensive, two-level pixel-based QA information. Therefore, the pixel reliability (PR) summary is simplified, five-level information on pixel quality (Table 1).
In the dataset, only pixels labelled as good (0), marginal (1), or cloudy (3) occurred. The 250 m 16-day VI QA is a detailed QA system represented by a 16-bit string, with nine flags (Didan et al., 2015) . To evaluate the influence of pixel quality on classification performance, three flags were used for filtering the TS in subsets: 1) VI quality (bits 0-1), 2) VI usefulness (bits 2-5), and 3) land-water mask (bits 11-13) (Table 2).
Several spectral indices have been proposed for flood mapping, including such VIs as the NDVI and EVI (Peng et al., 2005; Kwak et al., 2015) and indices specifically developed for water surface mapping, such as the normalized difference water index (McFeeters, 1996; Gao, 1996) , the LSWI (Chandrasekar et al., 2010) , the modified normalized difference water index (Xu, 2006; Ordoyne & Friedl, 2008) , and the open water likelihood (OWL) algorithm (Guerschman et al., 2011) .
Despite its proven performance (Guerschman et al., 2011) , it was decided that the OWL algorithm would not be considered, which makes use of third-party data products (elevation model). The use of external covariates would complicate the interpretation in comparisons between classification performances resulting from pixel quality, vegetation cover, and hydrological period. Because of the large number of environmental, climatic, and hydrological conditions, data
![]()
Table 2. VI quality QA of the MOD13Q1 product, applied for data subsetting in the assessment of influence of pixel quality on flood-mapping accuracy (adapted from Solano et al., 2010 ).
products, and the classification and validation approaches used in previous studies, it was not possible to identify a possible most effective spectral water index. Instead, the LSWI was applied, which is calculated from the original near-infrared and midinfrared bands and has been widely recommended for MODIS inundation mapping (Gao, 1996; Chandrasekar et al., 2010) —see (1).
(1)
Several flood-mapping studies have recommended the conjunctive use of at least one water index and one VI (Ordoyne & Friedl, 2008; Guerschman et al., 2011) . Therefore, the predictive power of the LSWI was tested with the MOD13Q1 NDVI and EVI as covariates (Figure 2).
Land use and cover types have been identified to modulate spectral indices, in part independently from biophysical parameters captured by VIs (Liu & Kafatos, 2005) . Therefore, the suitability of the MCD12Q1 product to improve
![]()
Figure 2. Data-processing pipeline for predictor selection and the assessment of the influences of validation technique, pixel quality, vegetation cover and hydrologic period on the flood-mapping performance.
flood mapping was evaluated. The annual 500 m resolution Land Cover Type 2 (UMD) bands were resampled to the MOD13Q1 250 m resolution. The MOD13Q1 and MCD12Q1 products of the Northern Pantanal were from the h12v10 tile and were obtained as HDF format with the sinusoidal projection from the NASA Land Processes Distributed Active Archive Center (Didan et al., 2015) . They were then re-projected to the UTM (Datum WGS1984, zone 21S) using the MODIS Reprojection Tool (MRT) (NASA/USGS, 2011) and subset to the 8300 km2 study area.
2.2.1. Classification Algorithm, Data Processing, and Validation
Flood stages during the inundation period were visited with a weekly frequency on average (Fantin-Cruz et al., 2010, 2011; Girard et al., 2010) . Thus, the dates of gauge level measurements do not exactly coincide with the pixel acquisition dates of the MOD13Q1 temporal composites, which also vary inside a scene, principally during the wet season. To associate the pixel values in time with the inundation measurements at the corresponding gauges, a two-fold preprocessing was applied. Water levels between two dates with observed inundation were linearly interpolated. To estimate the date of ending inundation, the last nonzero measurement was extrapolated, using the maximum rate of flood reduction by day TL observed at each stage, as in (2).
(2)
Here, TL is the variation of inundation during a time interval tn+1 − tn, and Ln+1 − Ln is the difference of inundation levels during this interval. The function min selects the lowest value for each stage and each hydrological year.
After the dry falling of all gauges in both study sites, field measurements of the respective hydrological years were interrupted. Because the Pantanal flooding cycle is largely monomodal, the assumption that all sites remain dry after the first dry falling can be considered valid (Fantin-Cruz et al., 2011) .
To minimize the LR model bias caused by unbalanced sampling and avoid biasing of accuracy metrics by using an excessive number of nonflooded pixels (Chen et al., 2013) , samples were maintained from this dry period until an equilibrium of flooded and nonflooded observations was reached. According to this precondition, MOD13Q1 multiband pixelwise values and associated PR and VI QA levels and MOD12 values were extracted and linked, according to the acquisition date (DOY), with the respective observed or estimated in situ flood estimates (Figure 2). The entire tabular dataset with binarized flood-monitoring results and multiband and multiproduct pixel values consisted of 1444 ground truth measurements.
To identify the most-suitable predictors for flood mapping and evaluate the influence of validation technique on accuracy metrics, the dataset was filtered for PR values = 0. Then a generalized linear model with a logit link function was adjusted to the pixels of the limnimetric stages of the entire TS to assess the predictive power and significance of a complete model with the four tested independent variables (LSWI, NDVI, EVI, land cover). The Wald statistic and its significance levels were used to assess the predictors’ contributions to the model. Diagnostics of multicollinearity of independent variables were conducted using the variance inflation factor (VIF). Variable selection of the final multivariate models used in the further analyses was then conducted by the AIKE information criteria.
For further analysis, the complete dataset was then split into subsets. The training subset was extracted from 66 MODIS scenes—35 of them for the SESC transect site from two hydrological years (2005/2006 and 2006/2007) and 31 for the Pirizal LTSS site for the hydrological years 2004/2005 and 2005/2006. The subset for external validation was composed of 16 scenes for the SESC transect site of the hydrological year 2006/2007 and from 18 scenes for the Pirizal LTSS site of 2008/2009. Because many previous studies on flood monitoring used no or insufficient field data for external validation (Klein et al., 2014, 2017; Huang et al., 2014; Sun et al., 2012; Chen et al., 2013) , true error rates (external accuracy) were also compared with the training accuracy (internal accuracy) and a k-fold (k = 10) cross validation (Royston & Altman, 2010) , both based on the training subset.
2.2.2. Influence of Pixel Quality on Classification Performance
To evaluate the influence of pixel quality, including TS gap filling, on classification performance, 14 data subsets were generated by applying filters according to the PR summary and the detailed VI quality QA (Table 1). VI quality values of 00 and 01 were filtered according to the 15-level VI usefulness (VI_U). Inland waters identified in the land-water mask were classified as flooded in the output maps. Because none of the in situ flood stages fall in this category, they did not influence the validation outcomes. Furthermore, the entire TS has been visually inspected for known consistency issues between the VI quality and usefulness flags (https://lpdaac.usgs.gov/about/news_archive/).
To guarantee the major possible sampling amount, this assessment was used for the entire dataset with the field observations available (n = 1444). Two sets based on the PR band considered only pixels with values of 0 and 1, respectively. Twelve further sets included sets of decreasing VI_U. A gapless 15th set was built using the SG filter as implemented in the TIMESAT software (Jönsson & Eklundh, 2002, 2004) , reconstructing the spectral indices of the TS. After extensive empirical tests using VI_U and PR levels as TS gap filling and smoothing parameters, a dataset was selected for further analysis, reconstructed by the following model weights: PR_0 = 1.0, PR_1 = 0.8, and PR_3 = 0. The model parameters were: window size = 4, number of envelope iterations = 3, and adaptation strength = 2.
Pixel sample sizes at the field-monitoring stations of unique QA levels were, in some cases, insufficient for statistical inference. Therefore, VI_U strata had to be partially aggregated. According to predictor selection from the overall flood-mapping performance assessment, a reduced LR model was than developed for each aggregated subset. The subsets were validated separately and cumulatively.
2.2.3. Influence of Vegetation Cover and Hydrological Period on Classification Performance
For evaluation of the influence of vegetation cover and land use, the comparisons of mapping accuracy at the two field sites—two pan-sharpened multispectral 2.5 m SPOT 5 imageries from the 2007 dry season—were visually interpreted, defining three land cover classes: forest, Cerrado savanna, and grasslands. To minimize interference from pixel quality, without excessively reducing the sampling amount, the entire PR_0 dataset was used for this assessment without splitting. Because both gauging sites in the floodplain are under conservation, there was no abrupt land cover change (e.g., deforestation) during the period of the evaluated MODIS TS.
Furthermore, a three-level grouping of hydrological phases was conducted separately for each hydrological year of the TS, differentiating the periods based on when all limnimetric stages were flooded (flood phase), when some of the stages were flooded (transition phase: rising and receding floodwaters), and when none of the stages was flooded (dry phase). OA, Cohens K, and commission—false-positive rate (FPR)—and omission—false-negative rate (FNR)—errors were calculated for vegetation classes and the three hydrological periods.
Model building and validation were conducted with the R software (R Core Team, 2014) , and extensive scripts were coded for automatic data extraction and aggregation.
3. Results
3.1. Variable Selection and Model Performance for the Unfiltered Dataset
For the selection procedure of independent variables, the spectral indices LSWI, EVI, and NDVI of the entire PR_0 dataset (n = 1249) were tested for their predictive power in flood mapping. Furthermore, the resampled MCD12Q1 product (LAND) was included in the full model, in which four classes were differentiated in the study area: deciduous broadleaf forests (LandDB), evergreen broadleaf forests (LandEB), savanna (LandS), and wood savanna (LandWS) (Table 3). Both VIs were highly correlated (r = 0.89) and presented high VIFs of 5.2 and 3.5 for the EVI and NDVI, respectively. If used as unique independents, none of the two VIs are significant predictors of flooding. Despite having a higher Wald coefficient than the categorical class LandEB (p > 0.001), the NDVI is not significant (p > 0.260), and its VIF is above the maximum inclusion threshold of 5, as suggested by Logan (2010) . After elimination of the NDVI, a more parsimonious model was obtained, with a slightly better Akaike information criterion (AIC) of
![]()
Table 3. Full and reduced LR models for the flood prediction using the PR_0 calibration dataset (n = 1249) for the Pirizal LTSS and SESC transect water-level stages.
LandDB is the categorical variable reference class.
712.10 versus 712.84 for the complete model.
The inclusion of the MCD12Q1 product particularly improved flood classification in forest stands, which accounted for 24% of the study area. Principally, overestimates of flooded pixels in the forest stands during the transitional and dry seasons were reduced, as shown in an example from August 13, 2007 for the SESC transect (Figure 3). Although all forest stands in the Northern Pantanal felt dry during the peak of the dry season, the complete model shows few false positives (FPs), whereas FPs are frequent in the reduced model.
To minimize the nontraceable influence on model outcomes of the MCD12Q1 covariate was removed in all following evaluations on the influences of pixel quality, vegetation cover, and hydrological phase on classification performance, even though the reduced model had a lower predictive power and a reduced AIC of 755.64.
In the external validation, the reduced model showed an overall accuracy of 77.2% and a K of 0.525, with a tendency of underestimating the flood extent (FPR of 19.9% versus FNR of 27.2%) (Table 4). Both, k-fold and internal validation strongly overestimated real multiyear classification accuracy, with the OAs more than 10 p.p. higher and Ks approximately 0.25 higher than the same metrics obtained in external validation.
![]()
Figure 3. Flood classification of the SESC transect from August 13, 2007 (DOY 225), obtained by the (a) reduced (LSWI, EVI) and (b) complete (LSWI, EVI, MCD12Q1 Land Cover Type 2) models, and (c) MCD12Q1 Land Cover Type 2 (UMD): at this date, evergreen and deciduous broadleaf forests in the Northern Pantanal felt dry at least two months before.
![]()
Table 4. OA, percentage of FPRs and FNRs, and K for internal and k-fold (k = 10) validation of the complete dataset (PR_0 stratum) and external validation.
This strong reduction in the K is not only caused by the number of detected misclassifications, but also by the disparity between FPR (19.9%) and FNR (27.2%).
3.2. Model Performance as a Function of Pixel Quality
The semihumid climate of the Northern Pantanal is reflected in a strong cyclical distribution of PR levels. The percentage of lowest-quality PR values increases from the beginning of the rainy season in October (DOY 273), reaching the highest median values of approximately 30% in early December (DOY 337). During the entire 13-year MOD13Q1 TS, the highest rate of pixel degradation occurred at the beginning of January (DOY001) of 2003, when 89.7% of pixels had a PR value of 3. The highest rainfalls in the region occurred in February and March; during these two months (DOY 33 through DOY 81), 79% of the pixels had the highest PR_0 levels, and the PR_3 levels were mostly less than 7% (Figure 4).
From the second week of April (DOY 97) until the end of the dry season in the second week of September (DOY 225), almost all 16-day composites had more than 90% of PR_0 level pixels and less than 0.5% of the PR_3 level. Until the second half of October (DOY 257), PR_1 percentages increased up to approximately 49%, but with almost no PR_3 pixels.
To obtain sample sizes of ≥85 for validation, the 15 original QA strata were aggregated into seven datasets. OA and K are almost equivalent for the subsets QA_0100*-PR_0 through QA_0110*-PR_0, as sample sizes only increase in 12 and 16 pixels, respectively, if QA_0101 and QA_0110 pixels are added to the QA_0100 strata (Table 5). The performances of the PR_3 and the PR_3-QA_110 (pixels with quality better than PR_3 and lower than QA_110) decrease to significantly lower OA and K levels of 81.0 and 0.58 (PR_3 - PR_0) and 67.1 and −0.004 (PR_3 - QA_0110), respectively. The negative K means that classification agreement is less than would be expected by chance. This is caused by the misclassification of the only five nonflooded candidate pixels as flooded, resulting in the high rate of FPs (FP = 60%).
The 195 candidate pixels of the TIMESAT-PR_0 subset are classified accurately, with an OA of 89.7%, a K of 0.74, and almost equal FP and FN rates of
![]()
Figure 4. Variation of PR (PR = 0, PR = 1, and PR = 3) values during the hydrological years from 09/29/2000 (DOY 273) through 09/14/2012 (DOY 257), discharge at the Porto Cercado gauge and TRMM 3B43 precipitation averaged for the study area.
![]()
Table 5. LR model performance (internal validation of the training dataset) according to pixel quality: *includes pixels of QA_0100/QA_0101/QA_0110 level or higher without PR_0; **pixels of PR_3 without PR_0 and PR_1.
approximately 10%. OA accuracy remains high for the TIMESAT-QA_0110 subset (85.9%). Once again, K turns negative because of a 100% FP rate.
For the test sites and the period with field observations, the PR_3 and TIMESAT subsets contain 1444 candidate pixels for the aggregated strata. These two subsets include all pixel qualities and/or reconstructed values (Table 6).
As formed by more than 86% by high-quality pixels of the PR_0 subset, the
![]()
Table 6. Model accuracies (internal validation of the training dataset) as a function of pixel quality given by MODIS VI PR and VI usefulness (VI_U): the 12 classes with 0 ≤ PR ≤ 1 and VI_U between 0000 and 1100 were aggregated in six—OA, overall accuracy; K, kappa index; FP, false positives; FN, false negatives.
full PR_3 dataset has good performance metrics in the internal validation, with an OA of 86.9 and K of 0.738. It is largely outperformed by the TIMESAT set, which has a performance at least as good (88.6%/0.773) as that of the QA_0100 subset (88.3%/0.765). In all subsets, numbers of FPs and FNs are stable, ranging between 9.5% and 14.2%. The importance of TS filtering and gap filling for detection of spatial flood patterns in the study area becomes more visible in a comparison for the flood maps obtained from the classification of the PR_0 versus the TIMESAT subsets of the hydrological year 2006/2007 (Figure 5). Although unclassified pixels outside permanent water bodies at the beginning (DOY 273, 289) and end (DOY 241, 257) of the hydrological year can be considered mostly nonflooded areas and could thus be interactively filled in, no assumption about the flooding state of unclassified pixels from December through beginning February (DOY 337-033) can be made.
Despite of the good performance in the SG dataset, with less FPR in forest stands during the dry season and with OC and K values only approximately 0.6 p.p. and 0.009 lower than obtained for the PR_0 subset classification (Table 6), in the SG subset flooding is underestimated at the beginning of the wet season.
3.3. Vegetation Cover and Flooding Cycle versus Classification Performance
As shown for an internal validation of the reduced classification model based on the LSWI and EVI for the PR_0 dataset which only considers pixels with minimal atmospheric interference (n = 1249), classification performance varies by vegetation type and even more by flooding period (Figure 6).
Here, validation pixels were grouped in three major physiognomic classes according to the visual interpretation of SPOT imagery: 1) grasslands with tree cover of less than 5% had an average percentage error (PE, 100% OC) of 11.7, 2) savannah formations with 14.0, and 3) woodland/forests with 14.6. The average PE was 2.9% during the dry season (all test sites dry), 15.0% in the wet season (all test sites flooded), and 25.3% during the transition periods (>0% and <100% of test sites flooded), when, on average, approximately 52% of the validation sites
![]()
Figure 5. Flood maps as obtained from the classification of (a) the PR_0 and (b) the SG-reconstructed subsets of the hydrological year 2006/2007.
were flooded.
All vegetation types show a similar parabolic relation between flooding extent and classification error. The highest PEs of all vegetation types occur during the hydrological transition periods (25.3%), followed by the flood season (15.0%),
![]()
Figure 6. Relative error ratio (1-OA/100) and its average by flooding and predominant vegetation formation (PR_0 strata): point symbols are scaled by frequency and colored by predominant error type average. PEs by flood phase are 2.9% (dry season), 25.3% (transition), and 15.0% (wet season).
and the dry season (2.9%). For the adjusted loess model, the estimated maximum PEs are 38% when 52% of the grasslands test sites are flooded, 54% when 43% the Cerrado savanna sites are flooded, and 60% when 36% of the forest sites are flooded. Adjusted loess curves show furthermore that classification performances of both woody formations, but principally in the case of broadleaf forests, are superior during the hydrological dry season compared with the transition and wet seasons, when field stages were partially or completely flooded.
Relations between FPR and FNR are highly dependent on the flooding period, and vary little by vegetation class. During the dry season, only FPRs were observed (2.1%), and they were predominant during the transition periods, with 13.0% FPR compared with 7.3% FNR. During the wet season, only FNs (9.7%) occur.
4. Discussion
4.1. Predictor Selection and Classification Performance
The exclusive use of both tested VIs resulted in mapping models similar to a null model. The NDVI slightly outperformed the EVI. Better results obtained by the use of the NDVI (Fayne et al., 2017; Bergé-Nguyen & Crétaux, 2015) or EVI (Almeida et al., 2015; Padovani, 2010) as a single inundation predictor were probably because the study areas included large, deeply flooded floodplains and lakes with no or little vegetation cover.
In the Pantanal, such areas occur principally along the Paraguay river outside the study area. Under such conditions, performance in flood mapping is generally high if spectral indices are used that include a near-infrared band (Huang et al., 2014; Ordoyne & Friedl, 2008) . For unvegetated or little-vegetated water bodies, spectral VI generally shows an unambiguous water signal. The gauged sites of this study were flooded only between 96 and 204 days, with average depths at dates with flooding between 3.5 cm and 42.8 cm. None of the gauged sites represented open water throughout the area of an entire MOD13Q1 pixel.
If used in conjunction with the LSWI, both VIs provide significant predictors and improve classification results, with a slightly better performance for the EVI. The information from both VIs is redundant, shown by an increasing AIC if both VIs are included. As found in previous investigations (Xiao et al., 2005; Sakamoto et al., 2007) , the application of a spectral surface water index (LSWI) in conjunction with a VI results in satisfying inundation mapping models.
A previous study with multitemporal Landsat TM imagery in the Pantanal showed that spectral water indices not only vary by vegetation density and green biomass, as mainly captured by spectral VIs, but also by the structure and physiognomies of vegetation types (Zeilhofer, 2006) . Therefore, the inclusion of the MCD12Q1 LUC product as covariate shows an improvement in model adjustment (AIC = 712.1 versus 758.4).
The expressive differences in performance metrics obtained in a comparison of validation methods emphasize the importance of obtaining an independent validation subset to estimate real classification performance and to consider the applied validation approach in comparisons between different flood-mapping studies. In the approach, an equilibrated number of flooded and nonflooded pixels was used to avoid a positive bias in the validation. Specifically, imagery from August through October, the peak of the dry season in the Northern Pantanal, was not considered. According to the classifications of the driest imagery of each hydrological year, about 97% of the study area and all validation gauges felt dry. As shown by Chen et al. (2013) , a large number of dry areas biases the OA-measured classification performance of MODIS composites positively.
In this case, internal and 10-fold cross validation performed almost equally, with a difference of only approximately 1% in the OA and K. In the external validation, OA dropped 12 p.p., and 0.25 in the K, showing strong decreases in the accuracy metrics, previously reported as well for multiclass classification problems (Oldeland et al., 2010) . The differences in decrease between the OA and K metrics are attributed to an increasing disequilibrium between FPR and FNR because of the reduced sampling amount in the validation dataset, which does not influence OA but tends to reduce K, which represents the level of agreement of two datasets corrected by chance (Viera & Garrett, 2005) . The FPR/FNR ratio was 0.79 in internal validation, 0.80 in the k-fold cross validation, and 0.73 in external validation. Even with a considerable amount of independent validation data, principally K of the external validation was biased.
Considering that large open water bodies were not used as validation sites, the results are in congruence with previous studies and seem to outperform slightly the flood map accuracy obtained from MODIS eight-day composites with a 500 m spatial resolution. In a flood mapping from the Everglades, Florida, using six predictors from the MODIS MCD43B4 and MYD11A2 products, Ordoyne & Friedl (2008) obtained an OC of 83.7% for a cross-validation approach, which was slightly inferior to the cross-validation results for the gap-filled subset in this study. Using the OWL algorithm (Guerschman et al., 2011) based on the eight-day, 500 m resolution MOD09A1 product, Huang et al. (2014) achieved, for internal validations of two study areas, OAs of 97.0% and 94.8% in Australian wetlands, which were characterized by a high number of validation sites inside open water bodies. This remarkably high OA was, however, accompanied by a much lower K of approximately 0.6, indicating reduced detection accuracy. This was not explored systematically in this study, but numerous empirical tests showed that model performance could be improved by approximately 1% - 5% in the K, if the global mapping model were substituted by a local model fitting for each image (Chen et al., 2013; Heimhuber et al., 2018) . Previous studies have shown that MODIS classifications tend to underestimate flooded areas, principally in highly structured floodplains with straight channels, where flooded patches below the spatial sensor resolution are underdetected (Sakamoto et al., 2007; Islam et al., 2010) . If one opts for such an individual classification model adjustment, expert knowledge on flooding dynamics in the wetland should be included (for example, at which DOY a certain vegetation type tends to become dry) to define adequate cutoff values for each mapping model.
4.2. Intra-Annual Variability of MOD13Q1 Product Quality and Its Influence on Classification Performance and TS Smoothing
Under the semihumid climate in the Pantanal, the cyclicity of which is similar to that in other large tropical floodplains, the pixel quality of the MOD13Q1 product strongly varies throughout the hydrological year. In the northern Pantanal floodplain, deteriorated pixel reliabilities (PR_1) in the 16-day composite surge at the end of the dry season (Figure 3), end of August and beginning of September (DOY > 241), and the lowest qualities (PR_3) occur at the end of September and beginning October (DOY > 273). Lowest-quality pixels first peak at the beginning of December and remain high until the middle of January (DOY = 17). Composites begin to be built up mostly by high-quality pixels from the middle of April.
This average cycle means that the highest number of low pixel qualities occurs during the filling phase of the northern Pantanal and not during the season of highest rainfalls, which peaks in March (Figure 3(a)). Rainfall precipitation in the austral summer of the Upper Paraguay River Basin is governed by the South America Summer Monsoon (SASM) (Bergier et al., 2018; Marengo et al., 2015) , where rainfall events are of a convective nature and develop mainly in the afternoon, whereas the atmosphere may be little clouded and have low aerosol optical depth during satellite passage. On the other hand, aerosol optical depth in the Pantanal peaks in September, caused by allochthone and autochthone biomass burning in the region (Sena et al., 2013) .
The QA_0001 VI Usefulness (VI_U) level was the most restrictive filter, followed by the PR of 0 (PR_0) and QA_0010. The PR and VI_U from the two QA systems are not unambiguously related: the PR_0 subset (87% of all training and validation pixels) includes pixels tagged by the 0001 and 0010 VI_U levels. Comparisons showed that down to a VI_U level QA_0110, a subset that comprises 94% of the total dataset, OA and K remain more than 90.0% and 0.79, respectively, for the individualized subsets, and 88.2% and 0.76% for the accumulated datasets, respectively. Performance metrics show a steep drop if PR_3 pixels or pixels with a VI_U lower than 0110 are included.
Based on the validation results of pixel quality strata, recommendations for smoothing and fulfilling of degraded candidate pixels were elaborated for adaptive SG filtering, a standard TS analysis method known to perform well for overall noise reduction and for a large range of seasonal patterns (Chen et al., 2004; Kandasamy et al., 2013) .
Liu et al. (2017) showed that SG performed well for the VI cycle in a tropical monsoon forest, where vegetation greenness recovers rapidly at the beginning the rainy season, but gradually decreases as the dry season comes. The flooding extent and flooding likelihoods in the study area have a very similar periodicity. As shown here for the PR_0 and SG reconstructed series during the hydrological cycle of 2007/2008 at the Pirizal test site, in situ gauges were first flooded at DOY 337 (December 3) (Figure 7). When, approximately two weeks before, all gauges were still dry, the flooding extent jumped to approximately 87% of field gauges water covered. The flooding remained almost stable until the composite of 81 DOY (end of March), but the receding was gradual, and the Pirizal LTSS test site only became completely dry at DOY 177, about three months later. This periodicity is widely reproduced by the flooding likelihoods in the PR_0 and the SG-reconstructed datasets. In some cases, the probabilities of the smoothed dataset fit the p.p. of flooded sites even better than the PR_0 set, principally for the composite of DOY 17 and during the receding of flood water (for example, DOY 129).
Weighted moving average filtering in the SG is designed to fit the upper envelope of a TS (Li et al., 2018) . Contrary to Hird & McDermind (2009) , no tendency to overestimate high values was detected. The probability of flooding, as estimated by the SG and PR_0 TS, is almost identical during the high flood period (2007/337 through 2008/097).
The SG filter was found to perform particularly well in less noisy TS (Jönsson & Eklundh, 2002; Kandasamy et al., 2013) and to outperform other methods in
![]()
Figure 7. Probability of flooding estimated for the SG-reconstructed and PR_0 strata TS (boxes) and percentages of gauges flooded (red line) during the hydrological cycle 2007/2008 at the Pirizal test site: gauges were flooded for the first time at DOY 337 (December 3).
overall noise reduction (Hird & McDemind, 2009) . The high classification performance for the SG TS achieved here is in congruence with the strength of the method, because inundation in the Pantanal has a widely monomodal regime, which is well represented through the classification of the denoised LSWI and EVI.
Because differences in the major performance were marginal for the QA_0110 or PR_3 strata, which were identified as appropriate quality thresholds to substitute original pixels with estimated values, PR levels were used for parameter definition in the SG TS reconstruction.
4.3. Classification Performance as a Function of Flooding Cycle and Vegetation Cover
The influence of vegetation cover on classification errors in MODIS-derived inundation maps has been widely reported (Sakamoto et al., 2007; Ordoyne & Friedl, 2008) , but without quantitative assessment of classification accuracies throughout the individual vegetation formations in tropical floodplains.
Even based on an internal validation, the present study showed that the OA values averaged for vegetation types differ in a range of approximately 5% only. Because in situ measurements in single imagery are occasionally less than 10 for the PR_0 strata, no meaningful K values could be determined for these comparisons. Because this analysis was based on the PR_0 dataset, the results are not biased by pixel quality or by the validation procedure, as the relative differences in multiyear-averaged OAs between vegetation classes almost do not change between internal, k-fold, or external validations. In this study, classification performance decreases according to vegetation formations characterized by increasing tree cover or vegetation density, but less than expected according to previous studies. Islam et al. (2010) for example depicted an elevated difficulty to detect flooding under dense vegetation in the monsoon-driven floodplain system of Bangladesh.
It is believed that these results favor a combination of different factors, such as the phenological cycle of Pantanal tree species, its time coupling to precipitation/flood dynamics, and specificities of spectral surface water indices, such the LSWI. In the northern Pantanal, the highest flood stages are little delayed (see Figure 1) to precipitation maxima. Most tree species in the Pantanal are flood tolerant (Scremin-Dias et al., 2011) , and their tree canopy is fully developed during the rainy/flooding season, whereas seasonal forest and woodlands drop leaves during the dry season. Surface water indices are sensible to leaf humidity (Gao, 1996; Xiao et al., 2005) . Therefore, the true positive rate in classification may not only be caused by the real detection of flooded ground, but also by the high leaf water content. In addition, a high true negative rate in forest stands occurred, because of the dry season leaf-drop, in a period when these stands are never flooded. Thus, it is supposed that classification performance in dense forest stands in the southern Pantanal may be reduced, because the highest floodwaters occur in the beginning of the dry season with reduced foliage of the riverine broadleaf forests and woody savannahs, which have a fully developed canopy during the beginning of the rainy season, when their stands are not flooded yet.
Furthermore, the classification performance in stands without or with reduced tree cover is reduced because of macrophyte cover (Huang et al., 2014) . During the shallow flooding of grasslands on alluvial lobes of the northern wetland (mostly < 1 m), a dense herb layer of amphibian and emerging macrophytes develops (Pott et al., 2011) , which decay during the dry season. On deeper water bodies, fixed or free-floating water Eichornia hyacinths prevail, which can develop throughout the year and which can show NDVI and EVI values as high as forest stands.
The described interactions further explain the improved classification performance, if the MODIS 12 product is included in the LR model, despite its reduced spatial resolution. As shown, cycles of vegetation greenness and phenological development in the wetland differ strongly between vegetation types, and these variations are not sufficiently captured by VIs alone.
The assessment showed, furthermore, that flood-mapping performance not only depends on vegetation types, but also varies during the hydrological year, and that flooding seasonality affects classification performance more that the proper vegetation cover. Whereas the highest vegetation densities in the northern Pantanal coincide with the climatic wet season during flooding, the lowest classification accuracies were observed during the transition phases, with approximately two-times-higher error rates than in the wet season. Comparative studies showed that the flood-mapping performance of MODIS imagery is reduced because of spatial resolution (Chen et al., 2013; Heimhuber et al., 2018; Kuenzer et al., 2019) , because one pixel may include flooded and non-flooded areas. In the study area, misclassifications during high floodwaters and the transition periods occur because of the assemblage of small stream channels and forest-covered levees, principally at the braided alluvial fans, represented by the limnimetric gauges of the SESC transect (Fantin-Cruz et al., 2011) .
Previous studies have also attributed misclassifications to the shallowness of flooding (Dao & Lou, 2015) or sediment content (Gibbs et al., 2016; Garcia et al., 2015) , which can be elevated by wind interaction in shallow waters. Ordoyne & Friedl (2008) have applied MODIS TS classification not only to map flooding but for estimates of flood depths. They observed higher flood probabilities in their LR at more-deeply-flooded areas. Therefore, the steep performance loss during the transition period to subpixel variations is attributed to the low flood heights over higher-lying lobes near the margins of the wetland, partially inundated by stagnant rainwater. Such areas are represented by the Pirizal test site, where flooding depths have been approximately 60% less than 20 cm. Unpublished simulations using a high-resolution digital elevation model showed that if a gauge is flooded, parts of the corresponding pixel may be dry because of the microrelief, even inside a same vegetation unit. According to these findings, it is recommended that future studies should evaluate whether parameter definition in TS reconstitutions can be done not only according to pixel quality, but also by a priori information on hydrological period and vegetation type in the developing multidimensional convolution coefficients (Chen et al., 2004) .
Classification performance during the transition period is further diminished if the high percentages of low-quality pixels are maintained, e.g. in classifications without prior application of a TS reconstruction technique. In application studies, which use classifications as inputs for further wetland assessments, this should be considered. Important ecological processes occur principally during this period, such as the steep oxygen depletion phenomena, regionally known as “Decoada,” when, by flooding, organic material begins accumulating from the dry season and decomposes, causing high fish mortality (Calheiros & Hamilton, 1998) .
5. Conclusion
To improve our knowledge of how pixel quality, vegetation cover, flooding stage, and validation procedures affect inundation classification performance in large tropical floodplains, multiyear in situ water-level data were compared with mapping results obtained for the corresponding MOD13Q1 TS, and the following conclusions were drawn.
1) Under the monomodal climatic and hydrologic regimes in the Pantanal, the application of TS reconstitution techniques is highly recommended to obtain gap-free flood maps, which can be achieved without significantly decreasing mapping performance. Despite its strong simplification, in large tropical floodplains with a wet/dry climate, pixel QA given by PR should be, in most cases, sufficient to parameterize TS reconstitution algorithms, such as the SG used here. Thus, PR = 3 pixels should not be considered, and PR = 1 should be weighted for the definition of convolution coefficients.
2) Both vegetation cover and, even more importantly, hydrological period affect flood classification accuracy. The least accurate results were obtained for areas of dense forests/woodlands during the hydrological phases of filling and drying, which were mainly attributed to shallow flooding and/or the partial flooding inside pixels that occurred during these hydrological phases in large parts of the study area. Therefore, it is recommended for future studies to include information on vegetation cover and hydrological period (or expected flooding depth) in TS reconstruction efforts for gap-free flood classification in large tropical wetlands.
3) Validation metrics are equal for k-fold, 10-fold, and internal validation and heavily overestimate true error rates in more than 10 p.p. (OA) and 0.25 (K).
4) Despite its reduced spatial resolution, the use of the MCD12Q1 improves flood-mapping performance because it captures vegetation attributes that are not sufficiently represented by VIs.
5) The major limitation of this research was the delineation of the study area to include only the Northern Pantanal, centered around the available field validation stations. The middle and southern Pantanal has additional vegetation formations and its flooding regime can be delayed by more than three months in comparison with the precipitation peaks.
6) Due to the large differences in outcomes as a function of validation method, which showed a decrease of about an quarter in OA, an operational network of in situ flood measurements within the wetland should be implemented. Newly available satellite systems with a simultaneous high spatial and temporal resolution should be tested for their performance of flood monitoring in the Pantanal.
Acknowledgements
We gratefully acknowledge the support of National Council for the Improvement of Higher Education (CAPES).