Identification of Global Extended Pseudo Invariant Calibration Sites (EPICS) and Their Validation Using Radiometric Calibration Network (RadCalNet)

Abstract

This study introduces a global land cover clustering using an unsupervised algorithm, incorporating the novel step of filtering data to retain only temporally stable pixels before applying K-means clustering. Unlike previous approaches that did not assess the pixel-level temporal stability, this method provides more reliable clustering results. The K-means identified 160 distinct clusters, with Cluster 13 Global Temporally Stable (Cluster 13-GTS) showing significant improvements in temporal stability. Compared to Cluster 13 Global (Cluster 13-G) from earlier research, Cluster 13-GTS reduced the coefficient of variation by up to 1% and increased the number of calibration locations from 23 to over 50. This study also validated these clusters using TOA reflectance from ground-truth measurements collected at the Radiometric Calibration Network (RadCalNet) Gobabeb (RCN-GONA) site, incorporating data from Landsat 8, Landsat 9, Sentinel-2A, and Sentinel-2B. The GONA Extended Pseudo Invariant Calibration Sites (EPICS) GONA-EPICS cluster used for the validation provided statistically comparable mean TOA reflectance to RCN-GONA, with a reduced chi-square test indicating minimal differences within the cluster’s uncertainty range. Notably, the difference in reflectance between RCN-GONA and GONA-EPICS was less than 0.023 units across all the bands. Although GONA-EPICS exhibited slightly higher uncertainty (6.4% to 10.3%) compared to RCN-GONA site (<5%), it offered advantages such as 80 potential calibration points per Landsat cycle and reduced temporal instability, and it provided alternatives to reduce the reliance on single sites like traditional PICS or RCN-GONA, making it a valuable tool for calibration efforts. These findings highlight the potential of the newly developed EPICS for radiometric calibration and stability monitoring of optical satellite sensors. Distributed across diverse regions, these global targets increase the number of calibration points available for any sensor in any orbital cycle, reducing the reliance on traditional PICS and offering more robust targets for radiometric calibration efforts.

1. Introduction

The availability of satellite data has played an essential role in advancing the under-standing and monitoring of changes on the Earth’s surface. However, a critical requirement for ensuring the reliability of these changes identified through satellite imagery is that the sensors used for image collection are radiometrically calibrated in compliance with international standards [1,2,3].

Even though the sensors on board satellites for Earth observation undergo rigorous testing before their launch, their radiometric calibration may degrade over time. Several techniques have been developed over the years to continuously assess the radiometric performance of sensors once they are in orbit [4,5]. One such technique involves the utilization of homogeneous, temporally stable regions on the Earth’s surface, referred to as pseudo invariant calibration sites (PICS) [6,7,8,9,10]. These calibration areas, predominantly situated in North Africa, have been the subject of numerous studies that have demonstrated their suitability as calibration targets. However, these locations also present certain limitations, including cloud cover, potential temporal changes that can lead to the misidentification of calibration drift specific to the site, and the necessity of several months or even years to accumulate a sufficiently extensive dataset for drift assessment and calibration purposes [11,12,13].

To overcome some of the limitations of traditional PICS, Shrestha et al. [14] conducted a pixel classification analysis across North Africa, specifically focusing on pixels characterized by a temporal uncertainty of less than 5%. This temporal uncertainty was calculated as the temporal standard deviation divided by the temporal mean. The classification process involved the use of a K-means clustering algorithm, leading to the identification of 19 distinct clusters. Notably, within this cluster set, Cluster 13 emerged as a strong candidate for potential selection as an EPICS. Cluster 13 was considered due to its substantial pixel count and the aggregation of pixels within contiguous regions across the continent providing near-daily observations depending on the sensor.

Hasan et al. [15] conducted a comprehensive investigation into the potential of EPICS as a target for radiometric calibration and stability monitoring for various satellite platforms, including Landsat 7 ETM+, Landsat 8 Operational Land Imager (OLI), Sentinel-2A, and Sentinel-2B. Their study estimated the temporal stability of Cluster 13 to be approximately 3%, with temporal uncertainties potentially up to 2% higher than those observed in PICS, depending on the spectral bands. Additionally, the spatial uncertainty associated with EPICS is approximately 5%, primarily due to their extensive coverage.

Nevertheless, the authors assert that the advantages of achieving near-daily calibration and increased sensitivity for change detection outweigh this spatial variability. That study demonstrated the capability of EPICS to enhance temporal resolution, resulting in multiple observations per sensor cycle, with an average frequency of approximately one observation every 1.4 days for Landsat 8 OLI. This represents a substantial enhancement compared to the traditional PICS approach, which typically provides observations every 18 to 20 days for the same sensor. Additionally, the authors revealed that a single randomly selected location within Cluster 13 can effectively represent the entire cluster with a 2% uncertainty, and this uncertainty decreases exponentially with the inclusion of more locations.

Fajardo et al. [16] expanded upon Shrestha and Hasan’s work by conducting a global-scale evaluation, improving the temporal resolution to an average of one observation every 1.02 days. However, during this global classification, regions displaying temporal instability became apparent. Consequently, the authors developed a methodology to evaluate each location within this classification to guarantee the exclusion of temporally unstable locations from the analysis, as well as locations with potentially different spectral characteristics. Despite the rigorous work that had to be performed to ensure that each location included in the global EPICS was temporally stable and spectrally similar, this global classification showed great potential to establish a global-scale target that could provide multiple observations per day of the Landsat 8 cycle using worldwide locations and removed the dependency on North African locations only. In addition, the global EPICS identified in that analysis demonstrated the capability of these global-scale targets to achieve a large dataset to perform stability analysis quicker than when using traditional PICS (months vs. years).

Given the challenges faced in the global classification process performed in the previously outlined analysis, particularly concerning the assurance of temporal stability for pixels allocated to distinct clusters and the integration of temporally unstable pixels, Fajardo et al. [17] developed a global mosaic of temporally stable pixels to support the radiometric calibration of optical satellite sensors. The authors focused on eliminating regions with temporal instability that could affect the accurate monitoring of sensor stability. The methodology involved two stages: first, the generation of data cubes containing time series data for specific regions of interest (ROIs) using Landsat 8 TOA reflectance data from 2013 to 2022, and second, the selection of a statistical test or combination of tests to perform a global analysis. Three ROIs were selected: a temporally stable region (Libya 4-CNES), a known unstable region in the Middle East, and a region in Brazil with a mix of stable and unstable pixels. Six statistical tests, including the Mann–Kendall test, Spearman’s rho test, and Pettitt’s test, were applied to detect long-term trends and change points. Cloud filtering and BRDF normalization were applied to minimize seasonal effects. Binary masks were created based on the results of each statistical test and their combinations. After comparing the masks to visually inspected samples, the test showing the highest agreement was selected as the reference mask.

In the second stage, global data cubes covering 1-degree latitude by 1-degree longitude were generated using Google Earth Engine (GEE), spanning latitudes −43 to 43 and longitudes −180 to 180. Two representative points per year (summer and winter) were selected for each pixel, with seasonal effects normalized using summer reflectance. After analyzing each pixel’s time series with combinations of statistical tests, the Spearman’s rho and Pettitt’s tests proved most effective in identifying temporally stable pixels while removing unstable ones. The global analysis produced binary masks, classifying pixels as stable or unstable. That study identified specific geographic areas suitable for radiometric calibration, although it did not classify pixels into land cover types or spectral clusters. The map generated from that work can be used for identifying new PICS or as a pre-filter for unsupervised clustering of stable pixels for calibration efforts, among others.

In this study, the global mosaic consisting of temporally stable pixels found in the previously mentioned study was used as the foundation for a global clustering process employing a K-means unsupervised clustering algorithm. It is essential to emphasize that this clustering was specifically limited to the regions previously identified as temporally stable. The global mosaic of temporally stable pixels obtained by the authors served as a data pre-filter for the K-means algorithm, thereby ensuring that only pixels demonstrating temporal stability were clustered into groups of pixels with similar temporal and spectral characteristics. This approach stands in contrast to the global classification conducted by Fajardo et al. [16], which lacked a temporal stability analysis and consequently resulted in the inclusion of unstable pixels across different clusters.

