*Article* **MODIS Aqua Reflective Solar Band Calibration for NASA's R2018 Ocean Color Products**

**Shihyan Lee 1,2,\*, Gerhard Meister <sup>2</sup> and Bryan Franz <sup>2</sup>**


Received: 9 August 2019; Accepted: 12 September 2019; Published: 20 September 2019

**Abstract:** Remote-sensing ocean color products have stringent requirements on radiometric calibration stability. To address a calibration deficiency in Moderate Resolution Imaging Spectroradiometer (MODIS) Aqua in recent years, the NASA Ocean Biology Processing Group (OBPG) developed a new calibration for reflective solar bands. Prior to the reprocessing of NASA's ocean color products for 2018 (R2018), the OBPG MODIS products had been based on calibration provided by the MODIS Calibration Support Team (MCST). Several modifications were made to the MCST calibration approach to improve the calibration accuracy for ocean color products. These include (1) applying 936-nm detector normalization to solar diffuser stability monitor (SDSM) data to reduce coherent noise; (2) modeling solar diffuser (SD) degradation wavelength dependency to determine SD degradation in near-infrared and shortwave infrared wavelengths; (3) computing detector gains using SD screen-closed data to better match ocean radiance levels in all bands; (4) performing a simple atmospheric correction to reduce bidirectional reflectance distribution function (BRDF) effects in desert trends; (5) estimating and using modulated relative spectral response (RSR) impact on ocean data to adjust the calibration coefficients; (6) using smoothing to characterize the temporal change in calibration; and characterizing response versus scan angle (RVS) changes using 2nd-order polynomials to improve spatial/temporal calibration stability. Relative to the previous R2014 ocean color products, the R2018 calibration removed the suspect late-mission global trends in blue-band water-leaving reflectance and some anomalously large short-term variability (spikes) in the temporal trend of chlorophyll concentration. This paper will describe the OBPG calibration with a focus on the differences between the MCST and OBPG approaches.

**Keywords:** satellite; calibration; solar diffusor; SDSM; desert trend; lunar calibration; RVS; MODIS; Aqua; ocean color

#### **1. Introduction**

The Moderate Resolution Imaging Spectroradiometer (MODIS) on the Earth Observing System (EOS) platform has 36 spectral bands to provide near-global observations every 2 days [1]. Of the 36 spectral bands, 20 are reflective solar bands (RSB) with a spectral range of 0.41 to 2.1 μm. Within the RSB, bands 8–16 are optimized to observe ocean biological processes (Table 1) and the rest of the RSB are designed for land and atmosphere applications but can also be used in ocean science application with reduced accuracy.

The MODIS is a key instrument aboard the Terra and Aqua Satellites. The two EOS spacecrafts, Terra and Aqua, were launched on 18 December 1999 and 4 May 2002, respectively, into sun-synchronous polar orbits at an altitude of 705 km. Terra is on a descending node with an equator crossing time of 10:30 AM (local time) and Aqua is on an ascending node with a 1:30 PM equator crossing time. Since they began operations, both instruments have continuously provided global earth observation data to study land, ocean, and atmospheric processes. However, due to the changing polarization sensitivities in MODIS Terra (MODIST) [2], the Terra ocean products have been relying on Ocean Biology Processing Group (OBPG) crosscalibration to correct the temporal calibration biases [3]. The crosscalibration uses MODIS Aqua (MODISA) global 7-day composites of water-leaving radiance as a truth field, combined with modeled atmospheric path radiance and surface contributions, to estimate the temporal change in gain and polarization as a function of time and scan angle.



Because of the MODIST-MODISA crosscalibration, the Terra ocean products are not considered an independent science data product and the accuracy of Terra products is dependent on the quality of Aqua products. Therefore, the calibration accuracy of MODISA is crucial, as the entire suite of MODIS ocean color products rely on the quality of Aqua products. Historically, the Aqua ocean color products are produced using the MODIS Calibration Support Team (MCST) radiometric calibration and are fine-tuned by OBPG crosscalibration [4] and vicarious calibration [5]. The Aqua ocean color products have been maintained at science quality until recent years [6]. As it is getting increasingly challenging to calibrate the aging MODISA sensor, a decision was made by OBPG to develop an independent radiometric calibration specifically for the ocean color products.

The main challenge of calibrating remote-sensing ocean color products is its sensitivity to calibration uncertainties. In general, calibration errors are magnified by anywhere from 5 to 100 times in ocean color products [7]. To produce meaningful climate data records, a long-term radiometric calibration stability of about 0.1% is often needed [8]. To achieve this stringent requirement, there are several advantages for an independent, end-to-end calibration solution. First, because it is all "in house", the radiometric calibration process can be efficiently integrated into the data production life cycle (Figure 1). For each version of calibration look-up-table (LUT), corresponding vicarious coefficients and crosscalibration LUT [4,5] are regenerated and a global, life-of-mission test processing is performed to assess the impacts. The uncertainties in radiometric calibration are significantly amplified during ocean color product retrieval [7]. This is further complicated by the subsequent changes in vicarious and crosscalibration in response to given changes made in radiometric calibration. Therefore, it is essential to integrate the testing of a radiometric calibration LUT into the data production life cycle in order to accurately evaluate a radiometric calibration LUT on ocean color products.

**Figure 1.** Ocean Biological Processing Group (OBPG) ocean color production process: L1A is the formatted raw instrument data. L1B is the calibrated radiance. L2/L3 are the level 2 and 3 ocean color products [8].

For each version of a calibration LUT, a global mission test reprocessing is performed that includes the first 4-day period of each month over the life of the mission. Each 4-day period is binned into a 9.2-km sinusoidal projection to form a 4-day mean global distribution of ocean color products for each month, which is then separated regionally (e.g., geographically and by water type) and trended with time. Temporal anomaly trends are also computed by subtracting the mean seasonal cycle from the global or regional time series. The ocean color products are reviewed by OBPG staff based on the known ocean biological processes to identify any product deficiency that could have originated from calibration issues. This process is vitally important to achieve the stringent radiometric calibration stability required to maintain ocean color products at science quality [9]. Maintaining radiometric calibration stability is especially difficult for Aqua MODIS as its radiometric response has changed significantly during more than 17 years of operation [10].

An ocean color specific calibration LUT can be optimized to the spectral properties of ocean scenes. For example, the solar calibration for RSB are performed at 2 different solar diffuser (SD) radiance levels (with and without SD screen attenuation). Normally, the land bands are calibrated at high SD radiance and ocean bands are calibrated using low SD radiances to maximize the signal-to-noise ratio (SNR). Since the spectral radiance of the ocean is closer to the low SD radiances, the OBPG solar calibration computes all detectors gains at the low SD radiance levels. Because the MODIS calibration algorithm uses a linear gain, calibrating detector gains at radiance closer to ocean radiance will reduce the impact of nonlinearity in detector responses on ocean color products. This approach, however, will increase nonlinearity impacts on high radiance scenes.

In addition, some radiometric characteristics can be better refined in subsequent OBPG crosscalibration and vicarious calibration and the radiometric calibration strategy can be simplified to improve calibration stability. For example, a previous study showed that the response versus scan angle (RVS) of band 8 (412nm) can be better characterized by a 4th-order fit [11]. However, our preliminary analysis showed that the fit residual reduction by the higher-order fit is insignificant. Compared to using the 2nd-order fit (same as in instrument prelaunch calibration), the 4th-order RVS fit is less stable temporally. Since the OBPG crosscalibration already used 4th-order polynomials to fine-tune the residual RVS uncertainties on ocean color products, using 2nd-order fit for RVS can improve temporal calibration stability while achieving similar RVS calibration accuracies on ocean color products.

The flow chart in Figure 2 summarizes the steps used to generate the OBPG RSB LUT. First, we used the yaw maneuver data to characterize the solar diffuser stability monitor (SDSM)/SD sunscreen transmission functions (tau) and Bidirectional Reflectance Function (BRF) functions and to store them in an internal LUT. In each subsequent SD calibration event, the SDSM calibration computes the relative SD reflectivity (H factor) using SD door open data. After the H factor is determined, the SD calibration computes the detector gain (m1) from the detector response and the SD radiance estimate from H factor and solar incident radiance. The computed m1 is then adjusted by modulated relative spectral response (mRSR) impact on ocean color products based on the typical ocean spectral radiance [12]. The mRSR impact is computed by optical degradation estimated by each detector's radiometric gain (m1) and electronic gain determined by electronic self-calibration (Ecal).

To estimate temporal change in RVS, desert and lunar observations are used to track instrument response changes at various scan angles. At the desert site, the time series of calibrated instrument reflectance, after correcting the atmospheric contributions, are computed as desert trends, which are used to track the RVS change at the 16 repeated observation angles. The calibrated lunar irradiance is computed from the monthly lunar observation. The calibrated lunar irradiances are compared to RObotic Lunar Observatory (ROLO) model predictions to produce lunar trend, which is used to track the RVS at the lunar observation angle.

