Intercalibration of MERIS, MODIS, and OLCI Satellite Imagers for Construction of Past, Present, and Future Cyanobacterial Biomass Time Series

: Satellite imagery has been used to monitor and assess Harmful Algal Blooms (HABs), speciﬁcally, cyanobacterial blooms in Lake Erie (the USA and Canada) for over twelve years. In recent years, imagery has been applied to the other Great Lakes as well as other U.S. lakes. The key algorithm used in this monitoring system is the cyanobacterial index (CI), a measure of the chlorophyll found in cyanobacterial blooms. The CI is a “spectral shape” (or curvature) algorithm, which is a form of the second derivative around the 681 nm (MERIS/OLCI) or 678 nm (MODIS) band, which is robust and implicitly includes an atmospheric correction, allowing reliable use for many more scenes than analytical algorithms. Monitoring of cyanobacterial blooms with the CI began with the European Space Agency’s (ESA) Medium Resolution Imaging Spectrometer (MERIS) sensor (2002–2012). With the loss of data from MERIS in the spring of 2012, the monitoring system shifted to using NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS has bands that allow computation of a CI product, which was intercalibrated with MERIS at the time to establish a conversion of MODIS CI to MERIS CI. In 2016, ESA launched the Ocean and Land Color Imager (OLCI), the replacement for MERIS, on the Sentinel-3 spacecraft. MODIS can serve two purposes. It can provide a critical data set for the blooms of 2012–2015, and it offers a bridge from MERIS to OLCI. We propose a basin-wide integrated technique for intercalibrating the CI algorithm from MODIS to both MERIS and OLCI. This method allowed us to intercalibrate OLCI CI to MERIS CI, which would then allow the production of a 20-year and ongoing record of cyanobacterial bloom activity. This approach also allows updates as sensor calibrations change or new sensors are launched, and it could be readily applied to spectral shape algorithms.


Introduction
Blooms of cyanobacteria are generally considered a harbinger for anthropogenic eutrophication [1]. Potentially massive blooms can have long durations and detrimental environmental and human health impacts [2]. Cyanobacterial blooms are often toxic and can contribute to hypoxia when the blooms senesce and increase the biological oxygen demand. These blooms have a particular affinity for warm stratified water, and as such, blooms of cyanobacteria may become more prevalent in the context of a warming climate [3]. Cyanobacteria can produce various toxins (such as microcystins, anatoxins, and saxitoxin) that pose a health risk and cause mortalities in domestic and wild animal populations [4]. Additionally, several species have compounds such as geosmin that can cause taste and odor issues in drinking water. As a result of these detrimental influences, cyanobacterial blooms are classified as Harmful Algal Blooms (HAB) and should be monitored to reduce potential deleterious impacts.
Perhaps the most effective monitoring, detection, and assessment tool for describing cyanobacterial bloom dynamics is synoptic remote sensing [5,6]. In order to detect any significant increase in cyanobacterial bloom concentration and frequency, there must be an effective and robust algorithm to estimate the extent and severity of these blooms. After an effective algorithm is determined, there needs to be a way to merge the detection technique across different sensors as new ones are deployed and older ones cease operations. These satellite intercalibrations work under the basic premise that two instruments should make identical observations over the same place and time with consistent spatial and spectral responses and identical view angles [7]. Obviously, this is impossible in a real-world situation. However, it is possible to compensate for any error effect by applying a correction if biases are small and the relationship between the two remotely sensed parameters is linear. Many intercalibration studies have been conducted using consistent, bright targets as the intermediary. For example, Yu et al. [8] use the Sonoran Desert, a consistently bright target that is visible to both the GOES-East and GOES-West satellites, to perform an intercalibration between the two satellites based on getting the same reflectance from each sensor. Bouvet et al. [9] describe an important new calibration network (The Radiometric Calibration Network; RadCalNet) that can be used for vicarious calibration based on a set of reference sites of relatively stable surface reflectance in different parts of the world. This network is an important advancement in calibration for land targets. The aquatic environment is far too dynamic for a similar correction technique (and also, the MODIS ocean color bands saturate over bright land targets). A more typical approach for aquatic applications is a point-by-point or pixel-to-pixel matchup [10,11]. The pixel-topixel matchup has some caveats. Spatial heterogeneity in water bodies can be substantial and subject to change. Instrument characteristics (ground field-of-view, sensor response functions and others) can all introduce mismatches between sensors even when viewing the same point on the water. The result is noise and biases in the resulting comparison. Here, we propose a technique using pixel integration across the cyanobacterial blooms detected within a basin as a means of reducing analytical error. We compare this technique with the more traditional pixel-to-pixel intercalibration analysis.