Moreover, this study validated the efficacy of the EPICS obtained through unsupervised classification for calibration tasks using the Radiometric Calibration Network (RadCalNet). This marks the first instance of EPICS validation using TOA reflectance data derived from ground-truth measurements. This validation highlights the potential use of these global EPICS for radiometric calibration and stability monitoring of optical satellite sensors. This approach eliminates the reliance on a singular optimal site and, instead, offers global coverage with multiple calibration opportunities worldwide. Detailed explanations regarding the classification of temporally stable pixels and the validation of global EPICS using RadCalNet data will be provided in the forthcoming sections.

2. Materials and Methods

This section will describe the classification methodology as well as the validation of an EPICS intersecting with the Gobabeb RadCalNet site. In addition, a detailed explanation of the validation steps taken will be provided.

2.1. Classification of Pixels on a Global Scale

2.1.1. Data Cube Generation and Clustering Process

In order to perform a global classification of the temporally stable pixels identified by Fajardo et al. [17] into different classes with similar spectral characteristics, TOA reflectance was used. This was used to cluster regions with not only similar cover type and spectral characteristics but also a similar response considering their atmosphere. To perform a global analysis, Landsat 8 OLI was chosen due to its reliable radiometric calibration better than 3% and its global data acquisition [5,18]. Due to the amount of data needed to perform this study, the world was split into a grid of data cubes of 1° latitude by 1° longitude from −43 to 43 latitude and 180 to 180 longitude. These data cubes were generated using GEE, a cloud-based platform developed by Google for the analysis of geospatial data and satellite imagery [19]. Each data cube contained a per-pixel temporal mean TOA reflectance spanning from 2013 to 2022 for each pixel for each spectral band in the visible and near-infrared (VNIR) and short wavelength infrared (SWIR) channels, as shown in Figure 1. For this study, a total of 9238 data cubes were generated.

Using the generated data cubes, a global mosaic was created and intersected with the global mosaic containing temporally table pixels only (filter mask), as identified by Fajardo et al. [17] in the work mentioned previously. This resulted in a global mosaic of 1° latitude by 1° longitude data cubes containing temporal mean TOA reflectance from 2013 to 2022 of temporally stable pixels only. An example of this is shown in Figure 2. Figure 2a shows the data cube containing the mean TOA reflectance created as described in Section 2.1 and shown in Figure 1. Figure 2b shows the data cube binary mask containing temporally stable pixels only as identified in Fajardo et al. [17] for a location in Central Africa. Lastly, Figure 2c shows the 1° by 1° temporally filtered data cube containing the mean TOA reflectance of temporally stable pixels only and input to the classification algorithm.

2.1.2. Classification Using a K-Means Clustering Algorithm

An unsupervised K-means clustering algorithm was selected to classify each pixel in the previously generated global mosaic. This choice was based on its computational simplicity, which is crucial when handling large datasets. Additionally, studies at both the continental and global scales have demonstrated the capability of K-means to cluster pixels with similar spectral characteristics globally. For instance, Cluster 13 achieved variability of less than 3% in the continental study and less than 5% in the global study, while providing multiple locations worldwide. This indicates that K-means can effectively group pixels from different regions, even across continents, that share the same spectral characteristics. These studies collectively highlight the potential of K-means to identify EPICS, yielding consistent results across various input datasets and initial conditions. This makes K-means a suitable choice for identifying EPICS, particularly given the large data volume and its adaptability for parallel processing.

The K-means algorithm used in this study is a standard implementation, configured to identify 160 clusters. It initializes randomly, uses the Euclidean distance as the metric, and incorporates all the layers of each data cube in the analysis. Details of this K-means clustering were described in the initial search of EPICS performed by Shrestha et al. [14].

The global mosaic of temporally stable pixels obtained by Fajardo et al. [17] served as a pre-classification filter. As described in Section 2.1.1, each data cube containing TOA reflectance was filtered using its corresponding temporal mask from this global mosaic before entering the K-means clustering algorithm, as shown in Figure 2. The main difference of this work compared to previous studies lies in the classification of pixels exhibiting temporal stability on a global scale, unlike the prior North African and global classifications, which lacked a per-pixel filter ensuring temporal stability.

2.2. Validation of EPICS

Even though the clusters identified using a K-means clustering algorithm have demonstrated great potential in previous studies, a comparison with TOA reflectance derived from ground-truth measurements has never been performed. This work presents, for the first time, a validation of one of the clusters identified in the classification performed in this work using the Radiometric Calibration Network (RadCalNet). In the upcoming subsections, a detailed explanation of the RadCalNet and Cluster comparisons will be provided.

2.2.1. Study Areas

RadCalNet Gobabeb Site

The RadCalNet, introduced to the public in 2018, comprises automated sites that collect essential measurements for teams overseeing the in-flight monitoring of spaceborne sensors. These sites ensure that radiometric calibrations are traceable to the International System of Units (SI), promoting data consistency across different sensors and agencies. The data collected at these automated sites are valuable for vicarious calibration and radiometric monitoring, undergoing quality control and processing to obtain TOA reflectance. These sites provide measurements at 30 min intervals for nadir-viewing sensors, with 10 nm spectral resolution at wavelengths between 400 and 2400 nm and 400 and 1000 nm depending on the site [20].

As of 2023, the RadCalNet comprises five automated sites situated in Railroad Valley Playa in the United States (RVUS), La Crau in France (LCFR), Gobabeb in Namibia (GONA), and artificial (BTCN) and sand targets (BSCN) in Baotou, China. For the purposes of this work, the Gobabeb site in Namibia (RCN-GONA) was chosen as the RadCalNet site for EPICS validation. The primary rationale for selecting the Gobabeb site was the presence of pixels from the RCN-GONA recommended ROI in the global clustering performed in this study. Sites such as RVUS and LCFR were excluded from the temporal stability filtering process; consequently, these RadCalNet sites were absent from the clustering process. Despite the presence of BSCN in the identified clusters, the Gobabeb site was preferred for validation, given its spectral range covering 400 to 2400 nm [21].

The RCN-GONA is situated in Gobabeb, Namibia, with the Gobabeb Namibia Research Institute located 10 km from the site and within the Namib-Naukluft National Park. The site is periodically monitored with an ASD FieldSpec spectroradiometer, collecting data across the spectral range from 400 nm to 2400 nm at 10 nm spectral intervals. Various atmospheric parameters, including the Angstrom coefficient, 550 nm aerosol optical depth (AOD), column water vapor, atmospheric pressure, and column ozone, among others, are collected at the site. The TOA reflectance product is representative of a 30 m radius centered in latitude 23.6002S and longitude 15.11956E [22,23]. Figure 3 shows an area of a Landsat 8 OLI image over WRS-2 Path 179/Row 76 where the red cursor and yellow circle show the center point of the RCN-GONA ROI.

GONA-EPICS

A cluster resulting from the K-means clustering algorithm described in Section 2.1.2, which contains pixels intersecting with the RCN-GONA site, was identified and selected as GONA-EPICS. This cluster was chosen because the pixels within the RCN-GONA site were included as part of it, suggesting that they exhibit similar spectral characteristics.

Figure 4 shows WRS-2 Path 176 Row 79, where the RCN-GONA site is located (left), along with a close-up of the ROI and surrounding pixels. The red cursor marks latitude 23.6002S and longitude 15.11956E, highlighting the RCN-GONA ROI and indicating the GONA-EPICS pixels within it. In both images, the magenta pixels represent the GONA-EPICS area. These figures illustrate that several regions within this WRS-2 Path/Row share spectral characteristics similar to those of the RCN-GONA site, leading to their inclusion in the same cluster. Although Figure 4 focuses on the EPICS pixels around the RCN-GONA site, this EPICS cluster is also present in various locations worldwide and will be discussed further in Section 3.2.

The mean TOA reflectance for the GONA-EPICS cluster was extracted following the methodology presented by Fajardo et al. [16]. This technique involves using zonal masks to retrieve the mean TOA reflectance for each location within the EPICS clusters, specifically for the GONA-EPICS, indicated by the magenta pixels shown in Figure 4. The masks were created and reprojected to UTM to align with the sensor’s projection system, thereby facilitating the retrieval of the TOA reflectance.

2.2.2. Sensors Used for Validation

Landsat 8 and Landsat 9