The raw m1 and RVS values are estimated at the temporal interval of the data collection frequency, i.e., solar, lunar, and desert calibration. A temporal smoothing is applied to generate the time series of m1 and RVS. The smoothed m1 and RVS time series are used to populate the RSB LUT at the predefined time intervals.

**Figure 2.** Aqua RSB calibration procedures: The calibration process produces an RSB look-up-table (LUT), which contains temporal gain (m1) and response versus scan angle (RVS) parameters at time steps of m1t and RVSt respectively. The grey boxes are the source data; the yellow boxes are the calibration processes; the intermediate calibration parameters are outlined by red lines; and the red box is the final product of RSB calibration LUT. Ecal: electronic self-calibration; Yaw: spacecraft yaw maneuver.

In this paper, we will describe the calibration methodology and procedures applied by the OBPG in producing the RSB calibration LUT for the R2018 ocean color reprocessing and forward stream processing. The paper will focus on describing the differences between OBPG MODISA calibration and the approach used in the MCST calibration. The principle mechanisms of MODIS RSB calibrations will not be described in detail, as they have been discussed extensively in past literatures during the 19+ years of MODIS operation.

#### **2. MODIS L1B Calibration Algorithm**

The MODIS L1B calibration algorithm uses a nominal gain (m1) and response versus scan (RVS) angle function to convert the detector response of each earth view pixel to top-of-atmosphere (TOA) reflectance (Equation (1)) [13].

$$\text{Ref} = \text{dn}^\*(b, \text{ms}, d, p) \text{ \* } \text{m1}(b, \text{ms}, d, t) \text{ / } \text{rvs}(b, \text{ms}, p, t) \tag{1}$$

where Ref is reflectance, b is band, ms is mirror side, d is detector, t is time, and p is pixel number. Notice the scan angle in RVS is determined by pixel number and is independent of the detector. Both m1 and RVS are time dependent as both parameters are changing over time. Operationally, the instantaneous m1 and RVS values are interpolated from the calibration LUT that stores m1 and RVS parameters at predefined time intervals.

#### **3. Detector Gain**

The MODIS RSB radiometric gain is primarily determined by SD and SDSM calibration assembly (Figure 3). Table 1 lists the center wavelengths of the SDSM detectors and Aqua RSB. During solar calibration events, the SDSM calibration estimates the relative SD reflectivity by the ratio of the measured SDSM detector's sun and SD view responses. The SD radiance during solar calibration can be computed by the solar incident angle on the SD and SDSM-estimated SD reflectivity. Finally, the SD calibration determines the detector gain as the ratio of instrument response and estimated SD radiance [11,13].

**Figure 3.** Schematic of solar diffuser (SD)/SDSM calibration assembly. The SD sunscreen has ~7.5% transmission and SDSM sunscreen has ~1.4% transmission.

#### **4. Method**

#### *4.1. SDSM Calibration*

The MODIS SDSM is designed to track the change of SD reflectivity during on-orbit operation. The knowledge of SD reflectivity is the key to accurately computing detector gain, as the reflectivity of the SD is known to degrade on-orbit due to UV exposure. The MODIS SD degradation is described by the H factor, which is the ratio of instantaneous SD reflectivity to the reflectivity of the pristine diffuser. The OBPG method of estimating SD degradation was described previously in Reference [14]. The notable differences between our method and the one used in the MCST calibration are (1) use of detector 9 normalization to remove the wavelength-dependent coherent noise; (2) use of a wavelength model to estimate SDSM detector 9 degradation; and (3) use of the wavelength model to interpolate and extrapolate SD degradation measured at SDSM detector wavelengths to RSB wavelengths. Figure 4 shows the SDSM H factors and model-estimated SD degradation for bands 5–7.

**Figure 4.** (**a**) SD reflectivity (H) for SDSM detectors 1 to 9 and (**b**) SD wavelength model-estimated SD degradation at MODIS Aqua (MODISA) bands 5 (1240 nm), 6 (1640 nm), and 7 (2130 nm).

#### 4.1.1. Detector Gain: m1

The primary MODIS L1B product for RSB is the calibrated TOA reflectance factor [13]. Operationally, a LUT storing a time series of calibration coefficients is used to convert the instrument responses to TOA reflectance factors (see Equation (1)). The calibration coefficients *m*<sup>1</sup> for band (b), detector (d), and mirror side (ms) are computed as follows:

$$m\_1(b,d,ms) = \frac{BRF\_{SD}(b\_{\prime\prime}) \* H(b) \* \cos\theta}{dn\_{SD}^\*(b,d,ms) \* d\_{ES}^2} \tag{2}$$

where b is band; ms is mirror side; d is detector; *dn*∗ *SD* is the background-subtracted, temperature-corrected digital count from the SD measurement; dES is the earth–sun distance; θ, γ are the solar zenith and azimuth angles on SD; H is the SD degradation; and *BRFSD* is the BRF of the SD characterized by yaw maneuver data. For data collected with an SD screen, *BRFSD* is the combined effects of SD screen transmission function (tau) and SD BRF [14].

The m1 coefficient is computed at each solar calibration event. During each solar calibration event, the earth–sun distance, instrument temperatures, and solar zenith and azimuth angles are computed from spacecraft telemetry. The SD BRF is computed using the solar zenith and azimuth angles, and the SD BRF LUT is derived from yaw maneuver data [14]. The RSB H factor is interpolated from the H factors at SDSM detector wavelengths based on the SD wavelength model [14]. dn\*SD is the detector response at SD view subtracted by the space view detector response and then adjusted by temperature effects based on the measured instrument temperatures.

The SD calibration is performed both with and without SD attenuation screen to provide 2 calibration radiance levels. The SD attenuation screen has a transmissivity of ~7.5%, resulting in screen open and closed SD radiance ratios of ~13. The high and low SD radiances are designed to calibrate land and ocean bands, respectively, to maximize the SNR within the band's dynamic range. Figure 5 shows that the typical ocean radiance is between the high and low SD radiances for bands below 600 nm. Above 600 nm, the typical ocean radiance is lower than the screen-closed SD radiance. Since the MODIS L1B calibration assumes a linear function, it is desirable to have the calibrated radiance be similar to the scene radiance to reduce biases due to detector response nonlinearity.

**Figure 5.** Comparison of typical ocean radiance (blue curve) and average Aqua SD radiance with (SAA\_close) and without (SAA\_open) SD attenuation screen assembly.

The ocean band detectors saturate at SD radiances without SD screen attenuation. Figure 6 shows MODISA band 9 (443 nm) dn\*SD during solar calibration events. Band 9 detectors are saturated when the SD screen is open (Figure 6a). The screen-closed dn\*SD shows a decreasing trend and an annual cycle due to seasonal variation in earth–sun distances. Using the screen-closed dn\*SD, we compute m1, as shown in Figure 6b.

**Figure 6.** MODISA band 9 SD calibration: (**a**) detector 1 background-subtracted, temperature-corrected dn at SD view for screen open (black) and screen closed (red) and (**b**) computed m1 for detector 1 (black) and detector 10 (red).

For MODIS land bands, the m1 can be computed by either SD screen open or closed data, as both high and low SD radiances are within the detectors' dynamic range. Figure 7 shows that the band 2 detector gains computed using high (screen open) and low (screen closed) SD radiances have up to 1.5% differences. The ratios of m1 computed from high/low SD radiances also changed over time and are detector dependent. Since m1 is computed as the linear fit, the differences in high/low radiance m1 indicate the detector gain is not linear over the dynamic range. The temporal trend in high/low radiance m1 ratios indicates the detector response nonlinearity is changing over time. Figure 7b shows that the temporal change in detector response nonlinearity is detector dependent. Similar detector response nonlinearity behaviors are observed for most land band detectors, and the temporal change of detector response nonlinearities are most prominent for bands 1 and 2. The detector nonlinearity behavior is likely caused by the readout electronics as the detector gain of MODIS high-resolution bands (bands 1–4) have been shown to be subframe and radiance dependent [15]. To optimize calibration for ocean color products, we used low SD radiance data to compute m1 to better match the typical ocean surface radiance.

**Figure 7.** MODISA band 2 SD calibration: (**a**) detector 1 m1 computed from screen open (black) and closed (red) data and (**b**) ratios of m1 computed using SD screen open and closed data for detectors 1 (black) and 10 (red).

#### 4.1.1.1. m1 Time Series

MODISA performs solar calibrations every 3 weeks after the first few years of operation [16]. The computed m1 (Figures 6 and 7) has noise of up to 0.5 %. As the change in detector gain is expected to be gradual, a boxcar smoothing function is used to compute the running average of m1. The smoothing interval is set at 3 months for green and blue bands and 1 year for red and Near-InfraRed (NIR) bands. The smoothing intervals are selected to maximize the temporal m1 smoothness while retaining its temporal trend. In Figure 8, MODISA band 8 shows sharp changes in detector gains between 2011 and 2012. The same features are shown in the other blue and green bands with reduced magnitudes. The shorter smoothing time interval for blue/green bands is used to retain the temporal features. For red/NIR bands, the m1 changes are gradual and do not have distinct features (Figure 8b), and the longer interval is used to generate smoother m1 time series.

