Hyperspectral Empirical Absolute Calibration Model Using Libya 4 Pseudo Invariant Calibration Site

: The objective of this paper is to ﬁnd an empirical hyperspectral absolute calibration model using Libya 4 pseudo invariant calibration site (PICS). The approach involves using the Landsat 8 (L8) Operational Land Imager (OLI) as the reference radiometer and using Earth Observing One (EO-1) Hyperion, with a spectral resolution of 10 nm as a hyperspectral source. This model utilizes data from a region of interest (ROI) in an “optimal region” of 3% temporal, spatial, and spectral stability within the Libya 4 PICS. It uses an improved, simple, empirical, hyperspectral Bidirectional Reﬂectance Distribution function (BRDF) model accounting for four angles: solar zenith and azimuth, and view zenith and azimuth angles. This model can perform absolute calibration in 1 nm spectral resolution by predicting TOA reﬂectance in all existing spectral bands of the sensors. The resultant model was validated with image data acquired from satellite sensors such as Landsat 7, Sentinel 2A, and Sentinel 2B, Terra MODIS, Aqua MODIS, from their launch date to 2020. These satellite sensors differ in terms of the width of their spectral bandpass, overpass time, off-nadir viewing capabilities, spatial resolution, and temporal revisit time, etc. The result demonstrates the efﬁcacy of the proposed model has an accuracy of the order of 3% with a precision of about 3% for the nadir viewing sensors (with view zenith angle up to 5 ◦ ) used in the study. For the off-nadir viewing satellites with view zenith angle up to 20 ◦ , it can have an estimated accuracy of 6% and precision of 4%. predict TOA reﬂectance in these bands adequately. The plots illustrate that the model has predicted TOA reﬂectance sufﬁciently in 11 Sentinel solar reﬂective bands as both the datasets accumulated along the unity line.


Introduction
The absolute radiometric calibration refers to converting an image into digital numbers (DN) or the recorded voltage by a satellite sensor into physical quantities. Those physical units can be of at-sensor spectral radiance (W m −2 s r −1 µm −1 ) or apparent top of the atmosphere reflectance [1,2].

Absolute Calibration
Satellite instruments are usually calibrated before launch. Unfortunately, no matter how well these instruments are, their performance may degrade with time once in space. This degradation could be due to thermal, mechanical, slow deterioration of the electronic system or UV radiation exposure. So, it is necessary to characterize and assess the sensor's performance throughout its lifetime. Absolute calibration can be a primary tool to evaluate the sensor's performance from the very pre-launch stage to throughout its on-orbit operations. Continuous absolute radiometric calibration of the sensors will improve the data interpretability and quality of remotely sensed data. It will ensure that the satellite sensor system's image data is acceptable and radiometrically accurate to its intended user community. Apart from absolute calibration, pre-launch and post-launch on-orbit calibrations are used for image data analysis and monitor their radiometric response over time. To validate mission requirement performance and follow SI traceability, series of components are calibrated and characterized during prelaunch calibration. Post-launch calibration includes assessing onboard calibrators, vicarious calibration, cross-calibration between satellite sensors, etc. [3][4][5][6][7]. In case of Vicarious calibration, the limited frequency of data collections had been improved by RadCalNet but it had come at the expense of field equipment at calibration sites and can be labor-intensive [8]. On the other hand, Pseudo-Invariant Calibration Sites (PICS) have been widely adopted to perform cross-calibration and monitoring sensor stability [1][2][3][4]9] overcoming the limitations of some calibration techniques.

PICS Stability and Identification of Stable Areas
Pseudo invariant calibration sites have been proven as the least expensive method of on-orbit calibration. Over the last 14 years, a significant expansion in PICS use has been observed to monitor the long-term top-of-atmosphere (TOA) reflectance trends from different sensors [9][10][11]. These PICS sites are known for their spatial uniformity by manifesting stable spectral characteristics over time, higher reflectance, and negligible atmospheric impact on upward radiance [12]. Cosnefroy et al. chose twenty desert sites that had 3% or better spatial uniformity. After eliminating the directional effects, temporal variability of 1% to 2% was found [13]. The Committee on Earth Observation Satellites (CEOS) entitled six North African PICS: Libya 4, Mauritania 1, Mauritania 2, Algeria 3, Libya 1, and Algeria 5, and they were exhibiting 3% or less temporal variability throughout all the bands [14]. In 2010, Helder et al. concluded that Libya 4, Libya 1, and Algeria 3, Arabia 2, Egypt 2, and Egypt 1 were competent enough to monitor the long-term trends and had variability of less than 3% [15]. Later, South Dakota State University Image Processing Laboratory (SDSU IP Lab) did the stability analysis to find optimal regions exhibiting 3% or less temporal, spatial, and spectral variability. Optimal regions within Libya 4, Niger 1, Sudan 1, Niger 2, Egypt 1, and Libya 1 contemplated less temporal uncertainty. Among them, Libya 4 demonstrated as the most temporally stable site with a temporal variation of less than 3% and could be broadly used in radiometric calibration work [1,3,10,16].

PICS for Absolute Calibration Model Development
In 2004, Govaerts et al. developed an absolute calibration method using bright desert calibration sites [17]. In 2012, the model was improved by including atmospheric polarization, consideration of spheroidal aerosol particles, and the surface bidirectional reflectance function. This model had an accuracy of 3% over Libya 4 with the largest error in the blue channel [18]. Helder et al. developed an absolute calibration model over Libya 4 (PICS) in 2010 employing Terra MODIS and EO-1 Hyperion sensors. The model demonstrated accuracy within 3% in the visible and 6% in the shortwave infrared (SWIR) region [3]. A desert daily exoatmospheric radiance model (DERM) founded on a well-calibrated geostationary Earth orbit (GEO) sensor over a single PICS was developed by Bhatt et al. in 2013 [4]. The reference Meteosat-9 DERM and ray-matched calibration consistency were within 0.4% for Meteosat-8 and 1.9% for Meteosat-7. Likewise, GOES-10 and GOES-15 were calibrated using the GOES-11 DERM followed by a consistency within 1% and 3%, respectively [9].
In 2014, Mishra et al. expanded the work of Helder et al. on the absolute calibration model considering the BRDF effect caused by the off-nadir measurements of the sensor and seasonal variations due to solar position change. This empirical absolute calibration model had an accuracy of the order of 3% with an uncertainty of about 2% across all bands [2]. However, Terra MODIS (reference radiometer) collection 6 data products were used in this absolute calibration model. The research work of empirical absolute calibration model by Mishra et al. [2] was extended by Raut et al. [1] in 2019 for the other five additional Saharan Desert PICS. For Egypt 1, Libya 1, and Sudan 1, the model had an estimated accuracy of approximately 3% along with a precision of approximately 2% for the sensors used for the analysis. However, Raut et al. found Niger 1 and Niger 2 sites provided less accuracy with similar precision because of an insufficient amount of reliable Hyperion data over these PICS sites [1]. As the work of Raut et al. was an extension of the work of Mishra et al., they followed the same procedures to structure the model. Moreover, they used Terra MODIS (reference radiometer) collection 6.1 data products for the model development.