Landsat 8, launched on February 2013, has on board the OLI sensor and the Thermal Infrared Sensor (TIRS). The OLI sensor uses a pushbroom architecture with 14 individual focal plane modules that contain over 69,000 detectors measuring the visible, NIR, and short-wave infrared spectral regions. The OLI images the Earth in a Sun-synchronous orbit with a 185 km swath width and a field of view of 15 degrees. The OLI sensor images at a spatial resolution of 30 m for multispectral bands, 100 m for thermal bands and 15 m for the panchromatic band, and at a temporal resolution of 16 days [24,25,26].

Landsat 9, launched on 27 September 2021, is very similar to Landsat 8, with both satellites having nearly identical imaging systems to ensure continuity in data acquisition. Landsat 9 hosts the Operational Land Imager-2 (OLI-2) and the Thermal Infrared Sensor-2 (TIRS-2). The spectral bands and spatial resolution of OLI-2 and TIRS-2 are almost identical to those on Landsat 8, maintaining the same spectral ranges for the multispectral, panchromatic and thermal bands (30 m, 15 m and 100 m, respectively). A key difference is that OLI-2 on Landsat 9 provides a 14-bit data download for each pixel, compared to 12-bit for OLI on Landsat 8. These sensors have an absolute calibration within 3% [27,28].

Sentinel-2A and Sentinel-2B

The Sentinel-2 mission, part of the European Commission’s Copernicus program, includes the Sentinel-2A and Sentinel-2B satellites, both equipped with Multi-Spectral Instruments (MSI). Launched in June 2015 and early 2017, respectively, these satellites operate in a Sun-synchronous orbit at an altitude of 786 km, phased 180° relative to each other. They offer a revisit time of 10 days for each platform and 5 days when combined. The MSI sensor captures medium- to high-resolution imagery across 13 spectral bands, with spatial resolutions ranging from 10 to 60 m. Distributed across 12 separate modules, the MSI focal plane detectors enable the satellites to capture images with a 20.6-degree field of view and a swath of 290 km. These sensors are calibrated to within 5% accuracy [3,29,30].

RadCalNet

Reflectance at the RCN-GONA site is measured within a 30 m radius disk using a Cimel 12-filter BRDF Sun photometer and an ASD FieldSpec spectroradiometer. These instruments provide surface reflectance measurements with uncertainties ranging from 3% to 4.3% across wavelengths from 400 nm to 2400 nm in 10 nm steps. The site also monitors atmospheric conditions, such as aerosols, water vapor, and ozone content, using a Sun photometer and NASA data to produce TOA reflectance products. Additionally, an AllSky-340 sky camera installed 10 km away at the Gobabeb Namib Research Institute captures cloud cover images every 10 min. RadCalNet’s processing methodology for converting site-specific data into TOA reflectance includes using the MODTRAN 6.0.2r2 radiative transfer code, which uses inputs such as the latitude, longitude, altitude, and atmospheric conditions specific to each data collection time and site [31].

2.2.3. Data Processing for Validation

TOA Reflectance Retrieval

• Landsat 8 and Landsat 9

For the L8 and L9 sensors, the TOA reflectance retrieval was achieved by converting digital numbers (DNs) to TOA reflectance, as calculated using Equation (1) [32].

$$\rho \lambda ={\displaystyle \frac{M\rho \times Qcal+Ap}{cos\left(\theta SZ\right)}}$$

(1)

where Mρ and Aρ are the multiplicative and additive scaling factors, Qcal is the quantized and calibrated product for the pixel values (digital numbers or DNs), and θSZ is the per-pixel solar zenith angle.

• Sentinel-2A and Sentinel-2B

Due to variations in processing, the Sentinel-2A and Sentinel-2B sensors require two distinct methodologies for computing TOA reflectance. For the Multispectral Instrument (MSI) data collected prior to 25 January 2022, the TOA reflectance was calculated using Equation (2) [29]

$$\rho \lambda ={\displaystyle \frac{D{N}_{cal}}{Q}}$$

(2)

where $\rho \lambda$ represents the TOA reflectance for Sentinel-2 Level 1, $D{N}_{cal}$ is the calibrated and quantized product pixel value, and Q is the quantization value or reflectance scaling factor obtained from the metadata file. In addition, considering updates to the Sentinel-2 processing system, the TOA reflectance for data collected after January 25 was calculated using Equation (3).

$$\rho \lambda ={\displaystyle \frac{D{N}_{cal}+Offset}{Q}}$$

(3)

where $\rho \lambda and D{N}_{cal}$ are the Sentinel-2 TOA reflectance and the calibrated and quantized product pixel value, respectively, $Q$ is the quantization value or reflectance scaling factor, and $Offset$ is the radiometric offset obtained from the metadata file.

• RadCalNet

As mentioned in the RadCalNet Gobabeb Site section, the RadCalNet collects data every 30 min for nadir-viewing sensors and provides a TOA reflectance product from these collections. For this analysis, the RadCalNet collection time was closest to the sensor’s overpass time. Considering the hyperspectral nature of RadCalNet data, a comparison with the multispectral spaceborne sensors used in this study required the convolution of the RCN-GONA TOA reflectance spectra with the Spectral Response Function (RSR) of each space sensor. Figure 5 illustrates the hyperspectral RCN-GONA TOA reflectance data alongside the RSR for Landsat 8, which was used to obtain the multispectral RCN-GONA TOA reflectance for validation purposes. This was also achieved for all the other sensors used in this study.

Cloud Filtering

To eliminate contaminated pixels, a cloud filter was implemented for every sensor used. The first cloud filter was implemented in the data downloading process, where every scene marked by USGS as having more than 10% cloud cover was not downloaded into the Image Processing Laboratory (IP LAB) archive. In addition, for every scene downloaded and for any ROI, images were discarded if 50% or more of the pixels within the ROI were cloud-covered. Conversely, images with less than 50% cloud coverage were retained for further analysis, and a per-pixel cloud filter was applied. For the EPICS processing, this 50% threshold was applied on a Path/Row basis. Thus, if pixels in the GONA-EPICS of a specific Path/Row were 50% or more covered by clouds, the image was rejected to prevent cloud contamination; otherwise, it was accepted and subjected to per-pixel filtering.

For Landsat 8 and Landsat 9, the pixel quality assessment band from Level 1 Collection 2 was employed to create per-scene cloud filtering masks. This band contains specific Bits that are essential for filtering various atmospheric and surface elements. For cloud filtering in particular, the relevant Bits used are as follows: Bit 0 for fill values, Bit 1 for Dilated Clouds, Bit 2 for Cirrus, Bit 3 for Clouds, Bit 4 for Cloud Shadows, Bit 9 for Cloud Confidence, Bit 11 for Cloud Shadow Confidence, and Bit 15 for Cirrus Confidence. More information about the cloud filtering using the quality control band can be found at https://www.usgs.gov/media/files/landsat-8-9-olitirs-collection-2-level-1-data-format-control-book accessed on 10 February 2022. Figure 6 illustrates a Landsat 8 image and its generated cloud binary mask obtained for the WRS-2 Path 179/Row 76, where white pixels are clouds and black regions are cloud-free pixels. In addition, for Sentinel-2, the cloud mask product was applied for each resolution to remove contaminated cloud pixels following the threshold methodology.

BRDF Normalization

The mean TOA reflectance recorded by a satellite sensor can be affected by the atmospheric conditions, view geometry of acquisition, Sun’s position at the time of the collect and therefore seasonal effects. The impact of changes in the Sun’s position on the TOA reflectance of the target can be described by the Bidirectional Reflectance Distribution Function (BRDF). By modeling the BRDF, the seasonal impact on TOA reflectance measurements can be minimized. Kaewmanee developed a 4-angle BRDF quadratic model, refined by Hasan et al. by including all the interaction terms between solar and view angles in the quadratic model (IPLAB BRDF model), as shown in Equation (4) [6,33].

$$\begin{gathered} {\rho }_{predicted}={\beta }_0+{\beta }_1{X}_1+{\beta }_2{Y}_2+ {\beta }_3{X}_2+ {\beta }_4{Y}_2+ {\beta }_5{X}_1{Y}_1+{\beta }_6{X}_1{X}_2+{\beta }_7{X}_1{Y}_2+{\beta }_8{Y}_1{X}_2+{\beta }_9{Y}_1{Y}_2 \\ +{\beta }_{10}{X}_2{Y}_2+{\beta }_{11}{X_1}^2+{\beta }_{12}{Y_1}^2+{\beta }_{13}{X_2}^2+{\beta }_{14}{Y_2}^2 \end{gathered}$$