**Figure 8.** MODISA-normalized m1 computed from SD calibrations (symbols) and the smoothed m1 used for calibration (red curve) for mirror side 1 and detector 1 of (**a**) band 8 and (**b**) band 16.

#### *4.2. RVS*

The RVS is the ratio of the detector response at a particular scan angle to the reference scan angle. In MODIS, the RVS is normalized to the SD angle, where the nominal detector gain (m1) is determined. Because the MODIS scan mirror is not protected, the MODISA RVS has changed significantly over time due to nonuniform mirror degradation [17]. We track the temporal RVS change by trending the instrument responses over pseudo-invariant calibration targets, i.e., SD, moon, and desert, observed at different scan angles [17,18]. For MODISA, the SD and moon are always observed at fixed scan angles. The desert can be observed at 16 angles from the repeatable orbital tracks. The SD is observed about every 3 weeks during solar calibration events. The moon is observed monthly except for the summer months when the moon is too close to the horizon. The desert site can be observed daily at one of the 16-day repeat scan angles, if the scene is free of clouds.

By comparing the relative trend difference at different scan angles, we can estimate the temporal change in RVS. The temporal trends at each angle and calibration target are normalized to the first measurement, and the RVS at each scan angle is compared relatively. Therefore, this method determines the change in RVS relative to the first measurement, where the Aqua RVS can be approximated as prelaunch measured values.

One of the limiting factors of this method is that the target brightness needs to be within the dynamic range of the detectors. Among the calibration targets, desert sites are too bright for bands 10 to 16 and certain parts of the moon are also too bright for bands 13 to 16, but the saturated moon pixel radiance can be estimated by the band ratioing method [19]. Therefore, for bands 1–9, 17–19, and 26, desert, moon, and SD are used to track RVS and only moon and SD are used to track RVS for bands 10–16 (Table 2).


**Table 2.** Calibration targets used to track MODISA temporal RVS change and fitting function used for RVS characterization. The desert site is not used for bands 10–16 due to saturation.

Ideally, all calibration targets should have the same brightness to reduce the detector response nonlinearity impact. However, the brightness of SD, moon, and desert are significantly different, which adds additional uncertainty to RVS estimates using a combination of calibration targets.

The time-dependent RVS is computed by fitting 2nd-order polynomials over the temporal relative response vs. scan angles from all available calibration targets (Table 2). As described earlier, 2nd-order polynomials are more stable as they are less affected by noise in lunar and desert trends. The residual higher-order RVS behaviors are removed by the OBPG crosscalibration. This process optimizes the overall RVS correction for ocean color products when used in conjunction with OBPG crosscalibration. However, the temporal RVS by itself might not be the best possible RVS characterization for other disciplines without additional product specific corrections.

Since the RVS is normalized at the SD angle, the lunar and desert trends are computed using the time-dependent m1 (see Section 4.1), so the moon and desert trends are essentially the RVS trends at their respective observation angles. Using the time-dependent detector m1 for lunar and desert trends remove the potential uncertainties from temporal change in mirror side and detector-to-detector gain ratios. Using time-dependent detector m1 is especially important for lunar trends, as the moon is not a uniform target and cannot be evenly sampled among detectors. Because the moon image is an ellipsoid, more lunar pixels are sampled by middle detectors than the edge detectors. If the temporal gain ratios between middle and edge detectors changes, this will cause biases in the lunar trend.

#### 4.2.1. Desert Trend

We used the Libya 4 site to track temporal RVS change. Libya 4 is one of the widely used pseudo-invariant sites due to its long-term surface reflectance stability [20]. Despite being a stable target, the satellite observations showed significant bidirectional reflectance distribution function (BRDF) effects with respect to solar incident and instrument view angles. A previous study characterized such effects by simultaneously fitting the solar incident, instrument view angles, and RVS changes [17]. Although this method can reduce the BRDF effects in the observed data, the method is empirical and requires a good assumption of the basic shape of RVS trends.

Since satellite observations are made with different geometrics between the sun and the instrument, the changing atmospheric path radiance from different solar incident and instrument view angles can cause the BRDF effect to be seen in the desert observation. To remove atmospheric contributions, an atmospheric model, 6SV [21], is used to estimate Rayleigh scattering and gaseous absorption. Figure 9 shows the band 8 mirror side 1 RVS characterization process at a −42-degree scan angle. The instrument TOA reflectance (Figure 9a) shows a large annual cycle. Figure 9b shows that the 6SV model estimated atmospheric contribution based on the solar and instrument view geometry. Figure 9b shows that the annual cycle of atmospheric contribution matches well with the annual variations of the TOA reflectance (Figure 9a). Removing the atmospheric contribution from TOA reflectance (Figure 9c) significantly reduces the annual cycles in the desert trending. In Figure 9c, the solar incident power and atmospherically corrected desert trend has a residual noise of a few percent, a significant reduction

from the observed data. The residual noise in the desert trend could come from variations in aerosol, polarization effects, and possibly small short-term BRDF changes in the Libya 4 site itself.

**Figure 9.** Desert trend for Aqua band 8, mirror side 1 at a scan angle of−42 degrees: (**a**) top-of-atmosphere (TOA) reflectance; (**b**) atmospheric contribution in reflectance; and (**c**) TOA reflectance—atmospheric reflectance. The TOA reflectance is computed using Equation (1) with RVS set to 1.

To further demonstrate the importance of atmospheric correction, we compared the desert site RVS before and after atmospheric correction in January 2003 when the RVS should have minimal deviation from the prelaunch measured values. In Figure 10a, the sensor-observed TOA reflectance showed a 60% variation in response at different scan angles. After correcting for the atmospheric contribution, the RVS variation (Figure 10b) is significantly reduced. Since, the prelaunch measured RVS for band 8 is only about 2%, a 60% variation in RVS is unlikely after only 6 months of on-orbit operation. The results indicate the atmospheric correction removed the bulk of the instrument observation BRDF artifacts to improve RVS characterization accuracy.

**Figure 10.** Libya 4 desert observations for MODISA band 8, mirror side 1 during 2003 January: The plots show the normalized reflectance (ref) versus scan angle (RVS) before (**a**) and after (**b**) atmospheric correction. Notice in some angles, 2 clear sky observations were made from the 16-day repeat orbits in January 2003.

To reduce the residual noise in the atmospherically corrected desert trends, a 1-year running average was computed and used to produce smoothed temporal trends at each observation angle. Figure 11 shows examples of the computed desert trends and the smoothed RVS trends at −42 degrees for bands 8 and 9. The smoothed RVS desert trends will later be combined with lunar trends to characterize temporal RVS change.

**Figure 11.** RVS computed from desert trends of MODISA mirror side 1 for (**a**) band 8 and (**b**) band 9 at a −42-degree scan angle: The black symbols are the computed RVS at each desert observation, and the red curve is the temporal smoothed RVS.

#### 4.2.2. Lunar Trend

The moon is an ideal calibration target, as its surface is photometrically stable [22] and it can be observed by MODIS without atmospheric contamination. The geometrically corrected lunar irradiance can be obtained from the US Geological Survey's RObotic Lunar Observatory (ROLO) photometric model [23], using the observation time and the satellite position to determine the lunar phase and libration angles and sun–moon/satellite–moon distances. The lunar trend is computed as the ratio of the sensor-observed and model-predicted lunar irradiances. The temporal changes in the lunar trends provide the temporal RVS change at lunar observation angle [24].

The lunar trends (Figure 12) are normalized to the first observation to remove the residual biases between SD-calibrated and ROLO-predicted lunar irradiance. The normalized lunar trends are the temporal RVS changes at an Space View (SV) angle, which can be up to 20% in band 8 (412 nm). The temporal RVS changes are more prominent in the blue bands (8, 9, 3, and 10) and decrease at longer wavelengths. The temporal RVS trends resemble the m1 trends (Figure 8), where the trends in the blue bands reversed around 2011. For green to NIR bands, the temporal RVS changes are small with no clear reversing trends like the blue bands.

**Figure 12.** MODISA lunar trends for selected RSBs for (**a**) mirror side 1 and (**b**) mirror side 2: The symbols indicate each lunar calibration data point, and the curves are the fitted lunar trends.

The lunar trends show a quasi-annual cycle, which is likely due to the residual uncertainty in the ROLO model's libration correction [25]. To compute the RVS change at an SV angle, a 1-yr running average was computed to produce smoothed lunar trends, except for bands 8 and 9. A 2-piece fit was used for bands 8 and 9 before and after 2011 to better characterize the reversal in trends.

#### 4.2.3. RVS Characterization

To characterize the temporal RVS, the available desert trends are combined with lunar trends. Since the lunar and desert trends are computed with time-dependent SD m1, the lunar/desert trends are the RVS normalized to SD angles. For bands 1–9, 17–19, and 26, the temporal RVS is estimated as a 2nd-order polynomial fit between the scan angle and the smoothed RVS from lunar and desert trends (Table 2). For bands 10–16, the RVS change is small (<4%) and only lunar trends were used to characterize the temporal RVS change. Since the moon is always observed at the same scan angle, we used the prelaunch estimated RVS curve as the basic shape to characterize the 2nd-order RVS fit from the lunar observation angle alone. This method assumes that small changes in RVS do not significantly change the shape of the RVS.