Required Improvement in the Absolute Calibration Model
Few improvement features are needed to consider in the previous IPLab derived absolute calibration model like (1) using pixel-based angle information instead of scene center angle information, (2) utilizing four angles information: solar zenith and azimuth, view zenith and azimuth angles, (3) conversion of the angles from Spherical to Cartesian coordinate, (4) developing a hyperspectral BRDF model. Helder et al. (2013), Mishra et al. (2014), and Raut et al. (2019) developed empirical BRDF models that utilized scene center-specific Spherical angles. Precisely, only the solar zenith angles of Terra MODIS and view zenith angles of Hyperion data over PICS were used. Then the scene-center-specific TOA reflectances were plotted as a function of solar zenith and view zenith angles for 6 reflective bands of Terra MODIS. Finally, the model incorporated the TOA reflectances as a simple, continuous linear function of solar zenith angle and a quadratic function of view zenith angle, and hence the coefficients were calculated. However, these models were generalized in such a way that they could be used for all existing spectral bands [1][2][3]. However, their true capacity of working on the existing non-Landsat or MODIS equivalent spectral bands was never tested or validated.
In 2016, the Landsat archive was reorganized by the USGS into a tiered collection management structure titled Landsat Collection 1. It provided a reliable archive of known data quality to support pixel-level time-series analyses and data processing advancements [19]. However, the most current data processing level comes with angle files consists of per pixel four angles information: solar zenith and azimuth and view zenith and azimuth angles. To cope with the latest data type, an improvement in the absolute calibration model is required by using pixel-wise angles information instead of scene center-specific angles.
Farhad et al. proposed a new Bidirectional Reflectance Distribution Function (BRDF) model that considered all four angles: solar zenith and azimuth, view zenith, and azimuth angles. They found that the solar and view geometries were best defined in the Cartesian coordinates to preserve the nature of the data and to achieve a robust fit for the BRDF model. Moreover, the four-angle multi-linear interaction model provided the best BRDF model categorization, and after normalization, the estimated temporal stability was better than 3% over Libya 4 [20]. Considering the research of Farhad et al., conversion of the angles from Spherical to the Cartesian domain and developing an empirical BRDF model using four angles rather than just two angles can be contributed as improvement features in the absolute calibration model. Mainly, pixel-wise four angles BRDF information in cartesian coordinate can be a better representation of the angles.
Another important question arises on the previous SDSU IPLab derived absolute calibration models that whether the BRDF coefficients (derived as a simple continuous function of solar and view zenith angle from 7 spectral bands) are real or not. More importantly, what if there are bands in between the multispectral bands which have absorption features that are not as expected as the previous models suggested. This limitation can be overcome by utilizing a hyperspectral BRDF model which is the main objective of this study.

Landsat 8 as Reference Radiometer
When the absolute calibration models were developed firstly in SDSU IPLab, Terra MODIS was considered as one of the best-calibrated sensors in the reflective bands with an uncertainty of 2% in TOA reflectance [3]. However, Landsat 8 Operational Land Imager (OLI) sensor has comparable absolute radiometric accuracy as Terra MODIS and provides higher spatial resolution than Terra MODIS. Additionally, L8 OLI captures data with improved radiometric precision over a 12-bit dynamic range, which is the same as Terra MODIS and has an overall improved signal-to-noise ratio than Terra MODIS [21,22]. It is important to add that considering the L8 data availability in the SDSU IPLab image archive and as L8 OLI has similar or better radiometric accuracy as Terra MODIS, L8 OLI is employed as a reference radiometer for developing the new absolute calibration model.

Hyperspectral BRDF Model
Over the years, many empirical, semi-empirical, and physical BRDF models have been developed. Usually, this information is not offered at a hyperspectral resolution over a large range of wavelengths. Any well-calibrated sensor with hyperspectral imaging capacity can be a useful source to structure the BRDF effect over a wide range of wavelengths. A good source of the hyperspectral system: Hyperion EO-1 has 220 unique spectral channels varying from 0.357 to 2.576 micrometers along with a 10-nm bandwidth and also has an absolute calibration capacity of 5% [23]. Forming a BRDF model using Hyperion helps to overcome the limitation of interpolating BRDF coefficients from the multispectral level to the hyperspectral one and can provide true coefficients even in absorption bands. As a result, the model can truly be capable of predicting the TOA reflectance of sensors with non-Landsat equivalent bands such as Yellow (605.37 nm), Red-Edge (740.5 nm).

Objectives of the Work
This article presents a new empirical hyperspectral absolute calibration model for Libya 4 PICS (hyperspectral APICS) using L8 OLI and Hyperion EO-1. L8 OLI sensor has an absolute radiometric capacity with stated uncertainty of better than 3% [24]. Hence, it has been chosen as a reference radiometer. Hyperion EO-1 has been used as a hyperspectral source to develop a hyperspectral BRDF model. This hyperspectral APICS model is developed using an improved pixel-based four angles hyperspectral BRDF model. The four angles hyperspectral BRDF model exploits the solar and view geometries converted from Spherical to the Cartesian domain, unlike spherical angles used by Raut et al. and Mishra et al. This model accompanies the truly hyperspectral signature of BRDF coefficients interpolated from 10 nm to 1 nm spectral resolution. Thus, the hyperspectral APICS model is capable of performing absolute calibration in 1 nm spectral resolution and can predict TOA reflectance in all existing spectral bands of any sensor. The hyperspectral BRDF model will be validated using Landsat 7 ETM+ (L7), Landsat 8 OLI (L8), Sentinel 2A MSI (S2A), and Sentinel 2B MSI (S2B) (Red-Edge bands), Terra and Aqua MODIS. This paper is organized into four sections. The introduction section provides a brief review of absolute calibration and some of the earlier research work done using PICS to develop an absolute calibration model. Section 2 discusses the methodology used to develop the proposed hyperspectral absolute calibration model. Section 3 presents the model development and validation of results for the Libya 4 PICS site. Finally, Section 4 summarizes the article and provides potential routes for future work in this area.

Sensors Overview
The Landsat legacy has been supplying uninterrupted acquisition of high temporal resolution and multispectral data covering globally since 1972. Landsat 8, part of the Landsat series was launched on 11 February 2013 constitutes the OLI sensor and Thermal Infrared Sensor (TIRS). The OLI has been serving high-quality image data for Earth surveillance. The pre-launch calibration of the L8 OLI had an assessed uncertainty of about 3% in reflectance products. Moreover, post-launch calibrations of L8 OLI have been consistently verified uncertainties on the order of 2% or less [22]. Landsat 7 Enhanced Thematic Mapper Plus (ETM+) was launched on 15 April 1999 flies at a mean altitude of 705 km in a sun-synchronous orbit with an equatorial crossing time of 10.00 A.M. having 8 spectral bands. Before the launch of L8 OLI, the L7 ETM+ was considered the most stable of the Landsat series, with stated uncertainties of 5% [25]. In this study, Landsat Collection 1 data from their launch date to 2020 was used. The new features in Collection 1 data come with solar illumination and sensor viewing angle information, additional metadata, a quality assessment bands.
The pixel values in the region of interest (ROI) of the Libya 4 site were converted to TOA reflectance by using linear scaling factors specified in the associated product metadata for the L8 OLI, and L7 ETM+ sensors. The TOA reflectance value was calculated as below [21,26]: where, M ρ and A ρ are band-specific, reflectance-based multiplicative, and additive scaling factors, respectively, Q Cal is the calibrated DN pixel value and ρ λ is the estimated TOA reflectance. However, a cosine correction is essential as these coefficients do not comprise solar zenith angle (SZA) like in Equation (2): where α indicates the solar zenith angle.