CI Algorithm Usage
Since 2009, the National Oceanic and Atmospheric Administration (NOAA) has released short-term (<1 week) forecasts of cyanobacterial blooms in the western basin of Lake Erie [12]. The primary product used is a satellite proxy for cyanobacterial biomass called the cyanobacterial index (CI). The CI has been used extensively in various freshwater systems throughout the United States [13][14][15][16]. While the algorithm was initially developed by Wynne et al. [17,18] for the MERIS sensor on-board the Envisat spacecraft, it has since been successfully applied to MODIS [10,11]. This was necessary as the Envisat spacecraft ceased operation in April 2012. In early 2016, ESA launched the follow-on mission to MERIS, the Ocean and Land Color Imager (OLCI) on the Sentinel 3A spacecraft that is now operated by the European Organization for the Exploration of Meteorological Satellites (EUMETSAT). A second OLCI sensor was launched on the Sentinel 3B satellite in 2018 ( Figure 1). The MODIS sensor is deployed on two spacecraft: Terra (launched December 1999) and Aqua (launched May 2002). Both MODIS sensors are well past their designed mission life and could become inoperable at any time. The OLCI and MERIS sensors are less noisy and have a higher spatial resolution (for MERIS, high-resolution images were only available when specifically recorded or captured by direct broadcast [16]) and would generally be preferred to MODIS for monitoring cyanobacterial blooms [19]. As a result, it is advantageous to develop a way to calibrate the CI algorithm among the three sensors to create a seamless time-series blending the MODIS, MERIS, and OLCI sensors. Having a method for comparing the CI values between sensors is imperative as it could help direct monitoring of cyanobacterial blooms over decadal timescales and address the efficacy of management decisions designed to decrease the detrimental impacts of cyanobacterial blooms. Furthermore, having a technique to compare algorithms across remotely sensed platforms satisfactorily could be beneficial to intercalibrating other biogeophysical algorithms between sensors. Having interoperability between sensors allows for increased observational frequency, as well as an increase in spatial coverage. This will assist in overcoming data gaps due to clouds, orbital patterns, or sunglint. it could help direct monitoring of cyanobacterial blooms over decadal timescales and address the efficacy of management decisions designed to decrease the detrimental impacts of cyanobacterial blooms. Furthermore, having a technique to compare algorithms across remotely sensed platforms satisfactorily could be beneficial to intercalibrating other biogeophysical algorithms between sensors. Having interoperability between sensors allows for increased observational frequency, as well as an increase in spatial coverage. This will assist in overcoming data gaps due to clouds, orbital patterns, or sunglint. Figure 1. The time period when each sensor was actively collecting data. An arrow indicates that the sensor was still collecting data at the time of this writing.
The process described here has three steps: first, to intercalibrate the CI product between MERIS and MODIS; second, to intercalibrate MODIS-to-OLCI CI; third, to use MODIS as a bridge to establish a relationship between OLCI-to-MERIS CI ( Figure 1).

Characteristics of the MODIS, OLCI, and MERIS Sensors
The characteristics of the three sensors that are relevant to this study are described in Table 1. The OLCI, MERIS, and MODIS Terra sensors all have overpass times of approximately 10:30 a.m. local time. The MODIS Aqua sensor has an overflight time in the afternoon, approximately 3 h later. This time lag is sufficiently long to allow bloom locations to shift, introducing a significant bias when scenes from this sensor are compared to comparable ones collected during the earlier OLCI and MERIS overpasses. To ensure the greatest similarity in data collection in space and time [7], only the MODIS Terra data will be used in this study for intercalibrating the different satellite sensors. Consequently, any reference to the MODIS sensor onboard the Terra satellite will be referred hereafter as MODIST. The few references the MODIS instrument on the Aqua satellite will be denoted as MODIS-Aqua.
MODIST has an exact orbital repeat time of 16 days on a 705 km high orbit in descending node. The MERIS sensor has a repeat time of 35 days, which allows global coverage in ~3 days at an orbit height of 800 km. The OLCI sensor has a repeat time of 27 days at an orbital height of 815 km. OLCI has a spatial resolution of 300 m, as did MERIS when data was collected in the Full Resolution (FR) mode. FR collection was not routine over the USA and Canada until 2008 when the Canadian Center for Remote Sensing began acquiring direct broadcast data. The MERIS data set was routinely resampled onboard to a 1200m spatial resolution in the Reduced Resolution (RR) mode. The RR data was more reliably archived and will be used here, as the FR data was not reliably archived over much of North America. While the algorithm has been applied successfully to all three sensors (OLCI, MODIS, and MERIS) in a variety of water bodies, the validation of the algorithm Figure 1. The time period when each sensor was actively collecting data. An arrow indicates that the sensor was still collecting data at the time of this writing.
The process described here has three steps: first, to intercalibrate the CI product between MERIS and MODIS; second, to intercalibrate MODIS-to-OLCI CI; third, to use MODIS as a bridge to establish a relationship between OLCI-to-MERIS CI ( Figure 1).

Characteristics of the MODIS, OLCI, and MERIS Sensors
The characteristics of the three sensors that are relevant to this study are described in Table 1. The OLCI, MERIS, and MODIS Terra sensors all have overpass times of approximately 10:30 a.m. local time. The MODIS Aqua sensor has an overflight time in the afternoon, approximately 3 h later. This time lag is sufficiently long to allow bloom locations to shift, introducing a significant bias when scenes from this sensor are compared to comparable ones collected during the earlier OLCI and MERIS overpasses. To ensure the greatest similarity in data collection in space and time [7], only the MODIS Terra data will be used in this study for intercalibrating the different satellite sensors. Consequently, any reference to the MODIS sensor onboard the Terra satellite will be referred hereafter as MODIS T . The few references the MODIS instrument on the Aqua satellite will be denoted as MODIS-Aqua. MODIS T has an exact orbital repeat time of 16 days on a 705 km high orbit in descending node. The MERIS sensor has a repeat time of 35 days, which allows global coverage in~3 days at an orbit height of 800 km. The OLCI sensor has a repeat time of 27 days at an orbital height of 815 km. OLCI has a spatial resolution of 300 m, as did MERIS when data was collected in the Full Resolution (FR) mode. FR collection was not routine over the USA and Canada until 2008 when the Canadian Center for Remote Sensing began acquiring direct broadcast data. The MERIS data set was routinely resampled onboard to a 1200-m spatial resolution in the Reduced Resolution (RR) mode. The RR data was more reliably archived and will be used here, as the FR data was not reliably archived over much of North America. While the algorithm has been applied successfully to all three sensors (OLCI, MODIS, and MERIS) in a variety of water bodies, the validation of the algorithm is not the focus of this manuscript. For more detailed results of the algorithmic performance with these sensors, the reader is directed to Wynne et al. [19] for a detailed Remote Sens. 2021, 13, 2305 4 of 19 look into MODIS, Wynne et al. [18] for validation of MERIS, and Mishra et al. [16] for application to OLCI.