Figure 13 shows the characterized RVS for selected ocean bands and years. As mentioned earlier, band 8 has the largest temporal RVS change. In 2011, the RVS variation was the largest with the mirror response at the beginning of the scan (−55 degrees) ~30% higher than at the end of the scan (55 degrees). For the green to NIR bands, the RVS and their changes are considerably smaller.

**Figure 13.** MODISA RVS fitting at selected years for (**a**) band 8 (412 nm), (**b**) band 12 (547 nm), and (**c**) band 16 (869 nm), mirror side 1.

#### *4.3. Relative Spectral Response (RSR) Modulation*

The temporal radiometric calibration for Aqua MODIS can be computed from the time-dependent m1 and RVS estimated in Sections 4.1 and 4.2. The calibration coefficients convert the instrument response to TOA reflectance or radiance. However, the Aqua detector RSR has been shown to change over time due to the spectrally dependent optical degradation [12]. The temporal RSR change will cause bias in observed ocean radiances. Applying time-dependent RSR in ocean color product retrieval is not feasible due to complexity and computational time limitations. The current solution is to apply an adjustment in calibration coefficients to correct the estimate of the calibration bias in ocean color products due to RSR modulation.

The derivation and correction of the time-dependent modulated RSR impact on Aqua ocean color products is detailed in our previous study [12]. The calibration bias due to RSR modulation is caused by the spectral differences between calibration targets and ocean scenes (Figure 14a). The RSR modulation impact on ocean color products is the combined effects of biases in observed SD, desert, and ocean scene radiances due to RSR modulation. Figure 14b shows the biases in ocean color product radiances due to modulated RSR computed using the typical spectral radiances in Figure 14a. Figure 14b shows that band 8 is the only band significantly impacted by RSR modulation with up to

1% biases in 2011. The large band 8 RSR modulation is the result of its large out-of-band response convolved with the large spectral radiance difference between calibration targets and the typical ocean scene. For the remaining ocean bands, the impacts from modulated RSR are less than 0.1%.

**Figure 14.** (**a**) Typical radiance for ocean (red), SD (green), and the desert (black) and (**b**) mean calibration bias (L\_mRSR/L\_RSR0) due to modulated RSR for bands 8 to 16: L\_mRSR is the typical ocean radiance computed using modulated RSR (mRSR), and L\_RSR0 is the typical ocean radiance computed using prelaunch measured RSR (RSR0).

The biases of modulated RSR are estimated per band for all ocean bands (bands 8 to 16), where RSRs are measured for both in-band and out-of-band regions. Although the modulated RSR impact is dependent on scan angle due to the angular-dependent optical degradation [12], we apply the scan-angle mean impact for simplicity. The correction is meant to correct the bias in global radiometric trends. The angular dependency of modulated RSR impact is small and will be added to residual RVS characterization uncertainty and corrected by OBPG crosscalibration.

#### **5. Results**

#### *5.1. OBPG RSB Calibration LUT*

The main differences between OBPG and MCST calibration are summarized in Table 3. Additional differences could come from computation methodology and data selection. Figure 15 compares the temporal trends of band 8 m1 and RVS between OBPG and MCST LUT at selected scan angles. Figure 15 shows that the OBPG and MCST gain (m1/RVS) at an SD angle is within 1.5% different but can be up to 4% different at large scan angles. The RVS differences at an SD angle are the modulated RSR effects as described in Section 4.3. The scan-angle-dependent gain differences between OBPG and MCST LUT will change the global trends in ocean color products.

The OBPG RSB calibration LUTs are based on the MCST LUT with the time-dependent m1 and RVS coefficients replaced by the values derived in Sections 4.1 and 4.2. The correction of modulated RSR impact is applied to RVS to keep the m1 coefficients as a pure product of solar calibration for better calibration traceability. Applying the modulated RSR correction to either m1 or RVS produces the same results.

Since the OBPG LUT retains the MCST LUT format, it can be used in MODIS L1B processing outside of OBPG. The OBPG produces a monthly updated LUT, based on the monthly updated MCST LUT to extend m1, RVS, and modulated RSR correction. The crosscalibration LUT (see following section) is also updated monthly after each new RSB LUT is generated to keep the correction up-to-date.


**Table 3.** Calibration methodology differences between R2014 and R2018.

<sup>\*</sup> Based on MCST publications [11,17].

**Figure 15.** Time series of OBPG and MCST mirror side 1, detector-averaged normalized band 8 for (**a**) m1 and (**b**) RVS and (**c**,**d**) their differences at SD, space view (SV), Nadir, beginning of scan (BOS; −55 degrees), and end of scan (EOS; 55 degrees). Solid line: OBPG, dashed line: MCST. Colors indicate scan angles. The scan angle-specific detector gain is computed as m1/RVS, and m1\_t0 is the SD m1 at the beginning of the mission. Note: the legend in Figure 15a is applied to Figure 15b–d.

#### *5.2. Artifact Mitigation*

Based on the OBPG radiometric calibration LUT, vicarious and crosscalibrations were performed to correct absolute instrument calibration biases and to flat-field the residual RVS, mirror side, and detector-to-detector uncertainties [4,5]. The vicarious and crosscalibration were performed using Rrs products; therefore, any residual atmospheric correction errors will also be included. Without

vicarious and crosscalibration, a perfectly calibrated sensor could still produce biased and/or stripped Rrs images due to residual error from atmospheric correction.

The vicarious calibration produces a constant factor for each band, which has no impact on temporal trends [5]. However, the residual RVS correction can potentially impact global temporal trends if the correction is large [4]. Figure 16 compares the crosscalibration-computed residual RVS correction (M11) between R2014 and R2018 products. The magnitude of R2018 M11 correction is smaller than the R2014 correction after 2011 and similar to before 2011. The R2014 M11 correction varies significantly after 2011, indicating RVS calibration issues. The large M11 correction will change the global temporal trends, as the mean earth scene correction is not 1. This is the case in R2014 products as the scene-averaged M11 correction changes significantly after 2011 (Figure 17).

**Figure 16.** Band 8 (412 nm), mirror side 1, detector 1, crosscalibration correction factors (M11) for select years: (**a**) R2018 calibration and (**b**) R2014 calibration. M11 is the RVS correction (unitless).

**Figure 17.** Time series of band 8, mirror side 1, scan-averaged M11 correction for (**a**) R2018 and (**b**) R2014 products: The scan-angle-averaged M11\_avg indicates the correction on the global mean. Colors indicate detectors 1–10.

In R2018 products, the M11 is normalized to ensure no temporal trends are introduced by the crosscalibration. Compared to R2014, the R2018 M11 has larger detector spread. This is because the MCST calibration incorporates detector-to-detector gain correction using ocean sites [11] for the calculation of their m1, whereas in R2018, the striping correction is applied only during the crosscalibration. There is no impact of the M11 detector spread on image quality in R2014 vs. R2018 after applying the OBPG crosscalibration.

#### **6. Discussion**

To demonstrate the impact of the calibration approach described in this paper, the ocean color products from the R2018 reprocessing are compared with the previous R2014 reprocessing. The only changes between these two reprocessing events were the updated instrument calibration approach and the updated vicarious calibration [26]. Figure 18 shows a comparison of the two reprocessing versions for the time series of spectral water-leaving reflectances (Rrs) and derived chlorophyll concentration, as derived from monthly averages for the globally distributed deepest ocean waters (>1000-m depth). We restrict our analysis to this deep-water subset as it allows for an assessment over the vast majority of the ocean surface while avoiding complex coastal regions with highly variable terrestrial influences. In Figure 18a, the R2014 water-leaving reflectances (solid lines) for the shortest wavelengths (412 nm and 443 nm) show a strong downward trend after 2012, which is not associated with any known geophysical events. This late-mission drift in the blue was significantly reduced in the R2018 products. Similarly, the large seasonal variability in the R2014 chlorophyll product, especially after 2012, is much reduced in the R2018 products (Figure 18b). Such variability will arise if the calibration in the blue is erroneously low, as the observed radiance has a strong seasonal cycle due to variations in solar geometry and associated atmospheric path radiance, and the correction for that variable atmospheric contribution to the TOA radiance will result in an underestimation of water-leaving reflectance in the blue. The chlorophyll concentration is derived from the ratio of blue to green water-leaving reflectances, such that a lower blue–green ratio implies a higher chlorophyll concentration, thus giving rise to seasonal peaks in the R2014 chlorophyll products. There is also an overall bias (increase) in the blue wavelengths between the R2014 and R2018 results, which is caused by an update to the vicarious calibration [26] and which also contributes to reduced seasonal variability in the R2018 chlorophyll product over the full mission lifetime.

**Figure 18.** MODIS global deep-water temporal trends for (**a**) water-leaving reflectance (Rrs) and (**b**) chlorophyll. Solid line: R2018, dashed line: R2014 products.

