Using Drone-Captured Imagery and a Digital Elevation Model to Differentiate Eelgrass Species: Padilla Bay, Washington

In the Pacific Northwest, eelgrass offers critical habitat for salmon, crab, shellfish, and sea birds (Bulthuis, 1995; Thom et al., 2003) and provides a variety of ecosystem services, including stabilizing sediment, buffering storm surge, and combating coastal erosion (Nahirnick et al., 2019; Nordlund et al., 2016). In Washington State, Padilla Bay National Estuarine Research Reserve (PBNERR) is home to one of the largest beds of eelgrass, approximately 3240 hectares, in the contiguous United States (Bulthuis, 1995). PBNERR has been conserved for future estuarine research, education, and monitoring. The eelgrass population in Padilla Bay and the greater Pacific Northwest is primarily composed of two species: Zostera marina and Zostera japonica (Bulthuis, 1995). Z. marina, the native and most abundant species, is found throughout the bay’s mid- to subtidal regions (Bulthuis, 1995). Z. japonica, a nonnative eelgrass, is found throughout upper to middle intertidal ranges (Bulthuis, 1995; Kaldy, 2006; Kaldy, Shafer, and Magoun, 2015; Kim et al., 2016). Originally from Japan, Z. japonica is believed to have reached North American shores as an aquatic hitchhiker in the packing materials of clam exports (Bulthuis, 1995; Kaldy, 2006; Shafer, Kaldy, and Gaeckle, 2013).

Eelgrass, nearshore and offshore, can occur in expansive stands that are hard to access (Neckles et al., 2012). Monitoring large, multispecies eelgrass populations, like the ones at Padilla Bay, can be very labor intensive (Neckles et al., 2012). Neckles et al. (2012) suggested a multitiered approach for eelgrass research and monitoring. Teir 1 involves landscape approaches that target areas for more detailed (tier 2) monitoring. Tier 2 monitoring is performed using permanent plots for foot, boot, and diving surveys. Tier 3 surveys provide even more detail and capture a variety of measurements for eelgrass at a given site (Neckles et al., 2012). The tiers of monitoring give an accurate understanding of eelgrass presence, health, and systemwide responses to change on a large scale (Neckles et al., 2012). Remote-sensing analysis of intertidal and submerged eelgrass via overhead image capture is an effective tier 1 approach to identify areas for more detailed surveys (Nahirnick et al., 2019; Neckles et al., 2012; Young et al., 2008). Unoccupied aerial systems (UAS) have become a desirable tool used for tier 1 eelgrass monitoring to capture imagery at a lower cost and at a higher spatial and temporal resolution than is available from traditional aircraft or satellites (Duffy et al., 2018; Nahirnick et al., 2019). As such, UAS imagery can be used for remote-sensing analysis to answer questions regarding the spatial and temporal dynamics of eelgrass cover. Answering these questions through other approaches is challenging because detailed and frequent data cannot be collected as easily (Duffy et al., 2018; Nahirnick et al., 2019). There are multispecies stands of eelgrass throughout the world, and species differentiation in these stands using remote sensing has been challenging. As such, most UAS projects in temperate systems have not been successful in differentiating between Z. japonica and Z. marina and were limited to detecting eelgrass presence and mapping spatial extent (Duffy et al., 2018; Nahirnick et al., 2019). However, recent work by Hein (2022) was successful in using UAS imagery to differentiate between areas dominated by Z. marina or Z. japonica with an overall accuracy of 70% using a fusion of different image-processing and classification techniques. Hein’s (2022) work also builds on earlier research by Nahirnick et al. (2019) in contributing to the understanding of how UAS image quality and successful mapping of eelgrass are influenced by sun angle, cloud cover, wind speed, and tide stage. Careful consideration of each of these factors is critical for developing flight plans that will minimize glare from the water surface and yield usable imagery.

The distribution of Z. marina and Z. japonica is strongly influenced by elevation (Bulthuis et al., 2021; Ruesink et al., 2010). Z. marina has been primarily observed in Padilla Bay between 0.3 m above and −3 m below mean lower low water (MLLW) but can be found as high as 0.6 m above MLLW (Ruesink et al., 2010). Z. japonica is usually found in Padilla Bay between 0.3 and 0.8 m above MLLW (Bulthuis, 1995; Thom, 1990) and has been found at elevations as high as 1.5 m above MLLW (Ruseink et al., 2010). Furthermore, elevation data were used in past research by Young et al. (2008, 2015) as a factor to differentiate Z. marina and Z. japonica cover in aerial imagery acquired from Yaquina Estuary near Newport, Oregon.

Elevation is correlated with other environmental factors that influence Z. japonica or Z. marina presence, including temperature, light availability, and exposure to air (Kaldy, Shafer, and Magoun, 2015; Ruesink et al., 2010; Shafer, Wyllie-Echeverria, and Sherman, 2008; Shields, Moore, and Parrish, 2018). Z. marina in the Pacific Northwest has a narrow optimal temperature range of 5–8 °C and experiences high stress and a loss of productivity at temperatures of 15 °C and higher (Kaldy, Shafer, and Magoun, 2015; Thom et al., 2003). In contrast, Z. japonica can tolerate warmer temperatures, with optimal growth between 15 and 25 °C (Kaldy, Shafer, and Magoun, 2015). Since Z. japonica has a higher tolerance for warmer temperatures and desiccation stress, it occupies a habitat in the Pacific Northwest that was once barren mudflats (Ruesink et al., 2010; Shafer, Wyllie-Echeverria, and Sherman, 2008). Although Z. marina and Z. japonica have slightly different habitat ranges, there are numerous areas where these species overlap in Padilla Bay (Bulthuis, 1996). Ruesink et al. (2010) suggested that Z. marina and Z. japonica both have the most success in the middle elevations of their habitat ranges. However, in zones dominated by one species, subtle, centimeter-level, changes in elevation can be favorable for the other species, creating islands of Z. marina in depressions or Z. japonica on mounds throughout the intertidal area (Ruesink et al., 2010).

