ORIGINAL RESEARCH article

Front. Remote Sens., 10 October 2025

Sec. Multi- and Hyper-Spectral Imaging

Volume 6 - 2025 | https://doi.org/10.3389/frsen.2025.1670390

Practical aspects of providing pixel-level spectral Rrs error covariance in satellite ocean color products

  • 1. Science Application International Corp., McLean, VA, United States

  • 2. Ocean Ecology Laboratory, Goddard Space Flight Center, National Aeronautics and Space Administration, Greenbelt, MD, United States

  • 3. GESTAR II, University of Maryland Baltimore County, Baltimore, MD, United States

  • 4. Go2Q Pty Ltd., Sunshine Coast, QLD, Australia

Abstract

We previously established a derivative-based approach to generate a pixel-level spectral error covariance matrix in satellite-retrieved remote sensing reflectance, ∑Rrs. However, one practical issue is the delivery of the products without increasing the file size by an order of magnitude or more, considering that for N sensor spectral bands, there are N × (N+1)/2 covariance matrix elements to be specified at each pixel. The issue becomes more pertinent for hyperspectral imaging spectroradiometers such as the Ocean Color Instrument (OCI) on NASA’s Plankton, Aerosol, Cloud, ocean Ecosystem mission (PACE), which has 286 bands, resulting in ∼40,000 unique elements in ∑Rrs per pixel that would lead to a ∼60 GB Level-2 file for one 5-min granule. As a first step to tackle the issue, we took OCI and Moderate Resolution Imaging Spectroradiometer (MODIS) data to explore the possibility of approximating ∑Rrs using a third-degree polynomial, thereby decreasing the memory overhead to 4×N numbers. We found that ∑Rrs derived from the polynomial fitting matches well with the original value, with the difference smaller than 5%. We then compared the relative uncertainty in two derived ocean color data products (chla and Kd(490)) calculated using the original fully computed ∑Rrs and then using the polynomial model approximation for ∑Rrs, finding the absolute difference between the two approaches to be smaller than 0.5%. These evaluations suggest the polynomial approximation of ∑Rrs is suitable without degrading the scientific quality. By including the coefficients derived from polynomial fitting instead of the full error covariance matrix, a typical 5-min Level-2 file for OCI decreases from ∼60 GB to a more practical ∼1.7 GB.

1 Introduction

Due to the capability of obtaining a long time series of data at a global scale, satellite remote sensing has been widely used to study Earth System science and global change (; ). As a product fundamental to deriving a host of bio-optical and biogeochemical properties of the water column, ocean remote sensing reflectance (Rrs, sr−1) has been routinely generated from satellite radiometry and distributed by NASA’s Ocean Biology Distributed Active Archive Center (OB.DAAC) for more than 3 decades. While Rrs and properties derived from it have been validated against ground truth, there is no routinely generated pixel-level uncertainty product. Although the uncertainty in water-leaving reflectance at the pixel level is provided for the Ocean and Land Colour Instrument (OLCI) onboard Sentinel-3 (), it is underestimated, as only sensor noise is included. It should be noted that the European Space Agency’s Ocean Color Climate Change Initiative program (OC_CCI) does provide pixel-level uncertainty in Rrs (), computed using the weighted average of the uncertainty within optical water classes associated with that pixel. The uncertainty of each class is derived from validation against ground truth, which has an issue with representing the pixel-level uncertainty, as described by . Furthermore, this approach could not capture spatiotemporal variations in some key drivers of uncertainty in Rrs, such as geometry and atmospheric turbidity. Without associated uncertainty, ocean color data products are incomplete measurements, limiting the interpretation of subtle changes in global ocean biology and biogeochemistry. Uncertainty in ocean color products has been recommended by the Group on Earth Observations (GEO) and the International Ocean Colour Coordinating Group (IOCCG) within the framework of quality assurance for earth observation and should be required by any satellite mission (; ).