Another very useful subset of the oceans for assessment of product temporal stability is the mean over the clearest ocean gyres, far from terrestrial influences, where chlorophyll concentrations are very low (typically <0.1 mg m-3) and homogenous, and weak variations in phytoplankton abundance or physiology, which modulate the absorption of blue light much more than green light, are the only significant influences on the ocean water optical properties. Temporal anomalies of water-leaving reflectance for the green ocean channel (547 nm) and green land channel (555 nm) were computed by subtracting the mean seasonal cycle from the water-leaving reflectance time series as derived over oligotrophic waters. In the absence of any major geophysical events, we can expect this time series to be relatively stable over time. The temporal anomalies for R2014 (Figure 19a,c) show a long-term trend and shorter-term variabilities that are larger than in the R2018 results (Figure 19b,d). Also notable is that the R2014 results for the two green bands with very similar center wavelengths (Figure 19a,c) are

less consistent than the R2018 products (Figure 19 b,d). These results indicate that the R2018 products are more stable and consistent relative to the R2014 products.

**Figure 19.** MODIS water-leaving reflectance (Rrs) temporal anomaly in global oligotrophic water: (**a**,**c**) R2014 products for 547 nm and 555 nm and (**b**,**d**) R2018 products for 547 nm and 555 nm.

#### **7. Conclusions**

In this paper, we described the independently developed MODIS Aqua radiometric calibration used in OBPG R2018 ocean color products. The OBPG calibration is developed based on the general MODIS calibration principles with a focus on ocean color products. The in-house developed calibration process is integrated in the ocean-color-product production cycle which enables fast evaluation. Timely product evaluation facilitates accurate assessment of the combined changes in radiometric calibration and the subsequent OBPG vicarious and crosscalibration. The integrated calibration development process is key to achieving the high radiometric calibration accuracy required for ocean color science.

Because the calibration is developed specifically for OBPG ocean color products, it is optimized for OBPG ocean color processing. For example, the nominal detector gains are calibrated with low SD radiance to better approximate ocean surface radiance. Applying the OBPG calibration for high-radiance targets could result in biases due to detector response nonlinearity. Furthermore, the RVS is characterized by simple 2nd-order fits to improve temporal stability. Applying the OBPG crosscalibration is needed to obtain the optimal RVS characterization.

Using the OBPG-developed radiometric calibration, the R2018 ocean color products show improved quality over the previous products (R2014). The global trends of the R2018 products are more stable and do not have late-mission anomalous behaviors shown in the R2014 products. The R2018 products are more consistent across wavelengths than R2014 products. Compared to the R2014 coefficients, the OBPG crosscalibration in R2018 is more consistent over time, indicating a more stable RVS characterization in R2018. The improvement in R2018 ocean-color-product quality and crosscalibration stability are the results of improvements in R2018 radiometric calibration.

**Author Contributions:** S.L. development the methodology, performed the analysis and wrote the manuscript. G.M. provided the advice, oversaw the progress and contributed to the writing and revising of the manuscript. B.F. contributed to the ocean color product validation processes and the writing and editing of the manuscript.

**Funding:** This research was funded by "Advancing the Quality and Continuity of Marine Remote Sensing Reflectance and Derived Ocean Color Products from MODIS to VIIRS" (17-TASNPP17-0065), submitted to the Science Mission Directorate's Earth Science Division, in response to NASA Research Announcement (NRA) NNH17ZDA001N-TASNPP: Research Opportunities in Space and Earth Science (ROSES-2017).

**Acknowledgments:** The authors would like to thank MCST for providing in-depth knowledge of MODIS instrument characteristics and past and current calibration issues and developments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Glossary**


#### **References**

1. Xiong, X.; Barnes, W.L.; Chiang, K.; Erives, H.; Che, N.; Sun, J.; Isaacman, A.T.; Salomonson, V.V. Status of Aqua MODIS on-orbit calibration and characterization. In Proceedings of the Sensors, Systems, and Next-Generation Satellites VIII, Maspalomas, Spain, 13–15 September 2004.


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

**Xiaoming Liu 1,2,\* and Menghua Wang <sup>1</sup>**


Received: 18 December 2018; Accepted: 16 January 2019; Published: 18 January 2019

**Abstract:** The Visible Infrared Imaging Radiometer Suite (VIIRS) on the Suomi National Polar-orbiting Partnership (SNPP) and National Oceanic and Atmospheric Administration (NOAA)-20 has been providing a large amount of global ocean color data, which are critical for monitoring and understanding of ocean optical, biological, and ecological processes and phenomena. However, VIIRS-derived daily ocean color images on either SNPP or NOAA-20 have some limitations in ocean coverage due to its swath width, high sensor-zenith angle, high sun glint, and cloud, etc. Merging VIIRS ocean color products derived from the SNPP and NOAA-20 significantly increases the spatial coverage of daily images. The two VIIRS sensors on the SNPP and NOAA-20 have similar sensor characteristics, and global ocean color products are generated using the same Multi-Sensor Level-1 to Level-2 (MSL12) ocean color data processing system. Therefore, the merged VIIRS ocean color data from the two sensors have high data quality with consistent statistical property and accuracy globally. Merging VIIRS SNPP and NOAA-20 ocean color data almost removes the gaps of missing pixels due to high sensor-zenith angles and high sun glint contamination, and also significantly reduces the gaps due to cloud cover. However, there are still gaps of missing pixels in the merged ocean color data. In this study, the Data Interpolating Empirical Orthogonal Functions (DINEOF) are applied on the merged VIIRS SNPP/NOAA-20 global Level-3 ocean color data to completely reconstruct the missing pixels. Specifically, DINEOF is applied to 30 days of daily merged global Level-3 chlorophyll-a (Chl-a) data of 9-km spatial resolution from 19 June to 18 July 2018. To quantitatively evaluate the accuracy of the DINEOF reconstructed data, a set of valid pixels are intentionally treated as "missing pixels", so that reconstructed data can be compared with the original data. Results show that mean ratios of the reconstructed/original are 1.012, 1.012, 1.015, and 0.997 for global ocean, oligotrophic waters, deep waters, and coastal and inland waters, respectively. The corresponding standard deviation (SD) of the ratios are 0.200, 0.164, 0.182, and 0.287, respectively. Gap-filled daily Chl-a images reveal many large-scale and meso-scale ocean features that are invisible in the original SNPP or NOAA-20 Chl-a images. It is also demonstrated that the gap-filled data based on the merged products show more details in the dynamic ocean features than those based on SNPP or NOAA-20 alone.

**Keywords:** VIIRS; SNPP; NOAA-20; DINEOF; ocean color data; data merging; gap-filling

#### **1. Introduction**

Ocean color data are critical for monitoring and understanding of water optical, biological, and ecological processes and phenomena, and it is also an important source of input data for physical and biogeochemical ocean models [1]. Since the launch of the Visible Infrared Imaging Radiometer Suite (VIIRS) on the Suomi National Polar-Orbiting Partnership (SNPP) [2,3] in October 2011, ocean color products derived from VIIRS-SNPP have been routinely produced globally [3,4]. On 18 November 2017, the follow-up VIIRS sensor housed in the National Oceanic and Atmospheric Administration (NOAA)-20 satellite was launched as the first of four sensors in the Joint Polar Satellite System (JPSS) satellite series, and global ocean color data from NOAA-20 are also being routinely produced. For both SNPP and NOAA-20, chlorophyll-a (Chl-a) concentration [5–7], normalized water-leaving radiance spectra *nLw*(*λ*) [8,9], including new *nLw*(*λ*) data using VIIRS imaging bands [10], and water diffuse attenuation coefficient at the wavelength of 490 nm *Kd*(490) and at the domain of photosynthetically available radiation (PAR) *Kd*(PAR) [11,12], are all generated as standard VIIRS ocean color products. Some experimental products such as inherent optical properties (IOPs) [13,14], and a newly added quality assurance (QA) score product for measuring data quality [15] are also included for evaluation. In addition, to improve ocean color data quality over turbid coastal and inland waters [16–18], the shortwave infrared (SWIR)-based and near-infrared (NIR)-SWIR combined atmospheric correction algorithms have been used to routinely produce global VIIRS ocean color products for both SNPP and NOAA-20 [19,20]. Furthermore, VIIRS ocean color data processing algorithms have been significantly improved, including improved sensor on-orbit calibration using both solar and lunar approaches [21,22]. However, for either VIIRS-SNPP or VIIRS-NOAA-20, there are always missing pixels in the VIIRS-measured ocean color data imageries due to cloud cover and various other reasons, e.g., strong sun glint contamination, dust storms, very large solar- and sensor-zenith angles, etc. It is certainly useful to fill the gap of missing pixels before being used as input for ocean models and for various other applications.