(4)

where $X_1$, $Y_1$, $X_2$, $Y_2$ are the Cartesian coordinates projected from the spherical coordinates. ${\beta }_0$ through ${\beta }_{14}$ represent the model coefficients, and ${\rho }_{BRDF predicted}$ denotes the predicted TOA reflectance. The Cartesian coordinates $X_1$, $Y_1$, $X_2$, $Y_2$ used for this model are as follows:

$$X_1=\sin \left(SZA\right)\ast \cos \left(SAA\right)$$

(5)

$$Y_1=\sin \left(SZA\right)\ast \sin \left(SAA\right)$$

(6)

$$X_2=\sin \left(VZA\right)\ast \cos \left(VAA\right)$$

(7)

$$Y_1=\sin \left(VZA\right)\ast \sin \left(VAA\right)$$

(8)

where SZA, SAA are the solar zenith and azimuth angles and VZA, VAA correspond to the view zenith and azimuth angles.

The BRDF normalized TOA reflectance was calculated using the following equation:

$${\rho }_{BRDF Normalized}={\displaystyle \frac{{\rho }_{observed}}{{\rho }_{predicted}}}\ast {\rho }_{reference}$$

(9)

In Equation (9), ${\rho }_{BRDF Normalized}$ is the BRDF normalized TOA reflectance, ${\rho }_{observed}$ is the TOA reflectance observed by each satellite, ${\rho }_{predicted}$ is the per-observation BRDF-predicted TOA reflectance obtained from the BRDF model, and ${\rho }_{reference}$ is the reference TOA reflectance estimated based on a reference geometry of acquisition. For this study, BRDF models were derived for both the GONA-EPICS cluster data and the RCN-GONA data. These BRDF models were also obtained for each sensor considered in this study (Landsat 8, Landsat 9, Sentinel-2A, and Sentinel-2B). Reference solar angles were selected as shown in Figure 7. The green data represent the solar angles for GONA-EPICS, while the magenta data represent the solar angles for the RCN-GONA site. Solar reference angles were chosen where both datasets intersect, ensuring that both RCN-GONA and GONA-EPICS were normalized to a similar geometry of acquisition. The view angles for RCN-GONA were set to zero, as the TOA reflectance for this site is given for the nadir only. Additionally, the view reference angles for GONA-EPICS were chosen as the view angles for the selected reference solar angles. The reference angles chosen for this study were SAA: 92, SZA: 25, VAA: 0.3, VZA: 145 for GONA-EPICS, and SAA: 92, SZA: 25, VAA: 0, VZA: 0 for the RCN-GONA site.

2.2.4. Validation Methodology

The validation of the EPICS was performed using the GONA-EPICS and RCN-GONA site data. Given that the RCN-GONA site provides TOA reflectance data derived from ground measurements, this validation considers the RCN-GONA TOA reflectance as the true measurement. Two statistical tests were considered: the reduced chi-square test and Welch’s t-test. The reduced chi-square test evaluates the agreement between the GONA-EPICS data and the model RCN-GONA, using the GONA-EPICS uncertainty to determine whether the observed values fall within the expected range of the true measurement. Welch’s t-test compares the means of RCN-GONA and GONA-EPICS while accounting for the uncertainties in both datasets. By applying these complementary statistical approaches, the analysis aims to confirm the reliability of the EPICS validation.

Statistical Analysis for Validation

• Reduced Chi Square Test ${{\chi }^2}_{red}$

To validate the GONA-EPICS, the reduced chi-square test (${{\chi }^2}_{red}$) was used. The value of ${{\chi }^2}_{red}$ indicates the disagreement between the observed values and the values expected under a statistical model, accounting for uncertainties. In this study, the statistical model was based on TOA reflectance measurements from the RCN-GONA site, while the observed data were the GONA-EPICS TOA reflectance. The reduced chi-square test is based on the chi-square test described in Equation (10) [34].

$${\chi }^2={\displaystyle \sum }_{i=1}^n{\left({\displaystyle \frac{y_i-f\left(x_i\right)}{{\sigma }_i}}\right)}^2$$

(10)

where f(x) is a fit function applied to n data points, $y_i$ is the measurement of the quantity $y$, usually obtained experimentally when the quantity $x$ is $x_i$, and ${\sigma }_i$ is the uncertainty of $y_i$.

Finally, the ${{\chi }^2}_{red}$ is the ratio of ${\chi }^2$ to the degrees of freedom, calculated as $df=n-p$, where n is the number of data points and p is the number of parameters in the model.

The chi-square statistic is typically used to evaluate the goodness of fit between a model and observed data. In this study, however, a non-traditional or “inverse” approach was applied. The model, based on RCN-GONA data, is assumed to represent the “true” behavior of the RCN-GONA site, while the GONA-EPICS TOA reflectance data are evaluated against this model to assess the alignment with the expected values. In this inverse approach, the model is treated as deterministic—meaning it does not possess inherent uncertainty—and the observed data (GONA-EPICS TOA reflectance) must include some level of uncertainty. The uncertainty of this dataset was estimated in the GONA Cluster Uncertainty section.

Bevington et al. [35] provide details on the interpretation of the reduced chi-square test (${{\chi }^2}_{red}$). A “good fit”, where the model represents the observed data within the uncertainties, is indicated by a ${{\chi }^2}_{red}$ value close to 1. Conversely, values greater than one suggest a “poor fit”, indicating that the model does not adequately represent the data, even when uncertainties are considered. Additionally, ${{\chi }^2}_{red}$ values lower than one may suggest that the uncertainties were overestimated. However, a low ${{\chi }^2}_{red }$ value could also imply that the model does not properly describe the observed data, despite accounting for uncertainties. For this analysis, histograms of both datasets were examined to determine if their distributions overlapped. This approach helped assess whether a low ${{\chi }^2}_{red}$ value signifies a potential overestimation of the uncertainties rather than a statistical difference between the datasets.

• Welch’s t-Test

The reduced chi-square test, as mentioned earlier, provides insight into the alignment between the observed GONA-EPICS data and the model RCN-GONA, using the GONA-EPICS uncertainty to assess whether the GONA-EPICS values fall within the expected range of RCN-GONA. Welch’s t-test was applied to compare the means of RCN-GONA and GONA-EPICS while accounting for the uncertainties in both datasets. Unlike the reduced chi-square test, which evaluates the consistency of the GONA-EPICS values with the model RCN-GONA based on the GONA-EPICS uncertainty, Welch’s t-test determines the statistical difference between the two means, incorporating uncertainties from both datasets. Using both tests demonstrates consistency across different statistical tests. In this analysis, the RCN-GONA uncertainty used was that provided by the RadCalNet with the TOA reflectance product used in this study.

Welch’s t-test was applied considering the GONA-EPICS and the RCN-GONA uncertainty. The hypotheses of this test are that [36]:

Null hypothesis:

$$H_0={\mu }_{RCN-GONA}={\mu }_{GONA-EPICS}$$

(11)

Alternative hypothesis:

$$H_1={\mu }_{RCN-GONA}\ne {\mu }_{GONA-EPICS}$$

(12)

Welch’s t-test compares the difference between the sample means to determine if it is statistically significant, accounting for the unequal variances between the samples. Unlike the traditional application of Welch’s t-test, which uses standard errors derived from the sample variances, this analysis incorporates combined uncertainties to account for both the GONA-EPICS and RCN-GONA uncertainties as follows:

$$t={\displaystyle \frac{{\mu }_{RCN-GONA}-{\mu }_{GONA-EPICS}}{\sqrt{{U^2}_{RCN-GONA}+{U^2}_{GONA-EPICS}}}}$$

(13)

The null hypothesis posits that the two population means are the same. If the p-value is higher than the significance level (0.05), it indicates insufficient evidence to reject the null hypothesis, suggesting that any observed differences may be due to random variation rather than a true effect.

• Modified Savitzky–Golay (MSG) Filter