Study Area
The Laurentian Great Lakes lie between the border of Canada and the United States and feature three catchments that are routinely affected by cyanobacterial blooms: western Lake Erie, Saginaw Bay, and Green Bay [20] (Figure 2). These are large water bodies that have been routinely experienced cyanobacterial blooms over the last 20 years and will be the regions of interest for this study. These lakes are used by millions of people to supply drinking water, and cyanobacterial blooms can be detrimental to human health. For example, in 2014, the metropolitan area of Toledo, OH (USA), issued a "Do Not Drink" order on its municipal water supply due to contamination from the biotoxin, microcystin, caused by cyanobacterial blooms [21]. Green Bay and Lake Erie are separated by over 500 km and do not necessarily appear in the same swath on MERIS (1150 km) or OLCI (1270 km). The CI algorithm has been routinely run on all three basins to monitor and detect cyanobacterial blooms [22]. Western Lake Erie has the largest blooms of the three catchments, and these blooms exhibit a large degree of interannual variability in size and severity (interannual variability is >20 fold) [19]. Saginaw Bay has intermediate-sized blooms with a small degree of interannual variability [19]. Green Bay exhibits the smallest blooms of the three basins, with a high degree of interannual variability in its cyanobacterial biomass [20]. These three regions will be used to make the three CI intercalibration factors (MERIS-to-MODIS T , MODIS T -to-OLCI, and MERIS-to-OLCI). Using three distinct basins in three different lakes ( Figure 2) will yield a more robust dataset as opposed to using only one single basin. All three basins are shallow, relatively warm, eutrophic environments that have the preponderance of their nutrients delivered by a single river. The Fox River supplies 60% of nutrients into Green Bay [23], the Saginaw River supplies approximately 90% of the nutrients into Saginaw Bay [24], and the Maumee River, along with the smaller Cuyahoga and Sandusky Rivers, supplies approximately half of the nutrients into western Lake Erie [25] (Figure 2).

Algorithm
The cyanobacteria index (CI) is the algorithm that will be evaluated here. It has been used for over a decade to monitor for cyanobacteria in freshwater bodies. The CI algorithm [17,18] is defined as: where SS is the spectral shape around the 2 band (about 680 nm) that corresponds to the strong chlorophyll-a absorption peak ( Table 1). The SS takes the form of a second derivative and is sometimes called a baseline algorithm: where ρ is the dimensionless Rayleigh-corrected surface reflectance determined from the instrument-observed top-of-atmosphere radiance by only removing Rayleigh radiances and molecular absorption transmission loss, both corrected for elevation. See Table 1 for band information for each of the sensors such as ρ, and CI is dimensionless. MERIS-reduced resolution L1B, OLCI-reduced resolution L1B, and Terra L0 imagery were downloaded and processed with NOAA's Satellite Automated Processing System (SAPS), with methods presented in Wynne et al. [26]. SAPS used the NASA l2gen "rho_s" option to generate the Rayleigh-corrected reflectance, ρ, used in Equation (2). These were mapped, using nearest neighbor, to the same Albers Equal Area protection with a pixel size of 1100 m [26] ( Table 1). The results allow consistent comparison of the three sensors.

How Image Pairs Were Selected
For calibration, high-quality image pairs with well-defined cyanobacterial blooms are needed to avoid any artifacts or biases resulting from noise such as sunglint, thin clouds, haze, etc. Image pairs were selected using a sequence of steps. Dates with sameday pairs during blooms were identified. The cyanobacterial bloom season in the Laurentian Great Lakes lasts from 1 June to 31 October and is generally highest in July-October [27]. All images shown here were chosen from the July-October time period from 2002

Algorithm
The cyanobacteria index (CI) is the algorithm that will be evaluated here. It has been used for over a decade to monitor for cyanobacteria in freshwater bodies. The CI algorithm [17,18] is defined as: where SS is the spectral shape around the λ 2 band (about 680 nm) that corresponds to the strong chlorophyll-a absorption peak ( Table 1). The SS takes the form of a second derivative and is sometimes called a baseline algorithm: where ρ is the dimensionless Rayleigh-corrected surface reflectance determined from the instrument-observed top-of-atmosphere radiance by only removing Rayleigh radiances and molecular absorption transmission loss, both corrected for elevation. See Table 1 for band information for each of the sensors such as ρ, and CI is dimensionless. MERIS-reduced resolution L1B, OLCI-reduced resolution L1B, and Terra L0 imagery were downloaded and processed with NOAA's Satellite Automated Processing System (SAPS), with methods presented in Wynne et al. [26]. SAPS used the NASA l2gen "rho_s" option to generate the Rayleigh-corrected reflectance, ρ, used in Equation (2). These were mapped, using nearest neighbor, to the same Albers Equal Area protection with a pixel size of 1100 m [26] ( Table 1). The results allow consistent comparison of the three sensors.