Hyperion EO-1
The Earth-Observing One (EO-1) satellite was a hyperspectral sensor that was designed just for a one-year mission in orbit. It was launched on 21 November 2000 as part of NASA's New Millennium Program [27]. The purpose of the push-broom hyperspectral sensor was to provide high-quality calibrated hyperspectral data. Hyperion imaged across 242 bands, among which 196 bands were onboard calibrated spectral ranging from 400-2500-nm, at a minimal 10-nm spectral resolution with 30-m spatial resolution. Hyperion was a pointing satellite, unlike Landsat 8 and MODIS. It was capable of imaging up to 25 degrees from nadir view. In February 2011, the satellite was wiped out of fuel essential to maintain orbit. This incident introduced an adjustment in precession rate that directed to gradually earlier equatorial crossing times throughout its last five years. This orbital drifting came along with the earlier overpass times. It changed from approximately 10:00 A.M. to approximately 7:00 A.M. local time resulting in increasing solar zenith angles in the acquired images [1,23]. However, the satellite was decommissioned on 20 March 2017. Franks et al. reported ALI or Hyperion had no pronounced trend throughout 2016 and atmospherically corrected reflectance products are within 5 to 10% of mean pre-drift products. As a result, Hyperion is still considered as a high-quality hyperspectral resource until the end of the mission [23]. For this study, Hyperion EO-1 data from 2000 to 2017 was used. The pixel values of the image of EO-1 Hyperion were converted using Equation (3).
where, DN cal is known as the calibrated DN pixel value, h is scaling factor, d 2 is the earth to sun distance in A.U. unit, E sun is the radiance calibration converted to reflectance based on ChKur solar spectrum (ESUN(ChKur)) [28], φ is the sun elevation angles, and θ is the sensor look angle. The main goal of this sensor is to provide better spatial resolution (10 to 60 m) image data. Barsi et al. [24] proved that OLI and MSI had shown stable radiometric calibration, with approximately~2.5% consistent spectral bands matching. For this study, Sentinels' data from their launch date to 2020 was used.
By using Equation (4) the MSI sensors pixel values were converted to TOA reflectance by dividing it with a constant scaling factor. This scaling factor considered the exoatmospheric irradiance, Earth-sun distance, and cosine correction: where DN cal is considered as the calibrated DN pixel value and the value of Q = 10,000.

MODIS
A key part of NASA's (National Aeronautics and Space Administration, Washington, DC, USA) Earth Observing System MODIS is an instrument on-board the Terra and Aqua satellites. Terra MODIS was launched on 18 December 1999, whereas Aqua MODIS was launched on 4 May 2002. MODIS has 36 bands the highest temporal resolution obtains data at three spatial resolutions: 250 m, 500 m, and 1 km, which are coarser than the other sensors used in the study. A ±49.5 • scanning pattern at the Earth Observing System (EOS) orbit of 705 km has a 2330-km swath and near-daily revisit acquisition capability. Terra MODIS TOA reflectance products have approximated calibration uncertainty of 2-3% [29]. MODIS Collection 6.1 data having view zenith angle up to 20 • from their launch date to 2020 was used for this analysis. Table 1 represents the angular variation along with the number of scenes of EO-1 Hyperion, L8 OLI, L7 ETM+, S2A and S2B MSI, and Terra and Aqua MODIS sensors over the Libya 4 site.

Study Area (PICS) and ROI
The investigated study area Libya 4 PICS was situated in the Saharan Desert of North Africa. The South Dakota State University Image Processing Lab (SDSU IPLab (Brookings, SD, USA)) developed the PICS normalization process (PNP) algorithm. They combined the L8 OLI observations of several PICS into a single time series with a larger temporal resolution for satellite calibration. The PNP algorithm was employed to the L8 OLI image data to choose "optimal" regions within the Libya 4 PICS demonstrating 3% or less temporal, spatial, and spectral variability [16]. For developing the model, a sub-region at the scene center was chosen in the optimal region where Hyperion and L8 images were overlapping with each other. Another sub-region at the edge of the scene was taken to validate the model performance. Figure 1a,b shows the image of Libya 4 with Landsat 8 and Hyperion EO-1. Figure 1c represents the PICS optimal regions (white pixels) identified for Libya 4 site and the red rectangle are the selected scene center ROI as a sub-region within the optimal region. In Figure 1d, the green rectangles are the selected ROI at the edge of the scene as a sub-region within the optimal region. Table 2 provides the location of the selected region of interest. The scene center ROI for L8, L7, EO-1 was used to build the hyperspectral APICS model. To validate the model performance other three ROIs were used. To spatially match target and reference data for inter-calibration pixel observations for each ROI were averaged and only averaged data from 4 ROIs were used.