As mentioned above, the MSG filter was used to temporally interpolate the RCN-GONA TOA reflectance for the GONA-EPICS validation using the reduced chi-square test. The Savitzky–Golay filter, initially introduced by Savitzky and Golay, is a method used for smoothing data in the time domain. It applies a low-pass filter constructed from local least-squares polynomial approximations. This methodology performs temporal interpolation while preserving the potential patterns in the original dataset [37,38]. The MSG filter fits a polynomial to the dataset using a specified window size. This produces an output representing the polynomial value at the central point of the window. The window then shifts by one day, and the fit is performed again, continuing this process throughout the time series. For this analysis, a window size of 120 days and a polynomial fit of order three was used, considering this window size proved to be effective in other studies where the MSG was used for temporal interpolation purposes using remote sensing data [39].

• GONA Cluster Uncertainty

To apply the reduced chi-square test ${{\chi }^2}_{red}$, the uncertainty of the observed data was required. Consequently, the GONA-EPICS uncertainty was calculated using the methodology outlined in the Guide to the Expression of Uncertainty in Measurement (ISO-GUM) [40], as detailed in Equation (11). For this analysis, $U_{target}^2$ represents the temporal coefficient of variation (CV%), which is the ratio of the temporal standard deviation of the mean TOA reflectance to the temporal mean TOA reflectance including all the Paths and Rows within the GONA-EPICS cluster. This CV captures both temporal and spatial variability within the cluster. Additionally, $U_{BRDF}^2$ refers to the BRDF error, calculated as the root mean square error (RMSE) between the observed TOA reflectance and the TOA reflectance predicted by the BRDF model. Lastly, $U_{sensor}^2$ denotes the sensor uncertainty, with absolute radiometric calibration uncertainties of 3% for Landsat 8 and 9, and 5% for Sentinel-2A and 2B, specified as the radiometric calibration requirement for each mission.

$${U_{total}}^2=\sqrt{U_{target}^2+U_{BRDF}^2+U_{sensor}^2}$$

(14)

3. Results and Analysis

3.1. Classification Results

This study developed a global land cover classification to identify targets for radiometric calibration and stability monitoring of optical satellite sensors using a K-means clustering algorithm. In contrast to previous studies, this research focused exclusively on pixels that were screened and identified as temporally stable on a global scale. Earlier EPICS searches, such as those conducted by Shrestha et al. [14] and Fajardo et al. [16], did not incorporate a per-pixel stability analysis to exclude pixels with known temporal instability.

The K-means clustering algorithm developed in this study identified 160 clusters, each with distinct temporal and spectral characteristics. Figure 8 displays the multispectral signatures of the identified clusters, which include the spectral signatures of targets such as soil (including sands), rocks, vegetation, and water across all the bands. However, some clusters exhibit high reflectance in the shorter wavelengths and lower reflectance in the longer wavelengths, which is typical of targets covered by snow or ice. These clusters could also represent potential remaining clouds that were not detected by the Landsat 8 quality control band, as their spectral response similarly exhibits high reflectance in the shorter wavelengths and low reflectance in the longer wavelengths. Such targets pose challenges for calibration and stability monitoring of optical satellite sensors due to their inherent properties and higher spatial variability, as reflected by the error bars representing the temporal standard deviation of each cluster.

Given the challenges and higher variability associated with these targets, a threshold based on the standard deviation of previously used EPICS was applied, as shown in Figure 8 (right). This threshold, three times the standard deviation of previously identified EPICS, was used to identify potential new EPICS with similar temporal characteristics to the existing EPICS [15,17,41]. This threshold identified 75 clusters that represent potential EPICS candidates distributed globally. While this filtering process reduced the diversity of the clusters from 160 to 75, it removed targets unsuitable for stability monitoring and calibration efforts. The remaining clusters, however, hold potential for these tasks. It is important to note that the dynamic range of the identified targets represents an improvement over those identified in Shrestha’s work in North Africa, as it includes a broader range of targets beyond the narrow spectral range of sand types identified in that study.

To illustrate the pixel distribution of some EPICS identified in this study, Figure 9 displays 10 EPICS randomly selected from the 75 clusters identified and shown in the right chart. The figure shows that many of these clusters are widely distributed across the globe, with several of them present in multiple continents. Notably, a substantial number of clusters are concentrated in North Africa, suggesting that many of these EPICS may represent different types of sand or rocks. Additionally, the cluster depicted in cyan color shows a different distribution, primarily in South America, Central Africa, and East Asia. Based on this geographical distribution and visual inspection, this cluster appears to represent dense vegetation, such as that found in the Amazon rainforest.

This map highlights the K-means clustering algorithm’s ability to group pixels with similar spectral properties on a global scale. Furthermore, the map illustrates the EPICSs’ potential to provide numerous calibration points with comparable spectral characteristics despite being located in different regions of the planet. These EPICS provide major advantages, allowing sensors on any orbital path to image multiple portions of particular EPICS throughout each sensor cycle. This increases the number of observations for the calibration and stability monitoring of optical sensors compared to those achieved with traditional PICS, which are typically constrained to small locations in North Africa and may be limited by cloud cover.

Cluster 13-GTS (Global Temporally Stable) vs. Cluster 13-G (Original Global Class)

Fajardo et al. [16] conducted a global land cover classification to identify temporally stable pixels worldwide using a K-means clustering algorithm. However, their study did not filter out pixels and classified every pixel within the latitude range of −45 to 45 and longitude range of −180 to 180. The authors reported that several locations contained temporally unstable pixels, necessitating the development of a methodology to exclude such locations and retain only temporally stable pixels for the EPICS evaluation. This approach proved time-consuming and limited the number of locations available for radiometric calibration and stability monitoring.

In contrast, this work utilized the global mosaic of temporally stable pixels identified by Fajardo et al. [17] as a data filter prior to the K-means clustering application. This approach ensured that only temporally stable pixels were included in the analysis, eliminating the need for additional evaluations. Consequently, the pixels within the clusters output of the K-means algorithm are temporally stable and can be used directly without further processing.

Fajardo et al. [16] evaluated one of the EPICS identified in that clustering analysis for radiometric calibration and stability monitoring tasks. To achieve this, the EPICS known as “Cluster 13” was used. In that study, 160 clusters were identified, and “Cluster 13” was chosen as an EPICS due to its large pixel count, the aggregation of pixels, and its inclusion in the Libya CEOS ROI. That study demonstrated that EPICSs can be used for radiometric calibration and ongoing stability monitoring of the sensor’s performance, in addition to outperforming traditional PICS, such as Libya 4, in cases where the sensor’s performance, particularly temporal stability (measured as drift per year), needs to be assessed over shorter time frames. The EPICS proved to be better in these cases due to its capability to provide multiple observations per day for some days of the sensor’s cycle, which in turn created richer datasets that allowed a better estimate of the slope for the stability evaluation.

This study compares the cluster equivalent to “Cluster 13” identified in this work with “Cluster 13” from Fajardo et al. [16] to evaluate the potential impact of the pre-filter data on the resulting clusters. The “Cluster 13” from Fajardo et al. [16] will be referred to as “Cluster 13-G”, while the cluster identified in this work will be named “Cluster 13-GTS,” as it includes only temporally stable pixels. Figure 10 illustrates the pixel distribution of Cluster 13-G and Cluster 13-GTS, shown in red and green colors, respectively, as depicted in Figure 10.

To compare the datasets Cluster 13-G from Fajardo et al. [16] and Cluster 13-GTS generated in this work, 21 WRS-2 Path/Row locations identified as optimal in Fajardo et al. [16] were used. However, WRS-2 Path 30/Row 38 in North America and WRS-2 Path 100/Row 82 in Australia were not found in Cluster 13-GTS and were excluded from the comparison. The temporal mean TOA reflectance for 19 locations for both datasets was retrieved and the BRDF normalized using the BRDF model described in the BRDF Normalization section. The reference geometry of the acquisition for both datasets was VZA: 3, VAA: 105, SZA: 30, and SAA: 130 degrees.

The temporal mean TOA reflectance, temporal standard deviation, and temporal coefficient of variation (CV) were used to evaluate the variability of these EPICS, as shown in Figure 11. While the two datasets are statistically similar across all the spectral bands, their temporal variability differs. The CV of Cluster 13-G ranges from 2.8% in the NIR band to 4.5% in the blue band, while Cluster 13-GTS shows a decreased CV across all the spectral bands, with the highest at 3.7% and the lowest at 1.8%. This indicates that the temporal filter improved the clustering results by removing unstable pixels, resulting in clusters that are more compact with lower variability and reduced uncertainty, as shown in Equation (11). It is important to mention that, although this comparison included only 19 locations based on the selected sites for Cluster 13-G, Cluster 13-GTS can provide 58 observations per Landsat 8 cycle collected over North Africa and the Middle East regions, with a CV ranging from 2.2% to 3.7%. This indicates the compactness of Cluster 13-GTS and demonstrates that, even with an increased number of locations, the CV remained stable.