Historically, the ocean color community has relied on comparison with collocated in situ data to estimate uncertainty in ocean color products via “matchup” validation analyses. While the validation generates reasonable confidence in those products, it has some limitations (). In addition, while an analytical framework for estimating standard uncertainties in derived biogeochemical data products has been developed (), implementing this operationally requires pixel-level knowledge of standard uncertainties in Rrs(λ) and the error covariance between Rrs at different bands. Because of these limitations, we developed a derivative-based approach to calculate the pixel-level standard uncertainty in Rrs(λ) retrieved from the Moderate-Resolution Imaging Spectroradiometer (MODIS) aboard Aqua using a multiple scattering epsilon (MSEPS) atmospheric correction algorithm (), which has been used recently to complete a multi-mission ocean color reprocessing for heritage sensors by the Ocean Biology Processing Group (OBPG) at NASA. The Rrs(λ) uncertainty products generated by and were demonstrated to be reliable through evaluation using the Monte Carlo method as well as closure analysis with the statistics derived from validation.

While uncertainty can be readily generated, one practical issue is the delivery of all spectral uncertainty products without increasing data file sizes by an order of magnitude or more. A spectral error covariance matrix, ∑Rrs, is necessary. This matrix contains the variance in Rrsi) at the ith band as the diagonal elements and the covariance between Rrs at the ith and jth band, denoted as u(Rrsi), Rrsj)), as the off-diagonal elements. We note the off-diagonal terms in ∑Rrs are not negligible () and must be considered when more than one band is used to derive data products, for example, the band ratio algorithms for chlorophyll-a pigment concentration (chla, mg/m3) or the diffuse attenuation coefficient for downwelling irradiance at 490 nm (Kd(490), m−1).

Taking MODIS, for example, with N = 10 sensor spectral bands spanning 412 nm −678 nm, ∑Rrs has a dimension of 10 × 10. Considering that ∑Rrs is symmetric, a total of 55 (N × (N+1)/2) numbers must be stored at each pixel. The issue becomes more pertinent as we move forward to hyperspectral sensors with hundreds of bands, such as the Ocean Color Instrument (OCI), with 286 bands ranging from ultraviolet to shortwave infrared wavelengths on NASA’s Plankton, Aerosol, Cloud, ocean Ecosystem mission (PACE) (https://pace.oceansciences.org/about.htm), which launched in February 2024. It is unfeasible to deliver tens of thousands of ∑Rrs values at each pixel, given the broad (∼2,600 km) swath and moderate (1.2 km at the sub-satellite point) spatial footprint. As a first step toward delivering uncertainty products without unmanageably increasing the data file size, we investigate in this study the spectral variation of ∑Rrs retrieved from MODIS and OCI data using the MSEPS algorithm, with the goal of decreasing the dimensionality needed to reproduce it well enough for downstream scientific use.

2 Satellite data

Level-1A (uncalibrated) data of MODIS aboard the Aqua satellite and Level-1B (calibrated and geolocated) data of OCI from the PACE mission were downloaded from the Ocean Biology Distributed Active Archive Center (OB.DAAC). Level-1A MODIS data were processed to Level-1B using the SeaDAS software package (https://seadas.gsfc.nasa.gov) and the latest instrument calibration coefficients, which are also distributed by OB.DAAC. Rrs(λ) was derived by applying the MSEPS atmospheric correction algorithm to L1B data, and the uncertainty in Rrs(λ) was then calculated using the derivative-based approach established by . Here, we used MODIS data (i) over the Marine Optical BuoY (MOBY) () during 2002–2019; (ii) over the U.S. East Coast on 28 January 2012; and (iii) globally during 3–10 December 2019. We used OCI data over the U.S. East Coast on 13 May 2024.

3 Methodology

3.1 MSEPS atmospheric correction

We used MSEPS to calculate Rrs(λ) from spectral top-of-atmosphere (TOA) radiances Lt(λ), which can be expressed as Equation 1 (the spectral term, λ, is omitted for brevity):where tg(λ) is two-way gas transmittance, Lr(λ) (mW.cm−2.μm−1.sr−1) is the radiance from scattering by air molecules in the absence of aerosol, La(λ) (mW.cm−2.μm−1.sr−1) is the radiance from scattering by aerosols that also accounts for interactions between air molecules and aerosol scattering, Lf(λ) (mW.cm−2.μm−1.sr−1) is the radiance from scattering by surface whitecaps, Lg(λ) (mW.cm−2.μm−1.sr−1) is the sun glint, and tvr(λ) and tv(λ) are the diffuse transmittance from surface to sensor with the former only including molecular scattering and the latter including molecular and aerosol scattering. T(λ) is the beam transmittance for the view path. Lf (λ) is calculated using wind speed based on an empirical model (). A sun glint coefficient is calculated using a statistical model (), from which Lg(λ) can be derived (). The calculation of La(λ) is detailed by . After removing Lr(λ), La(λ), Lf(λ), and T(λ) Lg(λ) from Lt(λ), water-leaving radiance Lw(λ) (mW.cm−2.μm−1.sr−1) can be derived, and Rrs(λ) is then calculated by normalizing Lw(λ) to downwelling irradiance (the spectral term, λ, is omitted for brevity):where fb is the bidirectional reflectance correction calculated from the model presented by using chla, ts(λ) is the spectral diffuse transmittance from the Sun to the surface, θs is the solar zenith angle, and F0(λ) (mW.cm−2.μm−1) is the spectral extraterrestrial solar irradiance accounting for the Earth–Sun distance. F0 for OCI is calculated from the solar irradiance measurements from the Total and Spectral Solar Irradiance Sensor-1 (TSIS-1) (). F0 for MODIS is from .

3.2 Error covariance in Rrs

While the approach for calculating the error covariance matrix terms u(Rrsi), Rrsj)) is detailed by , we briefly recap here, taking MODIS as an example. It is based on a derivative approach. Two variables y1 and y2 are function of variables x, z, expressed by Equation 3:error covariance between y1 and y2, denoted by u(y1, y2), following the convention of (), can be derived from Equation 4:where u(xi, zj) represents the error covariance between variables xi and zj, and are the partial derivatives of y1 with respect to xi and of y2 with respect to zj, and T is the transpose operator. u(y1, y2) turns to variance when f1 and f2 are same, denoted by u2(y1). u(y1) is the standard uncertainty (). Based on Equations 2,4, u(Rrsi), Rrsj)) can be calculated from Equations 57:whereandXi is a vector of the uncertainty variables that are used to propagate the uncertainty through the atmospheric correction procedure, which is defined as Equation 8:where NIR corresponds to the two near-infrared bands used for MSEPS, is the aerosol optical thickness, and rh is the relative humidity. Lrfc is defined as Equation 9 (the spectral term, λ, is omitted for brevity):

