An Empirical Radiometric Intercomparison Methodology Based on Global Simultaneous Nadir Overpasses Applied to Landsat 8 and Sentinel-2

Abstract

The Simultaneous Nadir Overpass (SNO) method was developed by the NOAA/NESDIS to improve the consistency and quality of climate data acquired by different meteorological satellites. Taking advantage of the reduced impact induced by the Bidirectional Reflectance Distribution Function (BRDF), atmospheric effects, illumination and viewing geometries during an SNO, we created a sensor comparison methodology for all spectral targets. The method is illustrated by applying it to the assessment of data acquired by the Landsat 8 (L8), Sentinel-2A (S2A), and Sentinel-2B (S2B) optical sensors. Multiple SNOs were identified and selected without the need for orbit propagators. Then, by locating spatially homogeneous areas, it was possible to assess, for a wide range of Top-of-Atmosphere reflectance values, the relationship between the L8 bands and the corresponding ones of S2A and S2B. The results yield high coefficients of determination for S2 A/B with respect to L8. All are higher than 0.980 for S2A and 0.984 for S2B. If the S2 band 8 (wide near-infrared, NIR) is excluded then the lowest coefficients of determination become 0.997 and 0.999 from S2A and S2B, respectively. This methodology can be complementary to those based on Pseudo-Invariant Calibration Sites (PICS) due to its simplicity, highly correlated results and the wide range of compared reflectances and spectral targets.

1. Introduction

The growing number of Earth Observation (EO) satellites reflects the societal demand of products and services based on remote sensing data [1]. To address the questions of how and to what extent different data sources can be used for a given product or service, and how to perform its harmonization, it is necessary to execute an assessment of data interoperability [2]. Calibration and cross-calibration processes are essential to ensure data quality, functionality and interoperability [3,4]. Post-launch calibration and validation activities are important for all optical sensors and essential for those which lack on-board calibration systems [5].

Pseudo-Invariant Calibration Sites (PICS) based methodologies have been widely used for cross-calibration [6]. These methodologies require Bidirectional Reflectance Distribution Function (BRDF) modelization [7,8] and are influenced by the atmosphere [9,10], while their results are restricted to a limited range of reflectances [11,12].

One of the most used techniques for cross-calibration and inter-calibration is the Simultaneous Nadir Overpass (SNO) methodology [13,14,15,16,17,18] developed by the National Environmental Satellite, Data, and Information Service (NESDIS) of the National Oceanic and Atmospheric Administration (NOAA) more than a decade ago [19].

SNO-based methodologies consist on the direct inter-comparison and subsequent adjustment of measurements obtained by sensors that fly on-board satellites when these sensors observe the same area of the Earth’s surface in a very short time interval between acquisitions (Inter-Acquisition Time Gap, IATG) [3,14,20]. The key assumption behind the SNO is that when two different optical sensors acquire the same area nadir pointing and simultaneously, the BRDF, illumination geometry, viewing geometry and atmospheric influence on the measurements are equal for both observations [3,14,20], allowing the direct study of the compared radiometric performance of both sensors [3,20].

With the purpose of minimizing the aforementioned drawbacks in the PICS methodology and to simplify the process, we introduce a methodology based on multiple Simultaneous Nadir Overpasses (SNOs). In the frame of this work, we created a sensor comparison methodology for all spectral targets. We then used it to analyze the radiometric interoperability of data provided by the optical sensors on board the Sentinel-2A (S2A), Sentinel-2B (S2B) and Landsat 8 (L8) satellites. These sensors were chosen, due to their similarity in terms of spectral bands, spatial resolutions, and the interest of the earth-observation community in making them interoperable [4,21,22]. Firstly, this article presents a method to identify SNOs without using orbit propagators. Secondly, a technique is proposed to find spatially homogeneous areas to extract radiometric statistics. These data are analyzed and refined to discuss sensor interoperability and harmonization based on the compared Top-Of-Atmosphere (TOA) reflectances. The conclusions highlight the good radiometric correlation between L8 and S2 and the methodology consistency when compared with a PICS-based one.