Data Preprocessing
All the L7 ETM+, L8 OLI, S2A MSI, S2B MSI, EO-1 Hyperion, Terra MODIS, and Aqua MODIS image data of Libya 4 PICS used in this analysis were retrieved from the existing SDSU IPLab archive. L7 ETM+, L8 OLI as Landsat Collection 1 (LC1) products, and EO-1 Hyperion image datasets as Level 1 Terrain (L1T) products were previously downloaded to the SDSU archive through the United States Geological Survey (USGS, Reston, VA, USA) Earth Explorer (https://earthexplorer.usgs.gov/ (accessed on 27 March 2020)). It should be noted that the S2A MSI and S2B MSI image data were processed using Sentinel-2 Level 1C tiles retrieved from the Copernicus Open Access Hub (https://scihub.copernicus.eu/ (accessed on 2 October 2020)). For MODIS (Terra, Aqua) data, collection 6.1 image data products were accessed from MODAPS web service (https://modaps.modaps.eosdis. nasa.gov/ (accessed on 16 August 2020)). All the downloaded image products were preprocessed by each group to remove radiometric and geometric artifacts. The OLI, and MSI products were then scaled to 16-bit integer digital numbers whereas L7 ETM+ was delivered as 8-bit images. The L7 ETM+, L8 OLI, S2A MSI, S2B MSI, and EO-1 Hyperion image data were converted to TOA reflectance using the conversion coefficients listed in the associated product metadata and XML files. Additional details about the preprocessing steps can be found on the corresponding website. Figure 2 shows the flowchart of the procedures to determine the hyperspectral APICS model. Every step is discussed in the subsections.

Data Filtering
After calculating the TOA reflectance value from the specified ROI of each image, filtering was necessary to guarantee a cloud-free image. L7 ETM+, L8 OLI image data were filtered using the associated band quality assessment information. In the case of S2A MSI and S2A MSI, the cloud-mask product was used. The MSI cloud mask was available at a spatial resolution of 60 m, 10 m, and 20 m spatial resolution for each corresponding spectral band. For L8 OLI, Hyperion EO-1, Terra MODIS, and Aqua MODIS sensors, an empirical 2-sigma (±2σ) filtering approach (i.e., 2 standard deviations from the mean of the temporal TOA reflectance derived using all scenes) were applied for eliminating potential outliers. For the time series of S2A MSI, S2B MSI, L7 ETM+ sensors, ±1 temporal standard deviation filter was applied to detect the outliers. If any image's mean TOA reflectance derived from a predefined ROI, crossed the chosen threshold, a visual inspection of the image for all spectral bands was taken. Then if the visual inspection suggested it as clouds, shadows, or other artifacts that were not identified in the quality data, the entire scene (all spectral bands) was excluded from further analysis. Figure 3a is the mean hyperspectral reflectance profile of the target ROI over Libya 4 before applying the ±2 sigma filtering approach and Figure 3b represents the hyperspectral reflectance profile after using the ±2 sigma filtering approach. It is visible from Figure 3b that outliers are removed.

Drift Corrections to Hyperspectral Data
Because of mechanical stresses in the course of the launch, maneuvering in a harsh space environment, and aging of the sensor itself, the satellite sensors can show fluctuations in their radiometric response. Xin Jing et al. suggested a statistically significant drift in Hyperion EO-1 sensor response in the bands 8 to 16 based on the study of Libya 4 image data acquired from 2004 to final decommissioning in 2017 [5]. Therefore, the hyperspectral profiles from EO-1 over Libya 4 used in this analysis were corrected by accounting for the possible drift in the sensor response. Further details are described below.

Yearly Drift Parameter Estimation and Correction
The percentage change in drift was developed as a linear function of days since launch, was calculated using Equation (5): where, %Drift λ /year (reflectance per year) is the percent degradation per year in a band λ. Slope λ (reflectance per days since launch) and Intercept λ (reflectance) are the Slope λ and Intercept λ coefficients acquired from a least-squares linear regression of TOA reflectance in band λ as a function of days since launch. The reflectance in bands exhibiting potential drift was corrected as follows: where ρ λ,drift_corr is the Hyperion TOA reflectance after yearly drift correction, ρ λ is the TOA reflectance of Hyperion, yr represents the decimal year since launch and %Drift λ /year is the percentage yearly drift of any band λ exhibiting considerable degradation estimated using Equation (5) [5].

Four Angles BRDF Modeling
As we know most of the Earth's surface exhibits itself as a non-Lambertian target, the TOA reflectance of a given target can vary significantly with solar illumination and sensor viewing geometry. Usually, this effect can be modeled by the Bidirectional Reflectance Distribution Function (BRDF). However, sensors with a larger field of views such as the Advanced Very High-Resolution Radiometer (AVHRR) and MODIS (approximately ±49.5 • ) may exhibit significant BRDF effects. This phenomenon demands a BRDF model, in turn, which will affect the estimated TOA reflectance [20,30]. Hyperion EO-1 sensor's view to illumination angle observations were used to develop an empirical hyperspectral BRDF model. The solar zenith and azimuth angles and view zenith and azimuth angles were converted from a Spherical to a linear Cartesian coordinate. Then mirroring of the data to each quadrant was performed to achieve symmetry with respect to the scattering plane and to get a robust fit to the BRDF model. As a result, it had led the TOA reflectance as a continuous function of independent variables [20]. Equations (7) and (8) were used to convert the solar geometry and Equations (9) and (10) were used to convert the sensor view geometry from the Spherical domain to Cartesian coordinate.
where SZA and SAA are the solar zenith and azimuth angles, respectively, and VZA and VAA are the sensor viewing zenith and azimuth angles, respectively. SZA, SAA, VZA, and VAA were set of angles in the Spherical coordinate system and X 1 , Y 1 , X 2 , Y 2 were generated through conversion. Figure 4a shows the SWIR1 band TOA reflectance of Hyperion EO-1 plotted against the solar geometry in spherical coordinate after performing the mirroring of data into each quadrant. The SZA ranging from (−70 • to 70 • ) and SAA (−160 • to 160 • ). Mirroring assists to reduce the modeling error along with the edge effects due to the restricted operational variability in the view geometry [20]. Figure 4b corresponds to the view geometry in spherical coordinate having VZA ranging from (−20 • to 20 • ) and VAA (−280 • to 280 • ) after the mirroring. It can be observed in Figure 4a,b that the TOA reflectances are acting as four discrete data set. After converting the solar and view angles from Spherical to the Cartesian domain, in Figure 4c (solar geometry), d (view geometry) TOA reflectances appear as two continuous datasets with reduced edging effect. Farhad et al. recommended a multi-linear four-angles interaction model (15 coefficients model) [20] to study BRDF properties over PICS. To develop a simple empirical BRDF model for absolute calibration, an investigation to have appropriate BRDF parameters was performed. Table 3 presents the results of hypothesis tests on the full BRDF multilinear regression coefficients at the 95% significance level on Hyperion EO-1 band 185 (2284 nm) as a representative of all the available reflective bands. This statistical analysis was validated for all the Hyperion reflective bands except few bands at 942 nm (water vapor) and 1386 nm (cirrus).
The p-value for the t-statistic of the hypothesis tested that the corresponding coefficient was equal to zero or not. If the p-value was >0.05, that term became insignificant. It is visible from Table 3 that only the coefficients of X 1 2 , Y 1 2 , X 2 2 , Y 2 2 , X 1 . X 2 , Y 1 .Y 2 are the significant parameters in the model. However, the purpose of the task was to find a simple four-angle BRDF model from each scene based on a multilinear regression model. Hence, the model was derived using the solar zenith and azimuth angles, and sensor zenith and azimuth angles.  By assuming no interactions among these angles and based on the p-value, only highly significant four terms X 1 2 , Y 1 2 , X 2 2 , Y 2 2 should be taken into consideration. But from Table 3 we can see, the estimated standard error of the X 2 2 and Y 2 2 is −22.61 and 1.19 respectively, which indicates these parameters will provide significantly higher sensitivity in the model. Whereas the other two most significant parameters X 1 2 , Y 1 2 do not exhibit such higher standard errors. Generally, the BRDF is more driven by solar angles rather than view angles if the sensor has a narrow field of view [31]. Although Hyperion had the cross-track pointing capability to observe adjacent tracks, its limited view angles should have less impact on the BRDF model. In such a case if we consider the X 2 and Y 2 terms as quadratic components in the model, it may cause higher sensitivity and larger standard error to the model. However, considering higher standard error and huge sensitivity of the coefficients of X 2 and Y 2 terms were considered as linear components for the final model. Therefore, assuming no interactions between the angles the empirical four angles multiple linear least-squares regression model was derived as Equation (11): where ρ model is the model predicted TOA reflectance, C 0 (λ) is the derived intercept and C 1 (λ), C 2 (λ), C 3 (λ), C 4 (λ) are the BRDF model coefficients for X 1 2 , Y 1 2 , X 2 and Y 2 respectively. From Equation (11), it can be observed that the simple BRDF model has X 1 and Y 1 components as the quadratic terms and X 2 and Y 2 as the linear parameters.
This filter was settled on the underlying assumption that the actual coefficients dataset was not smooth and continuous over a long spectral width. Thus, the filter was employed to capture every distinct feature (including all the narrow bands) without distorting the signal tendency and increasing the precision. This modified Savitzky-Golay filter was used to subdivide a longer spectral window into a smaller one according to the chosen frame length (in terms of wavelength). Then the coefficients were interpolated according to step size for that small, specified window of chosen frame length. Finally, the interpolation was continued until the end of the spectral window. Other fitting tools like polynomial/linear interpolation/Fourier etc. requisite continuity of the dataset with a specific pattern over long spectral width for doing the interpolation. On the other hand, a modified Savitzky-Golay filter was sufficient to interpolate every distinct discontinuous reflectance and absorption feature in the dataset.
For this analysis "cubic interpolation" fit was chosen. It should be noted that the minimum size of frame length should be 10 as Hyperion EO-1 has 196 spectral channels ranging from 357 to 2576 nanometers with a 10 nm spectral resolution. However, an optimized frame length of 40 nm spectral resolution was considered for the analysis. Moreover, a step size equals to 1 was chosen as the coefficients were interpolated from 10 nm to 1 nm spectral resolution. Figure 5a,b represents the hyperspectral profile of X 1 2 coefficients, C 1 (λ) and Y 1 2 coefficients, C 2 (λ) respectively. They are plotted as a function of wavelength. The coefficients were interpolated using Modified Savitzky-Golay filter plotted as a function of each Hyperion EO-1 spectral bands in blue lines and red dots represents the actual X 1 2 coefficients and Y 1 2 coefficients value. A frame length of 40 nm spectral resolution, a 'cubic interpolation' fit, and step size of 1 nm (final output) were utilized.  Figure 6a,b represents the interpolated hyperspectral profile of X 2 coefficients, C 3 (λ) and Y 2 coefficients, C 4 (λ) respectively using the filter. The X 1 2 , Y 1 2 coefficients had shown more significance in the model manifesting higher values than, X 2 , Y 2 coefficients. As X 2 , Y 2 were structured as linear components, they had less significance in the model with a lower value. Figure 7 demonstrates the intercept, C 0 (λ) of the BRDF model which equals to the ρ h (λ) parameter. The ρ h (λ) was represented as the single hyperspectral profile of over Libya 4 for 196 bands.

Cross-Scale Factors
To link different sensors with diverging spatial, radiometric, and spectral resolutions into a common radiometric scale and to ensure data interoperability cross-calibration is necessary [31]. There are several techniques to do cross-calibration. One way is, analyzing the cloud-free imagery obtained from the coincident or near coincident scene pair approach from selected targets by two or more sensors. One of the two sensors is considered as a reference radiometer [32].
Hereby, for this analysis, cross-scale factors were calculated using Hyperion EO-1 and L8 OLI sensors. L8 OLI was used as a reference source for absolute calibration. To scale the Hyperion spectrum, an adjustment factor was calculated. After applying the adjustment factor, when the Hyperion TOA reflectance profile was integrated over the Landsat 8 spectral bandpass it had produced the same TOA reflectance estimate as Landsat 8. To do this, near-coincident scene pairs 6 days apart with solar and view zenith differences within ±5 • were considered. Because of that the scene overpassing time difference between the two sensors was not greater than 1 h 10 min. Therefore, this window had allowed less variation in the absolute cross-scale numbers. While following all the conditions only 14 near coincident events were found. For every pair of near coincident scene pair events, the Hyperion data were normalized using the 15 coefficients BRDF model [20] to make it angularly look like L8 OLI. Then to spectrally match Hyperion data to the seven L8 OLI solar reflective bands the band integration was performed in the Hyperion data. The cross-scale factors, K(λ) were determined as below: ρ Hyperion Banded to L8 (λ) = In Equation (12), the ρ L8 is the L8 TOA reflectance determined from a set of L8 and Hyperion near-coincident image pairs of 6 days apart. The parameter, ρ Hyperion Banded to L8 is the Hyperion simulated TOA reflectance from the near-coincident pair set that is normalized to angularly looks like L8 and then band-integrated or weighted according to L8 Spectral Response Function (SRF).
In Equation (13), the band-integrated TOA reflectance is estimated by integrating the SRF of the sensor with the Hyperion profile at each sampled wavelength, weighted by the corresponding SRF. Table 4 provides the mean K-scale factors with standard deviation derived for Libya 4 PICS calculated using Landsat 8 and Hyperion near coincident scene pairs. However, the cross-scale factor, K(λ) was treated as the function of wavelength. EO-1 has two different instruments with different fields of view operating in two different channels. Therefore, for the range of Hyperion VNIR bands corresponding to (426.82-925.41 nm, B008-B057) and for the range of the SWIR bands corresponding to (912.45-2395.50 nm, B077-B224) cubic interpolation was exploited separately to determine the K-scale values. As the hyperspectral APICS model was designed to facilitate absolute calibration in a 1 nm spectral resolution scale, the K-scale values were also interpolated into 1 nm spectral resolution.

Hyperspectral Absolute Calibration Model
The hyperspectral APICS model is represented as Here, ρ(λ, X 1 , Y 1 , X 2 , Y 2 ) are the predicted absolute TOA reflectance for a selected sensor for specified angles in Libya 4. X 1 , Y 1 , X 2 , Y 2 are the cartesian angles converted from spherical to cartesian coordinate. C 1 (λ) , C 2 (λ) are the quadratic BRDF coefficients for X 1 2 and Y 1 2 respectively. C 3 (λ), C 4 (λ) are considered as linear BRDF coefficients for X 2 and Y 2 respectively which is described in Section 2.7. The parameter ρ h (λ) is equal to C 0 (λ) specified as the derived hyperspectral BRDF intercept. Lastly, we have the hyperspectral cross-scale factor, K(λ) to normalize the Hyperion intercept spectral profile, ρ h (λ) into different sensor levels. It should be noted that the atmospheric parameter was not considered in the model. Barsi et al. found that image data with exact pixel-based angle information had a negligible impact on the atmospheric model upon the absolute calibration model [24].
For the assessment of this work, the model generates three metrics: (1) a percentage difference, Equation (15), determined as the ratio of the model-predicted and observed reflectance differences to the observed TOA reflectance; (2) normalized RMSE in Equation (16), estimated as a ratio of the root-mean-square of the model-predicted and observed reflectance differences to the observed TOA reflectance; and (3) a precision in Equation (17), calculated as the standard deviation of the model-predicted and observed reflectance differences to the mean of the observed TOA reflectance.

Validation with Landsat 8
Since, the hyperspectral cross-scale factor, K(λ) was calculated based on Landsat 8 and Hyperion observations, Landsat 8 data was used for the initial validation of the model. Figure 8 illustrates the comparison of the hyperspectral APICS model predicted TOA reflectance and the corresponding observed L8 bands. It is visible from the plot that the model predicted TOA reflectances are higher than the observed ones in all 7 bands. The cross-scale factors were contributing to scale the Hyperion spectrum into L8 OLI level. Apparently, the inherent bias was due to the used cross-scale number. Due to having insufficient Hyperion image data only a sample of 14 coincident scene pairs were used to calculate the cross-scale factors. This sample seemed insufficient to make the bias zero.  Table 5 shows the summarized results of mean absolute percentage difference, accuracy, and precision for different solar reflective L8 bands. The normalized RMSE (NRMSE), or accuracy between the model and L8 measured reflectance was under 3% in all bands, and precision was estimated well within 2.5% for Libya 4. Among all the bands, the model had exhibited a better result in the Red band. Figure 9a compares the hyperspectral APICS model-predicted and measured TOA reflectances in the L8 Red band (band 4) over Libya 4. The seasonal or angular variation was captured quite well by the model. As the best-case scenario, the normalized RMSE or bias found in this band was 1.95%, and the model predicted higher reflectance values comparing to the observed values. The precision describes the variation found in the respective band while measuring the same target repeatedly by the model and was estimated at 0.96%. However, in the Red band, the mean absolute percentage error between the model and measurement was 1.75%. It was also visible that in the middle of summer when the solar zenith angles were lower, the difference between the predicted and measured TOA reflectances got higher indicating strong solar angular dependency of the model.  Figure 9b represents the histogram of the corresponding percentage difference between model prediction and measurement in the Red band. The histogram plot tells that, the model works well within 3%. Most of the data, approximately 60%, lied within a 1-2.5% percentage difference from the measured TOA reflectances over Libya 4 with L8 OLI.

Validation with L7
For validating the model, the observations from different sets of well-calibrated satellite sensors were used. In this subsection, an assessment will be pictured between the Libya 4 model prediction and at-sensor reflectance derived from L7 ETM+ measurements. Figure 10 illustrates the comparison of the hyperspectral APICS model predicted TOA reflectance and the corresponding observed L7 bands. The plot illustrates that the model has predicted TOA reflectance adequately in all L7 bands as they are well distributed alongside the reference line. From the 1:1 line, it is visible that the model predicts a higher value than the observed reflectance in the Blue band. L7 and L8 sensors have their inherent dissimilarities. L7 has a coarser radiometric resolution which is 8-bit whereas L8 provides higher 12-bit image products with more precise geometry [21,26]. However, a significant scatter was observed in the L7 measured reflectances than in the L8 data which had caused higher precision numbers in the predicted model. This was probably due to the larger relative spectral response of L7 than L8. Table 6 exhibits the summarized results of mean absolute percentage difference, accuracy, and precision for six different L7 bands. The normalized RMSE (NRMSE), or accuracy between the model and L7 measured reflectance was about 3% or less in all bands. The estimated precision got higher comparing to L8 but was well within 3% for Libya 4. For the L7 sensor, the model had a better result in the Green band (accuracy 1.39% and precision 1.34%). But a detailed result of the NIR band will be explained to show how the model captures the water absorption feature. In Figure 11a the magenta circles show the absolute calibration prediction, and the black circles demonstrate ETM+ measurements for the NIR band. The model constantly followed the seasonal trend and predicts reflectance levels sufficiently. However, significant scatter, due to the inclusion of a water vapor absorption feature at 850 nm was observed in the measured reflectances in this band. This water absorption feature was not captured by the model contributing to a precision value of 2.20%. Besides, approximately 1-2% random variability was attributed mainly to the site's spectral behavior and atmospheric disagreements [33]. The assessed accuracy of the model was approximately 2.28%, well within the estimated 3% accuracy for Libya 4 and the calibration uncertainty of the ETM+ (approximately 5%) [25].  Figure 11b represents the histogram of the corresponding percentage difference between model prediction and measurement in the NIR band. The mean absolute percentage error between the model and measurement was 1.19%. The histogram plot for percentage differences was normally distributed and ranged from −6% to 6%. However, the model was good within ±3% as the majority, 81% of data tended to lie within ±3% percentage difference over Libya 4 with L7 ETM+.

Validation with S2A and S2B
In this subsection, an evaluation will be shown between the Libya 4 model prediction and at-sensor reflectance derived from S2A MSI and S2B MSI sensors for 11 bands. Figure  12a,b demonstrates the comparison of the hyperspectral APICS model predicted TOA reflectance and the corresponding observed S2A and S2B reflective bands, respectively. It should be highlighted that band 5 to band 8 are non-Landsat equivalent bands. No other APICS models developed at SDSU IPLab [1][2][3] were ever validated using the existing non-Landsat or MODIS equivalent spectral bands. But this new hyperspectral APICS model due to its design nature was able to predict TOA reflectance in these bands adequately. The plots illustrate that the model has predicted TOA reflectance sufficiently in 11 Sentinel solar reflective bands as both the datasets accumulated along the unity line.  Table 7 exhibits the summarized results of mean absolute percentage difference, accuracy, and precision for 11 S2A and S2B bands. The normalized RMSE (NRMSE), or accuracy between the model and measured reflectance was within 3.75% and 3% for S2A and S2B respectively. However, the inherent difference between S2A and S2B had caused the dissimilarity in the result of accuracy as well. S2A radiometry was slightly brighter than S2B. Moreover, inter-sensor calibration between S2A and S2B had provided an agreement of about 1-2% bias in the VNIR and none in the SWIR [34]. For both sensors, the precision was better than 2.5% for 11 solar reflective bands. The limitation of using L8 as the reference radiometer was that it does not have any Red-Edge bands and a wider NIR band like Sentinels. As a result, the interpolated cross-scale factors used in the model may not represent the true value. Thus it showed a higher bias in the few Red-Edge and wider NIR bands. Another important factor was that the model was generated over a desert site, not over any vegetative site having red edge spectra. As a result, the spectral response of Hyperion for the red edge band over a bright desert site was not a true representation of Red-Edge spectra.
S2A and S2B MSI narrow Red-Edge band located at 740 nm results was chosen as a representative of the non-Landsat spectral band. Figures 13a and 14a compare the modelpredicted (in red circle) and measured TOA reflectances (black circle) in the Red-Edge band over Libya 4 of S2A and S2B respectively. Though it was a non-Landsat equivalent band, the model was predicting TOA reflectances in this narrow band quite well maintaining overall seasonal variability for both the sensors. The corresponding estimated model accuracy and precision were better than 2% for both the sensors.  From Figure 14a it is visible that from 2019 to 2019.5 the predicted TOA reflectance is lower compare to measured values for S2B. At the same time frame, in Figure 13a, there is less variation in S2A data. It is because at the starting of each month radiometric calibrations are performed regularly in S2A data and the decontamination operations are usually scheduled once a year in September which ensures less variation in the radiometric response of S2A data. On the other hand, for S2B data the decontamination was performed only in November 2019 resulting in less variation of the radiometric response after November 2019 [35]. Figures 13b and 14b represent the histogram of the corresponding percentage difference between model prediction and measurement in the Red-edge band of S2A and S2B respectively. The mean absolute percentage error between the model and measurement was approximately 1.34% and 1.53% for S2A and S2B, respectively. The histogram plots for percentage differences were quite symmetric for both the sensors and varied up to ±4%. For both cases, the majority, 80% of data were within ±2% percentage difference, indicating the fact that the model worked well within ±2% over Libya 4 with Red-Edge band located at 740 nm.

Validation with Terra and Aqua MODIS
This subsection will present an assessment between the Libya 4 model prediction and at-sensor reflectance derived from Terra and Aqua MODIS sensors for all the available bands. Figure 15a,b illustrates the comparison of the hyperspectral APICS model predicted TOA reflectance and the corresponding observed Terra and Aqua MODIS bands, respectively. It is visible from the figures that Terra and Aqua MODIS data for different bands are not quite following the 1:1 reference line and they are poorly spread around the line. Moreover, the predicted TOA reflectance was lower in the shorter wavelength bands and higher in the SWIR channels. Both Terra and Aqua MODIS have a higher field of view (FOV) of ±49.5 • with a 2330 km wide viewing swath. Due to the nature of capturing images of Terra and Aqua MODIS, three different viewing geometries levels were observed in the measured TOA reflectance. Unfortunately, the model fails to capture these three different levels quite correctly and as a consequence, the data did not follow along the 1:1 line, especially in the shorter wavelengths.
However, a higher residual error was observed in Terra and Aqua MODIS data which had made bias and precision of the predicted model worsen. Unlike Terra MODIS, Aqua MODIS is an afternoon constellation. It has equator crossings in a northerly direction in the early afternoon at about 1:30 P.M. local solar time. On the other hand, EO-1 was a morning constellation, having equator crossings in the late morning around 10:30 A.M. Consequently, Aqua MODIS has higher solar azimuth angles comparing to Hyperion EO-1. The empirical BRDF model was built with Hyperion data and the solar azimuth angles of Aqua MODIS data were found different than the empirical BRDF model trained for. Hence, the derived empirical BRDF model may not be sufficient to model Aqua MODIS data accurately. Figure 16a,b illustrate the polar plots of Hyperion and Aqua MODIS solar and view geometries, respectively. Actual Hyperion, mirrored Hyperion, and actual Aqua MODIS angles are represented in black, yellow, and red dots, respectively. It is visible from the plots that there is no overlapping within the actual view and solar geometries of the Hyperion and Aqua MODIS. However, after mirroring, the view geometry of Aqua MODIS seems to overlap with mirrored Hyperion data. But for solar geometry, only a small portion of Aqua data is overlapping with the Hyperion's. It is worth remembering that the actual BRDF model was developed by mirroring the Hyperion data into four quadrants to reduce the edging effect. But for this scenario, when there is no actual overlapping between two sets of data and due to the absence of real data, the model coefficients become less useful to predict correctly.
For Terra MODIS, the angles were within the range with the range of the angles that were used to generate the BRDF model. Still, the model was unable to predict TOA reflectance distinctively for three different view geometries in Terra MODIS as shown in Figure 17a.  Table 8 shows the results of the mean absolute percentage difference, accuracy, and precision of both MODIS sensors. The normalized RMSE (NRMSE), or accuracy between the model and Terra MODIS measured reflectance was within 5% except for the Blue band and for Aqua MODIS it was within 4% in all bands. For Terra, the precision was found less than 4.5% and for Aqua, it was within 3%. A common channel for both sensors, the NIR band has been chosen for further illustration. Figure 17a compares the model-predicted and measured TOA reflectances in the Terra MODIS NIR band over Libya 4 from 2002 to 2020. Three sets of measured TOA reflectances were observed due to different ranges of view zenith angles. The observed TOA reflectance due to view zenith angle (1 • -5 • ), (7 • -11 • ), and (13 • -16 • ), data has been subdivided into three colors in green, black, and blue respectively. The model had predicted TOA reflectances preserving overall seasonal variability except for the three different levels of TOA reflectance and atmospheric scattering in the observed data. The estimated model accuracy and precision were approximately 2.34% and 2.17%, respectively. Figure 18a compares the model-predicted (magenta circle) and measured TOA reflectances (black circle) in the Aqua MODIS NIR band over Libya 4. The observed TOA reflectance due to view zenith angle (1 • -6 • ), (7 • -10 • ), and (16 • -17 • ), data has been subdivided into three colors in green, black, and blue, respectively. Again the model is not performing quite well to capture the three different levels of TOA reflectance in the observed data. The model had constantly predicted higher than the observed data in the NIR band. Aqua and Terra MODIS have almost the same spectral bandpass. The cross-scale adjustment factor used in the hyperspectral APICS model was also the same. But then for Aqua MODIS, the accuracy became 3.34% whereas Terra MODIS had an accuracy of 2.34% in the same band. However, the only difference was in the range of angles found for both the sensors. For Terra MODIS, the solar zenith angle was (16 • − 55 • ), the solar azimuth angle was (100 • − 167 • ), view zenith angle was (1.3 • − 17 • ) and view azimuth angle was (98 • − 292 • ). Whereas, in Aqua MODIS the solar azimuth angle ranged from (198 • − 260 • ) and other angles were quite similar to Terra. Conversely, there were no such larger solar azimuth angles found in the Terra data. Moreover, Aqua MODIS solar azimuth angles were outside the range of the angles that were used to generate the BRDF model as seen in Figure 16a. Hence, this scenario suggested an improvement in the BRDF model as well. Figures 17b and 18b characterize the histogram plot of the corresponding percentage difference between model prediction and measurement of the NIR band in Terra and Aqua MODIS respectively. The mean absolute percentage error between the model and measurement for Terra and Aqua was 1.82% and 2.89%, respectively. For Terra MODIS the histogram bin of percentage differences followed a normal distribution and the model performed well within ±5% over Libya 4 with NIR band. For Aqua MODIS the histogram plot of percentage differences had a distribution with larger residues varied up to ±7%. Overall, the model had operated well within ±4% over Libya 4 with NIR band.
However, the absolute calibration model developed by Helder et al. over Libya 4 demonstrated accuracy within 3% in the visible and 6% in the shortwave infrared region [3]. The empirical model developed by Mishra et al. had an accuracy of the order of 3% with an uncertainty of about 2% across all bands [2]. This new hyperspectral APICS model could predict the satellite measurements with a 3% or better systematic error and precision for nadir viewing satellites. On the other hand, for satellites with higher view geometries, the model was able to predict the measurements with 6% or better systematic error and precision. It can be said that for Nadir viewing satellites the new model has consistency with the prior absolute calibration model focused on Libya 4.

Uncertainty Analysis
Uncertainty is a parameter related to the result of a measurement, that depicts the dispersion of the values that could reasonably be attributed to the measurand [36]. For this uncertainty analysis, the cross-scale factor, intercept (ρ h ), and BRDF uncertainty were considered as primary sources of uncertainty.
Finding the correlation between variables is an important step in uncertainty analysis. For this analysis, an assumption of no correlation between variables was assumed. By assuming the identified uncertainty sources are statistically uncorrelated, following the International Standards Organization Guide to the Expression of Uncertainty in Measurements (ISO-GUM) method [36], the final uncertainty of the model can be stated as Equation (18). It is to be noted that the combined/total uncertainty is considered as one standard deviation or a coverage factor of k = 1.

Cross Scale Factor Uncertainty
There was a difference between Landsat 8 and Hyperion sensor revisit times. This difference created a variation/uncertainty in solar and sensor view geometries and atmospheric conditions. For finding the uncertainty of the cross-scale factor, geo-registration error, site non-uniformity, BRDF coefficients uncertainty, and BRDF model error uncertainty were considered as major sources of cross-scale factor uncertainty.

•
The Geo-Registration Error Uncertainty The geo-registration uncertainty was associated with the uncertainty if different ROIs were used each time. The geo-registration error uncertainty was stated as 0.026% [20].

•
Site Non-Uniformity The uncertainty of the Libya 4 site non-uniformity was calculated from each cloud-free scene using Equation (19) [36]. The number of observations was found 4000 which was the total number of pixels in the ROI.
where, s(q) = Variance of the mean, s(q k ) = Experimental standard deviation of the mean, n = Number of observations. However, the uncertainty of the Libya 4 site non-uniformity was less than 0.05% in all reflective bands.

• BRDF Coefficient Uncertainty
Uncertainty in the BRDF coefficients for each band was also calculated as the difference between the predicted TOA reflectances of 14 Hyperion coincident scenes multiplied with cross-scale gains and the associated L8 coincident scene pairs. This source of uncertainty considered the variability of the cross-scale gain and the band integration uncertainty (as an implicit factor). From Table 9 we can see that the estimated uncertainty ranged from 0.8% to 1.41%. CA band had the largest uncertainty value. A better result could be expected if there were consistent, reliable Hyperion image data available. Both the sensors had the same spatial resolution of 30 m so that no spatial resolution mismatch uncertainty was considered.

• BRDF Model Uncertainty
To find out the cross-scale factors, 15 coefficients BRDF [20] model were used to normalize the Hyperion data. This BRDF model uncertainty was associated with the 15 coefficients BRDF model. It had encompassed sensor overpass time difference and atmospheric variation. The BRDF model error uncertainty was calculated using the RMSE as a magnitude of the error. It can be written as Equation (20).
where, X obs = Observed TOA reflectance and X model = Predicted TOA reflectance.
The BRDF error uncertainty appeared as the most significant source of uncertainty contributing to the cross-scale factor uncertainty. The worst-case scenario was found on the SWIR2 channel estimated as 3.22%. The BRDF effect was more prominent in the higher wavelength bands which had explained the higher uncertainties in the SWIR1 and SWIR2 bands. Again, the shorter wavelength bands were more prone to atmospheric change such as aerosol optical depth and water vapor content. Resulting in higher uncertainty observed in CA (2.92%) and Blue (2.89%) bands.
Overall, the final cross-scale factor uncertainty was found approximately within 2.0% to 3.5%. The highest amount of cross-scale uncertainty was estimated in the SWIR2 band which was approximately 3.47%, details are in Table 9.

Intercept Uncertainty
As discussed in the methodology section the intercept was derived from the multilinear regression BRDF model. The associated intercept uncertainty was calculated as the standard error of the intercept coefficients. The standard error measures how precisely the model estimates the unknown coefficient values. The standard error of the estimated coefficients was calculated as the square root of the diagonal elements of the matrix. A detailed explanation will be found here [37]. From Table 10 we can see that the intercept uncertainty is less than 0.5% in all bands.  (11). This source had included uncertainties coming from solar and sensor position change, atmospheric condition change. From Table 10, it can be seen that BRDF model uncertainty has ranged from 2.70% to 4.84%. Again, due to the presence of aerosol content in the atmosphere, the CA and Blue bands showed larger uncertainties of 4.84% and 4.43%, respectively.

Sensor Uncertainty
The uncertainty coupled with the individual sensor calibrations was the final factor in this uncertainty analysis. The currently recognized uncertainties in the OLI and Hyperion calibrations are approximately 2% and 5%, respectively [5,22]. The total sensor uncertainty was 5.38%. Table 10 represents the estimated uncertainties for each source of uncertainty. By assuming no correlation between the variables, after propagating all the sources and by using Equation 18, the final estimated uncertainty was calculated. The final uncertainty ranged between 6% to 8%.

Model Specification and Limitations
The main drawback of exploiting an empirical model was that a large number of data samples were needed to shape the statistical model. The empirical hyperspectral APICS model was built with 363 cloud-free Hyperion scenes. This model was designed/trained with the solar zenith angle (20 • to 70 • ), solar azimuth angle (80 • to 160 • ), viewing azimuth angle (100 • to 280 • ), and viewing zenith angle of up to 20 • obtained from 363 Hyperion acquisitions. Therefore, it was expected that the model was not able to perform well beyond those trained angles.
However, to find out the specification of the model, the model was validated by different nadir viewing sensors choosing two different ROIs, one at the scene center and one at the edge of the scene. Two different ROIs were used for Sentinels and Landsats as Sentinels had smaller footprints than Landsats. These ROIs locations were mentioned in Table 2. All results were tabulated in Tables 11 and 12. As for Terra and Aqua MODIS, the performance of the model was categorized into three different sets of view zenith angles, as seen in Tables 13 and 14.   Table 14 represents the residual error between Aqua and Terra MODIS observed and the predicted data for three different levels of view geometries as shown in Figure 16b. As the worst-case scenario, in Terra MODIS, for view zenith angle ranged from (14 • -17 • ), the residual error was from (−11% to 6%). Again, in Aqua data with similar view zenith angles the residual error was (−3% to 7%). However, it strongly suggests a further study to improve the hyperspectral APICS model to accommodate sensors with a larger field of view.
Another important limitation is the use of multispectral sensor L8 to get the actual cross-scale factors for 7 spectral bands which were further interpolated into hyperspectral cross-scale factors. It does not have any Red-Edge band and a wider NIR band like S2A/S2B MSI does. It had led the interpolated hyperspectral cross-scale factors to become less stable on those bands. Incorporating the cross-calibration step using non-Landsat equivalent spectral bands may improve the overall model performance.
Moreover, while cross-calibrating between L8 and Hyperion, near coincident scene pairs of 6 days apart with solar and view zenith differences within ±5 • were chosen. L8 satellite launched in 2013. Hyperion launched in 2000 and was decommissioned in 2017. In 2011 Hyperion started drifting as the spacecraft ran out of fuel and introduced deorbiting and resulting in an earlier equatorial crossing time from 10.00 A.M. to about 7.00 P.M. [23]. Whereas the equatorial crossing time of L8 is about 10.00 A.M. However, a cross-calibration requires two sensors looking at the same site and same time to lower the BRDF effect and atmospheric variability in the acquired image. Therefore, by maintaining all the conditions, finding quality image pairs within these 4 years windows of these two sensors became a difficult task. Consequently, higher uncertainties and random errors were observed. From Table 9 it is clear that the cross-scale factors induced uncertainty is ranged from 2.1% to 3.5%. A larger number of cross-calibration opportunities with a quality product will be a helpful parameter to reduce the uncertainty in the model. With a recent study, in 2019, Shrestha et al. [38] generated a total of 185 hyperspectral profiles of the northern African region that had a similar spectral profile as Libya 4 (Clus-ter13). These findings could help to find near coincident scene pairs between Hyperion and Landsat 8 to perform cross-calibration.

Conclusions
Among all the Saharan PICS Libya 4 was established as the most stable site with a temporal variation of less than 3%. It is considered a reliable source to do the satellite's calibration assessment including the absolute calibration. This article demonstrates a new empirical hyperspectral absolute calibration model development using Libya 4 PICS and inspects if it can be used for absolute radiometric calibration of satellite sensors. Any hyperspectral instrument like CLARREO which will be launched in 2023 can use this model to calibrate the sensors. Moreover, even satellites with smaller footprints can use this model to calibrate the sensors over Libya 4 stable region.
This new hyperspectral absolute calibration model has used image data from the Landsat 8 OLI and EO-1 Hyperion sensors over Libya 4. Landsat 8 OLI was recognized as a well-calibrated reference radiometer, and Hyperion image data was used as a hyperspectral source. This model provides few improved features such as utilizing pixel-based angle information in the Spherical domain to obtain actual solar and view geometries over an ROI, an improved four angles empirical hyperspectral BRDF model, and exploitation of finer resolution BRDF coefficients (interpolated from 10 nm to 1 nm spectral resolution). Hence, the hyperspectral APICS model has equipped with the ability to perform absolute calibration applicable to any sensor with 1 nm spectral resolution over the stable region of Libya 4 PICS.
The model has been validated with L7 ETM+, L8 OLI, Sentinel 2A and 2B MSI (including three Red-Edge bands and the wider NIR band), and Terra and Aqua MODIS. For nadir viewing satellites, for solar zenith angle (15 • -60 • ), solar azimuth angle (95 • -160 • ), view zenith angle (0 • -5 • ), and view azimuth angle (55 • -280 • ), results showed that the model could predict the satellite measurements with better than 3% systematic error except for one red edge band of S2A. The precision was generally 3% or better across the spectrum from visible to Red-Edge and SWIR regions. On the other hand, for satellites like Terra and Aqua MODIS with larger view geometries of (1 • -17 • ), the model was able to predict the measurements better than 6% systematic error and 4% precision.
Moreover, a thorough analysis was performed to estimate the total uncertainty of the model. The total uncertainty ranged from approximately 6% to 8%. Nevertheless, improvements in hyperspectral APICS model accuracy and precision and reductions in uncertainty are possible. The low accuracy of the models is most likely driven by calibration differences between Landsat 8 and the other sensors. Random errors are most likely due to atmospheric conditions at the sensor overpass time. In general, a good quality hyperspectral image to do the cross-calibration can be a key factor to improve the accuracy.
The major drawback of the model is, it is restricted to view geometries up to 20 • and no atmospheric scattering like Rayleigh scattering, aerosol optical load, and gas absorption features are accounted for before modeling. There is a scope to improve the model by utilizing a great number of data samples specially focused on view geometry to shape the empirical BRDF model and it will lead to a more robust correction. Moreover, using a full atmospheric model developed through meteorological observations and climatological data can be performed to improve the model.  Acknowledgments: Thanks to Dennis Helder and Mahesh Shrestha for providing invaluable suggestions to clarify few concepts. At last, we enthusiastically acknowledge all the reviewers for their invaluable comments and suggestion for improving the manuscript.

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