Although elevation has been shown to have a major influence on the distribution of these two species, it can be challenging to obtain elevation data with sufficient detail. Bathymetric (blue-green) LIDAR has been used to measure benthic elevation because these sensors can penetrate the water’s surface if the water is clear (Webster et al., 2014). However, there is error in LIDAR elevation measurements in dense eelgrass stands, where the resulting measurements will represent the elevation of the top of the eelgrass mat, but the eelgrass mat may be floating some distance from the bottom (Webster et al., 2014). Further, the overall vertical error of this measurement strategy typically ranges between 15 and 20 cm (Cain and Hensel, 2018; Webster et al., 2014). Blue-green LIDAR was used to assess elevation in Padilla Bay in 2014 and resulted in a resolution of 2 m and vertical error of approximately 25 cm (ERDC, 2015). Another approach involves using global navigation satellite system (GNSS) equipment. GNSS equipment can provide measurements with a vertical error of 2 mm to 4 cm (Chang, Ge, and Rizos, 2004; Pingel, Saavedra, and Cobo, 2021), but it can be challenging to obtain measurements with a density that is sufficient to generate a detailed digital elevation model (DEM). For instance, previous elevation mapping conducted by the U.S. Geological Survey in Padilla Bay (Thorne et al., 2015) included bathymetry GNSS measurements every 100 m in the southern region of Padilla Bay, but these did not have the detail needed to generate a highly detailed DEM.

This project’s objective was to increase the accuracy of detecting Z. japonica and Z. marina cover in Padilla Bay, further quantifying the spatial distribution and seasonal dynamics of each species in this region. Successful mapping of eelgrass provides a baseline to monitor population health, understand distribution seasonally and interannually, ensure protection of eelgrass meadows, and assess impacts to eelgrass related to climate change. Differentiation of eelgrass species adds further insight into species-specific distribution, intraspecies competition, stress response, and distribution dynamics. Further, accurate mapping of Z. japonica and Z. marina can provide more accurate projections of change for eelgrass populations, target research and monitoring, and improve eelgrass conservation tactics and legislation (Neckles et al., 2012).

Analysis was based on the use of multispectral imagery captured using UAS (with methods outlined by Hein, 2022) and a highly detailed DEM with a resolution of 1 m that was created through the collection of hundreds of measurements using a GNSS probe. Predictions of eelgrass cover were made using multispectral imagery alone and by incorporating elevation data from the DEM. The best result was used to quantify variation in the spatial and temporal distributions of these two species of eelgrass in the study area for the months of April–September 2022. Data and results of this project will be used to direct species-specific research and monitoring efforts, implement UAS data collection, and better understand seasonal distribution of species.

Fieldwork was conducted at PBNERR during the months of April to September. Fieldwork included UAS flights, ground truth vegetation surveys, and elevation measurements. UAS imagery was then processed to perform object-based image analysis (OBIA) with ground truth vegetation surveys. Elevation measurements were used to generate a DEM, which then incorporated elevation into the OBIA with UAS imagery (Figure 1).

Fieldwork was conducted in the northern region of Padilla Bay. This area includes three permanent biomonitoring transects that were established by PBNERR in 2011 (Figure 2), stretching out 4 km perpendicular to the shore from east to west (Stevens et al., 2016). Transects were broken into six distinct zones. Zone 1, closest to the shore, is primarily barren mudflat. Zone 2 transitions from bare mud to Z. japonica dominance. Zone 3 is a mixed area with both Z. japonica and Z. marina species present. Zones 4 and 5 are still middle intertidal but are primarily dominated by Z. marina. Zone 6 is subtidal and also dominated by Z. marina. The three transects were designated as northern, central, and southern, denoted 1, 2, and 3, respectively (Figure 2). The spacing between the transects was approximately 100 m. Each zone on each transect included seven permanent biomonitoring plots placed 50–100 m apart, and plot spacing varied due to the variable length of each zone. Biomonitoring plots were marked with small pins that showed the outer boundary of a 0.5 × 2 m plot where annual surveys were conducted. Polyvinyl chloride (PVC) pipes (2.5 cm diameter and approximately 0.3 m tall) were located 1 m south of the biomonitoring plot sites to mark their location.

All research was conducted during periods of low tides, ranging from −0.17 to −1.24 m MLLW, between the months of April and September. This period was chosen to take advantage of daytime low tides and to capture eelgrass transition to peak biomass, which typically occurs during the months of June and July (Bulthuis, 1996).

The study area for this project included the first 2.5 km of biomonitoring transects, zones 1–4 and half of zone 5 (Figure 2), accounting for an area of approximately 0.75 km². The remainder of zone 5 and zone 6 rapidly declines in elevation from −0.35 m MLLW at the western edge of the study area to approximately −3 m MLLW in zone 6. PBNERR estimates elevation in some areas at the end of zone 5 to be around −0.9 m MLLW (Bulthuis et al., 2021). Due to the time limitations of tide windows, rapid elevation loss, and an increased likelihood of significant water inundation despite low tides, zones 5 and 6 were not included in this project.

UAS flights were conducted using a DJI Matrice 210 quadcopter carrying a Micasense RedEdge MX 10 band camera (Table 1) and were in full compliance with relevant Federal Aviation Administration Section 107 regulations. Optimal UAS imagery was obtained when wind speeds were less than 5 m s⁻¹, sun angles were between 30° and 52°, tide stage was between 0.48 and 0.76 m (−1.6 and −2.5 ft), and there was less than 10% cloud cover (Hein, 2022; Nahirnick et al., 2019). Windows of time providing these conditions occurred monthly through the spring and summer; however, high winds and rain did preclude flights during some of these times.

Predictions of sun angle and tide stage were critical for flight planning. Tide predictions for the La Conner, Swinomish Slough (NOAA, 2021) gauge were used to schedule flights. This tide gauge is 17 km SW of the study area, so 50 minutes were added to the La Conner predictions to determine the time of each day’s low tide at the study area (based on PBNERR estimates). Sun angle ranges for each day were assessed using NOAA’s (2022) sun angle calculator.

Ground-control points allowed for appropriate georeferencing of flight imagery (Duffy et al., 2018; Hein, 2022; Nahirnick et al., 2019). Prior to each flight, 33 highly visible ground-control point markers were placed on a subset of the biomonitoring plot PVC pipes. PVC pipes had recorded GNSS locations and were the same for each flight.