As a follow-up VIIRS instrument, VIIRS-NOAA-20 is essentially built the same as VIIRS-SNPP. Therefore, the sensor characteristics of the two instruments are very similar. Both SNPP and NOAA-20 operate at the 824-km sun-synchronous polar orbit, which crosses the equator at about 13:30 local time. There is about a 50-minute delay between the paths of NOAA-20 and SNPP, which makes the NOAA-20's path running through the middle of two adjacent SNPP paths, and vice versa. The overlap of the spatial coverages of the two sensors automatically fills each other's gaps caused by high sensor-zenith angles and high sun glint contamination, and it significantly reduces the missing pixels in the merged images. In addition, ocean color products from SNPP and NOAA-20 have the same spatial and temporal resolution, and are processed with the same software package, i.e., the Multi-Sensor Level-1 to Level-2 (MSL12) ocean color data processing system [3]. Therefore, the statistics of their ocean color products are very close to each other, and in fact the data can be directly merged without adjustment for matching up each other's statistical properties. However, there are still lots of missing data in the merged VIIRS SNPP/NOAA-20 products.

To completely fill the gaps of missing pixels in the merged VIIRS SNPP/NOAA-20 ocean color data, the Data Interpolating Empirical Orthogonal Functions (DINEOF) [23,24] method is used in this study to reconstruct the missing data in the ocean color images. The DINEOF exploits the spatio-temporal coherency of the data to infer a value at the missing location and has been successfully adopted in various applications using ocean remote sensing data [25–30]. With more and more availability and usage of ocean color data in recent years, the DINEOF method has also been applied to ocean color data from various sensors including the Moderate Resolution Imaging Spectroradiometer (MODIS) [31–33], the Spinning Enhanced Visible and Infrared Imager (SEVIRI) onboard Meteosat Second Generation 2 [34], and the Korean Geostationary Ocean Color Imager (GOCI) [35]. Most recently, Liu and Wang [36] used the DINEOF to fill the gaps of missing pixels in VIIRS-SNPP global ocean color data. In that study, the DINEOF is applied to VIIRS-SNPP global Level-3 binned ocean color data of 9-km spatial resolution and the DINEOF reconstructed ocean color

data are used to fill the gap of missing data. In particular, daily, 8-day, and monthly VIIRS global Level-3 binned ocean color data, including Chl-a concentration, *Kd*(490), as well as *nLw*(*λ*) at the five VIIRS visible bands are tested and evaluated. Results show that the DINEOF method can successfully reconstruct and gap-fill meso-scale and large-scale spatial ocean features in the global VIIRS Level-3 images, as well as capture the temporal variations of these features.

In this study, the DINEOF method is used to reconstruct and gap-fill a merged VIIRS SNPP/NOAA-20 global daily ocean color product, specifically the Chl-a data. With the reconstructed (gap-filled) VIIRS global daily Chl-a images, ocean features can now be well identified and observed both spatially and temporally. Some examples that demonstrate the advantages and usefulness of the gap-filled VIIRS SNPP/NOAA-20 merged global ocean color products are provided. The DINEOF is also applied to ocean color data from a single sensor, VIIRS-SNPP or VIIRS-NOAA-20, and the gap-filled data based on SNPP or NOAA-20 are compared with those based on the two-sensor merged Chl-a product.

#### **2. Data and Methods**

#### *2.1. VIIRS SNPP and NOAA-20 Ocean Color Level-2 and Global Level-3 Data*

The MSL12 is the official NOAA VIIRS ocean color data processing system for both SNPP and NOAA-20 (and all other follow-on VIIRS sensors), and it has been used for processing VIIRS data from Sensor Data Records (SDR or Level-1B data) to the Environmental Data Records (EDR or Level-2 data) [3,4]. MSL12 was developed for producing consistent ocean color products using the same ocean color data processing system for multiple satellite ocean color sensors [37–39]. In addition to the standard NIR-based atmospheric correction algorithm [8], the current version of MSL12 has been improved to include the SWIR- and NIR-SWIR-based atmospheric correction algorithms [16,19,20], including a recent algorithm improvement using the information from the short blue band, i.e., M1 in Table 1 [40], for improved satellite ocean color data over coastal and inland waters. Indeed, advantages of the SWIR-based ocean color data processing over global highly turbid coastal and inland waters for various applications have been demonstrated in several previous studies [41–43].


**Table 1.** The VIIRS reflective solar band (RSB) spectral band nominal center wavelength (in nm) for SNPP and NOAA-20.

The VIIRS SNPP and NOAA-20 global Level-3 ocean color product data are generated using the spatial and temporal binning from the corresponding Level-2 data [44]. In this study, VIIRS SNPP and NOAA-20 global daily Level-3 binned data of 9-km spatial resolution from 19 June to 18 July 2018 are produced. In the Level-3 data processing, the grid elements (9 × 9 km2 bins) are arranged in rows beginning at the South Pole. Each row begins at 180◦ longitude and circumscribes the Earth at a given latitude. For 9-km spatial resolution Level-3 data, there are only three bins in the first row near the South Pole, and the number of bins in each row gradually increases with latitude in the southern hemisphere. The maximum number of bins in a row is reached at the equator, and then the number of bins decreases with latitude in the northern hemisphere. Pixels containing valid Level-2 data are mapped to bins of 9 × 9 km2. Within each bin, statistics of mean or median are accumulated for periods of one day for daily ocean color products. Before the binning process, several Level-2 masks and data quality flags from MSL12 (e.g., land, cloud [45], stray-light [46], high sun glint [47], high sensor-zenith angle, high solar-zenith angle, etc.) are applied to VIIRS ocean color Level-2 data. As examples, Figure 1a,b shows the global Level-3 daily Chl-a images of 21 June 2018 for SNPP and NOAA-20, respectively. Obviously, there are many missing pixels in the original Level-3 VIIRS Chl-a

images. In fact, there are ~70% missing pixels in these global daily images due to cloud cover, high sun glint contamination, high solar- or sensor-zenith angles, and various other reasons.

**Figure 1.** The Visible Infrared Imaging Radiometer Suite (VIIRS)-measured global Level-3 Chl-a image on 21 June 2018 from (**a**) Suomi National Polar-Orbiting Partnership (SNPP), (**b**) National Oceanic and Atmospheric Administration (NOAA)-20, (**c**) merged VIIRS SNPP/NOAA-20 Chl-a image on the same day, and (**d**) gap-filled VIIRS SNPP/NOAA-20 merged Chl-a image on the same day.

#### *2.2. Merging VIIRS SNPP and NOAA-20 Global Level-3 Ocean Color Data*

As noted previously, SNPP and NOAA-20 operate at the same sun-synchronous polar orbit of 824-km. However, the positions of the two satellites in the orbit are intentionally arranged so that NOAA-20 leads SNPP by a half orbit, or ~50 minutes. The 50-minute difference between their paths makes the NOAA-20's nadir points always running through SNPP's gaps due to high sensor-zenith angle, and vice versa. This special arrangement of satellite positions in the orbit significantly benefits the observation of the ocean color from the space: the overlap of the spatial coverages of the two sensors fills most of each other's gaps in images due to high sensor-zenith angles and high sun glint contamination, and significantly reduces the missing pixels in the merged images. In addition, ocean biological features are assumed to have little change in ~50 minutes, so that the ocean color data from the two sensors can be directly merged.

The VIIRS SNPP and NOAA-20 instruments both have 14 reflective solar bands (RSBs) with three image bands (I1–I3) and 11 moderate bands (M1–M11) operating in the spectral region of 0.41–2.25 μm. There are slight differences between VIIRS-SNPP and VIIRS-NOAA-20 in the nominal center wavelengths of each band. Table 1 lists the VIIRS RSB spectral band nominal center wavelength for VIIRS-SNPP and VIIRS-NOAA-20. In the ocean color data processing, the VIIRS-NOAA-20 *nLw*(λ) of five bands (M1–M5) are adjusted to match up with those of VIIRS-SNPP, so that Chl-a, *Kd*(490), and *Kd*(PAR) products are equivalent from the two sensors, and can be merged. In this study, only Chl-a products are merged as a preliminary test.

The VIIRS SNPP and NOAA-20 have the same swath width of 3040 km, and horizontal resolution of 750 m for M bands, and 375 m for I bands at nadir. With the same spatial and temporal resolutions, and their ocean color products processed using the same MSL12 software package, Chl-a data from

the two instruments are statistically equivalent. Figure 2 shows the scatter plots of VIIRS NOAA-20 versus SNPP Chl-a data on 21 June 2018. The mean ratios of VIIRS NOAA-20/SNPP are 1.069, 1.024, 1.057, and 1.127 in global ocean, global oligotrophic waters (Chl-a < 0.1 mg/m3), global deep waters (water depth > 1 km), and global coastal and inland waters (water depth ≤ 1 km), respectively; and the corresponding standard deviation (SD) values are 0.212, 0.162, 0.175, and 0.336, respectively. Therefore, the statistics of the ocean color products from the two sensors are very close to each other, and they can be directly merged. In this study, the merged VIIRS SNPP/NOAA-20 ocean color products are calculated as weighted averages, where the weight is the number of valid pixels in each bin. The input datasets of the merging process are 30 days of daily SNPP and NOAA-20 global 9-km resolution Level-3 binned Chl-a data from June 19 to July 18, 2018, and the output data are merged SNPP/NOAA-20 global 9-km resolution data in Level-3 binned format.