3.3 Polynomial fitting of spectral u(Rrsi), Rrs(λ))

Spectral u(Rrsi), Rrsj)) are analyzed by fitting u(Rrsi), Rrs(λ)) as a function of wavelength. Here, λi and λ represent a specific wavelength and all the wavelengths, respectively. For example, u(Rrs(412), Rrs(λ)) is fitted as a function of MODIS bands in the wavelength range of 412–678 nm. Only visible wavelengths are considered in this study. We noted previously that the covariance was often spectrally smooth and decayed as a function of wavelength difference (). As a result, we experimented using a third-degree polynomial:where λ is the central wavelength in units of μm. As the matrix is symmetric, we do not “re-use” wavelengths when fitting polynomials as this will lead to an asymmetry. There are, therefore, ten points used to fit the coefficients for 412 nm (λ = 412, 443, 469, 488, 531, 547, 555, 645, 667, and 678), nine for 443 nm, eight for 469 nm, and so on. This means that there is no need to fit coefficients for wavelengths of 555 nm or longer, considering that there are four or fewer covariances to use for this band, which is equal to or less than the four coefficients needed for the fitting (so it is more efficient to reproduce the matrix exactly). The same fitting method is also applied to OCI. The coefficients a, b, c, and d in Equation 10, derived by means of least-square fitting, are stored in the Level-2 file for each band λi.

3.4 Uncertainty in chla and Kd(490)

Variance in chla (denoted by u2(chla)), calculated from the NASA standard algorithm (; ), and in Kd(490), denoted by u2(Kd(490)), calculated from the algorithm presented by is used here to evaluate the effect from polynomial fitting. Specifically, u2(chla) and u2(Kd(490)) are calculated twice. First, they are calculated using the ∑Rrs generated from the derivative approach, and then they are calculated using the polynomial model for ∑Rrs. The difference in u2(chla) and u2(Kd(490)) calculated from these two approaches is used to evaluate the effect of polynomial fitting. Based on the derivative approach, u2(chla) and u2(Kd(490)) are calculated from Equation 11:

MODIS bands at 443 nm, 488 nm, 547 nm, and 667 nm are used in the standard chla algorithm. Bands at 488 nm and 547 nm are used to calculate Kd(490). umc and umk are the model fitting uncertainty for chla and Kd(490), calculated by multiplying the mean relative uncertainty by chla or Kd(490). A mean relative uncertainty of 13% was derived by averaging the relative difference between model-fitted and in situ chla over all the data points that are used to fit the model coefficients in . A mean relative uncertainty of 10% was used for Kd(490) in this study.

4 Results

4.1 Effect of the fitting on u(Rrsi), Rrs(λ))

Figure 1 shows an example using a MODIS Aqua scene of the North Atlantic Ocean adjacent to the U.S. East Coast. In the example, we compare u2(Rrs(412)) and u(Rrs(443), Rrs(547)) calculated from the derivative approach and that derived from third-degree polynomial fitting. Values are higher at the edge than at the center of the swath. These higher values result from the high TOA signal due to the longer path length, which leads to higher sensor and forward model uncertainty. As a blue-to-green ratio of Rrs(λ) is used to calculate chla through the traditional ocean chlorophyll-type (OCx) algorithm (), we chose to show the error covariance between Rrs(λ) at the blue and green bands to evaluate the effect of modeling u(Rrsj), Rrsi)) with the polynomial. Like Figures 1, 2 shows an example using OCI data over the same geographic domain. We can see from Figures 1, 2 that u(Rrsj), Rrsi)) derived from the fitting matches well with that derived from the derivative approach, with the ratio in the range of 0.95–1.05. They show a similar spatial variation pattern.

FIGURE 1

FIGURE 2

Figure 3 further shows the comparison between u(Rrsj), Rrsi)) calculated from the derivative approach and that derived from third-degree polynomial fitting at the global scale. u(Rrsj), Rrsi)) is relatively higher in coastal and high-latitude regions. Like Figure 2, the third-degree polynomial fitted u(Rrsj), Rrsi)) matches very well with the original one, showing a similar spatial pattern, with mean ratios (standard deviations) of 0.989 (0.002), 0.99 (0.015), and 0.83(1.866) for u(Rrs(412)), u(Rrs(443), Rrs(547)), and u(Rrs(443), Rrs(667)), respectively. While the distribution is broader for u(Rrs(443), Rrs(667)), the absolute values are small, so this is a ratio of very small numbers.

FIGURE 3

While Figures 13 show that fitted u(Rrsj), Rrsi)) has a similar spatial pattern as the original fully computed one, Figure 4 shows a time series u(Rrsj), Rrsi)) derived from the fitting compared with the original over MOBY during 2002–2010, which shows consistent temporal variations, although u(Rrs(443), Rrs(667)) exhibits some bias. Considering the small value, such bias should not affect the uncertainty in the downstream products too much.

FIGURE 4

4.2 Fitting of u(Rrsi), Rrs(λ)) as a function of wavelength

Figure 5 shows two examples of the fitting of u(Rrs(412), Rrs(λ)) and u(Rrs(531), Rrs(λ)) versus wavelength over two pixels (denoted by A and B in Figure 1, with pixel A over coastal waters near the Chesapeake Bay and pixel B over open ocean waters of the North Atlantic Ocean. Like Figures 5, 6 shows examples of the fitting of u(Rrs(413), Rrs(λ)) and u(Rrs(530), Rrs(λ)) for OCI over two pixels (denoted by A and B in Figure 2). The high R-squared value (>0.99) indicates the reliability of the fitting. Note the comparatively higher error covariance over coastal waters than over clear waters across visible wavelengths.

FIGURE 5

FIGURE 6

4.3 Effect on uncertainty in chla and Kd(490)

Rrs(λ) is fundamental to deriving a host of bio-optical and biogeochemical properties of the water column. Taking chla and Kd(490) as examples, the effect of using the polynomial approximation for error covariance in Rrs(λ) is evaluated on u2(chla) and u2(Kd(490)). Specifically, u2(chla) and u2(Kd(490)) calculated using u(Rrsj), Rrsi)) derived from the derivative approach are compared with those calculated using u(Rrsj), Rrsi)) derived from polynomial fitting. We can see from Figure 7 that the relative uncertainty in chla and Kd(490), denoted by δ(chla) (u(chla) × 100/chla) and δ(Kd(490)) (u(Kd(490)) × 100/Kd(490)), calculated using original and polynomial fitted error covariance show similar spatial variation, with the absolute difference smaller than 0.5%.