Each flight was programmed using the DJI Pilot application. Prior to each flight, a picture of the Micasense Calibrated Reflectance Panel was taken. The camera was programmed to take pictures every 1 s during the flights. All flights were conducted at an altitude of 120 m above ground level, resulting in a ground resolution of about 6 cm. Total time for each survey flight was approximately 60 minutes, including time for one to two battery changes in the middle of the flight. Volunteers who placed ground-control markers were in communication with the UAS pilot. These volunteers also maintained visual contact with the UAS and scanned for other aircraft during each flight. Two flights were captured each month, totaling 12 UAS flights during the field season, and the best image from each month was retained for further image processing.

Ground-based vegetation surveys are crucial for analysis (Nahirnick et al., 2019; Neckles et al., 2012). Data from 42 supplemental plots were collected using a stratified random sample across flight zone extents. According to the PBNERR sampling protocol, percent cover was measured within a 2 m × 0.5 m plot using point intercept at 25 equally spaced intersections within a 0.5 m × 0.5 m grid (Table 2; Bulthuis et al., 2021). Ground truth surveys of supplemental plots commenced in April, starting 2 m south of the PVC pipe marking a PBNERR biomonitoring location. Every month, the distance from the PVC pipe increased by 1 m southward to avoid remeasuring previously impacted areas (Hein, 2022). In areas dominated by Z. japonica and Z. marina (e.g., zone 2 and zone 5), the sample size was increased due to mixed cover being most predominant in the study area (e.g., zone 3; Hein, 2022). All permanent biomonitoring locations in the study area were surveyed by PBNERR on 12 July 2022. Consequently, supplemental plots were not surveyed in July, and data from the PBNERR effort were used instead.

Cluster analysis of all vegetation survey data was conducted in R (Version 4.2.2, Vegan and Rmisc packages). Standardized Euclidean distance was utilized to identify clusters of various eelgrass covers (Gotelli and Ellison, 2012), including bare and low vegetation, mixed, Z. marina dominant, and Z. japonica dominant.

Prior to each day of elevation sampling, a permanent control mark, established by PBNERR, located north of transect 1 in zone 1, was measured to ensure accuracy of NTRIP corrections and to assess measurement consistency among days. Samples were collected using a nested sampling design that first used systematic measurements every 50 m in and between transects 1, 2, and 3. Then, random nested measurements were taken at 10% of the systematic measurement locations (Fromin et al., 2013; Webster and Oliver, 2007). Additionally, opportunistic measurements were taken of significant features (e.g., visible channels retaining water during low tides or mounds above the water line) in the study area (Kinzel, 2008).

The lowest tides of each month were prioritized for elevation sampling to increase the efficiency of using the Emlid Reach RS2 probes in the field and visibility of features noted for opportunistic measurements. Since elevation data collection in the study area was time and labor intensive and required the lowest tides, sample collection spanned most of the fieldwork season from May to August.

Images from each flight were aligned and georeferenced to create an orthomosaic using AgiSoft MetaShape Professional Version 1.8.4 (64 bit). For each flight, a photo of the Micasense Reflectance Panel as well as values from the downwelling sensor were used to convert the raw images into a calibrated reflectance index (Micasense, 2022). The workflow, as adapted from the U.S. Geological Survey (2017) and Hein (2022), generated a 10 band orthomosaic with a resolution of 6 cm (Bergner, 2023). The pixel size was then resampled to 10 cm to more cleanly translate to the 2 m × 0.5 m vegetation plots (Hein, 2022).

This NDVI mosaic was then used to reduce multispectral variance by masking out unvegetated areas in the 10 band reflectance full season mosaic. Sun glare issues, which increase multispectral variance, were mostly eliminated by using Hein’s (2022) recommendations to fly only at lower sun angles. Nevertheless, wind can generate a dappled water surface that generates scattered flecks of sun glare. To address these issues, pixels that exceeded the 95th percentile of reflectance in one or more bands were omitted from subsequent analyses (Hein, 2022). Principal components analysis (PCA) was then conducted in ENVI using wavelengths for reflectance from the 10 band reflectance mosaic as variables to reduce dimensionality (Bergner, 2023; Hein, 2022). The first three principal components (PCs) were retained in a PCA mosaic and used in image classification.

Image segmentation was performed in eCognition Developer 10.2 by combining the PCA mosaic and NDVI mosaic as project layers. The NDVI mosaic was only used for masking pixels with NDVI values indicating a low presence of vegetation based on a calculated threshold (see “Image Processing” subsection in “Results” section; Bergner, 2023; Hein, 2022). Then, the PCA mosaic generated from multispectral imagery was segmented based on its attributes (Bergner, 2023; Hein, 2022).

Image classification was performed in R. For each image segment, the mean and standard deviation of PC 1, 2, and 3 were used to predict eelgrass cover categories. Ground truth vegetation data were categorized using the results from the cluster analysis, creating training samples. Ground truth surveys including greater than or equal to 80% total eelgrass cover were retained for OBIA (see “Ground Truth Vegetation Surveys” subsection in “Methods” section; Bergner, 2023; Hein, 2022). Random forest (RF) classification, with 2000 trees, was used for analysis since it can provide more refined classification results and rule-based parameters that can be adjusted to better fit the data set (Chutia et al., 2020). Also, an RF classifier was used for analysis since it is resistant to overfitting (Pal, 2005), and it also works well for images with a high amount of noise (Rodriguez-Galiano et al., 2012). The model’s accuracy was evaluated using repeated K-fold cross-validation, which repeatedly and randomly subsets data into folds for training and validation to yield a robust estimation of the model’s performance (Krstajic et al., 2014; Zhong, He, and Chalise, 2020). Classification accuracy can be quantified using several metrics (Story and Congalton, 1986). The overall classification accuracy is simply the percentage of training samples that were correctly classified across all cover types. Producer’s accuracy was calculated to measure the proportion of training samples that were correctly classified. Conversely, user’s accuracy was calculated to measure the model’s accuracy in predicting eelgrass cover based on the training samples, which assesses the reliability of the classification results for each cover class. As such, producer’s accuracy, user’s accuracy, and overall accuracy were used to assess OBIA.

