Comparing PlanetScope to Landsat-8 and Sentinel-2 for Sensing Water Quality in Reservoirs in Agricultural Watersheds

Abstract

Agricultural runoff transports sediments and nutrients that deteriorate water quality erratically, posing a challenge to ground-based monitoring. Satellites provide data at spatial-temporal scales that can be used for water quality monitoring. PlanetScope nanosatellites have spatial (3 m) and temporal (daily) resolutions that may help improve water quality monitoring compared to coarser-resolution satellites. This work compared PlanetScope to Landsat-8 and Sentinel-2 in their ability to detect key water quality parameters. Spectral bands of each satellite were regressed against chlorophyll a, turbidity, and Secchi depth data from 13 reservoirs in Oklahoma over three years (2017–2020). We developed significant regression models for each satellite. Landsat-8 and Sentinel-2 explained more variation in chlorophyll a than PlanetScope, likely because they have more spectral bands. PlanetScope and Sentinel-2 explained relatively similar amounts of variations in turbidity and Secchi Disk data, while Landsat-8 explained less variation in these parameters. Since PlanetScope is a commercial satellite, its application may be limited to cases where the application of coarser-resolution satellites is not feasible. We identified scenarios where PS may be more beneficial than Landsat-8 and Sentinel-2. These include measuring water quality parameters that vary daily, in small ponds and narrow coves of reservoirs, and at reservoir edges.

1. Introduction

Agricultural runoff transports sediments and nutrients to downstream ecosystems where they may cause major water quality impairments. Sediments in the runoff cause elevation in turbidity (Turb) and a cascade of water quality problems such as a reduction in light penetration, low ecological productivity, alteration of food chains, and bioavailability of associated pollutants [1]. Nutrients brought in by the runoff can stimulate phytoplankton growth and promote harmful algal blooms (HABs) in downstream lakes and reservoirs [2]. A major problem with changes in water clarity and phytoplankton growth is that they are both spatially and temporally erratic.

The erratic nature of chlorophyll a (Chl-a) is evident in surface algal blooms, which may occur or disappear within hours or days depending on the underlying processes in a water body [3]. Algal blooms may be more frequent when nutrient loading through sediment transport is persistent. Since agricultural watersheds transport nutrients to lakes and reservoirs, especially during high flows, daily monitoring of sediments and phytoplankton in reservoirs may be the best way to detect pollution at the early stage of impairment.

Successful water quality monitoring programs require the collection of data that capture spatial and temporal trends in many water bodies [4,5]. Achieving this is often challenging because of limited resources. Many water bodies are not frequently monitored, resulting in uncertainties in determining if they support their beneficial uses and the magnitude to which they are impacted by nutrients and sediments [6]. Water quality monitoring programs often utilize rotational monitoring schedules to address this issue. The problem with infrequent monitoring and determining which water bodies to include in a rotational monitoring program is the possibility of missing out on significant pollution events because of their high temporal and spatial variations [7].

The availability of remotely sensed data from satellites provides the opportunity to monitor water bodies at greater temporal and spatial resolutions than field-based monitoring alone [8]. Studies have developed various methods to use spectral bands for water quality sensing. Examples include band ratios and band indices [9,10], spectral derivative techniques [11], and color space transformation [12]. The resulting spectral data are paired with the temporally and spatially coincident water quality data to establish empirical or analytical relationships [13]. Three key water quality parameters that have been studied using satellite remote sensing are Chl-a, Turb, and Secchi Disk depth (SD).

Satellite remote sensing of Chl-a is an active area of research because Chl-a is associated with adverse ecological effects at high concentrations during bloom events [14], and it is optically active with spectral signatures of green (G) reflectance and absorbance in red (R) and blue (B) [15]. The optical properties of Chl-a have been used to develop several empirical models using spectral bands and band combinations. For example, a ratio of R to near-infrared (NIR) is a common band ratio used to estimate Chl-a concentrations [9,16,17,18]. Other examples of band combinations used to monitor Chl-a in water bodies include a ratio of B to G [19] and the normalized difference between NIR and R at specific wavelengths [10]. First and second derivatives have also been used to study Chl-a in waterbodies. Becker et al. (2005) defined the first and second derivatives as shown in Equations 1 and 2, respectively.

d^1st = (ρ_n+1 − ρ_n)/(λ_n+1 − λ_n)(1)

d^2nd = (d_n+1^1st − d_n^1st)/(λ_n+2 − λ_n)(2)

where d^1st and d^2nd are the first and second derivatives, n is the band number, ρ is the surface reflectance, and λ is the wavelength.

Similarly, studies have used both single bands, such as B and R, and band combinations, such as B/R and G/R, to study SD and Turb [20,21]. Polynomial terms of single bands, such as R³, have also been used [22]. Studies have also compared visible bands from different satellite sensors to estimate water transparency [23].

Remotely sensed data from satellites with moderate to coarse spatial resolution imagery such as Landsat-8 (with a spatial resolution of 30 m) and the Sentinel-2 constellation (with a spatial resolution ranging from 10 to 60 m) have shown promising results in monitoring water quality in inland water bodies [24,25]. However, their applications for remotely sensing water quality in small reservoirs, narrow coves of large reservoirs, or for parameters such as Chl-a that can vary daily within a system may be limited. PlanetScope (PS) is a relatively new satellite that provides greater temporal (daily repeat schedule) and spatial (3 m) resolution than these other satellites (Table 1). The temporal resolution of Landsat-8 is 16 days and for Sentinel-2 is 10 days for each of its two polar-orbiting satellites giving an average temporal resolution of 5 days. [26]. Table 1 compares some of the major characteristics of PS, Landsat-8, and Sentinel-2. As such, we posit that where adequate characterization of water quality for small reservoirs using moderate-coarse resolution data (e.g., Landsat-8, Sentinel-2) might not be feasible, PS data with finer spatial (3 m) and temporal (daily) resolutions may be a viable alternative.