**Figure 2.** Scatter plots of VIIRS NOAA-20 versus SNPP Chl-a in (**a**) global all regions, (**b**) global oligotrophic waters, (**c**) global deep waters, and (**d**) global coastal and inland waters on 21 June 2018.

#### *2.3. Gap-Filling of the Merged VIIRS SNPP/NOAA-20 Data*

The DINEOF method [23,24] is an Empirical Orthogonal Function (EOF)-based technique, which allows the extraction of the dominant spatial patterns observed in a data time series through an iterative approach, while simultaneously filling in the missing data. During the DINEOF process, the original dataset is first stored in a spatio-temporal matrix with *m* × *n* dimensions, where *m* is the number of grids in the spatial domain and *n* is the number of time steps in the time series. Initially, the temporal and spatial mean is removed from the data, and all missing values are set to zeroes. The first EOF mode is then calculated by using the singular value decomposition (SVD) technique, and the missing values are replaced with the initial guess by the data reconstruction using the spatial and temporal functions of only the first EOF mode. The first EOF mode is then recalculated iteratively using the previous best guess as the initial value of the missing data for the subsequent iteration until the process converges. Subsequently, the number of EOFs increases one by one and for each EOF mode, the whole reconstruction is operated again until convergence. Then, by using a cross-validation technique, the optimal number of EOF modes can be determined. Finally, the reconstruction procedure is performed again, based on the optimal EOF modes, until convergence is reached. The process to determine the optimum number of EOF modes in the final reconstruction is fully automatic. For example, if the error (reconstructed−original) of the validation pixels decreases

gradually from mode 1 to mode 10, but the error starts to increase gradually from mode 11 to mode 13, then the first 10 modes are considered as optimum.

In this study, the DINEOF is utilized to reconstruct (gap-filled) the missing pixels in the merged VIIRS SNPP/NOAA-20 Chl-a data. Specifically, DINEOF is applied to 30 days of daily merged global Level-3 Chl-a data of 9-km spatial resolution from 19 June to 18 July 2018. It is noted that DINEOF is applied directly to the VIIRS Level-3 bin data files, rather than the mapped data files. However, it is a very large dataset for 30 days of global Chl-a data with 9-km spatial resolution, which makes the performance of DINEOF analysis very inefficient. To speed up the DINEOF process, the same procedure used for gap-filling of SNPP ocean color data by Liu and Wang [36] is used. Basically, the global dataset is evenly divided into 16 10◦-zonal sections between 80◦S and 80◦N, and DINEOF is applied to each of the zonal sections separately. In fact, DINEOF analyses on the 16 zonal sections are performed in parallel to further improve data processing efficiency.

#### **3. Results**

#### *3.1. Merged VIIRS SNPP/NOAA-20 Products*

The VIIRS SNPP and NOAA-20 global daily Level-3 Chl-a data of 9-km spatial resolution, from 19 June to 18 July 2018, were merged. As an example, Figure 1c shows the merged SNPP/NOAA-20 Chl-a concentration on 21 June 2018. The spatial coverage of the merged ocean color data was significantly improved compared with either SNPP (Figure 1a) or NOAA-20 (Figure 1b). Particularly, in Figure 1a,b, the gaps of high sensor-zenith angle occurred in both the northern and southern hemisphere. The gaps of sun glint, however, occurred only to the north of the equator because it was summer (June–July) in the northern hemisphere. In Figure 1c, gaps of missing pixels due to the high sensor-zenith angle and high sun glint contamination from both VIIRS sensors were almost filled, thanks to the overlap of the swath of the two sensors. The remaining missing pixels were due to a very small part of high sensor-zenith angles from the one sensor and a small part of high sun glint contamination from the other. They were significantly reduced. It can be clearly seen in Figure 1c that the gaps of the missing pixels due to the high sensor-zenith angles and high sun glint contamination were much narrower than those in Figure 1a,b. In the southern hemisphere, however, it was completely free of gaps from sun glint contamination and high sensor-zenith angle. In addition, the missing pixels due to cloud cover were also reduced in the merged Chl-a image, because of the shifts of the cloud pixels during the 50-minute period. Overall, the merged VIIRS SNPP/NOAA-20 global Level-3 daily data had ~38% more valid pixels than those from the SNPP or NOAA-20 data alone.

The major gaps of missing pixels in the merged VIIRS SNPP/NOAA-20 data were due to cloud covers. Since there was only a 50-minute delay between the paths of the two satellites, cloud usually had little change in its shape, but could travel a short distance. The merging only reduced some limited pixels on the edges of the cloud cover, but the general pattern of gaps due to cloud cover was still the same in the merged data as in the SNPP or NOAA-20 data. Indeed, the merging could not improve spatial coverage in regions with large areas of thick clouds, such as in the Arabian Sea and Bay of Bengal, equatorial Atlantic Ocean, and high-latitude Pacific and Atlantic Oceans. The summer monsoon usually brings large amounts of moisture and rainfall to the Arabian Sea from June to July.

#### *3.2. Gap-Filled VIIRS SNPP/NOAA-20 Products*

The DINEOF method was applied to the merged VIIRS SNPP/NOAA-20 global daily Level-3 Chl-a data from 19 June to 18 July 2018. The DINEOF has two options to reconstruct missing data: fully reconstructed images and filled images. A fully reconstructed image was calculated from the retained EOF modes using the DINEOF method. In the fully reconstructed image, all ocean color data were reconstructed on every ocean pixel (including non-missing pixels). The reconstructed image had no-gaps spatially. However, there were some small differences between reconstructed and original data even for non-missing pixels due to truncated EOF modes in computing all data

values. The filled image was a combination of the original image and the reconstructed image, i.e., missing pixels were filled with reconstructed data and original data were kept for non-missing pixels. As an example, the fully reconstructed global Chl-a image of 21 June 2018 is shown in Figure 1d. All gaps of missing pixels due to the high sun glint and large sensor-zenith angles were filled with valid pixels. The gaps of cloud-covered pixels were also reconstructed, especially in the northern Atlantic and Pacific, Arabian Sea and Bay of Bengal, and the equatorial Atlantic Ocean as well. The gaps of missing pixels were all filled very smoothly in Figure 1d. However, the large area of missing pixels in the Southern Ocean close to Antarctica were not reconstructed, since there were no valid pixels in the entire area from 19 June to 18 July 2018 due to high solar-zenith angles during the northern hemisphere summer time.

To evaluate the accuracy of the value from reconstructed pixels, 5% of valid (non-missing) pixels in the original Level-3 data were randomly selected, and intentionally treated as "missing pixels." These pixels were reconstructed and compared with the original data for validation. Figure 3 shows the density scatter plots of the reconstructed data versus original data in the global ocean, oligotrophic waters, deep waters, and coastal and inland waters. In general, most points are close to the 1:1 line in the global ocean, but the oligotrophic waters show the best results (Figure 3b), whilst the highest level of scatter is found in the coastal and inland waters (Figure 3d). Quantitatively, the mean and SD of the reconstructed/original ratios in oligotrophic water are 1.012 and 0.164, respectively. The mean and SD in deep waters are 1.015 and 0.182, respectively, which are slightly higher than those in oligotrophic waters. In the coastal and inland waters, the mean and SD are 0.997 and 0.287, respectively. Overall, for all pixels in the global ocean, the mean and SD of the reconstructed/original ratios are 1.012 and 0.200, respectively. It is noted that the locations of the validation pixels are selected randomly using the random number generator in the Interactive Data Language (IDL). In addition, the DINEOF package includes a built-in function to specify the "size of cloud surface" for cross-validation as a parameter of input. We tested the cross-validation with different cloud patch sizes and found that the cross-validation result does not change significantly with the size of the cloud surface. The main reason is that DINEOF considers not only spatial coherence, but also temporal coherence. Since the location of the cloud patches are selected randomly, even with the large size of cloud patch, the temporal variations in the time-series can still be used to infer the missing values. The advantage of the DINEOF over traditional interpolation is that it utilizes major EOF modes to reconstructs the missing pixels, which capture the major signals of variations in both spatial and temporal domains simultaneously.

#### *3.3. Ocean Features Revealed in the Gap-Filled VIIRS SNPP/NOAA-20 Chl-a Data*

Gap-filled daily Chl-a images reveal many large-scale and meso-scale ocean features that are invisible in the merged Chl-a images, or the original VIIRS SNPP, NOAA-20 Chl-a images. As shown in Figure 1d, the most obvious large-scale ocean features are the oligotrophic waters (Chl-a < 0.1 mg/m3) in the center of the five subtropical ocean gyres, i.e., North Atlantic, South Atlantic, North Pacific, South Pacific, and South Indian Ocean, which cover major parts of the global oceans. Because of the easterly winds near the equator and the equatorial upwelling, there is more enhanced Chl-a concentration (~0.1–0.5 mg/m3) in the equatorial Pacific Ocean and equatorial Atlantic Ocean regions than in the subtropical gyres. However, no significant Chl-a enhancement is found in the equatorial Indian Ocean. In addition, the Gulf Stream in the North Atlantic Ocean and Kuroshio in the North Pacific Ocean mark a clear boundary in the high latitude ocean regions: high Chl-a concentrations are found to the north of the boundary, while low Chl-a to the south.