Elevation measurements were used for regression kriging analysis in R. Regression kriging accounted for the strong relationship between distance from the shore and elevation (Bergner, 2023; Bunn, personal communication, 2022). The DEM, exported as a TIFF with a North American Datum (NAD)1983 Universal Transverse Mercator (UTM) Zone 10 North projection, was used for further analysis.

Previous work has indicated that the elevation in Padilla Bay is declining by 0.15 ± 0.15 cm/y due to sediment loss (Kairis and Rybczyk, 2010; Poppe and Rybczyk, 2022). Based on this minimal annual change and time constraints, it was concluded there was negligible elevation change during the fieldwork season. As such, a single DEM was produced using elevation data collected between May and September.

After conducting OBIA using multispectral data alone, a second round of OBIA was conducted using both multispectral data and elevation data. The NDVI mosaic was used to mask pixels that represented low vegetation (see “Image Processing” subsection in “Results” section). Then, the PCA mosaic and DEM were segmented together. All data from the PCA and DEM mosaic were retained aside from pixels masked out with the NDVI mosaic. OBIA used vegetation survey data, the mean and standard deviation of elevation data, and the mean and standard deviation of PCs 1, 2, and 3 from the multispectral imagery to predict eelgrass cover. RF classification was used for OBIA, and repeated K-fold cross-validation was used to assess accuracy.

Cluster analysis of ground truth vegetation survey data identified three cover classes used in further image classification: Z. japonica dominant, mixed, and Z. marina dominant. Survey locations with less than 80% total eelgrass cover were removed from analysis due to spectral variability. OBIA with PCA and NDVI mosaics derived from multispectral imagery was sensitive to spectral variation of environmental conditions during image capture. Elevation measurements taken in the field were used to create a DEM with a vertical accuracy of 4.3 cm. The DEM was used with the mosaics derived from multispectral imagery for OBIA. The overall accuracy rose when elevation and multispectral imagery were combined for OBIA, and the result was less sensitive to spectral variation between images.

Four clusters were identified, representing bare, Z. japonica, mixed, and Z. marina cover (Figure 3). Eelgrass clusters, Z. japonica, Z. marina, and mixed were centered between 75% and 80% total eelgrass cover (Figure 3; Bergner, 2023). Bare and low vegetation clusters were centered at approximately 4.5% total eelgrass cover (Bergner, 2023). Based on cluster analysis and spectral variation from exposed mud (Hein, 2022), plots with greater than or equal to 80% total eelgrass cover were retained for further analyses. As such, 131 of the 303 total vegetation surveys that contained greater than or equal to 80% total eelgrass cover were utilized as training data. Categories for Z. japonica, Z. marina, and mixed eelgrass cover were determined from the centroids for each cluster and used for OBIA (Table 3; Bergner, 2023). Since a large portion of vegetated survey locations contained less than 80% total eelgrass cover, a less than 80% vegetation eelgrass category identified regions that ranged from partial eelgrass coverage to mudflat (Table 3). Survey data containing less than 80% total eelgrass cover were not used for classification of eelgrass cover categories. Instead, plots with less than 80% total eelgrass cover were used to determine an NDVI threshold to be excluded from OBIA. Portions of the NDVI mosaic below the NDVI threshold were classified as a less than 80% vegetation class.

Orthomosaics were generated from over 20,000 images per flight. One flight per month was retained for image processing. Imagery had minimal streaking and was captured between 10:45 am and 4:45 pm (Figure 4; Table 4). Acquisition of all imagery occurred at a tide stage around −0.76 m MLLW, with a range of −0.44 to −1 MLLW (Table 4).

After conducting PCA, PCs 1, 2, and 3 were retained for OBIA (Table 5). PC 1 explained 82.1% of the variance in the imagery and had all positive loadings, with the strongest ones on the red edge₇₄₀ (9) and the near-infrared₈₄₂ (10) wavelengths. PC 2 explained 15.4% of the variance in the imagery and contrasted the red₆₅₀ and red₆₆₈ (5 and 6) wavelengths with the near-infrared₈₄₂ (10). PC 3 explained only 1.2% of the variance and contrasted the red edge₇₀₅ and red₇₁₇ (7 and 8) wavelengths from the near-infrared₈₄₂ (10) wavelength. The near-infrared₈₄₂ (10) wavelength had a strong contribution across all PCs and captured subtle spectral variations between Z. japonica and Z. marina. The red edge₇₀₅, red₇₁₅, and red₇₄₀ (7, 8, and 9) and red₆₅₀ and red₆₆₈ (5 and 6) wavelengths, both contrasting and contributing to near-infrared₈₄₂ (10), reflect complex factors like chlorophyll content, plant health, and structural differences between Z. japonica and Z. marina.

The NDVI mosaic, generated from multispectral imagery, had an overall accuracy of 77% when assessing total eelgrass cover (Table 6). The NDVI mosaic yielded a user’s accuracy of 70% identifying areas with greater than or equal to 80% total eelgrass cover and 82% identifying areas with less than 80% eelgrass cover (Table 6). Since the NDVI mosaic was used to mask areas with less than 80% total eelgrass cover, the NDVI’s accuracy for identifying those areas was the most important factor. To further assess the NDVI, values were compared with ground truth data. The average NDVI for all the plots that had less than 80% vegetated cover was 0.25 (Figure 5). However, the NDVI could vary from environmental factors like cloud cover and phenological differences like plant health. For instance, some ground truth plots with greater than 80% vegetation cover had NDVI values less than 0.25, and some plots with NDVI greater than 0.25 had less than 80% vegetation cover. Using the 0.25 threshold, 36% of the ground truth plots with less 80% vegetation cover were retained, while 64% of plots were removed. Conversely, only 1% of ground truth plots with greater 80% vegetation cover were removed from analysis using this threshold. The NDVI threshold of 0.25 was retained for segmentation.

Using multispectral imagery alone, the overall RF classification accuracy was 75% (Table 7). The producer’s and user’s accuracies for Z. japonica–dominant cover were 56% and 75%, respectively (Table 7). Mixed cover and Z. marina–dominant cover had similar accuracies (Table 7). The producer’s and user’s accuracies of mixed cover were 80.5%, and 72.3%, respectively. Producer’s and user’s accuracies for the Z. marina–dominant class were 75.7% and 78.4%, respectively (Table 7).