The spatial and temporal resolutions of PS for water resource remote sensing are promising as shown in recent studies including bathymetric mapping [27] and change detection of water quality parameters such as total suspended matter [28]. PS imagery produced up to 50% overall accuracy for mapping benthic habitats with multiple species, including coral reefs, macroalgae, seagrass, bare substratum, and dead coral. Up to 74.31% accuracy was recorded for mapping several species of seagrass [29]. PS has also been used together with Landsat-8 to improve on bathymetric maps produced by Landsat-8 [30].

The advantage of PS for water quality sensing in comparison to Landsat-8 and Sentinel-2 may be limited by the number of spectral bands. Landsat-8 and Sentinel-2 have additional spectral bands that may give them an advantage over PS to remotely sense water quality. The shortwave infrared (SWIR) bands have been shown to be less susceptible to Turb than the NIR band and, hence, better detect Chl-a [31]. The Red-Edge band in Sentinel-2 has also been shown to better detect Chl-a than the NIR band [32].

The objective of this study was to determine if the finer spatial and temporal resolutions of PS provided better water quality monitoring potential compared to Landsat-8 and Sentinel-2. To achieve this, we determined empirical relationships between the surface reflectance of each of the three satellites and Chl-a, Turb, and SD data from 13 reservoirs in Oklahoma. We also present a specific case study of an algal bloom in an Oklahoma reservoir to compare the utility of the three satellite platforms in monitoring an algal bloom with respect to the spatial, temporal, and spectral properties of each satellite.

2. Materials and Methods

2.1. Description of the Study Area

The state of Oklahoma has a total surface area of 181,196 km². It has approximately 3629 km² of water in reservoirs and ponds and 89,554 km of shoreline [33]. There are over 200 reservoirs distributed across the state that are greater than 0.2 km² [33] and qualify for rotational water quality monitoring in the state of Oklahoma’s Beneficial Use Monitoring Program (BUMP). Beneficial use refers to defined water uses (primary body contact recreation, public and private water supply, agriculture, aesthetics, and fish and wildlife propagation) that must be supported by water quality parameters according to prescribed standards [34]. Oklahoma also has thousands of reservoirs with surface areas less than 0.2 km² that include ponds ranging between 0.001 and 0.04 km². Most of these ponds were constructed to support livestock watering and recreation [35].

Water quality data were obtained from 13 reservoirs: Lake Arcadia, Broken Bow Lake, Canton Lake, Lake Eucha, Fort Gibson Lake, Foss Reservoir, Lake Hefner, Kaw Lake, Lake McMurtry, Oologah Lake, Lake Thunderbird, Waurika Lake, and a single cove in Grand Lake O’ the Cherokees (Grand Lake) called Horse Creek Cove. Grand Lake is the third largest reservoir in Oklahoma with a surface area of 188.2 km². We included Horse Creek Cove because Grand Lake has experienced HABs during the summer coinciding with the summer recreation season. Most of the HABs in Grand Lake start in Horse Creek Cove from nutrient-rich runoff from agricultural production in its watershed [36]. Figure 1 shows a map of Oklahoma and the state’s BUMP reservoirs (grey) including those studied here (colors). Table 2 shows some characteristics of the study reservoirs.

2.2. Water Quality Data

Water quality data for 12 of the 13 reservoirs were obtained from the BUMP, which samples 62 lakes ranging in size from 0.2 to 2 km², and 68 lakes greater than 2 km². The BUMP was designed to sample 130 reservoirs in Oklahoma quarterly, with a five-year rotational schedule [38]. The data for Horse Creek Cove were obtained from the Grand River Dam Authority (GRDA) water quality monitoring program.

We used Chl-a, Turb, and SD data collected between 2017 and 2020. We selected sample data that corresponded to imagery acquired within ±2 days of each sample date to allow for enough paired satellite and water quality data for each of the three satellite platforms [25]. The assumption was that a satellite image that was acquired within a ±2-day window of ground-based sampling would be representative of the water quality conditions if no major events such as flooding, severe weather, and limnetic mixing occurred [39]. When more than one image was available within the ±2-day window, the image with an acquisition date closest to the date of ground-based sampling was selected.

Chl-a samples were collected below the surface at a depth of 0.5 m and placed on ice and returned to the laboratory for filtration, extraction, and measurement of Chl-a concentrations in µg/L [40]. Surface water samples were collected from each sample site to measure Turb using a Hach© 2100P PORTABLE Turbidometer in NTU [41]. A black and white Secchi disc was used to measure transparency from the shady side of the boat. The depths at which the disc disappeared and then reappeared were taken, and the average value was reported as the SD.

2.3. Satellite Imagery

This study used imagery downloaded from the Planet Explorer website, the USGS Earth Explorer website for Landsat-8, and the Copernicus Open Access Sentinel Hub. Satellite images with less than 10% cloud cover were considered for analysis. When the cloud cover was more than 10% but did not cover the sample sites, the image was considered for analysis.

The PS images used in this study were acquired by the PS2.SD Dove satellites. The PlanetScope Lab geometrically corrected the images using Ground Control Points (GCPs) and digital elevation models (DEM) based on the Universal Transverse Mercator World Geodetic System 1984 (UTM WGS1984). The images were atmospherically corrected to surface reflectance using the 6SV2.1 radiative transfer code. The PlanetScope Lab also did radiometric correction using sensor telemetry and sensor model. The processed images were made available as open source for educational purpose. The image frame size is approximately 24 km × 16 km. The PS2.SD constellation provides four-band imagery globally with a spatial resolution of 3 m, and they have a daily flyover schedule. Each image is stored with information on the date, time of acquisition, and satellite ID [42].