**Figure 3.** Density scatter plots of gap-filled versus merged Chl-a values for (**a**) all selected validation pixels, (**b**) oligotrophic waters, (**c**) deep waters, and (**d**) coastal and inland waters on 21 June 21 2018.

In addition to these large-scale features, the gap-filled images also reveal many meso-scale ocean features, such as eddies and filaments associated with western boundary currents or coastal currents. Figure 4 shows the transformation of meso-scale eddies in the Gulf of Mexico (GOM) from 19 June to 18 July 2018. Ocean circulation in the GOM is dominated by the loop current (LC), an extension of the western boundary current system of the North Atlantic Ocean that loops into the GOM. The LC enters the Gulf through the Yucatan Channel and exits through the Florida Straits. The LC occasionally extends further northward into the GOM, approaching the northern shelf break. This long extension of the loop will eventually separate to form an anticyclonic loop current eddy (LCE). Formation of an LCE by separation from the LC happens irregularly every several months, and there can be a number of LCEs in the GOM at one time. Figure 4 shows three LCEs in the GOM: two of them (A and B) were aged LCEs that already moved to the west side of GOM, and one (C) was a new LCE that just separated from the LC. Loop current eddy A and B are about 100 km in diameter, and C is about 200 km in diameter. During the one-month period from 19 June to 18 July 2018, A and B were still moving continuously to the west, and A was showing some changes in its shape while interacting with the LC. The LC was very flat till 7 July, when the LC started to bulge northward, getting ready to shed a new LCE. The details of the transformation of these meso-scale ocean features were not available in either SNPP, NOAA-20, or in the merged images, but only revealed in the gap-filled images. The transitions of these meso-scale ocean features were very smooth both temporally and spatially. To compare with the original VIIRS images without gap-fillings, Chl-a global daily images from SNPP, NOAA-20, or merged SNPP/NOAA-20 can be found in the NOAA Ocean Color Team website (https://www.star.nesdis.noaa.gov/sod/mecb/color/), in particular, using the powerful image display tool, i.e., the Ocean Color Viewer (OCView) [48]. With OCView, the user can zoom in to the GOM region (or any other region) with high spatial resolution (2 km) images, and select SNPP, NOAA-20 or the merged SNPP/NOAA-20 images.

**Figure 4.** Daily Chl-a images of the Gulf of Mexico (GOM) from 19 June to 13 July 2018 with three loop current eddies (**A**, **B**, and **C**) marked on the 19 June 2018 Chl-a image.

#### *3.4. Comparison with Gap-Filled Data Based on VIIRS SNPP or NOAA-20*

The DINEOF was also applied to Chl-a data from a single sensor, VIIRS-SNPP or VIIRS-NOAA-20, and the gap-filled data based on a single sensor were compared with gap-filled data based on the merged products. Figure 5 shows the gap-filled data based on VIIRS-SNPP, VIIRS-NOAA-20, and merged VIIRS SNPP/NOAA-20 products in the GOM on 21 June, 30 June, and 18 July 2018, as examples. It can be seen that, in general, the gap-filled data based on the merged products show the same pattern as those based on SNPP or NOAA-20 alone. However, there are more details in the dynamic features in the gap-filled data based on the merged product than those based on SNPP or NOAA-20 alone, and the interactions of the loop current with surrounding eddies are different in the merged products. On 21 June, the gap-filled data based on merged VIIRS SNPP/NOAA-20 products (Figure 5c) show more details in the coastal eddy feature "A" than VIIRS-SNPP (Figure 5a) or VIIRS-NOAA-20 (Figure 5b) alone. On 30 June 2018, the feature of the big eddy "B" are different in the gap-filled data based on merged products (Figure 5f). From 30 June to 18 July, the eddy "B" was pushed northward by the loop current, and transformed from a north–south oriented shape to a west–east oriented shape. In the animation of the daily images (not shown here), the interactions of the eddy "B" with the loop current were more clearly shown in the merged products. On 18 July 2018, the eddy "C" was seen more clearly in Figure 5i than that in Figure 5g or Figure 5h. Generally, the merged data have more spatial coverage than that from either SNPP or NOAA-20 alone. Indeed, some ocean dynamic features could be partially blocked by pixels with high sensor-zenith angle, sun glint, or cloud in VIIRS SNPP or NOAA-20 images. But, merging data from the two sensors can significantly increase valid pixels in the coverage. Consequently, the gap-filled ocean color products based on the merged data set can capture more details of the ocean dynamic features.

**Figure 5.** Comparison of the VIIRS-derived ocean features in the GOM in the gap-filled images based on SNPP (left column), NOAA-20 (middle column), and SNPP/NOAA-20 merged (right column) data in 2018 on 21 June (top row) (**a**–**c**), 30 June (middle row) (**d**–**f**), and 18 July (bottom row) (**g**–**i**).

#### **4. Discussion and Conclusions**

In a recent study [36], we used the DINEOF to fill the gaps of missing pixels in daily VIIRS-SNPP ocean color data, and it was demonstrated that the DINEOF method can successfully reconstruct and reveal meso-scale and large-scale spatial ocean features in the global VIIRS Level-3 images, as well as capture the temporal variations of these features. Based on the same methodology, in this study, the DINEOF was applied on the daily merged VIIRS SNPP/NOAA-20 global Level-3 ocean Chl-a data. The VIIRS-SNPP and VIIRS-NOAA-20 have very similar sensor characteristics, and the statistics of the ocean color data derived from the two sensors are very close to each other. In addition, global VIIRS SNPP and NOAA-20 ocean color data are derived using the same MSL12 ocean color data processing system. As shown in examples, merging VIIRS-SNPP and VIIRS-NOAA-20 data almost removes the gaps of missing pixels due to high sensor-zenith angles and high sun glint contamination, and also significantly reduces the gaps due to cloud cover. The remaining regular gaps in the merged images are due to small gaps from high sensor-zenith angles in one sensor and sun glint contamination in another. Overall, there are ~38% more valid ocean color pixels in the merged VIIRS SNPP/NOAA 20 global Level-3 daily data than those in the SNPP or NOAA-20 data alone.

In the gap-filled Chl-a data, all missing pixels due to the high sensor-zenith angle and high sun glint were mostly reconstructed. The gaps of cloud-covered pixels were also reconstructed, especially in regions with large areas of thick clouds, such as in the northern Atlantic and Pacific, Arabian Sea and Bay of Bengal, and the equatorial Atlantic Ocean as well. Quantitatively, the mean ratios of the reconstructed/original data are 1.012, 1.012, 1.015, and 0.997 for the global ocean, global oligotrophic waters, global deep waters, and global coastal and inland waters, respectively; and the corresponding SD values are 0.200, 0.164, 0.182, and 0.287, respectively. Gap-filled daily Chl-a images reveal many large-scale and meso-scale ocean features and their variations, which are invisible in the original VIIRS-SNPP or VIIRS-NOAA-20 Chl-a images.

The gap-filled Chl-a data based on two-sensor merged products were compared with gap-filled data based on single sensor VIIRS-SNPP or VIIRS-NOAA-20 alone. It is demonstrated that the gap-filled data based on the merged products show more details in the ocean features than those based on a single sensor alone, and there are differences in the dynamics of ocean features. This is an important message that adding more sensors into the merged products will significantly improve the quality of gap-filled ocean color data. Previous missions such as the Sea-viewing Wide Field-of-view Sensor (SeaWiFS) (1997–2010), MODIS on the Terra (1999–present) and Aqua (2002–present) satellites, and the Medium-Resolution Imaging Spectrometer (MERIS) on the ENVISAT (2002–2012) had provided long-term high-quality global ocean color data to the community. Currently, in addition to the VIIRS on SNPP and NOAA-20, the Ocean and Land Color Instrument (OLCI) on the Sentinel-3A (2016–present) and Sentinel-3B (2018–present) satellites, and the Second-Generation Global Imager (SGLI) on the Global Change Observation Mission-Climate (GCOM-C) (2017–present) are all providing an unprecedented view of ocean optical, biological, and biogeochemical properties on a global scale simultaneously. Adding these sensors in the data merging process will significantly improve the spatial coverage in the merged global ocean color data, and the quality of gap-filled ocean color data as well.

**Author Contributions:** X.L. carried out the main research work for developing algorithm, obtaining the results, and analyzing the data. M.W. suggested the topic and contributed to the algorithm development and evaluation.

**Funding:** This work was supported by the Joint Polar Satellite System (JPSS) funding and NOAA Product Development, Readiness, and Application (PDRA)/Ocean Remote Sensing (ORS) Program funding.

**Acknowledgments:** We thank two anonymous reviewers for their useful comments. The views, opinions, and findings contained in this paper are those of the authors and should not be construed as an official NOAA or U.S. Government position, policy, or decision.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
