Classification of Wetland Vegetation Based on NDVI Time Series from the HLS Dataset

Abstract

Natural wetlands are intrinsically heterogeneous and typically composed of a mosaic of ecosystem patches with different vegetation types. Hydrological and biogeochemical processes in wetlands vary strongly among these ecosystem patches. To date, most remote sensing classification approaches for wetland vegetation either rely on coarse images that cannot capture the spatial variability of wetland vegetation or rely on very-high-resolution multi-spectral images that are detailed but very sporadic in time (less than once per year). This study aimed to use NDVI time series, generated from NASA’s HLS dataset, to classify vegetation patches. We demonstrate our approach at a temperate, coastal lake, estuarine marsh. To classify vegetation patches, a standard time series library of the four land-cover patch types was built from referencing specific locations that were identified as “pure” pixels. These were identified using a single-time high-resolution image. We calculated the distance between the HLS-NDVI time series at each pixel and the “pure”-pixel standards for each land-cover type. The resulting true-positive classified rate was >73% for all patch types other than water lily. The classification accuracy was higher in pixels of a more uniform composition. A set of vegetation maps was created for the years 2016 to 2020 at our research site to identify the vegetation changes at the site as it is affected by rapid water elevation increases in Lake Erie. Our results reveal how changes in water elevation have changed the patch distribution in significant ways, leading to the local extinction of cattail by 2019 and a continuous increase in the area cover of water lily patches.

1. Introduction

Coastal wetlands provide critical ecosystem services, such as wildlife habitat, fisheries support, carbon sequestration, flood protection, and water quality improvement. Specifically, around the US Great Lakes, as Harmful Algal Bloom (HAB) events are intensifying in extent and frequency, coastal wetlands are becoming increasingly important for mitigating nutrient runoff from agriculturally dominated watersheds to the lakes [1]. Natural wetlands are intrinsically heterogeneous and are typically composed of a mosaic of ecosystem patches with different plant types and hydrological conditions. The adaptation of these plant communities to a water-dominated environment is the basis for their use in improving the water quality in constructed wetlands [2]. Understanding the effects of wetland vegetation on the environmental conditions and ecological function of the wetland is key to determining which plants to grow in a constructed wetland with regards to the specific goals of the wetland-construction project (e.g., maximizing nutrient removal), and to correctly model the global contribution of wetlands to global greenhouse gas emissions and carbon sequestration. Wetland vegetation generates drag on the flow and therefore can influence the water movement through and around a vegetated patch [3]. The plant density and life forms affect the drag and thus control the flow velocity and residence time of the water in different parts of the wetlands, as well as the rate of deposition of suspended solids [4]. Furthermore, emergent plants with high transpiration rates can lower the water level [5]. The rates of nitrogen processing also vary among wetland vegetation communities [6]. Emergent vegetation shows a higher nitrogen retention rate and is more efficient at removing nitrogen than submerged macrophytes [7]. Identifying patch types is particularly important in predicting their contribution to nutrient uptake and sedimentation in coastal wetlands in agricultural watersheds, where the wetlands provide a last “line of defense” to improve runoff water quality. Vegetation uptakes carbon by photosynthesis, and different plants have different characteristic leaf areas and photosynthesis rates. Furthermore, vegetation plays an important role in methane transport, and, thus, the vegetation patch types have a strong impact on carbon sequestration and greenhouse gas flux rates in wetlands [8,9,10,11].

Traditionally, mapping the vegetation communities in a wetland requires intensive fieldwork, and some wetland areas are not easy to access. Remote sensing provides an effective, scalable, and economical alternative by which to monitor and map plant communities in wetlands. Many studies have been conducted to discern wetland vegetation communities using remote sensing imagery. Multispectral remote sensing data, for example, Landsat imagery and hyperspectral remote sensing imagery, have been used to discriminate and map wetland vegetation types [12,13,14,15,16,17]. However, classifying wetland plants comes with additional challenges in comparison to classification of dryland terrestrial plants [18]. Specifically, the boundaries between vegetation patches within a wetland are difficult to identify because of the intrinsically small scale of the wetlands and patches within them, combined with the short ecotones and sharp demarcation between vegetation patches, which is caused by the small scale of the spatial variability of wetland vegetation [17,19,20]. Thus, moderate-resolution remote sensing imagery, for example, MODIS, which has 250 m spatial resolution, is usually not sufficient for differentiating wetland vegetation patches.