Image processing 2.2.1. How Image Pairs Were Selected
For calibration, high-quality image pairs with well-defined cyanobacterial blooms are needed to avoid any artifacts or biases resulting from noise such as sunglint, thin clouds, haze, etc. Image pairs were selected using a sequence of steps. Dates with same-day pairs during blooms were identified. The cyanobacterial bloom season in the Laurentian Great Lakes lasts from 1 June to 31 October and is generally highest in July-October [27]. All images shown here were chosen from the July If one of the images was substantially off-nadir or experienced adjacency effects [28,29], the image pair was rejected. Each image pair was also rejected if the wind speed during acquisition time from either sensor exceeded 8 m s −1 to avoid perceived changes in surface cyanobacterial biomass due to vertical mixing [18]. The remaining image pairs were screened by visually inspecting the RGB true color composites (which are all processed with the same enhancement). Evidence of extensive clouds, cloud shadows, or sun glint in the bloom area in either image led to a rejection of that image pair. Each image pair used for the MODIS T -to-MERIS intercalibration is listed in Supplementary Table S1 (n = 42), and each image pair used for the MODIS T -to-OLCI intercalibration is listed in Supplementary Table S2 (n = 19). Example image pairs that were rejected based on a visual inspection are described in Supplementary Table S3.
The time differences of the final image pairs were considered acceptable. The largest time difference in the overflight times between the two different spacecraft was 78 min ( Figure 3). The average time difference between MERIS and MODIS T 's data acquisition was 26 min, with MODIS T having an earlier overpass 19 times and MERIS having an earlier overpass 23 times. The OLCI and MODIS T image pairs were similar with an average time difference between image acquisition of 25 min, with MODIS T having an earlier overpass 8 times and OLCI having an earlier overpass 11 times.
ble S1 (n = 42), and each image pair used for the MODIST-to-OLCI intercalibration is listed in Supplementary Table S2 (n = 19). Example image pairs that were rejected based on a visual inspection are described in Supplementary Table S3.
The time differences of the final image pairs were considered acceptable. The largest time difference in the overflight times between the two different spacecraft was 78 min ( Figure 3). The average time difference between MERIS and MODIST's data acquisition was 26 min, with MODIST having an earlier overpass 19 times and MERIS having an earlier overpass 23 times. The OLCI and MODIST image pairs were similar with an average time difference between image acquisition of 25 min, with MODIST having an earlier overpass 8 times and OLCI having an earlier overpass 11 times.

Image Quality Assurance
Radiance signal from the adjacent area surrounding a pixel (the object pixel) can propagate into the sensor field-of-view due to subsequent scattering and cause perturbations in the spectral signature at the pixel in question. This process is usually termed as adjacency effect [29]. In open ocean waters, the adjacency effect may not be an issue due to the spatial homogeneity. However, in coastal waters, it can affect the radiometric quality of the image data. In this study, adjacency issues affecting pixel quality were dealt with as follows. Initially, a 3 × 3 pixel neighborhood surrounding each pixel in the image was identified. Next, if at least one invalid pixel was identified within the 3 × 3 pixel neighborhood as being invalid due to clouds, sun glint, adjacency to land, etc., that pixel was flagged as unusable. When any two images were subsequently being compared as part of the satellite sensor intercalibration analysis, only the pixel pairs that did not suffer for adjacency effects were used. This assured only the highest quality data set was used in the subsequent analyses.

Image Quality Assurance
Radiance signal from the adjacent area surrounding a pixel (the object pixel) can propagate into the sensor field-of-view due to subsequent scattering and cause perturbations in the spectral signature at the pixel in question. This process is usually termed as adjacency effect [29]. In open ocean waters, the adjacency effect may not be an issue due to the spatial homogeneity. However, in coastal waters, it can affect the radiometric quality of the image data. In this study, adjacency issues affecting pixel quality were dealt with as follows. Initially, a 3 × 3 pixel neighborhood surrounding each pixel in the image was identified. Next, if at least one invalid pixel was identified within the 3 × 3 pixel neighborhood as being invalid due to clouds, sun glint, adjacency to land, etc., that pixel was flagged as unusable. When any two images were subsequently being compared as part of the satellite sensor intercalibration analysis, only the pixel pairs that did not suffer for adjacency effects were used. This assured only the highest quality data set was used in the subsequent analyses.

Pixel to Pixel Technique
The most common technique used for satellite intercalibration is the pixel-to-pixel technique. In this method, the value of corresponding pixels in regions of interest from Remote Sens. 2021, 13, 2305 7 of 19 paired images collected by two different sensors are extracted. In this study, CI values from paired pixels were extracted from the images (MODIS T /MERIS or MODIS T /OCLI), and least-squares regression analyses were performed. The regression line was forced through the origin, and if this was not so, the low values on the low end might become less than the slope intercept, rendering them meaningless. The intercalibration correction factor (assuming that the relationship of the data pairs passes through the origin) was the slope of the best fit linear regression line of the CI values for each sensor pair. This allowed determining intercalibration factors, which allowed the CI values between sensors to be interconverted [10,11,30]. Henceforth, this approach will be referred to as the pixel technique.

Integrated Pixel Technique
An alternative to cross-validating satellite imagery on a pixel-to-pixel basis is instead to compare the integrated pixel values across the same regions from paired MERIS, OLCI, or MODIS T images. This approach involves first identifying within paired satellite images all pairs of pixels that both have a positive CI value and have no image quality issues in either image (Section 2.2.2). These pixels were then summed to produce an aggregate, or integrated, CI value for each image. This is not an averaging but, instead, a sum of all of the CI positive values from each image. These pairs of summed CI values were then used as the basis for least-square regression analyses to determine the CI relationships between the MODIS T , MERIS, and OLCI sensors. As a result, each image pair produces a single pair of values, unlike the pixel technique, which may produce many pairs of values for each pair of images. Henceforth, this approach will be referred to as the integrated technique. As in the pixel technique, the regression line was forced through the origin, and if this was not so, the low values on the low end might become less than the slope intercept, rendering them meaningless.
The integrated technique has a major advantage over the pixel technique. It inherently reduces one key source of error, the mismatch from the spatiotemporal change in bloom distribution between the two sensor overpass times. Cyanobacterial blooms often are patchy with considerable spatial heterogeneity. The dynamic nature of the coastal environment can cause significant changes between overpasses due to shifts in these patches. These changes can cause mismatches in pixel-to-pixel inter-comparisons. However, an integrated technique would be expected to reduce the error from the spatiotemporal variability of the bloom by aggregating across the patchiness in the CI field.