3.2. Validation of EPICS Results

3.2.1. GONA-EPICS

In this work, as previously mentioned, 160 different clusters were identified. In addition to the previously discussed “Cluster 13”, another cluster of interest is “Cluster 160”, referred to as the GONA-EPICS cluster. This cluster was selected considering that it includes pixels from the RCN-GONA site. Consequently, the GONA-EPICS was chosen to validate these EPICS using the RCN-ROI TOA reflectance derived from ground measurements, as mentioned earlier. Figure 12 displays the pixel distribution of this GONA-EPICS worldwide. It is evident that this EPICS includes pixels in most continents, with a significant presence in Africa, the Middle East, and Asia, as well as several locations in Australia, North America, and South America.

Figure 12 also shows Landsat 8 images of all the WRS-2 Path/Row combinations intersecting this cluster that are archived in the IP LAB and used for this analysis.

Table 1. WRS-2 Path/Row used in this study to retrieve the TOA reflectance over the GONA-EPICS for each day of the Landsat 8 acquisition cycle.

Day of the Cycle	Path	Rows	Day of the Cycle	Path	Rows
1	37	38, 41	10	173	47
	190	43, 44		189	46, 48, 48
2	181	73		205	43
3	172	39	11	43	33
	188	45, 47		100	80, 82
	204	43, 46, 47		132	32
4	99	79		164	42
	163	37, 43		180	75
	179	42, 45, 76		196	41, 43
5	33	38	12	171	41
	170	39, 42		187	43
	202	44, 46		203	43, 46
6	40	33	13	162	48
	177	39		178	39, 41, 79
	193	37		194	41
7	152	43	14	185	43
	168	51		201	43, 47
	184	43		233	77
	200	47, 48
8	38	37, 38	15	39	37
	127	32		128	31, 32
	159	40, 46		160	47
	175	46		176	40
	191	37, 42		192	37, 42, 44
9	166	41	16	30	38
	182	46		167	40, 51
	198	42		183	46, 48

lists these WRS-2 Path/Row combinations, showing the Landsat 8 acquisition cycle days and the corresponding WRS-2 Path/Row imaged on each day. The table indicates that, for each day of the Landsat 8 acquisition cycle, there is at least one observation, with some days having up to three collections. The ability to obtain multiple collections per day facilitates the rapid and accurate assessment of sensor temporal stability, allowing for timely corrections if any changes are detected.

Figure 13 displays the BRDF-normalized mean TOA reflectance of the GONA-EPICS for all the Landsat 8 VNIR and SWIR bands. Outlier rejection for this dataset was performed using a 3-sigma approach, where scenes exceeding this threshold were identified as potential outliers and removed if a visual inspection confirmed the presence of undetected clouds that were missed by the cloud binary mask described in the Cloud Filtering section.

The GONA-EPICS exhibits a CV% of less than 6% for all the bands except the blue and SWIR2 bands, which have CVs of 6.3% and 6.1%, respectively, as shown in

Table 2. GONA-EPICS temporal mean, temporal standard deviation and temporal coefficient of variation.

	CA	Blue	Green	Red	NIR	SWIR1	SWIR2
Mean TOA reflectance	0.20	0.20	0.23	0.30	0.36	0.42	0.35
Temporal std.	0.01	0.01	0.01	0.01	0.01	0.02	0.02
CV (%)	5.7	6.3	5.3	4.8	3.9	4.4	6.1

. Notably, this EPICS offers 80 potential calibration points for every 16 days of the Landsat 8 cycle, in contrast to the single calibration point provided by traditional approaches such as PICS or the RCN-GONA site. A total of 80 locations were used based on their availability in the IP LAB archive; however, as shown on the map, many additional locations with EPICS pixels could potentially be added and utilized.

3.2.2. RCN TOA Reflectance

To validate the GONA-EPICS presented above, the RCN-GONA TOA reflectance derived from the surface reflectance measurements was used, as previously discussed. Figure 14 and

Table 3. RCN-GONA mean TOA reflectance, temporal standard deviation and temporal coefficient of variation.

	CA	Blue	Green	Red	NIR	SWIR1	SWIR2
Mean TOA reflectance	0.19	0.19	0.23	0.31	0.35	0.42	0.35
Temporal std.	0.009	0.009	0.012	0.015	0.016	0.017	0.017
CV(%)	4.8	4.7	5.1	5.0	4.4	4.1	4.7

show the RCN-GONA TOA reflectance for a 30 m disk centered at latitude 23.6002°S and longitude 15.11956°E, with the shaded area representing the reported RCN-GONA uncertainty. The figure indicates that the RCN-GONA site exhibited substantial variability and instability from 2017 to 2020, with stabilization occurring after mid-2020. Consequently, only data from after 2020 were used for this validation.

3.2.3. Statistical Analysis Results

GONA-EPICS vs. RadCalNet Uncertainty

As mentioned above, the GONA-EPICS mean TOA reflectance was validated using the RCN-GONA site data collected between May 2020 and May 2023. Figure 15 shows the mean TOA reflectance for the GONA-EPICS and the RCN-GONA site during this period. The figure demonstrates that, for all the bands, the mean TOA reflectance of the RCN-GONA site aligns closely with the GONA-EPICS TOA reflectance. However, the red band shows a slightly higher TOA reflectance compared to the GONA-EPICS.

The GONA-EPICS was validated using the reduced chi-square test (${{\chi }^2}_{red}$), which, as mentioned earlier, requires the uncertainty of the observed data. In this study, the uncertainty associated with the GONA-EPICS site was assessed. The total uncertainty was estimated as detailed in the GONA Cluster uncertainty section, with the resulting values for all four sensors evaluated presented in

Table 4. Sources of uncertainty and total percent uncertainty for all the sensors considered in this study for all the spectral bands.

	Sensor	CA	Blue	Green	Red	NIR	SWIR1	SWIR2
Temporal std.	L8	5.4	5.9	4.9	4.2	3.5	3.9	5.5
	L9	5.5	5.9	5.0	4.3	3.6	4.1	5.6
	S2A	6.2	6.4	5.4	4.7	3.7	4.3	5.3
	S2B	5.9	5.9	4.5	4.2	2.9	2.6	4.4
BRDF	L8	5.3	5.8	4.8	4.3	3.4	4.0	5.6
	L9	5.3	5.8	4.9	4.4	3.5	4.2	5.8
	S2A	6.1	6.3	5.4	4.9	3.7	4.5	5.4
	S2B	5.6	5.7	4.5	4.4	2.9	2.6	4.6
Sensor’s unc.	L8	3.0	3.0	3.0	3.0	3.0	3.0	3.0
	L9	3.0	3.0	3.0	3.0	3.0	3.0	3.0
	S2A	5.0	5.0	5.0	5.0	5.0	5.0	5.0
	S2B	5.0	5.0	5.0	5.0	5.0	5.0	5.0
Total uncertainty	L8	8.1	8.8	7.5	6.7	5.7	6.4	8.4
	L9	8.2	8.8	7.7	6.9	5.9	6.6	8.6
	S2A	10.0	10.3	9.1	8.4	7.3	8.0	9.1
	S2B	9.5	9.6	8.0	7.8	6.5	6.2	8.1

. This table displays each source of uncertainty considered, including the temporal standard deviation, BRDF error, and sensor-specific uncertainty.

The Landsat 8 and Landsat 9 sensors generally exhibit lower temporal standard deviation values across all the bands, ranging from 3.5% to 5.9%. In contrast, the Sentinel-2A and Sentinel-2B sensors show slightly higher values, particularly in the CA and blue bands, with Sentinel-2A reaching up to 6.4%. Sentinel-2B demonstrates better performance in the SWIR bands, with the lowest temporal standard deviation of 2.6% in the SWIR1 band.