The other challenge in classifying wetland vegetation is that the reflectance of different types of wetland canopies is often very similar, and these indistinct reflectance signatures are combined with the reflectance of the underlying soil and water surface [21,22,23]. To overcome this difficulty, high spatial resolution imagery has been used to perform the classification to reduce the spectral mixing effect of different wetland vegetations and underlying soil or hydrological regimes within single pixel [24,25,26]. However, high spatial resolution imagery is typically only available for a given wetland (a small location) with a very sparse return time (at best, 2–3 cloud-free images per year, usually less, and, in most years and locations, none). Therefore, it is not guaranteed that the imagery covers the wetland of interest during the growing season when the vegetation is visible. Alternatively, the spectral difference of wetland vegetation at different times of the year has been demonstrated to be effective in discriminating dryland vegetation types globally [27] and, specifically, in arid/semiarid regions’ dryland vegetation [28]. This approach has been demonstrated in wetlands for the identification of individual target wetland plant species in salt marshes [28,29,30,31]. Similar to our work, Sun et al. [32] used a time series of the Normalized Difference Vegetation Index (NDVI) derived from HJ-1 at a 30 m resolution to classify and map multiple vegetations at a salt marsh in China, but did not continue to monitor vegetation changes over a long (multiple years) period. Sun et al. [32] further indicated the Landsat data are suitable for similar analyses. Not every remote sensing dataset is suitable for constructing time series in terms of wetland vegetation classification. Moderate spatial resolution data have more frequent and regularly repeated images of the same location, which is important for constructing a time series, but its spatial resolution is not high enough to discern the patches of different wetland vegetation types. Moreover, the frequent cloud cover in coastal region results in significant data gaps in remote sensing data [32]. The harmonization of the Landsat and Sentinel-2 (HLS) data by NASA has made it possible to acquire imagery at a moderate-to-high spatial resolution and at a high temporal resolution, which has created the opportunity to build time series of wetland vegetation characteristics at sufficient spatial and temporal resolutions. NDVI is one of the most widely used indices to assess vegetation from remote sensing imagery, and the construction of the NDVI time series has enabled the quantification of the penological dynamics of different wetland vegetations.

The overall goal of this study was to develop a classification approach for wetland vegetation that utilizes vegetation-type characteristic plant phenology, as observed by the NDVI time series from NASA’s HLS dataset. The classification approach that we demonstrate and test here could be used to classify wetland vegetation on a large scale.

2. Materials and Methods

2.1. Study Area

The study was focused on Old Woman Creek (OWC) National Estuarine Research Reserve (NERR) (Figure 1a). It is a temperate, mineral soil, estuarine marsh, located at the south-central shore of Lake Erie, longitude: −82°30′3″, Latitude: 41°22′39″. It is an Ohio State nature preserve and is managed as a cooperative partnership between NOAA NERR and the Ohio Department of Natural Resources (ODNR). It experiences rapid fluctuations in hydrology and water levels as the river mouth is periodically blocked by a sand barrier that is deposited at the lake shore and affected by additional long-term (interannual) water level fluctuations due to the recent water elevation increase of Lake Erie (~1 m over the last decade). OWC was selected for this study because it represents many coastal wetlands of the Great Lakes region. It includes 61 Ha of natural wetland area, currently consisting of four characteristic vegetation patch types. These patch types are defined through their vegetation status, i.e., (1) open water; (2) mud flats; or subclassified among types of patches with emergent vegetation according to the most common (by area) emergent vegetation species. These vegetated types are: (3) American water lily (Nymphaea odorata), (4) yellow lotus (Nelumbo lutea), (5) cattail (several species and hybrids of Typha spp.), and (6) reed (several species of Phragmites spp.). During the past 10 years, patches of reed as well as mud flats have been present, although they occupied a very small (<5%) fraction of the area. However, due to the rising water levels of Lake Erie and the consequent deepening of the wetland water, the mud flats disappeared by 2016, reed in 2018, and cattail in 2020. As our analysis centers on 2017–2020, the only patch types we considered were: open water, water lily, lotus, and cattail. Long-term meteorological, hydrological, and water quality data for OWC are available through the NERR data interface https://cdmo.baruch.sc.edu/ (31 December 2021); long-term radiation, energy, and carbon flux data are available through the Ameriflux database, site ID US-OWC [33].