Model Evaluation
Leave-one-out cross-validation, often used in data-scarce scenarios, is an n-fold crossvalidation approach where n is the total number of samples in the dataset. In each fold, the process leaves one sample out for model validation and uses the remaining samples for model calibration, making it similar to a jack-knife estimation [31]. Thus, in each iteration, the model is calibrated with (n−1) data points and validated on the left-out single data point, allowing a validation for each data point. Similar to that idea, here, we followed the leave-one-lake-out validation approach for cross-validation. Least-square regression models for sensor intercalibration were adjusted, using data from two study regions (e.g., Western Lake Erie and Saginaw Bay) and validated using the third region (e.g., Green Bay) as an independent dataset, giving a "leave-one-lake-out" approach. The model (M) was then applied to the appropriate sensor (typically MODIS T ) for the third lake, and the MERIS or OLCI values for that lake were used as the observed (O) reference data. Thus, the models were calibrated three times, using data from pairs of lakes in combination, as well as validated three times. This "leave-one-lake-out" approach demonstrated the robustness of the regression coefficients or correction factors with a geographically independent data set. Additionally, instead of a single slope, this validation approach provided an uncertainty range around the mean slopes. The arithmetic mean of the regression coefficients from the three sets of calibration coefficients relating the MERIS-to-MODIS T CI and OLCI-to-MODIS T CI were finally used to convert OLCI CI observations to MERIS CI-equivalent observations, thereby connecting the MERIS and OLCI time-series. Calibrating for combinations of each of the three regions separately allows an assessment of calibration consistency and validation with the third region.

Error Metrics
We used bias and Mean Absolute Error (MAE) as error metrics [32] for evaluating the calibrations. Bias quantifies the systematic error or the difference between the modeled value (M) (in this case, the MERIS/OLCI CI modeled from MODIS T CI) and the observed or reference value (O) (actual MERIS/OLCI CI). As the analysis used log-transformed data, the mean multiplicative bias was determined: where M is the modeled CI; O is the observed; n is the sample size (number of pixels for the pixel technique or number of image pairs for the integrated image technique). The closer the multiplicative bias is to a value of 1, the less bias the comparisons have and the better the model accuracy. For example, a bias of 1.2 indicates a model that overestimates by about 20%, and a mean bias of 0.8 indicates a model that underestimates by about 20% (there is not an exact match between multiplicative error and percentage error [32]. The MAE captures the error magnitude, highlighting the absolute error between the observed and modeled quantities. The MAE is defined as: Equation (4) uses the same variables as Equation (3). The MAE in Equation (4) is multiplicative also. An MAE value of 1.2 indicates that the error is about 1.2× the observed value. Multiplicative error metrics are most appropriate for distributions that have errors that are proportional to magnitude.

Pixel Technique MERIS and MODIS T
The pixel-to-pixel matchups for all MERIS-MODIS T CI pairs in Supplementary Table S1 produced a slope of 2.76 with an R 2 of 0.82 (Figure 4 (left)). The residuals (or difference between pairs) in Figure 4 (right) suggests a skewed error distribution, with a high density of pixels below the zero residual line. The correction factor using the pixel technique to convert from MODIS T to MERIS is The value of 2.76 (vs. unity) is expected, as the CI algorithm uses different wavelengths for the λ 3 term in Equation (2): 748 nm band for MODIS T vs. 709 nm band for both MERIS and OLCI. Because water absorbs about 3-fold more at 748 nm than 709 nm, the ρ (water reflectance) at 709 nm will be about 3× that for ρ at 748 nm, resulting in a similarly higher value for the MERIS CI compared to that for MODIS T . The error metrics are presented in Table 2.

Integrated Technique MERIS and MODIST
The integrated technique uses the same 42 image pairs from Section 3.1.1 (see Supplementary Table S2). The intercalibration correlation between the two returned a slope of 2.92 (compared to 2.76 for the pixel technique) and an R 2 of 0.97 ( Figure 5). Therefore, the correction factor to convert from MODIST to MERIS using the integrated technique is: As noted before, the value of 2.92 is consistent with the difference in the λ3 term, with error metrics in Table 2. Both the integrated technique and the pixel technique had a negligible mean bias (1.009 and 1.003, respectively). However, the pixel technique had a much larger error, with an MAE 44% greater than the integrated technique. The technique also generated a slightly smaller slope (5.5%) than the integrated technique (2.76 vs. 2.92). The integrated technique also showed a uniform distribution of values close to the regression line, indicating linearity and consistency across the range. The difference between the two techniques lies primarily in the challenge of using the pixel technique to successfully match the pixels because discrepancies caused by any bloom movement or differences in ground field-of-view between sensors can manifest as errors due to local biases. Supplementary Table S3 shows examples of image pairs that were not selected as they were deemed flawed, along with the rationale for rejection.  Table S2). The intercalibration correlation between the two returned a slope of 2.92 (compared to 2.76 for the pixel technique) and an R 2 of 0.97 ( Figure 5). Therefore, the correction factor to convert from MODIS T to MERIS using the integrated technique is: As noted before, the value of 2.92 is consistent with the difference in the λ 3 term, with error metrics in Table 2. Both the integrated technique and the pixel technique had a negligible mean bias (1.009 and 1.003, respectively). However, the pixel technique had a much larger error, with an MAE 44% greater than the integrated technique. The technique also generated a slightly smaller slope (5.5%) than the integrated technique (2.76 vs. 2.92). The integrated technique also showed a uniform distribution of values close to the regression line, indicating linearity and consistency across the range. The difference between the two techniques lies primarily in the challenge of using the pixel technique to successfully match the pixels because discrepancies caused by any bloom movement or differences in ground field-of-view between sensors can manifest as errors due to local biases. Supplementary Table S3 shows examples of image pairs that were not selected as they were deemed flawed, along with the rationale for rejection.  Supplementary  Table S1. The symbol shapes represent the three different basins used with triangles from western Lake Erie, circles from Green Bay, and squares from Saginaw Bay.