For the BRDF errors, the results from the Landsat 8 data collected over the GONA-EPICS dataset show values ranging from 3.4% in the NIR band to 5.8% in the blue band, indicating consistent performance across bands with slightly higher errors at shorter wavelengths. Similarly, the Landsat 9 data over GONA-EPICS exhibit BRDF errors ranging from 3.5% to 5.8%, with the lowest error in the NIR band and the highest in the blue band. Sentinel-2A data present a broader range of BRDF errors, from 3.7% in the NIR band to 6.3% in the blue band, showing comparable results to the Landsat sensors. Conversely, the Sentinel-2B results over the GONA-EPICS show the lowest overall BRDF, ranging from 2.6% in the SWIR1 band to 5.7% in the blue band, indicating improved performance in terms of BRDF normalization.

The total uncertainty includes not only the temporal variability of the site and the BRDF normalization errors but also the sensor’s total uncertainty. Landsat 8 shows a total uncertainty ranging from 5.7% to 8.8%, with the lowest reflectance at 5.7% and the highest at 8.8%. Landsat 9 exhibits a similar range, from 5.9% to 8.8%, with slightly higher total uncertainty in most bands compared to Landsat 8. Sentinel-2A records the highest total uncertainty among the sensors, with a range from 7.3% to 10.3%, indicating consistently higher uncertainty. Sentinel-2B, while similar to Sentinel-2A, shows a slightly lower range, from 6.2% to 9.6%. It is important to remember that these uncertainties include the absolute radiometric calibration requirement for each of the sensors, 3% for Landsat 8 and Landsat 9 and 5% for Sentinel-2A and Sentinel-2B. However, the actual absolute radiometric calibration might be lower for these sensors; therefore, the GONA-EPICS can also be lower. The total percent uncertainty presented in Table 4 was used to validate the GONA-EPICS with the RCN-GONA data through the reduced chi-square test. The results of this validation will be discussed in the following section.

Reduced Chi-Square Results

In this work, the validation of GONA-EPICS using the RCN-GONA site involved four different sensors. The mean TOA reflectance data from these sensors collected over the GONA-EPICS were compared to the convolved RCN-GONA mean TOA reflectance for each sensor. The bands evaluated included CA, blue, green, red, NIR, SWIR1, and SWIR2 for all the sensors.

Figure 16 illustrates the histograms of the mean TOA reflectance for Landsat 8 OLI over GONA-EPICS (depicted in green) and the mean TOA reflectance for the RCN-GONA site (depicted in magenta) for all the bands. Additionally, this figure shows the results of the reduced chi-square test applied to both datasets. The reduced chi-square test, as mentioned earlier, evaluates how well the GONA-EPICS TOA reflectance aligns with the expected values derived from the RCN-GONA dataset. All the reduced chi-square values (${{\chi }^2}_{red}$) are below 1, except for the red band.

The interpretation of the ${{\chi }^2}_{red}$ results, as detailed in the Reduced chi square test ${{\chi }^2}_{red}$ section, demonstrates the level of agreement between the observed data and the expected values obtained from the RCN-GONA dataset. This is also evident from the histograms, where the TOA reflectance of GONA-EPICS and the RCN-GONA site overlaps across all the spectral bands. For the red band, despite a larger difference between the RCN-GONA data and the GONA-EPICS data, the TOA reflectance of GONA-EPICS encompasses the mean TOA reflectance of the RCN-GONA data, and both datasets are statistically comparable within the uncertainty, as shown by the ${{\chi }^2}_{red}$ results. Additionally, the fact that the ${{\chi }^2}_{red}$ values are substantially lower than the one for most bands suggests that the uncertainties of the observed data might be overestimated.

Figure 17 summarizes the reduced chi-square (${{\chi }^2}_{red}$) results obtained for all four sensors used in this study. The ${{\chi }^2}_{red}$ values are substantially lower than one for all the bands, except for the red band in the Landsat sensors. Although the ${{\chi }^2}_{red}$ value for the red band is slightly above 1 (1.14), it still indicates that the observed data (GONA-EPICS) and the model (RCN-GONA) are statistically the same within the GONA-EPICS uncertainty.

Additionally, it is important to note that the potential overestimation of the uncertainty for the other bands may be attributed to the fact that the total uncertainty used to perform the reduced chi-square test includes the absolute radiometric calibration uncertainty of the sensors. For this analysis, the calibration uncertainties were set according to the specifications for each sensor (3% for Landsat 8 and Landsat 9, and 5% for Sentinel-2A and Sentinel-2B). The low ${{\chi }^2}_{red}$ values suggest that the actual absolute radiometric calibration of these sensors may be better than the specified requirements.

Welch’s t-Test Results

The p-values obtained after the application of Welch’s t-test to all four sensors considered in this study are shown in

Table 5. The p-values obtained after the application of Welch’s t-test to the Landsat 8, Landsat 9, Sentinel-2A and Sentinel 2B GONA-EPICS TOA reflectance vs. RCN-GONA TOA reflectance.

Sensor	p-Value
Band	CA	Blue	Green	Red	NIR	SWIR1	SWIR2
L8	0.987	0.925	0.764	0.467	0.992	0.926	0.968
L9	0.986	0.953	0.770	0.589	0.904	0.979	0.822
S2A	0.856	0.874	0.836	0.734	0.993	0.926	0.783
S2B	0.842	0.955	0.744	0.688	0.961	0.975	0.817

.

The results show a consistent failure to reject the null hypothesis across all the sensors (Landsat 8, Landsat 9, Sentinel-2A, Sentinel-2B), indicating that the mean differences between GONA-EPICS and RCN-GONA are statistically insignificant, meaning they are statistically the same. The high p-values obtained in all cases suggest that the observed differences are likely due to random variation rather than actual differences between the datasets. This consistency not only reinforces the statistical similarity of the datasets but also highlights the reliability of GONA-EPICS as potential calibration targets aligned with RCN-GONA. Therefore, the findings support the use of GONA-EPICS as calibration sites, confirming their suitability for radiometric calibration across multiple sensors.

Difference Between the RCN-GONA and the GONA-EPICS

To quantify the differences between the RCN-GONA TOA reflectance and the GONA-EPICS, three metrics were employed: mean error (ME), mean absolute error (MAE), and RMSE. The mean error (ME) measures the average error between the RCN-GONA TOA reflectance and the GONA-EPICS TOA reflectance, providing insight into the direction of the error (positive or negative). In addition, the mean absolute error (MAE) represents the average magnitude of the errors between the RCN-GONA and the GONA-EPICS, regardless of their direction. Finally, the RMSE determines the average size of the errors, giving more weight to larger differences and highlighting the potential impact of the remaining outliers across all the spectral bands. The ME, MAE, and RMSE were calculated using Equations (15)–(17).

$$ME={\displaystyle \frac{1}{n}} {\displaystyle \sum }_{i=1}^n(RC{N}_i-EPIC{S}_i)$$

(15)

$$MAE={\displaystyle \frac{1}{n}} {\displaystyle \sum }_{i=1}^n\left|RC{N}_i-EPIC{S}_i\right|$$

(16)

$$RMSE={\displaystyle \frac{1}{n}} {\displaystyle \sum }_{i=1}^n{(RC{N}_i-EPIC{S}_i)}^2$$

(17)

where n is the number of observations, RCN is the MSG-interpolated RCN-GONA data and EPICS is the mean TOA reflectance of the GONA-EPICS cluster.

Figure 18 and

Table 6. The ME, MAE, and RMSE between RCN-GONA and GONA-EPICS for all the sensors and all the spectral bands.

	Sensor	CA	Blue	Green	Red	NIR	SWIR1	SWIR2
ME	L8	−0.0003	0.001	0.006	0.017	0.001	0.004	−0.0005
	L9	0.001	0.003	0.008	0.017	0.003	0.005	−0.005
	S2A	−0.004	−0.004	0.005	0.010	−0.0002	−0.002	−0.008
	S2B	0.002	−0.0001	0.006	0.010	−0.003	−0.001	−0.006
MAE	L8	0.009	0.010	0.011	0.018	0.011	0.014	0.018
	L9	0.009	0.010	0.012	0.018	0.012	0.014	0.019
	S2A	0.011	0.011	0.011	0.014	0.011	0.016	0.018
	S2B	0.009	0.009	0.009	0.013	0.009	0.011	0.015