April had the lowest vegetation coverage. Since areas with low vegetation and exposed mud increase multispectral variance (Hein, 2022; Svane et al., 2022), the largest number of misclassifications happened in April. For example, there was a clear overprediction of Z. japonica–dominant cover in April. As a result, mixed and Z. marina–dominant cover were underpredicted in April (Figure 6).

May had the second lowest vegetation coverage of the field season and thus the second highest number of misclassifications of eelgrass cover. Also, there was a dark streak spanning across zones 1 and 2 that was classified as Z. japonica–dominant cover but that contained little to no eelgrass. The streaking probably resulted from a passing cloud during image capture (Table 4). Z. marina–dominant cover was underpredicted in May but was not in an area with obvious streaking, whereas mixed cover was accurately predicted.

June imagery was captured at a different tide stage than the other flights (Table 4). Turbulence, foam, and debris carried by the incoming tide in the southwest corner of the study area (transect 3, zone 5) necessitated removing this portion of the imagery from subsequent segmentation and classification. This area was manually classified as “no data,” but it was likely dominated by Z. marina cover, as predicted for other months throughout OBIA. By June, the spatial distribution of predicted Z. japonica–dominant cover was more accurate. Z. marina–dominant cover was overpredicted in zone 3, which also resulted in an underprediction of mixed cover.

July and August had the highest image quality, with only moderate streaking. Predictions for all cover classes were the most accurate during these months.

In September, a large algae bloom of Cladophora columbiana occurred before the start of UAS flights and was noted during ground truth data collection. This bloom also caused subsequent predictions to indicate Z. japonica–dominant cover in areas where algae accumulated, primarily the upper portion of zone 1. For other months, the upper portion of zone 1 also included eelgrass predictions caused by algae and detritus deposits. The predictions of eelgrass in this region were manually changed to an algae/detritus class for all months, with the largest area occurring in September. All eelgrass cover categories were predicted accurately for September once the algae/detritus class was manually created.

Since there were no significant differences in elevation measurements of the control mark (p = 0.265; Bergner 2023), all sample days were used for DEM generation. Elevation measurements had associated root mean squared errors (RMSE). Various RMSE thresholds were used to generate DEMs. The result with the highest accuracy and best visual representation of the study area was used. The DEM that was generated using only measurements with an RMSE less than 1.8 cm, using 879 of the 961 recorded elevation measurements, was ultimately selected (Bergner, 2023). Regression kriging was used to generate the DEM to account for the relationship between elevation and distance from the shore (Meng, Liu, and Borders, 2013). The experimental variogram from regression kriging looked appropriate for further analysis (Figure 7). Data were best fit with an exponential variogram (Figure 7). The model was validated using K-fold (10 folds) cross-validation, resulting in a mean error of −0.096 cm, correlation of 0.99, and RMSE of 4.3 cm (Table 8). The elevation predictions ranged from −0.35 to 1.53 m MLLW, with an average elevation in the study area of 0.43 m MLLW (Figure 8). The DEM showed distinct channels and mounds throughout the study area that corresponded to observations in the field (Figure 8).

Overall classification accuracy using elevation data and multispectral imagery was 88.9% (Table 9). The producer’s and user’s accuracies for Z. japonica–dominant cover were 67.7% and 91.7%, respectively (Table 9). Z. japonica–dominant cover was slightly overpredicted in May (Figure 9). The overprediction in May was most notable in the streak spanning across transect 2 to transect 3 in zones 1 and 2 that was caused by variable cloud cover (Table 4; Figure 4). Some of the regions of this streak primarily had less than 80% vegetated cover and were misclassified as Z. japonica–dominant (Figure 9). Also, there were mounds located in zone 5 of transect 3 that had sparse Z. japonica cover. At one of these mounds, Z. japonica–dominant cover was predicted in August (Figure 9). Although Z. japonica was present in this area, the ground truth vegetation surveys noted that cover was not enough to be considered dominant. The misclassification could have been caused by algae present in the area. With the incorporation of elevation data, mixed cover was predicted more uniformly across zone 3 compared to results using multispectral imagery alone (cf. Figures 6 and 9). Elevation in zone 3 was the most uniform compared to the other zones in the study area. As a result, mixed cover primarily contained mixed species stands instead of mounds of Z. japonica and depressions of Z. marina. Additionally, most ground truth vegetation survey plots in zone 3 had both species interspersed with one another. The producer’s and user’s accuracies for mixed cover were 92.2% and 85.5%, respectively (Table 9). Z. marina–dominant producer’s and user’s accuracies were both 92% (Table 9). Spatially, Z. marina–dominant cover stayed consistent across all the months, and it was the highest performing class using this classification technique (Table 9). Incorporating elevation data also showed the location of deep channels in the study area where Z. marina was present. For example, there was Z. marina–dominant cover in a channel located in zone 3 of transect 1, which was predicted for every month (Figure 9).

Based on these results, Z. japonica–dominant cover was lowest in April and increased slightly in May. Between June and July, Z. japonica increased by approximately 0.01 km² (Figure 10). Mixed cover increased the most dramatically during the field season by an area of 0.07 km² between the months of April and May and stayed stable through the rest of the months (Figure 10). The area of Z. marina–dominant cover did not change significantly between April, May, July, and August (Figure 10). June likely had a similar area of Z. marina–dominant cover as every other month, but the section of no data in transect 3 of zone 5 lowered the total area of this cover class (Figure 10). Results suggest that the algal bloom accounted for 0.04 km² in September (Figure 10).

The multispectral mosaic had the greatest variance between zones 1 and 2, where the most exposed mud was located. The contribution of multispectral variance to misclassifications was lowered significantly with the addition of the elevation mosaic. The greatest variance of the elevation mosaic occurred at the northernmost and southernmost edges of the study area and at the end of zone 5 and did not contribute significantly to misclassifications. As such, pairing multispectral imagery with elevation for OBIA yielded a high accuracy of different eelgrass cover types. Compared to past research using UAS, the different eelgrass cover detection in this study was performed to similar accuracies of the past research that was confined to detecting presence and absence of eelgrass. Table 10 summarizes different research that used imagery for eelgrass mapping. The results of this study highlight the advances made in eelgrass detection through advancements in technology and incorporation of new data.