The Landsat-8 images used in this study were obtained by the Operational Land Imager (OLI) every 16 days, with spatial resolution 30 m, and included 11 bands of analysis ready surface reflectance data. The USGS Earth Resource Observation and Science (EROS) Center provides open-source Landsat Analysis Ready Data (ARD) for the United States with image frame size 5 km × 5 km. The images are processed to the highest scientific standards and level of processing required for direct use in assessing land cover change. The Landsat Surface Reflectance Code (LaSRC) uses the coastal aerosol band to perform aerosol inversion tests, auxiliary climate data from MODIS, and a unique radiative transfer model. The LaSRC hardcodes the view zenith angle to “0” and the solar zenith and view zenith angles are used to correct atmospheric scattering.

The ARD images are projected to the Albers Equal Area Conic (AEA) map projection and are processed directly from Level-1 AEA scenes through Level-2 products using the WGS84 datum. Level-1 processing includes radiometric calibration and orthorectifying using GCPs and DEM data, and geometric correction using the spacecraft ephemeris and DEM data. Level-2 processing includes atmospheric correction [43].

The Sentinel-2 Multi-Spectral Instrument (MSI) provides open access images available to users via the Copernicus Open Access Hub. The Level-2A products, which are orthorectified Bottom-Of-Atmosphere (BOA) reflectance in cartographic geometry, were generated through the Payload Data Ground Segment (PDGS) utilizing the Sen2Cor processor. Radiometric and geometric corrections, including orthorectification and spatial registration, are done at the Level 1 stage of image processing. The image is divided into 100 km tiles in UTM WGS84 projection. Sentinel-2 imagery is acquired every five days by a constellation of two polar-orbiting satellites placed in the same sun-synchronous orbit, phased at 180° to each other [26]. The Sentinel-2 instruments acquire images with 13 spectral bands and spatial resolutions 10-60 m. Resampling was done to allow a combination of 20 m and 10 m bands to compute and map the water quality parameters.

Band spectra from the three satellites were extracted using the extraction tool in the spatial analyst toolset in ArcGIS 10.7.1. The extracted values were multiplied by the scaling factor (0.0001) that was used for data storage. The surface reflectance values were paired with their corresponding Chl-a, Turb, and SD data for regression analysis. For reservoir polygons that required more than one tile of PS imagery, we ensured that they all had the same satellite ID before they were used to extract spectral data.

2.4. Band Combination and Band Selection

Multiple linear regression with stepwise selection of terms was used to determine empirical relationships between each of the three water quality parameters (the dependent variables) and their corresponding surface reflectance values from the PS, Landsat-8, and Sentinel-2 imagery. The independent variables included the spectral bands in each satellite, band ratios, band indices, and spectral derivatives. Table 3 highlights the properties of each satellite band that are relevant for water quality monitoring and were used in the study. These band ratios are empirical scenarios informed by the spectral properties discussed in Table 3. Transforming bands for the purpose of improving the predictive tendency of empirical models is a common practice in remote sensing of water quality [9,18,44].

2.5. Best-Fit Model Selection and Validation

The first step in the model building was power transformation of the dependent water quality variables to increase the chances of improving the regression fits. Each of the three water quality parameters was transformed into its natural logarithm (LN) and square root (SQRT), and all three data types (actual values, LN, and SQRT) were used in the generation and selection of the best fit model.

This was followed by generating several independent variables through different methods of band combinations to increase the chances of developing strong regression models. The first set of independent variables included single bands, band ratios, and band indices. The band ratios were obtained by dividing two bands (e.g., ρ_R*ρ_NIR⁻¹). The band indices were the square of the normalized difference between surface reflectance values in two or more bands. These were paired with corresponding Chla, Turb, and SD (actual, LN, and SQRT) values for the development of regression models. The resulting statistics were recorded and evaluated.

In the following step, the independent variables were the surface reflectance values divided by their respective center wavelengths and 1st and 2nd derivatives of the surface reflectance spectra [11]. The first derivatives were designated d_n with values ranging from d₁ = (ρ_band-3 − ρ_band-2) ∗ (λ_band-3 − λ_band-2)⁻¹ to d_{n−1 =} (ρ_band-(n+1) − ρ_band-n) ∗ (λ_band-(n+1) − λ_band-n)⁻¹. The 2nd derivatives were designated d_2n with values ranging from d₂₁ = (d₂ − d₁) ∗ (λ_band-3 − λ_band-1)⁻¹ to d_2(n−1) = (d_n+1 − d_n) ∗ (λ_band-(n+2) − λ_band-n)⁻¹. The value n is the band number, λ is the center wavelength, and ρ is the surface reflectance. These were combined with their corresponding water quality parameters for the development of regression models. Finally, we deployed all the independent variables generated above with their paired water quality parameters and ran them in the regression analysis.

The criteria considered in this study were the p-value, R², Root Mean Square Error (RMSE), and Variance Inflation Factor (VIF). The regression terms that gave the best satisfaction of these criteria were selected for model building. The criteria were set as follows:

• The relationships should be significant on a 95% confidence interval (α ≤ 0.05);

• A strong relationship between the dependent variable and the predictors (R² ≥ 0.5);

• Low standard deviation of the residuals (RMSE) relative to the range of values;

• Low correlation between the predictors (VIF < 10).

Selection of the best fit model was achieved using Minitab-19 statistical software [49]. The null hypothesis was that the surface reflectance from the bands, band combinations, and spectral derivatives did not predict variations in Chl-a, Turb, or SD, at a significance level (α) of 0.05. Each of the predictors was selected and evaluated using the Forward Information Criteria (FIC) in Minitab. The function is a stepwise selection that takes away and adds one variable at a time and determines selected statistics as it proceeds. This function is exhaustive and reshuffles variables until it determines that the best criteria are satisfied.