MERIS and MODIST Statistics on Intercalibration
The calibration approach using "leave-one-lake-out" gave quite similar results between regions (Table 3, Figure 6A-C). Figure 6D-F shows the observed MERIS CI to modeled MERIS CI (using MODIST CI) relationships for each of the three regions. All regions behave similarly, with Green Bay having a slightly higher bias relative to western Lake Erie and Saginaw Bay and western Lake Erie has the lowest MAE of the three regions. In the case of this intercalibration, the validation data was outside the calibration range; however, the intercalibration still performed sufficiently, showing linearity and consistency in slope. This consistency indicated extrapolation outside the calibration range is valid. Table 3. Least-square fit parameters from MERIS-MODIST intercalibration using different combinations of the regions as seen in Figure 6A-C. WLE, SB, and GB, respectively, represent Western Lake Erie, Green Bay, and Saginaw Bay. Percentages in the slopes indicate the difference of that pair from a calibration based on all regions.   Supplementary  Table S1. The symbol shapes represent the three different basins used with triangles from western Lake Erie, circles from Green Bay, and squares from Saginaw Bay.

MERIS and MODIS T Statistics on Intercalibration
The calibration approach using "leave-one-lake-out" gave quite similar results between regions (Table 3, Figure 6A-C). Figure 6D-F shows the observed MERIS CI to modeled MERIS CI (using MODIS T CI) relationships for each of the three regions. All regions behave similarly, with Green Bay having a slightly higher bias relative to western Lake Erie and Saginaw Bay and western Lake Erie has the lowest MAE of the three regions. In the case of this intercalibration, the validation data was outside the calibration range; however, the intercalibration still performed sufficiently, showing linearity and consistency in slope. This consistency indicated extrapolation outside the calibration range is valid. Table 3. Least-square fit parameters from MERIS-MODIS T intercalibration using different combinations of the regions as seen in Figure 6A-C. WLE, SB, and GB, respectively, represent Western Lake Erie, Green Bay, and Saginaw Bay. Percentages in the slopes indicate the difference of that pair from a calibration based on all regions.

Pixel Technique OLCI to MODIST
The pixel-to-pixel intercalibration between OLCI and MODIST image pairs (Supplemental Table S2) produced a slope of 2.75 and an R 2 of 0.86 ( Figure 7). As was the case in the MERIS-to-MODIST analysis, the residuals show a skewed error distribution, with a high density of pixels in the mid-range (0.00-0.001), showing a negative bias. The relationship to go from MODIST to OLCI using the pixel technique is seen in Equation (7).

Pixel Technique OLCI to MODIS T
The pixel-to-pixel intercalibration between OLCI and MODIS T image pairs (Supplemental Table S2) produced a slope of 2.75 and an R 2 of 0.86 ( Figure 7). As was the case in the MERIS-to-MODIS T analysis, the residuals show a skewed error distribution, with a high density of pixels in the mid-range (0.00-0.001), showing a negative bias. The relationship to go from MODIS T to OLCI using the pixel technique is seen in Equation (7).
Unlike the MERIS-MODIS T intercalibration, where both techniques showed negligible bias, the integrated technique for OLCI to MODIS T led to a significantly smaller bias (dowñ 24% to 0.94). Likewise, the MAE showed~29% less error using the integrated technique relative to the pixel technique.
Unlike the MERIS-MODIST intercalibration, where both techniques showed negligible bias, the integrated technique for OLCI to MODIST led to a significantly smaller bias (down ~24% to 0.94). Likewise, the MAE showed ~29% less error using the integrated technique relative to the pixel technique.  Table S2). The right panel shows the residual errors around the calibration shown in the left panel.

Integrated Technique OLCI to MODIST
The intercalibration using the integrated technique (Figure 8) uses the same image pairs as the pixel technique (Supplemental Table S2). The R 2 increases from 0.86 to 0.98 relative to the pixel technique, and the slope decreases slightly from 2.75 to 2.71. The equation to go from MODIST to OLCI using the integrated technique is shown in Equation (8).
No residual bias around the regression line is evident, unlike for the pixel-to-pixel intercalibration. Most values fall close to the regression line, and those at the lower end are equally distributed above and below the line, with two points well above (points 15 and 11) and three below (points 8, 17, and 12) (Figure 8; Supplementary Tables S1 and S2).   Table S2). The right panel shows the residual errors around the calibration shown in the left panel.

Integrated Technique OLCI to MODIS T
The intercalibration using the integrated technique (Figure 8) uses the same image pairs as the pixel technique (Supplemental Table S2). The R 2 increases from 0.86 to 0.98 relative to the pixel technique, and the slope decreases slightly from 2.75 to 2.71. The equation to go from MODIS T to OLCI using the integrated technique is shown in Equation (8).
Unlike the MERIS-MODIST intercalibration, where both techniques showed negligible bias, the integrated technique for OLCI to MODIST led to a significantly smaller bias (down ~24% to 0.94). Likewise, the MAE showed ~29% less error using the integrated technique relative to the pixel technique.  Table S2). The right panel shows the residual errors around the calibration shown in the left panel.