2. Materials and Methods

The study time frame ranged from 2015-07-04 to 2019-11-05 with a worldwide scope. Most of the metadata needed, like image geometries, acquisition times and ground tracks, were available on the United States Geological Survey (USGS) (https://earthexplorer.usgs.gov/) and Copernicus Open Access Hub (http://scihub.copernicus.eu). S2 A/B data were L1C TOA reflectance products, and for L8, the Collection 1, Tier 1, of calibrated TOA reflectance products (L1TP). Both product sets are found ortho-rectified to the Universal Transverse Mercator projection (UTM/WGS84) and available on the same websites mentioned above [23,24].

2.1. Satellite Sensors

L8 was launched on 11 February 2013, from the Vandenberg Air Force Base in a circular sun-synchronous orbit with 16 days of repeat cycle, an altitude of ~705 km and an orbit plane inclination of 98.2°. It carries the Operational Land Imager (OLI) and the Thermal Infrared Sensor (TIRS) scientific instruments which have both a 15° field of view (FOV), providing a ~185 km swath. The OLI payload has nine spectral bands in the visible to near-infrared (VNIR) and short-wave infrared (SWIR) range with a spatial resolution of 30 m for bands 1 to 7 and 9, while the resolution of the panchromatic band (band 8) is 15 m [25].

The S2 mission is part of the Copernicus Program of the European Commission for global land surface monitoring and coastal waters. Currently, the S2 program consists on a constellation of two satellites, S2A and S2B, which were launched on 23 June 2015 and 7 March 2017, respectively, from the Europe’s spaceport near Kourou in French Guiana. They are located in a sun-synchronous orbit, sharing the same orbit plane of 98.6° inclination and with a 180° orbital phasing difference at ~786 km height [26]. Both satellites are equipped with the MultiSpectral Instrument (MSI) payload which, with a 20.6° FOV provides a swath of ~290 km. The repeat cycle is 10 days for each satellite or 5 days for combined use [26]. The MSI has 13 spectral bands in the VNIR and SWIR spectral domains with spatial resolutions of 10 m for bands 2 to 4 and 8, 20 m for the bands 5, 6, 7, 8A, 11, 12, and 60 m for the bands 1, 9 and 10 [23,26].

2.2. Ground Trajectory Determination and SNOs Finding

Finding SNOs required the ground track determination of each satellite and the associated time for each point on the track to define the ground trajectories. A method was developed to estimate ground trajectories without using orbit propagators. For L8, only OLI C1 L1TP product metadata data [27] was needed, while S2A and S2B required the full-cycle 143-orbit tracks published by the ESA in addition to the L1C product metadata [28]. In order to avoid errors induced by the usage of projections, we performed the calculations on the three-dimensional WGS84 Earth model.

2.2.1. Landsat 8 Ground Trajectories

According to the L8 product metadata, more than 99.8% of the samples have a roll angle value lower than 2 × 10⁻³⁰, the rest have a maximum value of 4 × 10⁻³⁰ representing ~1.4 km and ~2.8 km of off-nadir distance on the ground, respectively [29]. We considered these distances negligible. Each scene center was, hence, deemed a point belonging to the ground track, which we built using geodesics. The scene center time was calculated as the middle time between the scene start and stop times. It was then interpolated along the geodesic to build the ground trajectory. An example is shown in Figure 1.

2.2.2. Sentinel-2 Ground Trajectories

S2A, S2B and L8 perform track maintenance [26,29,30,31], however, S2 products differed from L8 because they are tiled using the Military Grid Reference System (MGRS) [23,26]. This system has no relation to the orbits, therefore, we used the 143-orbit cycle ground tracks provided by ESA [28] for the calculation of the S2A and S2B ground trajectories. All S2 L1C products whose footprint intersected with the track of their corresponding cycle were then linked to their associated track segment. The granule sensing time metadata value of the product provided the time for the track segment. We obtained the ground trajectories for S2A and S2B by concatenating these time-labeled track segments and interpolating the time for each point on the track assuming that the granule sensing time corresponds to the first line of the granule. An S2A and S2B ground tracks example is shown in Figure 2.

2.2.3. SNOs Determination

We define the perfect SNO as the event in which one sensor is orbiting exactly over the other at the same time. Under these conditions, both sensors virtually see the same nadir point with the same viewing and illumination angles. There are, then, two requirements, one spatial (orbit intersection) and the other temporal (simultaneity).

In practice, it is necessary to add tolerance to these requirements. In this work, we began by being strict in the geometric one. The orbit intersections are our starting points, while the tolerance is provided by the IATG. Once both have been established, an additional spatial tolerance is added as a radius, creating a circular area of study centered in the orbit intersection. It is worth noting that the spatial tolerance has a stronger influence on the viewing angles, while the temporal tolerance has a more significant impact on the solar angles.

The established ground trajectories were used to determine intersection points and IATGs of L8 and S2. We initially constrained the IATG to 24 h, obtaining 1397 S2A and 953 S2B overpasses with L8. Since the data extraction of each overpass required a significant amount of computing resources, we filtered out overpasses with an IATG longer than 30 min or whose S2 sample had more than 5% cloudiness value in the metadata. In this way, we reduced the number of overpasses to 96 S2A and 60 S2B (Appendix B).

This methodology required cloud-free areas. We set the cloudiness threshold in the S2 scene because it is smaller than the L8 one, and therefore, more restrictive. The cloudiness threshold value is the result of a trade-off between the chances of having clouds in the study area and the number of SNOs. A cut-off value of 5% allows, approximately, three times more SNOs than the corresponding to 0%.

Once the intersection point and IATG were determined, the remaining element to complete our SNO identification was the radius. We used the radiometric data as they were obtained from the official providers without any kind of manipulation (e.g., resampling). For this reason, an SNO comprised a single L8 tile and a single S2 tile, therefore, the radius was calculated as the minimum distance from the intersection point to the polygon resulting from the intersection of both tiles. This radius was further reduced when necessary to avoid large water bodies, as shown in Figure 3. The study area is a circle centered in the ground track intersection point. This keeps the methodology compatible with satellites with opposed orbit directions.

Figure 4 and Figure 5 show the world distribution of the SNOs areas used in our study. The circle size has been overscaled to ease the visualization.

2.3. Statistics Extraction

The coefficient of variation (CV) is an operator used as an indicator of the radiometric spatial uniformity [32]. The CV is defined as the ratio between the standard deviation and the average of values (Equation (1)).

$$CV=\frac{\sigma }{\left|\overline{x}\right|}$$

(1)

where $\left|\overline{x}\right|$ is the average value and σ its standard deviation.

Using the CV, an algorithm searched for spatially homogeneous areas [32,33] on each S2 band, and extracted TOA reflectance statistics for each individual area. The same areas were used to extract statistics on the corresponding bands of the L8 dataset. The distance to nadir, distance to the intersection of the ground trajectories, area, sun elevation and sun azimuth on the area centroid were also noted.

The bands used were all S2 VNIR and SWIR provided by the MSI with the exception of the water vapor (band 9), which was excluded due to the lack of a corresponding L8 OLI band [34]. The S2 cirrus band (band 10) was also excluded because the band does not provide useful ground data [35].

Although it could have been possible to several mosaic datasets to obtain larger SNO areas, it would have required merging, and in some cases, resampling. Our approach was strict into avoiding raster data manipulation; therefore, each SNO comprises only a pair of datasets: One L8 L1TP and one S2A or S2B L1C.

For a given SNO, each pair of the selected bands was processed individually (Figure 6).

2.3.1. Homogeneous Areas Creation

The S2 band was first convolved using the CV operator to obtain an image representing its spatial variability. The CV value corresponding to the percentile 1 of the sample was set as a threshold to binarize the CV image previously calculated. The pixels with values lower than the threshold were set to 1, and the remaining were set to 0. The erode operator [36,37] was then applied to reduce the size of the pixel clusters and to remove the smallest. The remaining clusters partially recovered their surface after applying a dilate operator [36,37] with a smaller kernel. Finally, the clusters were vectorized creating what we defined as Homogeneous Areas (HAs). For each SNO and S2 band, a set of Has was created.

We assumed that a HA in a higher resolution image is also homogeneous in a lower resolution one. For this reason, we choose the S2 bands for the Has creation since, with the exception of the coastal aerosol, the S2 bands have higher resolution than the L8. The use of a threshold from a percentile grants data from all samples. We reached the percentile value by trial and error, probing between 0.1 and 5.0. It was observed that the results had low sensitivity to this parameter. We expected the Has to minimize the possible effects of the resolution differences between the S2 and L8 bands, as well as small geometric errors. For the latter reason, the erode is stronger than the dilate; it creates a geometric margin-of-error distance for the HA. An example of the generation of homogeneous zones is shown in Figure 7 and Figure 8.

2.3.2. Statistics Retrieval

We noted the maximum, minimum, average and standard deviation of each pixel Digital Number (DN) and HA of the corresponding S2 and L8 bands. Following the restriction of maintaining the original raster values, the HA polygons, and not the images, were reprojected when the pair of SNO images were found in different projections. HAs smaller than nine OLI MS pixels (900 m²) were removed from the study. Additionally, sun angles, pass-time, locations, centroid distance to nadir and intersection point were recorded. Each HA was an independent data point regardless of the SNO to which it belonged. They were kept separately, so the spectral differences stayed separate.

2.4. Data Analysis

The S2 TOA reflectance was calculated applying the rescaling factor provided in the metadata to the S2 DNs, since the solar elevation across the scene has been already taken into account [38]. Following the product guide, we obtained the L8 TOA reflectance by applying the rescaling factors included in the metadata to the L8 DNs. Subsequently, these values were divided by the sine of the solar elevation at the HA centroid location at the L8 scene center time [39].

The L8 TOA reflectances were plotted versus the S2A and S2B TOA reflectances individually for each pair of corresponding bands, and a linear regression was calculated. Each SNO was represented with a different color, allowing the location of outliers and the identification of SNOs requiring visual inspection (Figure 9). The IATGs were, in general, too long to assume that cloudy areas were radiometrically stable. Shadow areas are solely illuminated by the atmosphere. For these reasons, cloudy and shadow areas were not used in this study (Figure 10). Scenes containing water flows and those containing perceptible atmospheric effects were also radiometrically dynamic (Figure 11), causing the presence of outliers. Only data points with justified anomalies were removed. The list of SNOs with valid HAs can be found in Appendix B.

The number of generated HAs largely differed depending on the SNO area texture. In order to balance their weight, we limited the number of HAs per SNO selecting the ones with the lowest standard deviation. No improvement was observed; therefore, the limitation was not used.

In the SNO used to illustrate B and C of Figure 11, L8 crossed over the intersection point 136 s after S2A. In cases like these, where anomalies are caused by the atmosphere, the entire SNO (all bands) is removed from the study.

A general schematic of the workflow followed for this study is shown in Figure 12.

3. Results

3.1. Data Analysis Remarks

During the data analysis, we found that anomalies are very noticeable in the linear regression figures. As an example, Figure 13 shows how a smoke plume creates a salient outlier. Another interesting finding was the good behavior of the snow-covered slopes in different angles. Figure 14 shows HAs located on such surfaces whose reflectances stay aligned in the graph. It is worth noting that the reflectivity values are much higher than one, indicating the presence of the expected snow BRDF anisotropy.

3.2. Correlation with TOA Reflectances

A linear regression was applied to each band pair for S2 A/B and L8 TOA reflectance sets, obtaining the results presented in

Table 1. Linear regression values for each band pair for L8 and S2A.

L8 vs. S2A
S2 Band	L8 Band	Band Name	Slope	Intercept	Coefficient of Determination
1	1	Coastal Aerosol	0.99111	−9.307 × 10−4	0.99993
2	2	Blue	0.99267	3.248 × 10−4	0.99907
3	3	Green	1.01095	−1.441 × 10−3	0.99969
4	4	Red	0.97671	−2.016 × 10−3	0.99874
8A	5	Narrow NIR	0.99087	−3.732 × 10−4	0.99953
8	5	Wide NIR	1.02380	1.470 × 10−2	0.99115
11	6	SWIR 1	0.98953	−2.443 × 10−3	0.99956
12	7	SWIR 2	1.00041	−3.116 × 10−3	0.99964

and

Table 2. Linear regression values for each band pair for L8 and S2B.

L8 vs. S2B
S2 Band	L8 Band	Band Name	Slope	Intercept	Coefficient of Determination
1	1	Coastal Aerosol	1.00396	−9.729 × 10−4	0.99991
2	2	Blue	1.01424	−8.120 × 10−4	0.99974
3	3	Green	1.02565	−3.433 × 10−3	0.99973
4	4	Red	1.00034	−3.760 × 10−3	0.99941
8A	5	Narrow NIR	1.00989	−2.867 × 10−3	0.99942
8	5	Wide NIR	1.01025	1.909 × 10−2	0.99381
11	6	SWIR 1	0.99641	−3.151 × 10−3	0.99973
12	7	SWIR 2	0.98597	2.116 × 10−4	0.99972

. Appendix A shows the graphs of the linear regression for each band pair for S2 A/B and L8. The standard deviation of the residuals was the measure of the uncertainty used to create confidence intervals.

With the exception of the wide NIR band, all regressions produce offsets whose absolute values are lower than 0.38% with high coefficients of determination. We assumed then that all TOA reflectance products yield zero for zero signal, so we removed the degree of freedom corresponding to the intercept and repeated the process to obtain the final results presented in

Table 3. Linear regression values for each band pair (L8 and S2A) with zero intercept.

L8 vs. S2A
S2 Band	L8 Band	Band Name	Slope	Coefficient of Determination
1	1	Coastal Aerosol	0.98911	0.99985
2	2	Blue	0.99402	0.99814
3	3	Green	1.00540	0.99931
4	4	Red	0.97124	0.99741
8A	5	Narrow NIR	0.98999	0.99907
8	5	Wide NIR	1.06036	0.97962
11	6	SWIR 1	0.98441	0.99906
12	7	SWIR 2	0.99299	0.99914

and

Table 4. Linear regression values for each band pair (L8 and S2B) with zero intercept.

L8 vs. S2B
S2 Band	L8 Band	Band Name	Slope	Coefficient of Determination
1	1	Coastal Aerosol	1.00188	0.99980
2	2	Blue	1.01277	0.99948
3	3	Green	1.01872	0.99932
4	4	Red	0.99326	0.99867
8A	5	Narrow NIR	1.00428	0.99877
8	5	Wide NIR	1.04793	0.98434
11	6	SWIR 1	0.98985	0.99935
12	7	SWIR 2	0.98645	0.99943

. The corresponding graphs can be found in Appendix A.

Excluding again the wide NIR band we found coefficients of determination higher than 0.997 and 0.999 for S2A and S2B, respectively.

3.3. Dependence from Other Variables

The linear fitting allowed us to analyze the correlation of the residual of each HA with the following variables:

• Average reflectance;
• Reflectance standard deviation;
• Solar elevation at the HA centroid;
• Solar azimuth at the HA centroid;
• HA centroid distance to nadir;
• HA latitude.

The same analysis was performed using the aggregated residuals of complete SNOs with the acquisition date and IATGs. No correlation was found between the residuals and any of the listed variables. For this reason, no restriction was set on the radius or sun angles. The radius is solely limited by the restrictions imposed by the SNO definition itself (Section 2.2.3). Being OLI and MSI pushbroom, the HA to nadir distances account for the across-track viewing angle. The lack of correlation between residuals and distances is the reason why we left the distance (or across-track angle) unconstrained in our dataset. Four examples are shown in Figure 15.

Figure 16 shows the distribution of HA distance to nadir differences between L8 and S2 A/B for each HA. It is worth noting that it is centered, symmetric and bell-shaped. The average is ~2.1 m with a standard deviation of 199.2 m, approximately 4.5 × 10⁻²⁰ of across-track viewing angle difference error (±3σ), that we can consider negligible.

No HA was seen with an across-track angle larger than 3.9°. For the purpose of this work, the viewing zenith angles can be considered virtually the same as the across-track angles. Similarly, the viewing azimuth angle difference is virtually the difference between orbit plane inclinations, ~0.4° in this case, recalling that S2A and S2B share the same orbit plane. It is worth highlighting that the difference in viewing azimuth angles becomes irrelevant as the HA distance to nadir diminishes.

3.4. Ground Classes Distribution

Since SNOs can be located anywhere in the world and HAs are randomly distributed within them, it was considered appropriate to inspect the distribution of ground classes across HAs. The Copernicus Global Land Cover Layers collection 2 at 100 m resolution (CGLS-LC100) [40] was used for this purpose. The CGLS-LC100 has the following features [40]:

• Derived from PROBA-V satellite observations for the 2015 reference year;
• Discrete classification with 23 classes;
• 100 m spatial resolution;
• An overall 80% accuracy.

Table 5. Legend for the 23 discrete classes of CGLS-LC100.

Color	Class ID	Class Name
	22	Oceans, seas. Can be either fresh or salt-water bodies.
	21	Open forest, not matching any of the other definitions.
	20	Open forest, mixed.
	19	Open forest, deciduous broadleaf. Top layer—trees 15–70% and second layer—mixed of shrubs and grassland, consists of seasonal broadleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
	18	Open forest, deciduous needle leaf. Top layer—trees 15–70% and second layer—mixed of shrubs and grassland, consists of seasonal needle leaf tree communities with an annual cycle of leaf-on and leaf-off periods.
	17	Open forest, evergreen broadleaf. Top layer—trees 15–70% and second layer—mixed of shrubs and grassland, almost all broadleaf trees remain green year-round. Canopy is never without green foliage.
	16	Open forest, evergreen needle leaf. Top layer—trees 15–70% and second layer—mixed of shrubs and grassland, almost all needle leaf trees remain green all year. Canopy is never without green foliage.
	15	Closed forest, not matching any of the other definitions.
	14	Closed forest, mixed.
	13	Closed forest, deciduous broadleaf. Tree canopy > 70%, consists of seasonal broadleaf tree communities with an annual cycle of leaf-on and leaf-off periods.
	12	Closed forest, deciduous needle leaf. Tree canopy > 70%, consists of seasonal needle leaf tree communities with an annual cycle of leaf-on and leaf-off periods.
	11	Closed forest, evergreen broadleaf. Tree canopy > 70%, almost all broadleaf trees remain green year-round. Canopy is never without green foliage.
	10	Closed forest, evergreen needle leaf. Tree canopy > 70%, almost all needle leaf trees remain green all year. Canopy is never without green foliage.
	9	Moss and lichen.
	8	Herbaceous wetland. Lands with a permanent mixture of water and herbaceous or woody vegetation. The vegetation can be present in either salt, brackish, or freshwater.
	7	Permanent water bodies. Lakes, reservoirs, and rivers. Can be either fresh or salt-water bodies.
	6	Snow and ice. Lands under snow or ice cover throughout the year.
	5	Bare/sparse vegetation. Lands with exposed soil, sand, or rocks and never has more than 10% vegetated cover during any time of the year.
	4	Urban/built up. Land covered by buildings and other manufactured structures.
	3	Cultivated and managed vegetation/agriculture. Lands covered with temporary crops followed by harvest and a bare soil period (e.g., single and multiple cropping systems). Note that perennial woody crops will be classified as the appropriate forest or shrubland cover type.
	2	Herbaceous vegetation. Plants without persistent stem or shoots above ground and lacking definite firm structure. Tree and shrub cover is less than 10%.
	1	Shrubs. Woody perennial plants with persistent and woody stems and without any defined main stem being less than 5 m tall. The shrub foliage can be either evergreen or deciduous.
	0	Unknown

specifies each CGLS-LC100 class type used in this work, as well as its corresponding color and ID. Figure 17 shows the class distribution across HAs per S2 satellite and band. It is worth recalling that HAs were extracted from MSI data.

Table 6. Linear fitting slopes for S2 band 4 and L8 band 4 for the most frequent classes.

MSI Band 4, OLI Band 4 (Red). Most Frequent Classes
Satellite	Slope	Correlation Index	Number of HAs	Class	Class Description
S2A	0.97124	0.99741	2376	all	All classes
S2B	0.99444	0.99873	1702
S2A	0.96765	0.99310	1558	5	Bare/sparse vegetation
S2B	0.97898	0.99876	655
S2A	0.97106	0.99599	142	3	Cultivated and managed vegetation/agriculture
S2B	0.98698	0.99378	221
S2A	0.97171	0.99511	471	2	Herbaceous vegetation
S2B	0.99318	0.99932	152
S2A	0.97181	0.99857	33	1	Shrubs
S2B	0.97378	0.99667	158

and Figure 18 show, as an example, the regressions of the S2A/B band 4 when they are calculated using the most frequent classes separately. For convenience, the figure and data obtained with no class discrimination are also included.

Regression figures corresponding to Table 6 can be found in Appendix C.

4. Discussion

In this study, we identified SNOs across the globe using a ground-agnostic method. As a result, for each pair S2A/L8 and S2B/L8, diverse ground types were targeted (Figure 4 and Figure 5). Unlike PICS-based methodologies where the targets are in well-defined, characterized and fixed locations [7,9,41], the proposed method prioritizes the intersection point regardless of its location on the Earth, with a limited IATG of 30 min. The method allows a strongly empirical approach, with no need for atmospheric or BRDF correction. Despite the spectral diversity of the study area, the coefficients of determination between L8 and S2 bands (wide NIR band excluded) were higher than 0.997 (S2A) and 0.999 (S2B). The SNO plus Has technique provides a wide range of reflectances for study, often from less than 0.1 to values close to 1 (Appendix A), while PICS-based methodologies are restricted to a reduced set [11,12].

The ground classes breakdown (Figure 17) shows the presence of both arid and vegetated classes in the samples. The class that includes deserts, defined as “bare/sparse vegetation”, is the most frequent. It is followed by “herbaceous vegetation”, “cultivated and managed vegetation/agriculture”, “shrubs”, and “herbaceous wetland”. The frequency distribution of these vegetated classes differs for each satellite and band. “Herbaceous vegetation” is more strongly represented in S2A, while “herbaceous wetland” is its counterpart in S2B (Figure 17). These differences can be attributable to the world distribution of SNOs for each satellite. No significant statistic impact is expected to arise from this disparate distribution of similar vegetated classes using the SNO-HA method.

Helder et al. [4] collected the S2A/MSI and absolute calibration results obtained through different models and combined them in a metric that should serve to remove any systematic errors within the models and provide for a per-model comparison [4]. We added the results provided by our method to the figure published in the aforementioned work (Figure 19).

The Copernicus Sentinel-2 Mission Performance Centre (MPC) also published ratios as radiometric validation indicators using Rayleigh atmospheric backscattering over deep ocean sites, in-situ data, spectrally characterized PICS and comparison with other sensors (S2A/MSI vs. OLI, S2B/MSI vs. OLI and S2B/MSI vs. S2A MSI) over Libya-4 [42]. Thus, Figure 20 and Figure 21 show our results compared to the ones obtained in the latter study.

Using a methodology based on SBAF applied to the PICS Algeria-3, Algeria-5, Egypt-1 and Libya-4. Barsi et al. [43] also compared S2A/B MSI with OLI. Continuing with the per-model comparison philosophy, we also added our results to the latest data published by Barsi et al. [44] in which they compared S2A/B MSI with OLI data (Figure 22 and Figure 23).

Our methodology shows that S2A yields higher reflectance values than S2B for most of the studied bands. We attribute this effect to the small and systematic difference in radiometry of about 1% detected by ESA in the Sentinel-2 L1C Data Quality Report and which it is currently under study by the MPC [42]. This discrepancy has also been previously observed by other authors, such as Helder et al. [4].

The wide NIR band shows the lowest coefficients of determination with the narrow NIR of L8, 0.980 and 0.984 for S2 A/B, respectively. The effect can also be observed as a large confidence interval in Figure 20 and Figure 21. This is consistent with the fact that their RSRs are largely different. Apart from the wide NIR, the red band of S2A shows the largest deviation. For this reason, we chose this band to focus on the ground classes breakdown. Recalling that the ratios are inverted in Figure 18, we see that S2A yields higher values in all the represented classes, which is in agreement with the effect under study by the MPC [42]. We also observe that S2B ratios have larger variability across classes when compared with S2A.

The S2A/MSI and S2B/MSI are largely interoperable with L8/OLI for all studied bands with the exception of the wide NIR, which correlates notably worse. The slopes in Table 3 and Table 4 are the corrective factors that account for the different responses of each band pair. The satellites S2A and S2B are orbiting in the same plane and maintained in a phase shift of 180° [26,28]. This creates a special case where the IATGs are always significatively longer than one day. For this reason, their interoperability cannot be directly studied using this methodology. It could be inferred from our results, but it falls outside the scope of this work.

This methodology could be applied to other polar sun-synchronous satellites with similar characteristics, such as GaoFen-6/WFC, PROBA-V, Deimos-1 [45,46], or the future SeoSat (Ingenio) [47], even if they are in opposed orbit directions. With our methodology, it is possible to implement semi-automated quality assessment and quality control procedures because there is no need for BRDF calculations or atmospheric correction. Apart from the outlier inspection, the rest of the processes can be automated, including metadata retrieval, ground trajectory determination, SNO identification, data download, HAs creation, statistics extraction and representation.

5. Conclusions

In this work, we propose an empirical sensor comparison methodology for all spectral targets, based on SNOs, which does not need BRDF modelling or atmospheric correction to yield robust results. Due to its empirical approach, model-based techniques like the Spectral Band Adjustment Factors (SBAF) and atmospheric corrections are avoided. With illumination, ground and atmospheric effects minimized, the differences between intrinsic sensor characteristics, RSRs included, appear naturally, encompassed as slopes and uncertainties. The global scope of the method allows the acquisition of diverse ground types, with different spectral signatures and wide TOA reflectance ranges. Hence, enabling the creation of heterogeneous datasets for comprehensive cross-calibration analysis. The ground classes distribution will cause the slopes to have a systematic error that diminishes as diversity grows in the data set.

This procedure was used to assess the interoperability of L8 with S2A and L8 with S2B. In both cases, the band correlation calculations yielded high coefficients of variation. Showing that they are largely interoperable, with the exception of the wide NIR band, for which a research path was hinted. The methodology consistency was reassured by comparison with the results provided by third-party methods.

The presented technique can be extended. Since the L8 and S2 orbit planes are very close so their orbits run almost parallel for long distances [48] the SNO intersection point could be redefined as a segment, allowing larger SNO datasets. MSI tiles can be seamlessly mosaicked, without resampling, when they are in the same projection creating a bigger sample. Moreover, the lack of correlation of the residuals with the HAs distance to nadir and IATG means that both can be increased. Note that the method behaved well even with strongly anisotropic BRDFs (Figure 14). The number of valid SNOs could also increase raising the maximum cloudiness threshold. Another possibility to explore would be the usability of cloudy SNOs with extremely short IATGs. We have found that even minor anomalies are salient in the regressions, therefore, it would be possible to replace the visual outlier removal process with an automated one.

SNOs and HAs calculations are simple when compared with BRDF and atmospheric corrections. Therefore, this empirical method can be seamlessly extended for the assessment of other sensors, being the ground trajectory determination is the only element needing customization.