Model-derived values of each parameter were plotted against measured data that matched sampling periods and sites. At least 40 data points were randomly selected from the data and set aside for comparison to model derived values in scatter plots. The R² was used to determine the extent to which the model-derived data agreed with the measured data.

Each best fit model was further tested using the K-folds cross-validation (cv) in Minitab. The advantage of the K-folds cv method over a single training/testing split is that every data point is used in both the training and testing sets. The data is divided into a certain number of partitions (K-folds) and each partition is used as a testing set while the rest of the data is used as a training set [50]. The average R² and RMSE for 10 folds of cv were used to determine the validity of the models obtained. The procedure for the development of models, selection of the best fit model, and its validation are summarized in Figure 2.

2.6. Case Study Application—Algal Bloom in Lake McMurtry, Oklahoma

We used an algal bloom event that occurred in Lake McMurtry, northcentral Oklahoma, from 27 November to 3 December 2019, to compare the spatial-temporal properties of the three satellites. PS had 4 cloud-free images (27 November, 30 November, 1 December, and 3 December 2019) available during this bloom event. In contrast, Landsat-8 had only one image with a cloud-free portion in the northern part of the lake (27 November 2019). Sentinel-2 also had only one image (1 December 2019) with a cloud-free portion that included the lake. These images were used to compare PS to Landsat-8 and Sentinel-2 in terms of their temporal and spatial resolutions. The Spatial Analyst tool of ArcMap 10.7.1 was used to compute pixel values of Chl-a in the areas of interest.

3. Results

3.1. Range of Values of the Three Parameters

Chl-a concentrations in the study reservoirs ranged from 0.60 to 540 µg/L. Turb ranged from 0.00 to 966 NTU and SD ranged from 8.0 to 400 cm. PS had more data points to match with water quality than Landsat-8 and Sentinel-2 because of its daily flyover schedule. For PS, there were 207 data points for Chl-a, 225 data points for Turb, and 158 data points for SD. For Landsat-8, there were 100 data points for Chl-a, 75 data points for Turb, and 96 data points for SD. For Sentinel-2, there were 129 data points for Chl-a, 101 data points for Turb, and 113 data points for SD. Figure 3 shows the range of values of Chl-a, Turb, and SD.

3.2. Best Fit Models

Equations (3)–(11) present the best fit regression equations for the three satellites.

PS based equations:

SQRT(Chl-a) = 1.24 − 89.10ρ_B + 79.10ρ_NIR + 9.67ρ_{((R+B−NIR)/(R+B+NIR))}² + 13.67ρ_{((G−NIR)/(G+NIR))}²(3)

LN(Turb) = 1.65 + 2.30ρ_(R/G) − 0.035ρ_(1/NIR) − 0.07ρ_(1/R)(4)

LN(SD) = 4.22 − 0.90ρ_(R/B) + 0.02ρ_(1/NIR) + 0.07 ρ_(1/R)(5)

where ρ is the surface reflectance value at a wavelength within the bands.

Landsat-8 based equations:

SQRT(Chl-a) = 4.778 − 45.5d₄ − 123.03d₅ − 26.33ρ_(B/λB)(6)

where d₄ = (ρ_SWIR1 − ρ_NIR)/(λ_SWIR1 − λ_NIR), d₅ = (ρ_SWIR2 − ρ_SWIR1)/(λ_SWIR2 − λ_SWIR1), λ_B, λ_SWIR1 and λ_SWIR2 are the center wavelength of the respective bands.

SQRT(Turb) = 4.67 + 9.86d₂ + 41.70d₂₄(7)