FIGURE 7

5 Conclusion

We assessed how well a polynomial fitting can reproduce the full spectral error covariance in Rrs from MODIS and OCI data, with the goal of decreasing the Level-2 file size and therefore facilitating efficient delivery of those products. We found that a third-degree polynomial works well for data processed using the MSEPS algorithm, which is the standard algorithm used in Rrs retrieval for all global multispectral ocean color missions processed by NASA. We first evaluated the fitting by comparing the values numerically to the original error covariance source data, finding ratio close to 1 (typically from 0.95 to 1.05). We also found that using the fitted approximation of the covariance matrix to estimate the uncertainty of the downstream products chla and Kd(490) led to only a small absolute difference of the relative uncertainty that was less than 0.5%. Thus, we conclude that the accuracy of the approximation appears sufficient for most data user purposes.

By including the coefficients derived from polynomial fitting rather than the full ∑Rrs, a standard Level-2 file will decrease from ∼231 MB to ∼165 MB for a 5-min granule shown in Figure 1. The standard Level-2 file contains τa (869), Ångström exponent, Rrs(λ) at wavelengths from 412–678 nm, chla, Kd(490), particulate inorganic carbon (pic), particulate organic carbon (poc), instantaneous photosynthetically available radiation (ipar), normalized fluorescence line height (nflh), and photosynthetically available radiation (par). Taking the OCI 5-min granule in Figure 2 as an example, the Level-2 file would be ∼60 GB, including full ∑Rrs, which decreases to ∼1.7 GB by including the polynomial fitting coefficients.

Another benefit of using the polynomial model to estimate Rrs(λ) error covariance terms comes into play during operational data processing. As highlighted in this study, the full set of sensor bands is not always used in bio-optical algorithms (e.g., two for Kd(490) and four for chla). The approach described can thus facilitate dynamic computation of only the ∑Rrs terms required by an algorithm, thereby reducing system memory overheads and improving computational efficiency.

Statements

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

MZ: Conceptualization, Data curation, Investigation, Validation, Writing – original draft, Writing – review and editing. AI: Writing – review and editing. BF: Conceptualization, Funding acquisition, Supervision, Writing – review and editing. AS: Writing – review and editing. PW: Funding acquisition, Writing – review and editing. LM: Writing – review and editing.

Funding

The author(s) declare that financial support was received for the research and/or publication of this article. NASA Terra and Aqua Senior Review for MODIS algorithm maintenance and the NASA PACE Project.

Acknowledgments

We are grateful to NASA’s Ocean Biology Distributed Active Archive Center for making available all the in situ data and MODIS, OCI data used in our analysis.

Conflict of interest

Author MZ was employed by Science Application International Corp. Author LM was employed by Go2Q Pty Ltd.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Summary

Keywords

error covariance, ocean color, PACE, OCI, remote sensing reflectance, MODIS

Citation

Zhang M, Ibrahim A, Franz BA, Sayer AM, Werdell PJ and McKinna LI (2025) Practical aspects of providing pixel-level spectral Rrs error covariance in satellite ocean color products. Front. Remote Sens. 6:1670390. doi: 10.3389/frsen.2025.1670390

Received

21 July 2025

Accepted

23 September 2025

Published

10 October 2025

Volume

6 - 2025

Edited by

Zhigang Cao, Chinese Academy of Sciences (CAS), China

Reviewed by

Agnieszka Bialek, National Physical Laboratory, United Kingdom

Sundarabalan Balasubramanian, Geo-Sensing and Imaging Consultancy, India

Updates

Copyright

*Correspondence: Minwei Zhang,

Disclaimer

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.

Outline

Figures

Cite article

Copy to clipboard


Export citation file


Share article

Article metrics