UAS image capture demands focused planning, considering factors like sun angle, wind speed, cloud cover, and tide stage. Optimal conditions in the Pacific Northwest are rare, complicating efforts to capture consistent imagery across different times of the year. Vegetation surveys in areas without high eelgrass cover (greater than 80%) present classification challenges. Spectral variation between species allowed for differentiation between cover types due to seasonal changes, canopy structure, and different chlorophyll content. Integration of elevation data with multispectral imagery provided a comprehensive understanding of species distribution and environmental interactions in coastal ecosystems on a seasonal scale that can be applied for more frequent and directed monitoring.

UAS flight scheduling requires careful consideration of sun angle, wind speed, cloud cover, and tide stage. The perfect combination of all factors will rarely, if ever occur, especially if the goal is to obtain comparable imagery at different times during the year. Although flight planning greatly benefited from insights provided by previous studies (Hein, 2022; Nahirnick et al., 2019), unavoidable variations in conditions among months resulted in image anomalies that created challenges for image classification. When the sun angle is high, specular reflection from the water surface results in glare that can obscure large portions of the imagery (Hein, 2022; Nahirnick et al., 2019). During months of peak eelgrass biomass, reducing specular reflection by conducting flights at lower sun angles conflicts with flying during lower tide stages, which often occur when the sun angles are higher. Within a single flight, the sun angle can change considerably. This was particularly apparent in the July and August imagery (Table 4), resulting in light E-to-W streaking in the portion of the flight that occurred during higher sun angles (Figure 4). The cause of streaking is not entirely clear, but it may have been caused by variation in the camera viewing angle on alternate flight lines. The UAS flight plan included eight flight lines, with the first and last both outbound (to the west). In the July and August images, when the sun was in the southeast, the outbound flight lines were noticeably brighter, and the inbound flight lines were darker. This may have resulted from slight variation in the camera viewing angle. If so, this streaking might be eliminated when using a gimbal-mounted camera that can maintain a constant camera viewing angle.

High wind can result in a dappled water surface that can create scattered glare over much of the study area regardless of sun angle. Each of the flights occurring with wind speeds higher than 5 m s⁻¹ were excluded from analysis. The best imagery was collected when cloud cover was less than 10% percent (Nahirnick et al., 2019), and imagery collected during cloud cover greater than 10% resulted in shadowed areas of the image that could not be corrected. An exception was the September orthomosaic, which had 20% cloud cover, but the low sun angle and relative cloud position limited shadowed areas in the imagery (Table 4). Overcast conditions result in diffuse rather than direct illumination (Köppl et al., 2021). Diffuse illumination causes light to be coming in at all angles, including some from directly overhead. For each flight, a preflight calibration was done using the Micasense Calibrated Reflectance Panel, and imagery was converted into reflectance to correct for variation in illumination between dates. Yet, the June image, collected under 100% overcast conditions, had higher reflectance on average and was noticeably brighter than all the other images that were collected with less than 10% cloud cover (Figure 4). This was the direct result of the illumination under overcast conditions coming in from directly overhead, even with a low sun angle, and generating specular reflectance and uniform glare over the entire image. Therefore, if the purpose of research is to compare imagery from different dates, then all imagery should be collected with cloud cover less than 10%.

Ideally, all flights would occur at the same tide stage. This would ensure that each day of image capture would have eelgrass in a similar position in the water column, and similar regions would be drained or have water present. Capturing imagery at varying tide stages could alter cover predictions because eelgrass blades may be horizontally oriented across the bottom during a very low tide, while they would be vertical in the water column at the same location during a higher tide stage. Exploratory UAS flights of a new study area at different tide stages can provide insight on the optimum conditions for detecting eelgrass.

Vegetation survey data were gathered from the months of April to September, including areas that did not have a large amount of eelgrass present. Exposed mud introduced noise into the multispectral imagery that swamped the more subtle signal generated by the variation in vegetation cover and species composition. For this reason, image classification was restricted to areas with greater than 80% total eelgrass cover (Hein, 2022). Future efforts should focus on obtaining data from sites with high eelgrass cover and reduce surveys at sites with little or no cover. Furthermore, there were only 22 Z. japonica–dominant sites in the classification data set, mainly from the months of June–September. This resulted in Z. japonica–dominant cover being underpredicted in April and May and having the lowest accuracies in OBIA (Table 9). For future analysis, increasing samples for Z. japonica–dominated areas (zone 2 in the study area) to better match the sample size of other cover classes would improve the model’s prediction accuracy for this cover type (Lin et al., 2019; Ma et al., 2017). Given the larger training sample size for mixed and Z. marina–dominant sites, it is not surprising that these two classes had the highest classification accuracy in OBIA (Table 9).

During vegetation surveys, percent cover was the only data collected due to time constraints in the field. Surveys on biomass were not conducted, and results do not represent thickness of the mat. Instead, results provide information on species abundance at the surface of the eelgrass mat. Although Z. marina had relatively stable cover through all the months of fieldwork, there were changes in Z. marina biomass through these months that were not captured by percent cover data. In contrast, Z. japonica underwent considerable seasonal changes in both percent cover and biomass throughout the growing season.