where d₂ = (ρ_R − ρ_G)/(λ_R − λ_G), and d₂₄ = (d₅ − d₄)/( (λ_SWIR2 − λ_NIR).

LN(SD) = 4.1974 − 2.018d₂ + 1.227d₃ + 18.76d₄ − 18.22d₅(8)

where d₂ = (ρ_R − ρ_G)/(λ_R − λ_G), d₃ = (ρ_NIR − ρ_R)/(λ_NIR − λ_R)

Sentinel-2 Based Equations:

SQRT(Chl-a) = 2.986 + 6.936d₃ − 0.861d₄ − 91.4d₈(9)

where d₃ = (ρ_RE1 − ρ_R)/(λ_RE1 − λ_R), d₄ = (ρ_RE2 − ρ_RE1)/(λ_RE2 − λ_RE1); and d₈ = (ρ_SWIR2 − ρ_SWIR1)/(λ_SWIR2 − λ_SWIR1)

LN(Turb) = 0.076 − 29.80ρ_SWIR2 + 4.72ρ_(R/B) − 3.38ρ_{((R−NIR)/(R+NIR))}²(10)

LN(SD) = 6.66 + 28.81ρ_SWIR2 − 15.34ρ_RE3 − 2.86ρ_(R/G)(11)

All the regression models for the three parameters provided significant relationships between the spectral data from each platform and their respective parameters. The p-values for the coefficients and constants were also significant (p < 0.01). Table 4 presents results of the statistical analysis for the overall regressions, coefficients, and constants in Equations (3)–(11).

The Chl-a R² values for Sentinel-2 (R² = 0.85) and Landsat-8 (R² = 0.75) showed that the spectral data from these two satellites predict Chl-a better than PS (R² = 0.58). Their Chl-a RMSE values were also lower than that for PS. These results imply that Landsat-8 and Sentinel-2 are preferable for detecting Chl-a reflectance over PS. The low VIF values (Maximum VIF < 10) in all three platforms indicate minimal collinearity between the regression terms. This was the case for all the parameters and all three satellites.

For Turb, PS and Sentinel-2 had comparable R² values (0.79 and 0.78, respectively), while the R² value for Landsat-8 was lower (R² = 0.60). All the RMSE values were relatively low, indicating reliability in using all three platforms for sensing Turb in the study reservoirs. PS and Sentinel-2 are the preferred satellites for sensing Turb in the study reservoirs based on their R² values.

For SD, the R² in Sentinel-2 was the highest (R² = 0.80) followed by PS (R² = 0.76). The Landsat-8 value for SD was the lowest (R² = 0.58). The RMSEs were low relative to the range of values in all three satellites. Based on these results, Sentinel-2 and PS are the preferred platforms for sensing SD in the study reservoirs.

Figure 4 presents scatter plots showing relationships between the values of Chl-a (µg/L), Turb (NTU), and SD (cm) predicted by the PS, Landsat-8 (L8), and Sentinel-2 (S2) spectra and those measured in the study reservoirs. The R² (R-sq) values are displayed at the top of each graph along with the associated parameters and satellite platforms. The x-axis represents the model-derived values, and the y-axis represents the measured values. All three satellites showed strong relationships between measured and model-derived variables for the water quality parameters.

Validation of the Best Fit Models

Figure 5 shows the average R² values obtained in the 10-fold cv of the parameters with PS, Landsat-8, and Sentinel-2. For Chl-a, Sentinel-2 had the highest average value (R² = 0.83), followed by Landsat-8 (R² = 0.69) and then PS (R² = 0.55). PS had the highest average R² for Turb (R² = 0.78) followed by Sentinel-2 (R² = 0.75), and Landsat-8 (R² = 0.51). For SD, Sentinel-2 had the highest value (R² = 0.78), followed by PS (R² = 0.74), and Landsat-8 (R² = 0.53).

Considering R², the 10-folds cv showed that all three satellites are reliable for sensing all three parameters. Comparatively, Landsat-8 and Sentinel-2 were more reliable than PS for Chl-a sensing, while PS and Sentinel-2 were more reliable than Landsat-8 for Turb and SD sensing.

Figure 6 shows the average RMSE values obtained for the three parameters with the three satellites. The average RMSE values were below 5.0 for all three cases, which are considered low in reference to the range of values utilized in the study. The RMSE values for Turb and SD were comparable in the three satellites, while the PS RMSE for Chl-a (4.62 µg/L) was almost twice the amount in Landsat-8 (2.37 µg/L) and 3.5 times the Sentinel-2 amount (1.32 µg/L).

Similar to the observations made in the R² validation results, the 10-folds cross validation showed that all three satellites are reliable for sensing all three parameters. The RMSEs also revealed that Landsat-8 and Sentinel-2 were more reliable than PS for Chl-a sensing, while PS and Sentinel-2 were more reliable than Landsat-8 for Turb and SD sensing.

3.3. Case Study Application—Algal Bloom in Lake McMurtry, Oklahoma

The onset of an algal bloom was reported in Lake McMurtry on 27 November 2019. The bloom grew and dissipated between 27 November and 5 December 2019. The Landsat-8 image that was available during these blooming dates (27 November) had a cloud-free portion in one of the areas where the bloom occurred. Sentinel-2 had only one available image with a cloud-free portion that included Lake McMurtry on December 1. There were four cloud-free PS images available during the blooming period. This event allowed for temporal and spatial comparisons in Chl-a mapping using models derived for PS, Landsat-8, and Sentinel-2. Figure 7 shows the PS and Landsat-8 maps for Chl-a on 27 November, 2019. Note that Chl-a concentrations were more visible in the littoral and riverine coves of the reservoir just as the bloom was observed when these sites were visited.

The PS-based Chl-a map showed the same pattern of Chl-a distribution as the Landsat-8 map. This shows that the models developed for both satellites estimate Chl-a similarly. The key difference in Figure 7 lies in the sizes of the pixels for PS (3 m) and Landsat-8 (30 m). The PS image gave a higher definition (finer spatial resolution) of Chl-a distribution in the focus area of the reservoir compared to Landsat-8. This observation is important for water quality sensing at the edges of the reservoir and in the riverine coves. The coarser spatial resolution of Landsat-8 could have caused the lower concentration range (13.09–147.64 µg/L) in those edges and coves compared to the PS range of values (0.00–296.13 µg/L). Figure 7 shows that the small pixel size from PS is better at detecting Chl-a near the lake’s edge and in the coves of the reservoir, which is where the bloom was most concentrated.

Figure 8 shows the PS and Sentinel-2 maps for Chl-a on 1 December 2019. We used this figure to compare Chl-a from PS imagery to that from Sentinel-2 imagery for this single date during an algal bloom. The PS-based Chl-a models gave the same pattern of Chl-a distribution in both PS and Sentinel-2 maps, demonstrating that both models estimated Chl-a similarly. The PS image gave a higher definition (finer spatial resolution) of Chl-a distribution in the focus area of the reservoir compared to Sentinel-2. Unlike Landsat-8, the Sentinel-2 model gave a similar concentration range as the PS model, probably owing to its higher spatial resolution (10 m). Figure 8 again shows that on 1 December 2019, the small pixel size from PS was better at defining the distribution of Chl-a (than Sentinel-2) near the lake’s edge and in the riverine coves of the reservoir, which was where the bloom was most concentrated.

Figure 9 presents PS Chl-a maps in the northwest cove of the reservoir. The images were acquired on dates between 27 November and 3 December. The algae had dispersed in the cove by 3 December. The bloom peaked and started dissipating just before this time, and no bloom was observed during ground-based observation on 5 December. Figure 9 shows several PS images from the same bloom from 27 November to 3 December, 2019. These images show the dynamic nature of the algal bloom as it changed daily. We also note that during this period of the bloom, there was only a single image available for both Landsat-8 and Sentinel-2. Of the 3 satellite platforms, PS was the only one with enough available/usable imagery to follow the dynamic nature of the algal bloom during this period.

4. Discussion

4.1. PlanetScope Nanosatellites

PlanetScope (PS) nanosatellites were originally developed for terrestrial and aquatic applications [51]. Our research contributes to a growing body of research [28] showing that PS is also applicable for use in water quality management applications. The accuracy of the models developed from PS spectra (based on R² and RMSE values) for three key water quality parameters were relatively similar to those that have been developed for other water bodies using other satellites such as Landsat-8 and Sentinel-2 [22,25,44].

The independent variables of the Chl-a, Turb, and SD models for PS were the single bands, band ratios, and band indices. The spectral bands in the best fit PS model for Chl-a were B, NIR, ((R + B − NIR)/(R + B + NIR))², and ((G − NIR)/(G + NIR))². Chl-a reflects G and absorbs B and R during photosynthesis [15], which is probably why the B band has a negative linear relationship with Chl-a in the equation. We suggest that the square of the normalized difference between the visible bands (R + B) and the NIR band explains the impact of Turb on Chl-a reflectance, since both visible and NIR bands interact with both parameters [10,52]. The NIR band has also been used in previous satellite studies, especially at wavelengths close to the Red-Edge. NIR has been used mostly in combination with R [16]. Red may not be applicable in this case due to possible interference from Turb; agricultural runoff, shoreline erosion, and in-lake resuspension make Turb a primary source of reservoir impairment in the study area. The square of the normalized difference between G and NIR is not as common as that between R and NIR [10] at specific wavelengths. Our empirical results indicate that ((G − NIR)/(G + NIR))² contributes significantly to Chl-a detection.

The independent variables that were selected for Turb were R/G, 1/NIR, and 1/R. The positive linear relationship between R/G and the natural log of Turb likely explains the ferric oxide-laden sediment load in the study area. The second and third terms had a negative linear relationship with Turb, and their effect may be minimal when the exponential function is considered. The opposite is the case for SD in which the inverse values of R and NIR gave positive linear relationships while the R/B gave a negative linear relationship with the parameter. Oklahoma’s Turb is mostly reddish in color, and it has an inverse relationship with SD [47]. Additionally, since water absorbs NIR, clear waters (high SD) would translate to low NIR reflectance. Similarly, R is associated with Turb and minimizes light penetration. The reverse order of independent variables for Turb and SD is consistent with past studies showing an inverse relationship between Turb and SD [53].

4.2. PlanetScope Compared to Landsat-8 and Sentinel-2

Although PS has finer spatial and temporal resolutions than Landsat-8 and Sentinel-2, the latter two satellites have more spectral bands. These additional bands contributed to better models for Chl-a in terms of the coefficient of determination (R²) and the error values. Those additional bands are not as susceptible to Turb as the NIR band [31]. Additionally, the spectral derivatives showed strong relationships with Chl-a in Landsat-8 and Sentinel-2. They also showed strong relationships with Turb and SD in Landsat-8.

In addition to the visible bands, the two SWIR bands in Landsat-8 and Sentinel-2 contributed significantly to Chl-a reflectance in the study reservoirs. This observation is consistent with the results obtained in previous studies. The SWIR band (λ = 1500–2500) appears to be less susceptible to scattering than the NIR band because of its longer wavelengths. This allows minimal Turb interference to this band compared to NIR [31]. The Red-Edge band in Sentinel-2 has been shown to better detect Chl-a than Landsat-8 bands [32] because it is not as susceptible to absorption or scattering. Wavelengths in this region have a strong linear relationship with Chl-a reflectance [54].

The PS and Sentinel-2 models outperformed Landsat-8 for Turb and SD, probably owing to their finer pixel sizes (3 and 10 m, respectively). Even though Sentinel-2 has a coarser spatial resolution than PS, our study suggests that the 10 m pixel is fine enough for detecting Turb in the study reservoirs. The agricultural practices, along with streambank and shoreline erosion, and in-lake resuspension of sediment in the study area are common sources of continuous sediment transport, which elevates Turb levels in waters in the region. The same explanation applies to SD since its values trend inversely with Turb.

PS is a commercial satellite whereas Landsat-8 and Sentinel-2 images are open source. Therefore, it may not be economically feasible to use PS to remotely sense water quality for all water resource applications. The bloom event from Lake McMurtry highlights a scenario where the finer spatial and temporal resolutions of PS can be used to remotely sense Chl-a. The algal bloom developed over several days in the reservoir (Figure 9). During this period, there were several cloud-free PS images compared to only one for Landsat-8 and Sentinel-2 images during the blooming event. These return periods of PS allowed for the detection and monitoring of bloom dynamics that changed daily and could not be measured from the other satellite platforms. This temporal advantage of PS provides an opportunity for lake managers to make informed decisions regarding response efforts such as the duration of beach closures during HABs.

The finer spatial resolution of PS may allow for the remote sensing of water quality in smaller ponds and reservoirs compared to coarser-resolution satellites. Small (farm) ponds are important in agricultural watersheds and they often experience high nutrient and sediment loads. Oklahoma has thousands of ponds ranging in size between 0.001 km² and 0.04 km². Because of their small size, they are difficult to monitor with ground-based sampling or with Landsat-8 or Sentinel-2. For example, a 0.001-km² pond would be covered with only one Landsat-8 pixel or 10 Sentinel-2 pixels (considering Sentinel-2 data with 10-m pixel size). For the same 0.001-km² pond, a PS image would have approximately 111 pixels. As such, PS data with finer spatial resolution (3 m) could be used to assess more pixels in the ponds and better detect small-scale changes in water quality.

This study covered 13 reservoirs located in predominately agriculturally dominated watersheds that are dominated by agricultural practices. Additional studies are needed to determine if our models developed for agricultural landscapes could be used in other reservoirs in the region and around the world. Additional studies are also required to test and refine our models in areas with dominant human activities other than agriculture.

A limitation to using PS imagery is that it has fewer spectral bands compared to Landsat-8 and Sentinel-2. PS also has many satellites each providing image tiles with properties that may vary from each other. This may require normalization of those image tiles, which was beyond the scope of this study. We selected image tiles with the same satellite ID to minimize the likelihood of this source of error. However, a future application of PS imagery for water quality monitoring requires normalization of multiple image tiles to ensure comparability. Such image analysis, in addition to acquiring the imagery, makes it costly for PS use for water quality monitoring at a regional scale. Furthermore, the scalability of PS imagery for cost recovery in water quality monitoring is limited to areas with high frequency of cloud-free days. A key benefit of PS data is using its finer spatial and temporal resolutions to measure dynamic water quality parameters in small reservoirs and coves of larger reservoirs that are not feasible to monitor with moderate-coarser resolution satellites and are inaccessible for ground-based monitoring.

5. Conclusions

Traditional water quality monitoring is limited to specific sample locations and infrequent sample visits. These limitations create spatial and temporal gaps during which water quality impairment may occur and pose exposure risks. This study explored the utility of satellites to detect water quality in those spatial and temporal gaps in agricultural watersheds. We compared the spatial, temporal, and spectral resolutions of PS to those of Landsat-8 and Sentinel-2 for sensing three water quality parameters in Oklahoma reservoirs. We used three years (2017–2020) of surface reflectance spectra for PS, Landsat-8, and Sentinel-2 to develop empirical models for Chl-a, Turb, and SD in 13 reservoirs in Oklahoma. Multiple linear regression models were developed using surface reflectance from the bands, band ratios, and band indices of these water quality parameters.

Overall, the spectral data from each of the three satellite platforms (PS, Landsat-8, and Sentinel-2) had significant relationships with Chl-a, Turb, and SD. We conclude that the finer spatial and temporal resolutions in PS give it the capabilities that are comparable to Landsat-8 and Sentinel-2 for sensing key water quality parameters in the study reservoirs. However, this capability is moderated by the fewer spectral bands and the cost to acquire the PS imagery given that Landsat-8 and Sentinel-2 data are free. We identified scenarios where it may be more beneficial to use PS compared to Landsat-8 and Sentinel-2. PS is better suited for measuring water quality variables that vary daily, in small ponds and reservoirs where satellites with coarser spatial resolution do not provide enough pixels and at the edge of water bodies where variation in water quality requires fine-scale detection.

Future research should work to combine PS, Landsat-8, and Sentinel-2 to develop a predictive monitoring tool that utilizes their tradeoffs to increase our ability to understand water quality dynamics. Key areas of a combined satellite approach [55,56,57,58] could include targeted monitoring of hot spot areas, frequent monitoring of phytoplankton and sediments from nutrient-rich runoff, wider spatial and temporal coverage than ground-based monitoring programs can accomplish, and ultimately, effective use of resources for water quality monitoring.

Author Contributions

Conceptualization: all authors; methodology: A.S.M.; software: A.S.M. and M.E.M.; validation: A.S.M., H.G., A.R.D. and M.E.M.; formal analysis: all authors; investigation: A.S.M.; resources: K.L.W., S.H.S. and A.R.D.; data curation: A.S.M.; writing—original draft preparation: A.S.M.; writing—review and editing: all authors; supervision: K.L.W., S.H.S. and A.R.D.; project administration: K.L.W., S.H.S. and A.R.D.; funding acquisition: K.L.W., S.H.S. and A.R.D. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

The Grand River Dam Authority and the Water Quality Division of the Oklahoma Water Resources Board provided the data for this study. The Lake McMurtry Friends provided support for algal bloom observation from 27 November to 5 December 2019.

Conflicts of Interest

The authors declared no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Map of Oklahoma showing the BUMP reservoirs across the state (grey polygons and speckles) and the 13 reservoirs used in this study (in colors). The spatial data for this map were obtained from the Oklahoma Water Resource Board’s website [34].

Figure 2. Flow diagram showing the procedure for model development, model selection, and model validation.

Figure 3. Histograms of the Chl-a (µg/L), Turb (NTU), and SD (cm) data from the study reservoirs that correspond to their respective PS, Landsat-8, and Sentinel-2 images.

Figure 4. Scatter plots showing relationships between the model-derived concentrations and measured concentrations of Chl-a (µg/L), Turb (NTU), and SD (cm) with PS, Landsat-8 (L8), and Sentinel-2 (S2). The R2 (R-sq) values are displayed at the top of each graph along with the associated parameters and satellite platforms.

Figure 5. Average R2 values in the 10-fold cv of the parameters with PS, Landsat-8, and Sentinel-2.

Figure 6. Average RMSE values in the 10-fold cv of the parameters with PS, Landsat-8, and Sentinel-2.

Figure 7. PS (Left) and Landsat-8 (Right) maps of Chl-a in Lake McMurtry, northcentral Oklahoma. Both satellites acquired their images on 27 November 2019, during an active algal bloom event. The overview maps at the top show the location of Lake McMurtry in Oklahoma (Top Right spec in a red box) and Lake McMurtry (Top Left) showing the focus area delineated in a red box. The color bars represent concentration ranges as estimated by each of the two satellites.

Figure 8. PS (Left) and Sentinel-2 (Right) maps of Chl-a in Lake McMurtry, northcentral Oklahoma. Both satellites acquired their images on 1 December 2019, during an active algal bloom event. The overview maps at the top show the location of Lake McMurtry in Oklahoma (Top right spec in a red box) and Lake McMurtry (Top left) showing the focus area delineated by a red box. The color bars represent concentration ranges as estimated by each of the two satellites.

Figure 9. PlanetScope maps of Chl-a on four dates of image acquisition in Lake McMurtry, northcentral Oklahoma. The upper left panel is from an image acquired on 27 November 2019; the upper right panel is from an image acquired on 30 November 2019; the lower left panel is from an image acquired on 1 December 2019; the bottom right panel is from an image acquired on 3 December 2019. The color bars in the legends represent Chl-a concentration ranges on each date.

Table 1

Characteristics of Landsat-8, Sentinel-2, and PS.

Characteristics	PS	Landsat-8	Sentinel-2
Revisit time (temporal resolution)	Daily	16 days	10 days with each satellite (Sentinel 2A and 2B). Five days with combined satellites.
Spectral resolution	Four 3-m bands	Eight 30-m bands, two 100-m bands, one 15-m panchromatic band (11 bands)	Four 10-m bands, six 20-m bands, and three 60-m bands (13 bands)
Pixel size (spatial resolution)	More pixels in small areas/reservoirs	Few or no pixels in small areas/reservoirs (e.g., with area 0.001 km2 or less)	Few pixels in small areas/reservoirs (e.g., with area 0.001 km2 or less)
Bandwidth in nm (visible and NIR)	Blue: 465–517; Green: 547–595; Red: 650–682; NIR: 846–888	Blue: 435–512; Green: 533–590; Red: 636–673; NIR: 851–879, Shortwave IR1 (SWIR1): 1570–1650; Shortwave IR2 (SWIR2): 2110–2290	Blue: 458–523; Green: 543–578; Red: 650–680; Red-Edge (RE1): 698–713; Red-Edge (RE2): 733–748; Red-Edge (RE3): 773–793; NIR: 785–899; SWIR1: 1565–1655; SWIR2: 2100–2280
Availability of free imagery	10,000 km2 per month for education purpose	Unlimited	Unlimited

Table 2

Characteristics of the reservoirs used in this study. This information was obtained from the Oklahoma Water Resource Board’s website [37].

Reservoir	Surface Area (km2)	Trophic Status	Impairment Status
			Chl-a	Turb
Arcadia	7.40	Hypereutrophic	Impaired	Impaired
Broken Bow	57.50	Mesotrophic	Not impaired	Not impaired
Canton	32.00	Hypereutrophic	Insufficient data	Impaired
Eucha	11.60	Eutrophic	Impaired	Not impaired
Fort Gibson	60.30	Eutrophic	Insufficient data	Not impaired
Foss	35.61	Mesotrophic	Insufficient data	Impaired
Hefner	10.11	Hypereutrophic	Insufficient data	Not impaired
Grand	188.20	Eutrophic	Insufficient data	Not impaired
Kaw	68.96	Hypereutrophic	Insufficient data	Impaired
McMurtry	4.67	Eutrophic	Insufficient data	Impaired
Oologah	119.22	Mesotrophic	Insufficient data	Impaired
Thunderbird	24.60	Hypereutrophic	Impaired	Impaired
Waurika	40.87	Eutrophic	Impaired	Impaired

Table 3

Spectral bands and their spectral properties that allow for monitoring Chl-a, Turb, and SD.

Spectral Bands and Band Ratios	Wavelength Range, nm (λi–λn; i = 1)	Properties
Blue (B)	PS: λ465–λ517Landsat-8: λ435–λ512Sentinel-2: λ458–λ523	This is the region of deepest light penetration in clear waters. However, most of Oklahoma lakes are turbid. The B band is susceptible to scattering in the atmosphere and water [45].
Green (G)	PS: λ547–λ595Landsat-8: λ533–λ590Sentinel-2: λ543–λ578	The reflectance peak of different concentrations of Chl-a are at wavelengths in this region [46].
Red (R)	PS: λ650–λ682Landsat-8: λ636–λ674Sentinel-2: λ650–λ680	The Chl-a absorption peak is at λ660 [15], which falls within the R band. Ferric-rich soils in Oklahoma [47] end up in reservoirs through surface runoff, making the R band a crucial spectral signature for Turb (reflectance), and also for Chl-a and SD detection when used as a ratio to other bands.
Near-infrared (NIR)	PS: λ846–λ888;Landsat-8: λ851–λ879;Sentinel-2: λ785–λ899	This band is absorbed in water [15]. Its high reflectance will indicate the presence of substances other than water.
Red-Edge (RE)	RE1: λ698–λ713RE2: λ733–λ748RE3: λ773–λ793	The RE band transitions between the R and NIR bands, and it uniquely correlates with Chl-a [32]
Shortwave infrared (SWIR)	Landsat-8:SWIR1: λ1570–λ1650SWIR2: λ2110–λ2290Sentinel-2:SWIR1: λ1565–λ1655SWIR2: λ2100–λ2280	The longer wavelengths in the SWIR band give it the advantage of minimal scattering by mineral Turb in the water, making it suitable to detect algal pigments [31]. It is also useful to differentiate between algal pigments and those in aquatic macrophytes [48]

Table 4

Statistics of the best fit regression models for PS, Landsat-8 (L8), and Sentinel-2 (S2).

Parameter	R2			RMSE			Maximum VIF
	PS	L8	S2	PS	L8	S2	PS	L8	S2
Chl-a	0.58	0.75	0.85	4.41 µg/L	2.04 µg/L	1.19 µg/L	5.47	1.20	2.74
Turb	0.79	0.60	0.78	1.61 NTU	1.54 NTU	1.60 NTU	2.37	1.43	2.57
SD	0.76	0.58	0.80	1.54 cm	1.50 cm	1.35 cm	2.01	3.59	6.90