Integrated Technique OLCI to MODIST
The intercalibration using the integrated technique (Figure 8) uses the same image pairs as the pixel technique (Supplemental Table S2). The R 2 increases from 0.86 to 0.98 relative to the pixel technique, and the slope decreases slightly from 2.75 to 2.71. The equation to go from MODIST to OLCI using the integrated technique is shown in Equation (8).
No residual bias around the regression line is evident, unlike for the pixel-to-pixel intercalibration. Most values fall close to the regression line, and those at the lower end are equally distributed above and below the line, with two points well above (points 15 and 11) and three below (points 8, 17, and 12) (Figure 8; Supplementary Tables S1 and S2).   Supplementary Table S2. The symbol shapes represent the three different basins used with triangles from western Lake Erie, circles from Green Bay, and squares from Saginaw Bay.
No residual bias around the regression line is evident, unlike for the pixel-to-pixel intercalibration. Most values fall close to the regression line, and those at the lower end are equally distributed above and below the line, with two points well above (points 15 and 11) and three below (points 8, 17, and 12) (Figure 8; Supplementary Tables S1 and S2).
While both the pixel technique and the integrated technique produced equivalent slopes for OLCI (2.75 vs 2.71, respectively) (Table 4), the subsequent error metrics (bias and MAE) were much different. For OLCI, the integrated technique had a slight negative bias (0.94), while the pixel technique had a strong negative bias (0.77). This differs from the MODIS T -MERIS, where both techniques had a negligible bias. The pixel technique also had a much larger error (1.76) than the integrated technique (1.25), just as was seen with MERIS. In this case, the pixel technique produced substantially large errors that may lead to lower confidence in the resulting calibration.

MERIS and MODIS T Statistics on Intercalibration
Figure 9A-C shows the calibration fit parameters for the OLCI-MODIS T intercalibration for each combination of the three regions, with results summarized in Table 5. Figure 9D-F shows the observed OLCI to modeled OLCI (using MODIS T ) relationships for each of the three regions. There is more discrepancy within the regions using the OLCI to MODIS T intercalibration among the three geographical regions than there was between the MERIS-MODIS T intercalibration. Likewise, there was much more variability geographically in the MAE, with western Lake Erie and Saginaw Bay having quite similar values (>1% change) but with Green Bay having a much lower error with MAE some 12% lower. Table 5. Least-square fit parameters from the OLCI-MODIS T CI intercalibration as seen in Figure 9A-C based on the integrated technique. Percentages in the slopes indicate the difference of that pair from a calibration based on all regions.

MERIS-OLCI Intercalibration
The final intercalibration coefficients estimated from the leave-one-lake-out approach and based on the integrated technique are summarized in Table 6. These coefficients are the slopes of the regression lines between all three sensors for all three regions. To compare OLCI to MERIS, the relationships between MODIS T and MERIS/OLCI were used as an intermediary using the cross-regional averages ( Table 6). Using the calibration coefficients from Table 6 yields the following set of equations (Equations (9) and (10) represent the updated Equations (6) and (8), respectively, after replacing the slopes with the cross-regional average slopes.): OLCI integrated = 2.72 × MODIS T,integrated (10) MODIS T,integrated = 0.36 × OLCI integrated (11) Currently, the MERIS calibration is more stable than the OLCI calibration, likely because there has been sufficient time to properly calibrate the sensor, therefore converting OLCI to MERIS is the more appropriate choice as compared to converting MERIS to OLCI. Substituting the CI MODIST with the CI OLCI from Equation (11) into Equation (9) yields: Solving Equation (12) (11) Currently, the MERIS calibration is more stable than the OLCI calibration, likely because there has been sufficient time to properly calibrate the sensor, therefore converting OLCI to MERIS is the more appropriate choice as compared to converting MERIS to OLCI. Substituting the CIMODIS T with the CIOLCI from Equation (11) into Equation (9) yields: Solving Equation (12) will yield the final relationship between the CIMERIS and the CIOLCI expressed as:

Discussion
The need to merge cyanobacterial index (CI) data from multiple sensors for environmental analyses will become increasingly common [33]. This will be necessary to determine the efficacy of management strategies to control cyanobacterial bloom intensity, such as nutrient reductions on decadal time scales. If agricultural best management practices are amended, the need to synoptically differentiate any impacts through the changing of these practices must be taken into account [34]. The timing, severity, location, and size of cyanobacterial blooms can be used to test the efficacy of these changes. Further, if there are differences in the physical manifestation of the bloom (i.e., timing, severity, location, and size), this could be a harbinger for an ecological shift within a given ecosystem. This study examined both the traditional pixel-to-pixel technique and a new image-to-image integrated technique for intercalibrating the CI algorithm from the MODIS T , MERIS, and OLCI sensors. Being able to reliably compare the CI values between sensors will allow the production of a >20 year CI time-series in numerous geographically dispersed water bodies. These time-series will prove useful in retrospectively evaluating previous management decisions and in documenting temporal changes in cyanobacterial bloom intensity.
In this study, the image-to-image integrated technique exhibited much lower error and reduced bias than the pixel-to-pixel technique. The residual CI distributions produced by the pixel technique showed skewed bias in the residuals ( Figures 4B and 7B), which did not occur with the integrated technique ( Figures 5 and 8). Additionally, as the MODIS T CI increased, the more negative the MERIS or OLCI residuals became when using the pixel technique ( Figures 4B and 7B). The integrated technique did not exhibit such biases over the entire range of CI values, with the higher values converging on the 1:1 line (Figures 5 and 8). Pixel-to-pixel comparisons would be expected to have this problem because any shifts in the distribution of the variable being measured between overflight times or alignment issues can introduce noise confounding intercalibration ( Figures 4B and 7B for the pixel technique; 5 and 8 for the pixel technique for integrated technique). The integrated technique, in contrast, provides a more robust analysis option for evaluating the influence of individual image pairs on the calibration (e.g., Figure 6A-C and Figure 9A-C). In a pixel-topixel technique, the statistical influence of the pixels from individual scenes cannot be readily evaluated and would require more complex analysis. Any skewing of pixels from CI image pairs at the low or high end (those pairs with maximum statistical leverage) would bias the results of the comparison. Skewing could occur with even a one or two-pixel offset in the actual (or apparent) location of an intense bloom between two images. This problem is mostly removed by the integration technique, as the offset pixels are contained within the integrated region of the bloom.
The calibration coefficient between the MODIS T and MERIS for the pixel and integrated techniques differed by 5% (Table 4), while the coefficient between the pixel and integrated technique for the MODIS T -OLCI matchup differed by 1%. However, as noted above, the pixel technique for MODIS T to OLCI had a severe bias, which can lead to questions on the accuracy of the calibration. The integrated technique did not have this problem. OLCI appeared to miss the lowest CI values (Figure 7) that MODIS T and MERIS detect. MERIS has had much more time for calibration, and it is now at the fourth reprocessing. As the OLCI calibration is developed and improves, it is likely to see shifts in the calibration, potentially reducing the difference. Most calibration is conducted on a band-by-band basis.
Subtle inter-band differences within the uncertainties of the band calibration can affect a shape algorithm such as the CI [35]. A spectral shape method can detect these residual differences, providing a means of improving the inter-band calibration [35]. The integrated technique can provide a more reliable method for intercalibrating other biophysical variables across sensors as well. A similar approach has been used to calibrate MODIS [36,37]. Specifically, data sets were integrated over space and time to compare satellites for such issues as trends and calibrating for cross-scan polarization.
We carried out sensitivity analysis using bootstrap simulation to find out how stable the slope parameters were with respect to the effect from the points of leverage or influence. We repeated the simulation for 1000 runs for MERIS-MODIS analysis. OLCI-MODIS sensitivity analysis was repeated 171 times based on the number of unique permutations possible with a sample size of 17. At every step of the iteration, 39 out of 42 of the MERIS-MODIS samples and 17 out of 19 samples were randomly selected without replacement. By doing so, 5% and 10% of MERIS-MODIS and OLCI-MODIS pairs were left out to evaluate their influence on the regression slope. The median slope from the sensitivity analysis was identical to the reported slopes in Figures 5 and 7 (Table 7). The variability of the MERIS-OLCI slope is expected to be within 2.87-2.99, which is −1.7% and +2.4% of the median value of 2.92. Similarly, the OLCI-MODIS slope is likely to be within 2.61-2.81 (−3.7% or +3.5% around the median of 2.71). While more samples may lead to more robust statistics, special attention was given to image quality, and only the very best image pairs were selected for analysis. Table 7. Descriptive statistics of regression slopes from the bootstrap simulation. This study shows the CI intercalibration is consistent between the three basins in the Laurentian Great Lakes and that the integration technique is robust. The intercalibration coefficients presented here were consistent with no outliers compared to those derived with the pixel-to-pixel technique used in Wynne et al. [10]. We should note that the CI algorithm does not require an atmospheric correction [26]. It was based on the Maximum Chlorophyll Index (MCI) presented by Gower et al. [38], who used top-of-atmosphere reflectance. The calibration coefficient between MODIS and MERIS derived here was much higher than in Wynne et al. [10], primarily due to an issue with the measurement units. The original MERIS reflectance data sets were expressed as L2 files from the original ESA processing as having units of per steradian (sr −1 ) (although presented as dimensionless). This translates as having the initial L2 products being off by a factor of π sr. By multiplying the 1.33 correction factor derived by Wynne et al. [10] by π, the new factor becomes 4.17, which is considerably closer to the correction factor of 2.94 derived here. The remaining 30 % difference (between 4.17 and 2.94) could be explained by several reasons. ESA has recalibrated and reprocessed the MERIS data several times since the analysis of Wynne [10]. This study also used rigorous quality assurance for both the image pairs selected and the pixels used for analysis. Further, Wynne et al. [10] used only MODIS Aqua imagery, which is collected about 3 h later than MERIS and OLCI. By employing MODIS Terra imagery in this study, which has comparable overpass times to MERIS and OLCI, it was possible to reduce any impacts of wind events or biomass changes in surface waters unaccounted for in the original calibration study. This, in turn, reduced misfit, thereby improving the calibrations.

OLCI-MODIS T Slope Variability
While this study focuses on intercalibrating the CI algorithm among three large satellite missions in the Great Lakes, the techniques presented here could be adapted to other areas and sensors. The integrated technique could be applied to any remotely sensed image products from any number of sensors and algorithms. The integrated technique is not designed to better approximate an in situ measurement but instead to determine how well an algorithm from one sensor approximates the same algorithm on a different sensor. As a result, the integrated technique could also be applied to small satellites, such as CubeSats, which generally have poor absolute radiometric calibration, but show linearity. Two excellent review articles highlight water quality parameters that could be tested using the integrated technique for sensor intercalibration [39,40].

Conclusions
Using an integrated technique for intercalibration provides better metrics and more consistent results for calibrating an algorithm between the sensors than the pixel technique. The more typical pixel-to-pixel technique did produce equivalent calibration coefficients (slopes) for MODIS T to OLCI but with greater uncertainty. The comparison of three separate lakes separated by up to 500 km showed a robust and consistent result. This integrated technique allows us to calibrate across the three sensors, MERIS, OLCI, MODIS T , and by inference, also MODIS-Aqua (as NASA calibrates Terra to Aqua). The results further allow us to adjust OLCI data sets to match MERIS. For the calibration identified here, OLCI underestimates the CI relative to MERIS by 6% (Equation (13)). The results of this study allow the development of a consistent time-series across MERIS, MODIS, and OLCI for the CI product. The calibration should work in other lakes which are large enough for detection by these sensors. The correction can be updated using the same scenes with changes in calibration. The results will help in the analysis of time-series and in the comparison of MERIS and OLCI for climate changes in cyanobacterial blooms. The integrated technique can be applied to other algorithms and sensors to allow continuity in longer-term datasets.