2.2. Remote Sensing Data

NASA’s Harmonized Landsat and Sentinel-2 (HLS) dataset [34] was used in this study to generate an NDVI time series. This dataset provides consistent surface reflectance from the Landsat 8, Sentinel-2A, and Sentinel-2B satellites. The merged product was run though a set of algorithms that include atmospheric correction, cloud and cloud-shadow masking, spatial co-registration and common gridding, illumination and view angle normalization, and spectral band-pass adjustment. The HLS dataset contains two products of surface reflectance at a 30 m resolution: HLS S30 from the Sentinel-2A and Sentinel-2B satellites, and HLS L30 from the Landsat satellite. These two products are gridded into the same tiling system, the Military Grid Reference System (MGRS), so that they can be stacked together to generate time series with more frequent observations. OWC is fully covered by the tile at location code 17TLF. We accessed the data through https://hls.gsfc.nasa.gov/data/ (31 December 2021).

WorldView-3 high-resolution multispectral imagery was used to provide a ground-truth reference by supervised classification. We purchased all the WorldView-3 images that included OWC, where >80% were cloud free and from during the growing season (July–September). Only one image from 17 August 2017 complied with these standards (011281167020_01_P001_MUL). The image was purchased from Apollo Mapping. It contains 8 multispectral bands with 1.2 m resolutions.

2.3. NDVI Time Series Construction

We utilized the growing season of 2017 (days 134–296 of the year), when we had a high-quality WorldView-3 image of the site at peak growing season and high confidence in the results of its expert-supervised classification based on a detailed ground-based survey we conducted that year, and obtained all the available HLS images for that period. We preprocessed the data by extraction of the red and near-infrared (NIR) bands for calculation of the NDVI values and the quality assurance (QA) band for filtering of the cloud pixels. We then clipped each image in the dataset to the extent of OWC’s boundary. The NDVI for each pixel of each clipped image was calculated as follows:

NDVI = (NIR − Red)/(NIR + Red).

(1)

The cloud cover was filtered by indexing of the cloud pixels from the QA band. This resulted in a quality-controlled time series of the NDVI values per pixel at a weekly temporal resolution. We used linear interpolation to gap-fill these NDVI time series over gaps due to clouds.

2.4. Classification Standards from WorldView-3