The PCA highlighted the variation in the vegetated areas in the multispectral imagery (Jolliffe and Cadima, 2016). PC 1 loadings across all bands did not contrast with one another, and this suggests that PC 1 represents a weighted index to overall scene brightness, with strongest loadings on the red edge and near-infrared bands. The red edge₇₄₀ (9) and near-infrared₈₄₂ (10) bands are indicators of variation in vegetation cover. Overall, PC 1 loadings suggested spectral differences between Z. japonica and Z. marina based on seasonal distribution of each species. Z. japonica cover was sparse in April and reached peak biomass in June, while Z. marina cover remained more consistent throughout all months. PC 2 contrasts aligned with the NDVI, suggesting that this PC was a vegetation index. A staple in remote-sensing analysis, the NDVI depicts vegetation health and photosynthetic rate by assessing the difference between light absorbed (red) and light reflected (near-infrared) (Huang, Tang, and Hupy, 2020). Z. marina is a larger species than Z. japonica, and the large blades of Z. marina cover a greater area compared to those of Z. japonica. Furthermore, Z. japonica–dominated areas were closer to or mixed with exposed mud. Increased vegetated cover will raise the NDVI. As a result, zones 4 and 5, predominantly Z. marina–dominated areas, had higher NDVIs in the study area compared to Z. japonica–dominated areas. Mixed cover had varied NDVIs for April and May, with smaller values than Z. marina–dominated areas. By June, mixed areas had a thick mat of eelgrass that included both species and had the highest NDVIs in the study area. As such, PC 2 is a signifier for differences between canopy structure and coverage of species. PC 3 contrasts suggested a normalized difference red edge index (NDRE). The NDRE shows the effect of stress as well as chlorophyll concentration in vegetation (Barnes et al., 2000). Although PC 3 only accounted for 1.2% of the total variance (Table 5), it was used in analysis because it increased the predictive power in differentiating the species by 15% in OBIA with multispectral imagery (Bergner, 2023). The increased predictive power of PC 3 suggests that there may be different chlorophyll concentrations between Z. marina and Z. japonica, or that one species undergoes different or more severe stress from factors like heat and desiccation. Collectively, the loadings for the PCs captured different features of the spectral signature of the eelgrass cover. The unique signals contributed to success in discriminating between the two species of eelgrass (Table 9).

Although the goal was to focus on sites with greater than 80% vegetation cover, the selected NDVI threshold was only partially successful in isolating this subset of the study area (Figure 5; Table 6). Future studies can explore alternative approaches to isolate sparsely vegetated areas prior to classification of different eelgrass species.

Using multispectral imagery alone, OBIA inaccuracies in April and May were the result of more sparsely vegetated regions (less than 80% vegetation cover) that were not masked out with the NDVI threshold of 0.25 (Figure 5). In these areas, the spectral signature for vegetation was masked by the spectral variation of exposed mud (Hein, 2022). To address this issue, different segmentation parameters could generate smaller “pure” segments that include only eelgrass or mud, and a higher NDVI threshold could be used. Eelgrass cover in Padilla Bay typically reaches its peak between June and July (Thom, 1990). For the months of April and May, approximately 70% of vegetation surveys contained less than 80% total eelgrass cover. Yet, many of these survey sites were being classified as an eelgrass cover category in OBIA primarily in April and May. Most misclassifications occurred for the Z. japonica–dominant and mixed cover classes in zones 2 and 3. Misclassifications in April included the Z. marina–dominant cover class as well. By June, when more Z. japonica had grown in the study area, misclassifications were reduced. Results from OBIA with multispectral imagery were the most accurate for the months of July–September because of the high density of vegetation throughout the study area (Figure 6). Results using multispectral imagery alone suggested that there were spectral differences between eelgrass species that resulted in different cover class classification in OBIA, but the variation of other spectral factors in the study area, like exposed mud, did not produce results robust enough for detailed, species-specific distribution analysis.

Further sampling at Padilla Bay, especially in the defined channels and mounds throughout the study area, would further improve the vertical accuracy of the DEM. Previous studies have shown that the distribution of these two species is strongly influenced by elevation (Bergner, 2023; Bulthuis et al., 2021; Ruesink et al., 2010). In some cases, slope can also influence the distribution of these species, but slope for the study area only ranged from 0 to 0.1%.With such minimal variation in slope, it is not surprising that it did not contribute to improved prediction of eelgrass species cover (Bergner, 2023). However, in different study areas, slope may be a potential predictor to complement elevation.

Incorporation of elevation data into classification analysis requires careful consideration of a study area’s characteristics. For example, study areas near river confluences or deltas may experience seasonal fluctuations in elevation due to variation in sediment inputs and scouring. Under these conditions, performing elevation data collection over a several month period and then merging these data to generate a single DEM would be inappropriate.

Combining both elevation data and multispectral imagery for OBIA yielded the highest accuracies for each cover class, and elevation was deemed to be a strong predictor for eelgrass cover. Using the DEM alone for OBIA in preliminary analysis yielded a static layer that could not assess seasonal change, and the model performed poorly at classifying Z. japonica–dominant cover, yielding a user’s accuracy of approximately 60% (Bergner, 2023). Only through the combination of elevation and multispectral imagery was Z. japonica–dominant cover predicted to an accuracy of 90% (Table 9). As such, multispectral imagery provides insight into important variability between species that would otherwise be lost, and elevation data are critical for removing spectral variation.

Further inspection of the spatial and temporal patterns generated by this model provides another way to evaluate the results. Z. japonica is known to decrease dramatically in biomass during the winter months and grow to peak biomass between June and July. Based on ground truth vegetation surveys, Z. japonica–dominant cover was not present in zone 3 in April. Also, Z. japonica–dominant cover was not as widespread in zone 2 in May compared to July, when this species reaches peak biomass in Padilla Bay (Bulthuis, 1996). April misclassifications may be the result of algae that has a higher NDVI than the 0.25 threshold for segmentation, while the May misclassification was the result of cloud cover in zone 2 and the lower portion of zone 1. Although Z. marina does experience a decrease in biomass during the winter, Z. japonica experiences a large die-off in winter (Bulthuis, 1996). As a result, Z. marina varied less in percent cover across all months compared to Z. japonica. The changes in cover of Z. japonica over the field season caused higher errors in OBIA compared to Z. marina and mixed cover categories. For example, the sample size for Z. japonica was lower, and the amount of variance caused by exposed mud was captured most frequently in Z. japonica–dominated areas.