RMSE	L8	0.011	0.012	0.014	0.022	0.015	0.019	0.022
	L9	0.011	0.013	0.015	0.022	0.015	0.019	0.023
	SA	0.014	0.014	0.013	0.018	0.015	0.020	0.022
	S2B	0.012	0.011	0.011	0.016	0.012	0.013	0.019

show the ME, MAE and RMSE for all the sensors for all the spectral bands considered in this study. The ME values show distinct differences between the sensors when comparing GONA-EPICS to RCN-GONA. For Landsat 8, the ME ranges from −0.0003 in the CA band, indicating a very close agreement to RCN-GONA, to 0.017 in the red band, where GONA-EPICS reports slightly higher reflectance values compared to RCN-GONA. The smallest difference for L8 was seen in the CA band, and the highest difference is seen in the red band. Landsat 9 follows a similar pattern, with ME values ranging from 0.001 in the CA band to 0.017 in the red band. In contrast, Sentinel-2A shows both positive and negative differences, with ME values ranging from −0.008 in the SWIR2 band, where GONA-EPICS reflectance is higher than RCN-GONA, to 0.005 in the green band. The smallest difference for Sentinel-2A is in the blue band, where the ME is −0.0002, indicating very close agreement between GONA-EPICS and RCN-GONA. Sentinel-2B also presents mixed results, with ME values ranging from −0.006 in the SWIR2 band to 0.006 in the green band. The smallest difference for Sentinel-2B is in the blue band, where the ME is −0.0001, showing a very good agreement between the two datasets.

In addition, the MAE values reveal that Landsat 8 and Landsat 9 tend to have slightly higher errors in comparison to Sentinel-2A and Sentinel-2B. For Landsat 8, the MAE ranges from 0.009 in the CA band to 0.018 in the red and SWIR2 bands. Similarly, L9 exhibits MAE values from 0.009 in the CA and blue bands to 0.019 in the SWIR2 band, where the CA band exhibits the lowest difference and the red band the highest difference for both sensors. In addition, Sentinel-2A presents a consistent performance, with MAE values from 0.011 in the CA, blue, and green bands to 0.018 in the SWIR2 band. Sentinel-2B generally performs slightly better, with MAE values ranging from 0.009 in the CA, blue, green, and NIR bands to 0.015 in the SWIR2 band.

The RMSE values further highlight the differences between GONA-EPICS and RCN-GONA when using different sensors. For Landsat 8, the RMSE values range from 0.011 in the CA band to 0.022 in the red and SWIR2 bands, indicating larger differences in these bands when comparing GONA-EPICS to RCN-GONA. Landsat 9 shows similar results, with RMSE values ranging from 0.011 in the CA band to 0.023 in the SWIR2 band, demonstrating slightly higher discrepancies in the SWIR2 band compared to Landsat 8.

In contrast, Sentinel-2A exhibits generally lower RMSE values, ranging from 0.013 in the green and SWIR1 bands to 0.022 in the red and SWIR2 bands. This suggests a closer alignment between GONA-EPICS and RCN-GONA for the Sentinel-2A sensor, especially in the green and SWIR1 bands. Furthermore, Sentinel-2B demonstrates the best performance among the sensors, with RMSE values as low as 0.011 in the blue and green bands, and peaking at 0.019 in the SWIR2 band. These values are consistently lower than the maximum RMSE values observed in the Landsat sensors, indicating that GONA-EPICS aligns more closely with RCN-GONA for Sentinel-2B.

Even though the metrics presented show slight differences between RCN-GONA and GONA-EPICS for some bands, the highest difference observed is less than 0.023 reflectance units across all the bands. This difference is minimal, particularly when considering the significant advantages offered by GONA-EPICS. With over 50 calibration points per cycle for a sensor like Landsat 8, GONA-EPICS provides distinct benefits for radiometric calibration and stability monitoring, enabling faster stability assessments compared to traditional methods such as PICS. This makes the slight differences in reflectance between the two datasets less significant in the context of the broader advantages provided by GONA-EPICS, especially when rapid stability monitoring or calibration of newly launched sensors is required. Additionally, this work demonstrates through the validation of GONA-EPICS using RCN-GONA that the EPICS identified in this study can deliver reliable results. Although the TOA reflectance from these EPICS does not originate from ground measurements, the validation shows that GONA-EPICS provides statistically equivalent TOA reflectance to the RCN-GONA site within their uncertainties, as indicated by the reduced chi-square results, with the added benefit of global coverage.

4. Conclusions

This study conducted a global land cover clustering of pixels with similar spectral characteristics using an unsupervised clustering algorithm. A key difference between this work and previous continental and global analyses is that the data were filtered to include only temporally stable pixels before applying the K-means clustering algorithm. In contrast, other studies did not evaluate the per-pixel temporal stability before the clustering stage. This approach was used to identify clusters with similar temporal and spectral characteristics that can be used for radiometric calibration efforts.

This study identified 160 clusters with different spectral characteristics. Of these, two were evaluated in detail. The first, named Cluster 13-GTS, is equivalent to Cluster 13 identified in the global analysis by Fajardo et al. [16] (Cluster 13-G). Cluster 13-G was obtained from a global land cover clustering that did not filter out any pixels before the application of the K-means, resulting in the inclusion of temporally unstable pixels and a reduction in the number of useful locations. The authors had to develop a methodology to evaluate the temporal stability of each location individually after the clustering stage. In contrast, Cluster 13-GTS, identified in this work, was obtained by applying a per-pixel temporal stability filter to the data cubes input to the K-means clustering algorithm, ensuring that only temporally stable pixels were retained and avoiding the need for post-classification evaluations. A comparison between Cluster 13-G and Cluster 13-GTS indicated that the new Cluster 13-GTS could reduce the coefficient of variation (CV) by up to 1% compared to Cluster 13-G while providing more than 50 locations, compared to only 19 in Cluster 13-G. This demonstrates that the new EPICS obtained in this work offers a larger number of calibration points with lower temporal variability, allowing for faster evaluation of the temporal stability of optical sensors with reduced uncertainties, as the CV directly affects uncertainty estimations.

Additionally, this study validated one of these EPICS using the TOA reflectance derived from ground-truth measurements collected at the RadCalNet Gobabeb (RCN-GONA) site, using data from Landsat 8, Landsat 9, Sentinel-2A, and Sentinel-2B. The validation was performed using a cluster that includes the RadCalNet Gobabeb pixels, named GONA-EPICS. The results showed that GONA-EPICS could provide statistically comparable mean TOA reflectance to that measured by RCN-GONA, as evaluated using a reduced chi-square test. This test indicated that RCN-GONA is statistically indistinguishable from GONA-EPICS within the cluster’s uncertainty. Furthermore, the reduced chi-square test suggested a potential overestimation of the uncertainty of GONA-EPICS, which ranged between 6.4% and 10.3%, with higher uncertainty observed for the Sentinel sensors. These uncertainties could be lower, considering the absolute radiometric calibration uncertainty used in this study was 3% for the Landsat sensors and 5% for the Sentinel sensors, as specified by mission requirements. However, the potential overestimation of the uncertainty identified through the reduced chi-square test may indicate that the actual calibration uncertainty of these sensors is lower than the mission requirements for both Landsat and Sentinel-2 missions. Additionally, the root mean square error (RMSE) between RCN-GONA and GONA-EPICS indicated that the difference between these datasets is less than 0.023 units of reflectance across all the bands. The advantage of GONA-EPICS lies in its 80 potential calibration points per Landsat cycle, providing a rich dataset in a short period and enabling faster evaluation of optical sensor stability. Although GONA-EPICS has higher uncertainty than the RadCalNet, it remains a valuable option when RadCalNet Gobabeb is unavailable due to cloud cover or other factors. GONA-EPICS also demonstrates the advantage of aggregating multiple pixels, reducing the likelihood of temporal instability, unlike RadCalNet data, which may not be suitable for long-term temporal stability evaluation.

The EPICS shown in this work highlight the potential of these global targets for radiometric calibration and stability monitoring of optical satellite sensors. Many of these EPICS require further exploration, as they provide a wide dynamic range and are located in various regions worldwide. These newly developed global targets offer several advantages, such as increasing the number of calibration points for any sensor in any orbital cycle and reducing the dependency on single locations like traditional PICS. This approach offers a new and reliable method, as demonstrated in the validation presented in this work, to monitor optical sensors, ensuring the timely provision of reliable data for scientific applications.