WorldView-3 imagery was clipped to OWC’s extent and then classified into four categories by the QGIS using the “dzetsaka” classification tool plugin version 3.4.3 (https://github.com/nkarasiak/dzetsaka (31 December 2021), https://www.karasiak.net/dzetsaka-how-to-make-your-first-classification-in-qgis/ (31 December 2021), [35]) with the classifier set as a Gaussian Mixture Model (GMM), using all the multispectral bands contained in the WorldView-3 imagery with 40 training samples. We tested the other classifiers available through QGIS (Nearest Neighbor, Support Vector Machine, and Random Forest), but they resulted in lower accuracy, and thus we selected GMM. To evaluate the classification accuracy, we compared it with ground observations at 100 randomly selected points that were classified by visual interpretation and ground-based site knowledge. For the area that correspond to each HLS pixel (roughly 25 × 25 WorldView-3 pixels), we calculated the fractions of the area covered by each of the four land-cover types as classified by the WorldView-3 image. We identified HLS pixel areas that contained predominantly a single land-cover type (with a coverage fraction greater than 80%), i.e., “pure” pixel areas. Five “pure” pixel areas were selected at random for each land-cover type. The HLS pixels that corresponded to these WorldView-3 “pure”-pixel areas were used as a reference to construct a standard HLS-NDVI time series per land-cover type.

2.5. Classification Using Dynamic Time Warping

With these standard NDVI time series for each land-cover type, we used the Dynamic Time Warping (DTW) approach to calculate the similarity and distance between each HLS pixel in the dataset and the four standards. The Dynamic Time Warping algorithm was first introduced and explored in the field of speech recognition [36], but it was then applied in many other fields, including time series clustering specifically applied to land-use and land-cover mapping from remote sensing data [37]. DTW measures the similarity between time series and minimizes the effects of shifting and distortion over time. Particularly suitable to phenological time series, the approach by Maus et al. [37] can handle time series with different lengths. Considering that the length of the NDVI time series from different years can vary, DTW was selected in this study to measure the similarity between the time series to ensure that the standard NDVI time series built from one year can be applied to other years to conduct the classification. We classified the vegetation throughout the wetland during the growing season (May–September) of each year, 2016–2020. For each pixel, the land-cover that had the smallest distance was assigned to each HLS pixel as its dominant type. The complete process of classification using the NDVI time series from HLS is illustrated in Figure 2.

2.6. HLS Classification Validation and Partial Composition of Patch Types within HLS Pixels

We used the classified WorldView-3 image to validate the HLS classification result and the classified type in HLS pixels to the most common patch type among WorldView-3 pixels within the area of the HLS pixel. The true-positive rate was calculated as the count of HLS pixels classified as the same class as the majority of the co-located pixels classified from WorldView-3 and divided by the total amount of the HLS classification pixels that correspond to that class among all the HLS pixels. We computed the fraction of each patch type within each HLS pixel by counting the number of classified WorldView-3 pixels that belong to that particular patch type within the corresponding HLS pixel area and dividing it by the total number of WorldView-3 pixels within the same HLS pixel area.

3. Results

Classification of WorldView-3 Image for Identification of “Pure” Pixels

The classification results of the WorldView-3 image (Figure 1b) are shown in Figure 1c, and a confusion matrix (

Table 1. Confusion matrix of the WorldView-3 classification result for 2017. The classified result was compared to 100 randomly selected points and classified by visual interpretation and ground-based site knowledge.

		Classified As				Total	User Accuracy
		Open Water	Water Lily	Lotus	Cattail
True Patch Type	Open Water	16	0	0	0	16	100%
	Water Lily	3	36	3	14	56	64%
	Lotus	0	0	8	2	10	80%
	Cattail	0	0	0	18	18	100%
Total		19	36	11	34	100	78%
% True		84%	100%	73%	53.0%

) was calculated to evaluate the classification accuracy of the WorldView-3 imagery. The overall accuracy was 78%; the open water and cattail patch types had the highest precision of 100%, and the water lily patch type had the lowest precision of 64%. Most of the water lily patch type that was wrongly classified was classified as cattail. Water lily had the least contiguous distribution pattern compared to cattail and lotus (Figure 1c), which might be the reason for its misclassification. Many of the vegetation patches, especially the water lily patches, were very small, and the wetland included many transitions from one type to another even within the very high resolution of the WoldView-3 image, leading to a large number of mixed pixels that challenged the classification and led to an overall accuracy of 78%.

Using the HLS-NDVI time series of the selected pixels that were identified as “pure” pixels (>80% of the corresponding WoldView-3 pixels were of a single patch type), we defined a characteristic NDVI pattern for the four land-cover types in OWC (Figure 3). Water lily and cattail showed a very similar temporal pattern, but cattail had a consistently higher NDVI throughout the season. The lotus patch type was distinguishable from these water lily and cattail patch types as it started with a lower NDVI early in the spring and peaked at a higher NDVI early in the summer. The open water patch type showed a similar pattern to the other three vegetation types but peaked at a lower NDVI intensity and later in the summer (Figure 3). Using the temporal signature of each land-cover type NDVI time series, we classified the pixels in the weekly repeated HLS images into different types. The HLS classification map is presented in Figure 1d. Spatially, the distribution of the vegetation patches in the HLS classification agrees with that of the high-resolution WorldView-3 imagery. The cluster of lotus patches at the north part along the main channel of OWC, the open water and main channel, and the water lilies at the outflow of OWC were captured in the HLS result.

Many of the HLS pixels were heterogeneous and included a mixture of different patch types, with the classified patch type only contributing to a fraction of the full area of the pixel (Figure 4). In such cases, the classified patch type was assumed to represent the type that covered the largest fraction of the pixel area. Water lilies tended to grow in small patches, which resulted in highly mixed HLS pixels for the water lily patch type (Figure 4b). Furthermore, water lily leaves float on the surface of the water and are not truly emergent. As a result, the water lily classification accuracy, at 50–64%, was lower than that of the other patch types (Figure 5). Overall, for the three other patch types, when an HLS pixel covered an area that had more than 60% of a certain type of land cover, the overall accuracy was higher than 80% (Figure 5). The classification accuracy tended to increase when the fractional cover increased, which indicates that when a pixel is less mixed, it is more likely to be classified in the correct class.

We then classified the vegetation at OWC in 2020 (Figure 6a) using the DTW approach for the HLS data and the characteristic NDVI time series we built from the 2017 HLS pure pixels. There were no cloud-free WorldView-3 images available for OWC during the growing season of 2020. The 2020 HLS-classification confusion matrix (Figure 6b) was created by comparison to the ground-based sampling of locations throughout OWC, which was conducted on 20 July 2020, and included 30 points selected at random throughout the wetland. Lotus showed a higher classification precision (80%) than water lily (66.7%). The overall accuracy for classification in 2020, when no high-resolution images were available, was 75.8% (Figure 6b), which validates the potential of using NDVI time series from HLS data to classify wetland vegetation. We used our approach to provide a consistent and continuous annual classification of the land-cover types within OWC (Figure 7a) since 2016, when HLS data became available. Our classified maps clearly show the ecological changes that the wetland has undergone due to the changes in mean water elevation (Figure 7b,c). Most notable is that the patches of emergent cattail vegetation had completely disappeared from the wetland by 2019.

4. Discussion

The time series of NDVI that we observed in the pure pixels defines each patch type (Figure 3). The differences between patch types provide some ecological insight with regard to their seasonal dynamics. The cattail and water lily classes show similar NDVI trends, except that the overall mean NDVI of cattail was higher than that of water lily. Lotus started growing later than cattail and water lily, but its spring green-up was faster, such that by around day the year (DOY) 175, its NDVI value exceeded that of cattail and stayed the highest among all the land-cover types. Its fall senescence started at about DOY 210, as with all the other types, but it showed the fastest browning rate, and, by around DOY 280, its NDVI became lower than that of cattail. Open water always showed the lowest mean NDVI value and followed a similar temporal NDVI trend to that of the vegetated land-cover types. We hypothesize that the reason for this similarity is that water pixels are rarely truly pure and almost always include some vegetation within the pixel.

The reflectance of a mixed pixel in remote sensing imagery is a combination of the reflectance of the land-cover types that are present in the pixel. The linear spectral mixture model is a widely used technique to unmix remote sensing pixels [38], and it assumes that the reflectance is a linear mixture of the reflectance from all of the land cover in the pixel. Therefore, the more a certain type of land-cover is present in a pixel, the more dominant the effect it has on the mixed reflectance and thus the NDVI value. Mixed pixels could influence our confidence in the location of the identified “pure” pixels. However, there is a scaling issue to consider. The “pure” pixels correspond with pixels of the coarser HLS. They are identified as areas within the high-resolution world view images that are continuously and uniformly classified as belonging to a single patch type. By choosing “pure” pixel locations that are as far as possible from the transition between one World-View3-based classified patch type to another, we effectively minimized the effect of mixed pixels within the high-resolution image for “pure” pixel location detection.

Nonetheless, mixed pixels within the classified HLS images remain a concern and a source of error and potential bias. For water lily and cattail, the higher percentage pixels at the very north part of the wetland were classified as their corresponding classes, but for those pixels at the southern part of the wetland, the fractions of cattail and water lily were often around 50% to 60%, with some lower percentage lotus (Figure 1d). As a result, water lily showed the lowest precision of 54.5% (Figure 5). Future progress may consider a spectral unmixing approach [39], treating the NDVI time series of pure pixels, rather than a spectral reflectance series, as the signature of each end member and estimating the mixture fractions of each end member within each HLS pixel. We predict that this approach will be advantageous when the mixing is simple and includes only two end members (patch types) but will lead to larger errors in pixels of a more complex mixture. In cases when some patch types are intrinsically characterized by small or narrow patches with a typical size much smaller than the ~30 m resolution of HLS (for example, open water mostly in narrow channels), the classification of mixed pixels can lead to a bias against these patch types. Nonetheless, for many applications related to modeling of wetland processes, where patch-level processes are resolved at the sub-grid-scale without explicitly tracking the patches’ locations, the purpose of the patch-type classification is to inform the model of the aerial distribution of patch types within the full wetland. In such an application, the site-level statistics over all the mixed pixels will average out classification errors, as the probability of misclassification is proportional to the fraction of the partial area of the large patch type (Figure 5). Another potential cause of misclassification that needs to be taken into consideration is the cloud detection omissions of the Sentinel-2 cloud mask [34]. An outlier-detection-and-removal method is needed in the future to improve the accuracy of the classification.

The areas most affected by mixed types within the same pixel are those where the characteristic size of the patches is small, with transitions from one type to another at a short spatial scale. For example, the south part of OWC does not show a tight alignment of the two classification results. In this part of OWC in 2017, small patches of all three types of vegetation were growing next to each other. The HLS pixels in this area of the wetland, therefore, always represent a mixture of at least two land-cover types. Most of the mixtures were identified as cattail in the HLS classification, except some of the pixels with a low cattail contribution that were classified as water lily or lotus (Figure 1d and Figure 4). Part of the main channel was not recognized because the main channel is narrow, and therefore, in many parts of the main channel, the HLS pixel area that covers the main channel is mixed with vegetation that contributes to a high NDVI, which resulted in the pixels along the main channel being misclassified as vegetation (Figure 1d, Figure 4 and Figure 7). The pixels of the HLS imagery at the edge of the wetland are mixtures of water and forest (non-wetland vegetation types with a higher NDVI than the wetland vegetation), which challenged the classification of some pixels at the edge of the wetland.

OWC is a Lake Erie coastal wetland and is experiencing rapid ecological and hydrological changes due to changes in Lake Erie’s water elevation. The Great Lakes system has experienced rapid water elevation increases in the past decade. Lake Erie is an extreme case of the Great Lakes, and its water elevation has increased by about 1 m since 2011 (Figure 7c). Since the wetland is hydrologically attached to the lake, its corresponding water elevation also increases, although the wetland experiences much greater short-term variation than the lake due to the river and estuary hydrodynamics. This rapid water elevation increase drives a rapid shift in the patch type composition of the wetland. The consistent availability of HLS images from 2016 up to present allow the characterization and quantification of these changes (Figure 7a).

5. Conclusions

We demonstrated an approach used to classify vegetation patch types using repeated NDVI data from HLS in a small, heterogenous wetland. We demonstrated the accuracy of our approach over multiple years and its applicability for long-term monitoring of the vegetation patch type dynamics in wetlands. Advanced models of wetland function, and particularly mechanistic models, that predict nutrient processing, sedimentation, and greenhouse gas budgets at the resolution of plant functional types require knowledge of the areal distribution of different plant functional types within the wetland, their seasonal dynamics, and their longer-term (inter-annual) dynamics (e.g., Figure 7). Remote sensing-derived input for patch type distribution could improve the accuracy of the prediction of the ecological function of wetlands, specifically where they are experiencing ecological and hydrological changes and disturbances. While high-resolution multispectral images have been previously demonstrated by others and by ourselves as a reliable source of information for vegetation patch location and identity, high-resolution images may not be consistently available for a given wetland site in every year and growing season.

We demonstrated the applicability of our approach for long-term monitoring in Old Woman Creek (OWC), a small, coastal estuarine, mineral soil marsh of Lake Erie. As Lake Erie is experiencing a rapid water level rise (roughly 1 m over the last decade), OWC, which is hydrologically connected to the lake, is becoming deeper and experiencing rapid changes in its ecological patch location and distribution. It has switched from a cattail-dominated ecosystem to one absent of cattail and dominated by floating-leaf vegetation (water lily and lotus).

The approach proposed in this study, including the characteristic time series we developed, can be directly applied to other wetlands that grow the same type of plants using our characteristic NDVI time series (Figure 3). In the case that other wetlands of interest support patch types that are characterized by different emergent plant species than the ones growing in the current study site, the method we have demonstrated in this study could be still applicable to other types of vegetation, provided that the locations and times of pure pixels are available, from which a characteristic time series could be derived for each patch type. Such pure-pixel locations and times could be identified from limited ground surveys or from expert classification of patch types at one point in time from a single available high-resolution image.