Each cover category had spatial distribution that could be used to interpret seasonal variation when the multispectral imagery and elevation data were combined. Z. japonica–dominant cover grew from the south (transect 3) to the north (transect 1). Transect 3 was deeper than transect 1 (Figure 8), and these growth patterns suggest that transect 1 may be on the shallower end of the preferred elevation range of Z. japonica. Zone 3 contained a more uniform, mid-tide elevation and is an area where the species may undergo less desiccation stress (Ruesink et al., 2010). Z. japonica–dominant cover occupied the largest area in July and August (Figures 9 and 10) and began encroaching into areas previously dominated by mixed cover. Encroachment of Z. japonica–dominant cover into previously mixed areas suggests Z. japonica may have been able to outcompete Z. marina in this part of the study area during the longer days, higher temperatures, and daytime low tides associated with July and August. Although subtle, mixed cover expanded into previously Z. marina–dominated areas between July and August, which could also suggest Z. japonica was able to expand into deeper areas due to favorable conditions. Similarly, Z. marina–dominant cover changed slightly from month to month in the transition zone of mixed cover, with the largest change occurring between August and September (Figure 10). Converse to Z. japonica, the increase in Z. marina–dominant cover between August and September may have been influenced by cooler temperatures, tide changes, and shorter periods of sunlight. Patches of Z. marina–dominant cover in zone 3 predicted for each month corresponded to the location of the two most significant channels in the study area and were consistent in each month of imagery (Figure 9). Overall, Z. japonica changed dramatically between April and September and encroached on areas with Z. marina present. The results suggest that Z. japonica responded rapidly to midsummer conditions in Padilla Bay both in dominant coverage areas and in mixed regions.

It is likely there were small portions of the study area misclassified as eelgrass that were in fact algae. Unfortunately, most algae in the study area float in mats and can vary in location on a day-by-day basis. Although algae introduced inaccuracy in OBIA, algae mats were quite rare: They were present in only 58% of ground truth survey plots, with 71% of those plots having less than 10% algae present. For this reason, algae did not occur as a distinct group during cluster analysis. Furthermore, algae that was present in a plot on the day of a ground survey was probably gone on the next flight day. The shifting location of these mats made it problematic to collect useful ground truth data for this cover type. Future research can explore classification of algae using low-altitude UAS flights to create virtual ground truth data and include an algae cover category.

OBIA with multispectral imagery and elevation provides insightful results for tier 1 monitoring. Overall, tier 1 monitoring is critical for understanding eelgrass distribution. Now, differentiation of eelgrass cover types using a landscape-wide approach can be performed with high accuracy. These results will allow researchers to establish a cover-specific baseline for Z. marina–dominant, mixed, and Z. japonica–dominant cover on annual and seasonal scales. Additionally, the spectral differences and distribution of species cover suggest that Z. japonica is encroaching on Z. marina in the middle of the summer, when Z. marina may be undergoing stress. Areas where encroachment was observed can be targeted for more research to understand competition and physiological responses to changes in temperature, sunlight, and tidal fluctuations to predict how eelgrass in Padilla Bay will respond to climate change stressors. Eelgrass stands can vary widely between each year, and in Padilla Bay, species can vary in location seasonally and annually based on light availability, temperature, and water inundation during low tides (PBNERR, unpublished data). A comparison of seasonal and interannual differences in the elevation ranges of each cover type could provide a better understanding of environmental factors impacting the location of eelgrass species and increase the prediction accuracy for areas with eelgrass species presence in future years. As such, tier 1 monitoring using these methods can guide a species-specific approach to identifying critical habitats, targeting areas for tier 2 and 3 monitoring, and informing broader conservation efforts. For example, past research differentiating native and nonnative eelgrass beds using aerial imagery has not been able to detect spectral differences between species and instead used elevation and species zonation in Yaquina estuary (Young et al., 2015). Additionally, Z. japonica was not detected in the mid-tidal portions of the estuary, even though it was identified at those elevations in ground surveys (Young et al., 2008). Further, research by Young et al. (2008, 2015) implies that mixed stands of eelgrass cannot be detected. Alternatively, UAS and remote-sensing technology used in this project allowed for detection of eelgrass in mixed, native, and nonnative stands across varying elevations over the course of 6 months, yielding a higher accuracy and more comprehensive view of species distribution and areas occupied by both species.

Successful detection of Z. japonica–dominant, mixed, and Z. marina–dominant cover categories from UAS imagery can be generated with OBIA using multispectral imagery and a DEM of the study area that meets centimeter-level vertical precision. This is the first study to distinguish between Z. japonica and Z. marina cover with a high accuracy, reaching 89%. Obtaining useful UAS imagery is critically dependent upon careful selection of flight conditions (Hein, 2022; Nahirnick et al., 2019). Uniformity of conditions like cloud cover, wind speed, and tide stage between each day of image acquisition was imperative. Although satellite imagery can be beneficial for mapping distribution of eelgrass over wide areas, the type of detail captured in this research was a product of using the UAS, since spectral variability was captured at finer resolutions, and data were gathered monthly to assess seasonal dynamics.

Although stratified random sampling design for percent cover estimates provided a start to planning vegetation surveys, plots with less than 80% total eelgrass cover were not used for image classification because of increased noise from exposed mud. As such, it is useful to increase the sampling effort in areas with high vegetation cover. Also, an increase in ground truth survey plots for Z. japonica–dominant sites, to match the number of samples more closely for mixed and Z. marina–dominant plots, is needed to increase classification accuracy. Multispectral imagery alone detected seasonal changes for eelgrass cover types but with some lingering errors. The fusion of elevation data with multispectral imagery yielded the best results.

Eelgrass monitoring using traditional monitoring approaches may provide a detailed and accurate depiction of species presence and distribution for a specific area, but the approach highlighted in this research is a powerful way to understand species-specific cover distribution in a larger area and at more frequent intervals. Traditional monitoring methods may overlook specific changes identified in this research. Further, past methods for remote-sensing analysis also do not allow for the detailed classification of cover types and could entirely miss areas where Z. japonica is present in the mid-tide range or where mixed stands of species have formed (Young et al., 2008, 2015). Overall, this novel approach can be used to detect changes to eelgrass cover quickly and with a species-specific focus. Classifying native and nonnative eelgrass cover from multispectral imagery captured via UAS can provide useful information about the spatial distributions of species as well as seasonal and interannual changes in the cover of eelgrass species. As a landscape-wide (tier 1) approach, overhead image capture and OBIA analysis of eelgrass stands represent a time-efficient way to characterize the distribution and area of eelgrass species cover that will yield more species-specific detail in multitiered monitoring efforts (Neckles et al., 2012). UAS image capture and OBIA of multispecies eelgrass stands can be used worldwide to assess challenges to species-specific monitoring of eelgrass. In the Pacific Northwest, these methods provide timely tactics to further understand the dynamics and potential impacts of Z. japonica.