Feasibility of Ground-Based Sky-Camera HDR Imagery to Determine Solar Irradiance and Sky Radiance over Different Geometries and Sky Conditions

Abstract

We propose a methodological approach to provide the accurate and calibrated measurements of sky radiance and broadband solar irradiance using the High Dynamic Range (HDR) images of a sky-camera. This approach is based on a detailed instrumental characterization of a SONA sky-camera in terms of image acquisition and processing, as well as geometric and radiometric calibrations. As a result, a 1 min time resolution database of geometrically and radiometrically calibrated HDR images has been created and has been available since February 2020, with daily updates. An extensive validation of our radiometric retrievals has been performed in all sky conditions. Our results show a very good agreement with the independent measurements of the AERONET almucantar for sky radiance and pyranometers for broadband retrievals. The SONA sky radiance shows a difference of an RMBD < 10% while the broadband diffuse radiation shows differences of 2% and 5% over a horizontal plane and arbitrarily oriented surfaces, respectively. These results support the developed methodology and allow us to glimpse the great potential of sky-cameras to carry out accurate measurements of sky radiance and solar radiation components. Thus, the remote sensing techniques described here will undoubtedly be of great help for solar and atmospheric research.

1. Introduction

In recent years, interest in the accurate measurements of solar radiation on the Earth’s surface has resumed. It was mainly triggered by the large increase in solar energy plants that have been deployed all over the world as part of the most ambitious strategies defined by different international administrations to mitigate climate change [1,2]. Therefore, there is now a need to improve the solar energy projections of these plants for different time scales. It is especially important for photovoltaic and thermosolar plants, due to their increasing development and perspectives over the coming years [3,4,5]. The success of these projections highly depends on improving solar radiation measurements, both in accuracy and worldwide availability. For short time scales, the solar energy forecast is strongly related to the atmospheric conditions and especially with the number of clouds and their temporal evolution [6]. Consequently, different instruments and measurement techniques have been developed.

The measurements of solar radiation on the Earth’s surface usually rely on broadband radiometers, covering the entire solar spectrum (280–3000 nm). These radiometers are commonly installed on solar trackers and the different components of the solar radiation are measured with moderate–high temporal resolution (~1–15 min). This instrumental setup is widespread, especially in the most populated areas around the world [7,8]. In addition, several methods to estimate the cloud cover, optical depth and aerosol optical depth (AOD) have been developed from the combination of different components of broadband solar radiation [9,10,11]. Unfortunately, in most ground-based stations, only the global horizontal irradiance (GHI) component is measured and the other components may be only estimated considering other atmospheric factors (e.g., clearness index) and using different parametric, physical or statistical methods [12,13]. This induces additional limitations on the spatial resolution of accurate ground-based measurements of all solar radiation components.

Conversely, geostationary satellites have been used to overcome these limitations because they provide global coverage with enough spatial resolution (~2–5 km) and moderately high temporal resolution (~10–30 min) that allows the determining of some cloud properties (cloud cover, cloud type and top height), its daily evolution and, consequently, the solar radiation throughout the day. The retrievals of solar radiation components on the surface rely on broadband satellite reflectance measurements and radiative transfer modelling. On the other hand, the spectral measurements of solar radiation from the ground and satellites have been used to determine the columnar content and properties of different atmospheric components (aerosol, cloud and gases) that significantly affect solar radiation on the Earth’s surface. Most of these measurements are also used to constrain solar radiation databases and forecasts provided from physical and statistical models [14,15,16]. The satellites and modelling estimations of solar radiation agree with ground-based measurements within ±20%, depending on the component [16,17,18,19]. The best results have been obtained in clear sky situations and therefore, there is room for the improvement of these estimations, especially when different cloud layers are present.

Sky-cameras combine the characteristics of broadband instruments, since they are able to measure the solar incident light in three RGB channels with the bandwidth varying from 20 to 100 nm, depending on the sensors and filters used. These cameras use array sensors and fisheye lens that allows individual pixels to observe portions of the atmosphere with a narrow field of view (FOV). Thanks to their panoramic view of the sky-dome, around 180° FOV, they were originally designed to determine the cloud cover. However, much progress has been made recently in instrumental development and data processing techniques. Consequently, the different properties of clouds (cloud amount, type and optical depth) have been obtained [20,21,22], as well as calibrated radiances for each channel [23]. In addition, the aerosol optical properties have also been inferred using different techniques [24,25,26,27].

The hemispherical FOV and high angular resolution make sky-cameras an interesting instrument to use to obtain the solar radiation over different geometries, with moderate spectral resolution. The attempts to obtain and forecast the broadband solar radiation are usually based on statistical correlations and machine learning methods that frequently determine the GHI and direct normal irradiance (DNI) components [19,28,29,30]. For solar energy applications, especially in those photovoltaic plants with fixed tilted angles in which the direct horizontal irradiance (DHI) is also important, solar radiation components need to be additionally projected on the plane of the array.

Therefore, the main objective of this work is to provide the accurate and calibrated measurements of sky radiance and broadband solar irradiance using High Dynamic Range (HDR) images from a sky-camera and, by taking advantage of its high angular resolution and spectral characteristics, develop a methodology to derive the sky radiance and the solar DHI component under arbitrary geometries. These methodologies and retrievals may be of great interest to the solar energy community, since they allow the monitoring and forecasting of solar radiation incidents on an arbitrary plane with a single and affordable instrumentation under all sky conditions.

The paper is organized as follows: the site and instrumentation are detailed in Section 2; an exhaustive characterization of the sky-camera is tackled in Section 3; Section 4 describes the methodological approach to derive the solar radiation over arbitrary geometries; we discuss our results in Section 5; and the main outcomes of this work are summarized in the conclusion.

2. Site Measurements and Data

The data used in this work were provided by the instrumentation deployed in the Burjassot Atmospheric Station (BASS), which is maintained by the Solar Radiation Group of the University of Valencia. It is located northwest of the city of Valencia, within its metropolitan area, and around 10 km away from the Mediterranean Sea. The BASS is devoted to measuring and monitoring radiation, aerosols, clouds and some other atmospheric components by means of a large set of remote sensing and in situ instrumentation [31,32,33,34,35]. It also forms part of different international measurement networks, such as the aerosol robotic network (AERONET) [36], European Skynet Radiometers [36], e-profile and EARLINET, and belongs to the ACTRIS research infrastructure.

The main instrument used in this work was a SONA-201D sky-camera integrated by SIELTEC Canarias S.L. and originally designed to determine the cloud cover using all sky images. The system was composed of a body camera based on a color CMOS sensor and a fisheye lens with a 180° FOV that was encapsulated in an IP67 (dust-tight and water immersion resistant up to 1 m depth) enclosure with a borosilicate dome, and thermally stabilized to maintain the internal temperature at 45 °C, regardless of environmental conditions. The camera sensor was made by e2v (model EV76C570), had a 1600 × 1200 pixel resolution and mounted an RGB filter with an infrared protection window, which gave an exposure time range from 16 µs to 1 s. This sensor had several light blocked lines that were used to implement an adaptive dark current filter, conferring a sensor signal to noise ratio (SNR) of about 38 dB. Its sensitivity and SNR were adequate for obtaining night images with a low dark noise by using a long integration time while also being able to register day images with a short integration time without risk to the sensor integrity. In spite of the raw images having a 10-bit depth, they were forced into an 8-bit truncation. The images were scheduled and stored through a SIELTEC proprietary server connection, based on a PC platform, using an Ethernet network and GiGe protocol.

We used an automatic Cimel CE-318 sun photometer that used a double collimator with a 1.2° FOV to measure the direct sun and sky radiance in different bands, from UV to near infrared spectral regions. The bandwidth for the UV channels was 2 nm, and 10 nm for the others. The direct sun intensity was measured at 340, 380, 440, 500, 675, 870, 940, 1020 and 1640 nm, and was used to determine the spectral aerosol optical depth (AOD), Angström exponent (AE) and water vapor content using the 940 nm channel. Sky radiance was also measured in seven channels (380, 440, 500, 675, 870, 1020 and 1640 nm) over different geometries. The combination of both measurement schedules was used to determine the aerosol’s microphysical and intensive radiative properties [37,38].

The measurements of the broadband solar radiation were carried out by Kipp and Zonen radiometers mounted on a Kipp and Zonen Solys-2 solar tracker: two CMP21 pyranometers for global and diffuse horizontal components and a CHP1 pyrheliometer for the direct. Additionally, the measurements of the broadband down- and upward radiation were performed in solar and thermal spectral regions by CM11 pyranometers and CGR4 pyrgeometers, respectively. Another CM11 pyranometer measured the solar broadband radiation on a 30° tilted surface, oriented to the South.

In this work, we have used data from February 2020 to June 2021, in which our sky-camera was set up to provide HDR images every minute. Broadband radiation was also stored at a 1-min resolution with two Campbell CR-1000 data-loggers. The Cimel CE-318 was configured to perform the direct sun measurements with ~10-min frequency while the sky radiance schedule was more variable, with a typical temporal resolution of 30 min, depending on the solar zenith angle (SZA).

3. Characterization of the SONA Sky-Camera

3.1. Geometric Calibration

Our sky-camera was set horizontally in a still location, so each pixel always received light from the same portion of the sky-dome. Therefore, several empiric methods were applied to properly characterize the fisheye lens of the sky-camera, in order to establish the observation geometry. This implied a projection of the celestial coordinates on the images in order to obtain the polar coordinates of every pixel and the position of the optical center, and to study the possible distortion pattern of the fisheye lens.

To geometrically characterize the sky-camera lens, we developed a Celestial Bodies Tracker (CBT) software. The CBT automatically identified a multitude of celestial bodies with well-known trajectories visible under the FOV of our camera and tracked them during the several days of observation. This allowed us to obtain enough samples to establish the relationship between the sky coordinates of the celestial bodies, taken from an ephemerides database [39], and their relative coordinates in the sky-camera polar coordinate system. The software provided the projection type, optical center and focal length, as well as the coefficients that related to the polar and Cartesian coordinates of our sky-camera.

The pixel zenith angle (PZA, θ) varied with the radius, measured from the optical center to each pixel (Figure 1a), and several fisheye projections [40] were used to correlate them by a regression fit. We chose the equidistant projection since it showed the best correlation coefficient (

Table 1. Regression fit results for different fisheye projections, where r represents the radius from the pixel to optical center, f is the focal length and θ is the PZA.

			Centre
Fisheye Projection		f (µm)	X ± 2	Y ± 2	R2
Equidistant	$r=f·\theta$	0.84	577	584	1.0
Orthographic	$r=2·f·sin\left(\theta /2\right)$	0.01	594	572	0.8
Stereographic	$r=2·f·tan\left(\theta /2\right)$	0.02	576	589	0.9

). Therefore, the optical center of the image was located at the pixel coordinates (X = 577, Y = 584) and the focal length was (f = 0.84 µm). The Pixel Azimuth Angle (PAA, φ) increased counterclockwise, with the North of the image (φ = 0°) in the upper part. There was an indication of a misalignment of 5.8° towards the East (left side) with respect to Geographic North (Figure 1b).

We performed additional accuracy tests of the CBT projections, which resulted in sky-coordinates position uncertainties of ±2 pixels.

Once the PZA and PAA were calculated, the pixel solid angle (PSA, Ω) could be obtained from Equation (1). The integration limits were derived for every pixel of the image from the PZA and PAA matrices (Figure 1a,b).

$$\Omega =\frac{dS}{r^2}={{\displaystyle \int }}_{{\theta }_1}^{{\theta }_2}{{\displaystyle \int }}_{{\phi }_1}^{{\phi }_2}\sin \left(\theta \right)d\theta d\phi$$

(1)

The PSA slightly decreased with the radius and showed small values in the range of [6.5, 10] × 10⁻⁶ sr (Figure 1c). Thus, each pixel of the camera observed a different portion of the sky-dome with a narrow FOV.

3.2. Sensor Characterization

Our sky-camera sensor was a Focal Plane Array (FPA) based on CMOS technology, which meant that each cell of the matrix was independent from the surrounding cells. Under illumination, each cell stored a charge, which was proportional to the irradiance received and created an electric current that could be measured by an analog-to-digital converter.

The main sources of uncertainty in this kind of sensor were the non-uniformity effects, due to the manufacturing process: dead and hot pixels, which worsen with the aging of the sensor, and thermal noise and dark current due to the physics of the device [26]. To estimate these effects, a set of measurements was collected in a dark environment. As expected, the non-uniformity, thermal noise and dark current effects were negligible (results not shown) thanks to the thermal stabilization and the dark current correction loop built into the sensor. The non-uniformity effect was also negligible, mainly due to the intrinsic characteristics of the sensor and its production process. The sensor did not present dead or hot pixels at exposure times below 25 ms. Considering that the maximum time that we used the sensor during day was about 36 µs, this issue only became apparent during nighttime measurements, which were not used in this study.

3.3. Camera Response Function

The camera response function, or “characteristic curve”, represented the sensor response to variations in the exposure that were the product of exposure time and the irradiance [41].

We empirically obtained the response function, as part of the HDR constructing method, following the methodology proposed in [41], which was based on the physical property of reciprocity. This method used a sequence of images to obtain the relationship between the exposure and its corresponding digital counts (DCs). Assuming that this function was monotonically increasing and capturing different exposures times of a static scene, we should have obtained different fragments of the curve (Figure 2a). Finally, a linear optimization process based on mean squared error minimization was used to fit these fragments together (Figure 2b).

Using this technique, the precise knowledge of the exposure times was critical to obtain the response function, and especially to build a HDR image properly. In our sky-camera, the nominal exposure time differed from the real exposure time. Thus, we calculated this non-linear effect by comparing the nominal and real exposure time ratios of two images taken with different exposure times. By doing this iteratively and varying the exposure times, we finally obtained the function that related the nominal to the real exposure times [42,43]. To improve the performance of the response function under any conditions, we recalculated it as the average of the response functions in different sky conditions, according to their corresponding SZA and considering different exposition times in each sequence (not shown).

The original algorithm [41] proposed a triangular weighting function to minimize the non-linear effects of the sensor when low or high DCs were measured, emphasizing the smoothness and fitting terms toward the middle of the curve at around 128 DC (Figure 3). In practical terms, we did not observe a significant dependency of the resulting response function with different weighting functions. However, the use and the choice of the appropriate weighting function were decisive in constructing the HDR images accurately (Section 3.4).

While the triangular function enhanced the contrast of the image, giving more weight in the middle of dynamic range, we used a weighting function w(z) in Equation (2). It gave the same weight within the entire linear range of the sensor (50 < DC < 200), reducing the weight in the non-linear region (DC < 50 and DC > 200) (Figure 3).

$$w\left(z\right)=1-{\left(2z-2^b+1\right)}^a$$

(2)

In Equation (2), z referred to the measured DCs, while b represented the number of bits of each sample (8-bit in our case) and a was an even factor greater than 2. The a coefficient modified the width of the plateau and the slope of the non-linear region of the function. The appropriate values for our camera were between 12 and 24 (Figure 3). To use this function properly, it had to be normalized with respect to the maximum in absolute value and the offset be adjusted to finally obtain a function that was worth 1 in the plateau and 0 in the tails of the function.

In looking for a good compromise between computational complexity and smoothness of result, we selected 900 pixel locations per sequence, where each sequence had different pixel locations and consequently, had its own response function. The pixel locations were chosen randomly but with restrictions based on our experience to guarantee the convergence of the fitting algorithm, thereby avoiding ”anomalous” pixels. These restrictions were that the pixel value had to be monotonically increasing (with at least a 33% margin) and the distance between pixels had to be greater than 16 pixels. If not enough pixels were suitable, the sequence was discarded.

To obtain the response function, we used 1000 images, with exposure times of between 0.016 ms and 1s, over several days under different conditions of illumination, in terms of SZA, and environmental conditions. Then the response function for the three channels (RGB) of our sky-camera was the average of all of them (Figure 4). The resulting average response function presented a non-linear behavior in the lowest and highest DCs range, specifically from 0 to 40 DCs for all channels and above 227, 211 and 215 for blue, green and red channels, respectively. The shape of the sensor response function did not depend on the number of pixels chosen to construct it, nor on their location. On the contrary, the use of anomalous pixels, such hot or dead pixels, significantly affected the shape of the curve (bumpy curve) and the use of insufficient pixel locations resulted in a poorly defined curve.

3.4. Constructing the High Dynamic Range Imagery

The function response g(Z_i,j) of the sky-camera could be used to convert the pixel values (in DCs) of any image into relative irradiance values Er_i (in arbitrary units) by knowing its exposure time Δt_j and following Equation (3) [41].

$$\ln E{r}_i=g\left(Z_{i,j}\right)-\ln \Delta t_j$$

(3)

Therefore, the HDR image could be constructed by a superposition of several images taken at different exposure times. The relative irradiance Er_i could be obtained by Equation (4) [41], where the sub-index i represented the position of each pixel in the corresponding image j within the sequence, while w(Z_i,j) was the weighting function.

$$\ln E{r}_i=\frac{{{\displaystyle \sum }}_{j}w(Z_{i,j})·\left(g\left(Z_{i,j}\right)-\ln \Delta t_j\right)}{{{\displaystyle \sum }}_{j}w(Z_{i,j})}$$

(4)

For daylight scenes, we took six images with times of between 16 µs and 36 µs, using 8 bits of bit-depth per channel. In addition, we used our weighting function (Section 3.3) to guarantee the same weight for all pixel values within the linear range of every image of the sequence.

The main advantage of using our HDR images was to largely expand the dynamic range within the scene, and the extremely linear behavior achieved, thanks to our improvements. In addition, we also obtained a significant reduction in the background noise, improving the response and the signal to noise ratio.

3.5. Radiometric Calibration

The radiometric calibration of the sky-camera was completed by a comparison against the radiative transfer simulations carried out by libRadTran model [44] and the relative irradiance from our HDR imagery, following the methodology proposed in [23,26]. It was based on the monochromatic simulations of the sky radiance, using the effective wavelength of each channel of the camera [45]. The effective wavelength (λ_eff) was spectrally weighted by the illumination conditions (i.e., the atmospheric irradiance), following (5).

$${\lambda }_{eff}=\frac{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}\lambda ·E\left(\lambda \right) · S\left(\lambda \right) d\lambda }{{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}E\left(\lambda \right) · S\left(\lambda \right) d\lambda }$$

(5)

where λ was the wavelength, E(λ) was the spectral irradiance reaching the instrument and S(λ) was the spectral response of the channel (Figure 5). The integration interval was from λ₁ = 400 nm to λ₂ = 998 nm. To calculate λ_eff, a set of 200 libRadTran simulations for different spectral irradiance corresponding to different clear sky conditions were carried out by changing the aerosol Angström parameters and the SZA [23]. The Henyey–Greenstein phase function [46] and the fixed values of the asymmetry parameter, single scattering albedo (SSA) and total ozone column (TOC) were used in all simulations.

The average values of the effective wavelengths for each channel of the sky-camera are shown in

Table 2. Effective wavelengths calculated for each channel of our sky-camera. The ± 1 standard deviation is assumed as the uncertainty.

Wavelength (nm)	λeff ± 1 σ(nm)
Blue	480 ± 6
Green	541 ± 5
Red	615 ± 6

.

The radiometric calibration established a relationship between the relative radiance provided by the HDR images (in arbitrary units) for each pixel (i.e., observation geometry) and channel (RGB) and the corresponding radiance values obtained from the RTM simulations in radiometric units (W m⁻² sr⁻¹ nm⁻¹). To calibrate our HDR images, we used 800 clear sky images from February to October of 2020, accounting for a large number of cases under different atmospheric conditions and sun locations.

We used monochromatic simulations using the λ_eff (Table 2) with an angular resolution of 1° for the viewing geometry (PZA and PAA) at the geographical coordinates of the BASS. The range of the variation of the most important parameters used in the RTM simulations are provided in the Supplementary Material. A standard mid-latitude atmosphere and the surface albedo were set as constant in all simulations. The surface albedo did not show a significant seasonal variation since the BASS is within the metropolitan area and the vegetation does not change in the surroundings. Therefore, it was obtained as the 5-year average, in the period 2015–2020, from the AERONET data. The concentration of the atmospheric components (TOC and column water vapor, CWV) and aerosol properties (AOD, SSA, asymmetry parameter and the phase function) were obtained from the AERONET database, as well as solar geometry that was varied in the simulations in correspondence with the selected images. The AERONET derived aerosol phase function was converted to 200 Legendre moments with the pmom tool in the libRadtran package before running the simulations. The AERONET data for all spectral values were interpolated to the λ_eff of our sky-camera.

We obtained a calibration factor k_λⁿ for each channel (sub-index λ) and each image (super-index n) fitting the simulated radiance L_λⁿ (i, j) by linear regression to match the relative Irradiance Er_λⁿ (i, j) following the Equation (6), where σ (i, j) was the sky-camera PSA and i and j were the matrix indeces.

$$L_{\lambda }^n\left(i,j\right)=k_{\lambda }^n·\frac{E{r}_{\lambda }^n\left(i,j\right)}{\Omega \left(i,j\right)}$$

(6)

Figure 6 shows the probability density function of the clear sky imagery calibration factors resulting from each camera channel. Finally, the calibration coefficient (K_λ) was obtained as the average of the whole sample of every k_λⁿ and ± 1σ_λ was considered its associated uncertainty (

Table 3. Calibration coefficients (K_λ) of each channel of the sky-camera and its absolute and relative uncertainty.

λ (nm)	$K_{\lambda } \pm {\sigma }_{\lambda } \left(\mathsf{\mu }\mathrm{W} c{\mathrm{m}}^{-2} \mathrm{n}{\mathrm{m}}^{-1}\right)$	Relative Uncertainty (%)
480	0.26 ± 0.02	7%
541	0.34 ± 0.02	7%
615	0.50 ± 0.05	10%

). The relative uncertainty of the calibration coefficients for the blue and green channels was 7%, and increased to 10% in the red channel.

4. Methods to Obtain Solar Radiation

4.1. Determination of the Sky-Radiance

Once the camera was calibrated, we knew the irradiance reaching every pixel of an image for the three RGB channels E_λ(θ, φ) and it was straightforward to determine the sky radiance considering the solid angle by following Equation (7).

$$L_{\lambda }\left(\theta ,\phi \right)=\frac{E_{\lambda }\left(\theta ,\phi \right)}{\Omega \left(\theta ,\phi \right)}$$

(7)

The radiance from any point of the sky-dome could be extracted from an image by choosing the appropriated viewing angles (θ, φ) and translating them into the Cartesian coordinates (x, y) of the image.

Once the sun coordinates (θ_s and φ_s) in an image were known, the radiance in the almucantar plane of the sun could be obtained by just fixing θ = θ_s and varying φ. The radiance in the principal plane could be extracted by fixing the (φ = φ_s) or (φ = φ_s + π) for any θ. In addition, the radiance at different isolines of the scattering angles could also be extracted from the images, using the relationship between the pixel scattering angle (Θ) and PZA and PAA for any given sun position following Equation (8) [47,48].

$$\cos (\Theta )=\cos \left({\theta }_s\right)\cos \left(\theta \right)+sin\left({\theta }_s\right)\sin \left(\theta \right)\cos \left(\left|{\phi }_s-\phi \right|\right)$$

(8)

Figure 7 shows an example of the representation of these three sky geometries in a SONA image. Due to the viewing geometry described in Section 3, the almucantar is represented by a radial circle at a fixed θ = θ_s, while the principal plane is a straight line that passes through the optical center and through sun position at fixed φ = ±φ_s.

4.2. Broadband Diffuse Irradiance on a Horizontal Plane

Our sky-camera did not use any shadowing element to blind the sun. Therefore, the circumsolar area in most of the images remained saturated in some or all RGB channels, even when using the shortest available exposure time. The size of this saturated region of the images was highly dependent on the atmospheric conditions, with it being smaller in clear sky situations and then increasing with the amount and size of atmospheric scatters (aerosols and clouds). Consequently, only the diffuse radiation could be obtained using our sky-camera.

The broadband diffuse horizontal irradiance of each channel (DHI_ch) could be obtained from the spectral and angular integration of the product between the radiance projected from every viewing direction to the horizontal plane L_ch(θ, φ)·cosθ and the suitable spectral response of the sensor S_ch(λ), following Equation (9).

$$DH{I}_{ch}={{\displaystyle \int }}_0^{\frac{\pi }{2}}{{\displaystyle \int }}_0^{2\pi }{{\displaystyle \int }}_{{\lambda }_1}^{{\lambda }_2}P(\Theta ) S_{ch}\left(\lambda \right) L_{ch}\left(\theta ,\phi \right)\cos \theta sin\theta d\theta d\phi d\lambda$$

(9)

$$P(\Theta )=\{\begin{array}{c}0, \Theta \le 5°\\ \\ 1, \Theta >5°\end{array}$$

(10)

In Equation (9), ch referred to R, G and B channels, and λ₁ and λ₂ were the spectral integration limits of each channel. P(Θ) was the logical NOT of the penumbra function [49,50] reformulated in terms of the scattering angle (Equation (10)) to subtract the circumsolar area. We chose a limit angle of 5 degrees to mimic the FOV of the shadow ball used in the measurement of the horizontal diffuse irradiance by the pyranometer deployed in the solar tracker.

The numerical integration in Equation (9) was calculated by uncoupling and discretizing the spectral and angular terms Equations (11) and (12), and by considering the calibrated irradiance E_ch(i, j), PZA (i, j) and scattering angle (Θ_ij) matrices at pixel positions (i, j), where DHI_ch is the sum of all pixels of an irradiance calibrated sky-camera image with a scattering angle of more than 5° and a zenith angle below 90° as is shown in Figure 8.

$$DH{I}_{ch}={\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{Q}_{ch} P\left({\Theta }_{i,j}\right) E_{ch}\left(i,j\right) \cos \left({\theta }_{i,j}\right)$$

(11)

$$Q_{ch}={\displaystyle \sum }_{k=1}^{M-1}{S}_{ch}\left(k\right)·\left({\lambda }_{k+1}-{\lambda }_k\right)$$

(12)

$$DH{I}_{SONA}=DH{I}_R+DH{I}_G+DH{I}_B$$

(13)

Q_ch was the convolution of the spectral response S_ch(k) of each channel with 1-nm wavelength step, where k = [1, M − 1] corresponded to λ = [300, 1100] nm. The broadband diffuse horizontal radiation (DHI) was obtained by adding up the DHI_ch of each sky-camera channel following Equation (13).

In practical terms, P(Θ_ij) was a mask matrix, where pixels with scattering angle below 5° or with zenith angle above 90° were zero and the rest were one (Figure 8a). Multiplying this mask by the irradiance matrices (Figure 8b–d) and summing up the resulting product of all pixels, we finally obtained the DHI_ch and then the sky-camera DHI (DHI_SONA).

4.3. Broadband Diffuse Irradiance on Arbitrary Planes

A similar approach to that described in Section 4.2 could be used to determine the broadband diffuse tilted irradiance (DTI) on an arbitrary plane DTI_β,ϕs, with tilt angle (β) and orientation angle (ϕ_s), measured from the South (positive to the West) (Figure 9). In practice, it was reduced to construct a mask that only considered the light directions reaching the tilted plane, obtaining the equivalent pixel zenith angle (PZA’; θ’) and appropriately changing the limits of the angular integration in Equation (9).

In our PAA matrix, the North direction corresponded with (φ = 0°) and φ increased counterclockwise, so by naming the PAA matrix on the rotated plane as PAA’, φ’ is therefore the difference of the φ and ϕ_s angles (Equation (14)).

The PZA’ on the tilted plane depended on the tilted beta angle and varied with φ’ (Figure 9). For pixels over the horizontal plane, θ’ was |θ – γ| and for pixels under the horizontal plane θ’ was |θ + γ|, as we can see in Equations (14)–(16).

$$\phi ’=\phi -{\varphi }_s$$

(14)

$$r=\sqrt{x^2+y^2}; y’=r·\cos \left(\phi ’\right); \gamma =\arcsin \left(\frac{y’·\sin \left(\beta \right)}{r}\right)$$

(15)

$$\begin{gathered} if \phi ’ >90 \mathrm{and} \phi ’ <270 \to \tau =1, \mathrm{else} \to \tau =0 \\ \theta ’=\left|\theta +{\left(-1\right)}^{\tau }· \gamma \right| \end{gathered}$$

(16)

To obtain the DTI_β,ϕs, we projected the normal direction by multiplying by cos(θ’) and applying the same spectral correction as in the horizontal plane. Note that the G(θ’_{i, j}’) function limited the pixels that we had to consider to calculate the DTI_β,ϕs.

$$DT{I}_{\beta ,{\varphi }_s}={\displaystyle \sum }_{i=1}^N{\displaystyle \sum }_{j=1}^N{Q}_{ch} P\left({\Theta }_{i,j}\right) G\left({\theta ’}_{i,j}\right) E_{R,G,Bch}\left(i,j\right) \cos \left({\theta ’}_{i,j}\right)$$

(17)

$$G(\theta ’)=\{\begin{array}{c}1, \theta ’\le 90°\\ \\ 0, \theta ’>90°\end{array}$$

(18)

Figure 10 shows an example of a PZA’ matrix obtained with the method described above for a tilted plane with β = 30° and ϕ_s = 0° (South oriented).

Note that the matrix PZA’ represented the zenith viewing angles of the sensor pixels on a tilted plane. In Figure 10, the zenith angles between 90° and 70° in the North direction are shaded by the plane itself. In contrast, the South direction is not shaded. In the tilted plane, the pixel viewing angles depended on the γ angle (derived from the calculation of the Euler angles for the rotation of the coordinate axes). The heart shape effect obtained depended on the tilt angle of the plane (the steeper the plane, the greater the effect).

5. Results

Different sets of RGB sky radiance in the almucantar plane, DHI irradiance and DTI_β,ϕ were extracted from the SONA HDR images. These sets were used for testing the reliability of our results in comparison to independent and well-calibrated measurements.

5.1. Validation of Sky Radiance

The validation of the RGB sky radiance obtained from HDR imagery was carried out from February to October 2020 by comparing the sky-camera radiance results to the AERONET almucantar radiance dataset at 440, 500 and 675 nm (the closest wavelengths to SONA λ_eff). To ensure that both instruments observed the same portion of the sky, the SONA almucantar radiance was extracted from the images averaging within the FOV of the Cimel (1.5°). Note that we did not apply any spectral band shift or correction to account for the difference between the Cimel (10 nm spectral band at 440, 500 and 675 nm) and our camera (~300 nm spectral band at 480, 541 and 618).

For scattering angles below 10°, we observed a camera sensor saturation, different for each of the RGB channels, which increased with AOD and the presence of some types of clouds surrounding the sun. Furthermore, in clear sky conditions with an AOD of aerosols of below 0.1, we had no saturated pixels, but a glare around the sun appeared, perhaps due to undesirable optical issues. To avoid these issues, measurements with a scattering angle of below 10° or with a dirty dome were rejected for the validation.

This validation was carried out for all sky conditions, which were classified as clear, partially cloudy and overcast skies and attended to the symmetry with respect to the principal plane and the average blue to red ratio (BRR) of each image. To determine how symmetrical the image was, we used the mean square deviation (MSD) between the left and right sides of the image with respect to the principal plane. Partially cloudy skies did not present any symmetry, while clear and overcast skies did. Thus, the BRR was used as a second criteria to distinguish between them. We empirically determined that HDR images with an MSD < 30 guaranteed symmetry, which corresponded to clear sky conditions if the average BRR > 2.8 as well.

Note that AERONET only provided almucantar radiance measurements in clear sky and partially cloudy conditions when the sun was not covered. This implied additional restrictions to our validation dataset.

We used the Relative Mean Bias Deviation (RMBD) to evaluate the differences between the sky-camera and the AERONET measurements, and the confidence level (CL) of the sky-camera measurements was determined by the percentage of cases with an RMBD of below 10% (CL₁₀).

A significant agreement between AERONET and SONA was observed in all cases, with correlation coefficients of more than 0.90 in most of them (

Table 4. Comparison between AERONET almucantar and SONA estimated radiances from 12 February to 31 October 2020. N represents the number of samples for each scenario and R is the correlation factor.

	Clear Sky				Partially Cloudy Sky				All Scenarios
λ	RMBD	CL10	N	R	RMBD	CL10	N	R	RMBD	CL10	N	R
(nm)	(%)	(%)	(#)	(#)	(%)	(%)	(#)	(#)	(%)	(%)	(#)	(#)
675	9	91	12,680	0.96	28	42	36,279	0.81	23	51	48,959	0.82
500	6	85	8720	0.98	10	75	23,116	0.92	8	79	31,836	0.93
440	4	91	12,320	0.98	10	76	35,513	0.90	8	81	47,833	0.91

). However, this agreement was different with respect to sky conditions and spectral dependence.

In clear sky conditions, a very good agreement was observed in 25% of cases. The results indicate a high degree of correlation (R > 0.96) for all channels, with an RMBD < 10%. However, better performance was observed in the blue, green and red channels with the RMBD increasing to 4%, 6% and 9% at 440, 500 and 675 nm, respectively. In addition, most of the cases showed an RMBD smaller than 10%, with high values of CL₁₀ (91%, 85% and 91% for 440, 500 and 675 nm, respectively).

Under cloudy sky conditions, 75% of cases also showed high correlation coefficients, but in this case with a larger RMBD compared to the clear sky conditions. Notable RMBD differences were observed among the red (RMBD = 28%) channel. And the blue and green channels, both showing an RMBD = 10%. In addition, low CL₁₀ values, especially on the red channel, indicated a larger spread of the data with respect to the clear sky conditions. We suspect that at least part of these differences was due to a disparity in the base time between the AERONET and SONA measurements. The whole AERONET almucantar measurement sequence took a few minutes (~5–15 min, depending on the SZA), having only one timestamp for the entire almucantar. Opposingly, we extracted the complete almucantar from a single HDR image of the sky-camera, which were stored every minute. Note that under changing sky conditions, any given position in the sky could change from clear to cloudy, and vice versa, in just a few seconds. This could be the origin of the large differences between the results, and may be especially important in partially cloud conditions.

The results obtained for all scenarios sum up the characteristics of clear sky and partially cloudy sky conditions. For this reason, it behaved as a weighted average due to the proportion of cases and was closer to the results showed in partially cloudy conditions, since they represented 75% of the cases.

Our results for clear sky conditions were similar to [43] and demonstrated an improved performance relative to [23,26]. However, to our knowledge, this was the first time that sky radiance was measured and validated in all atmospheric conditions using a sky-camera. Moreover, despite all issues of the comparison discussed above, our results show a good agreement with those of AERONET.

5.2. Determination of Broadband Horizontal Diffuse Irradiance

We obtained the DHI from our radiometric HDR imagery dataset, from February 2020 to March 2021, following the methodology explained in Section 4.2 and applying a 5° scattering angle mask to compute the DHI_SONA. Hereafter, in this section, when speaking of saturated pixels, we will refer to those that have not already been considered within the mask. Therefore, the DHI_SONA was determined following two different approaches in terms of how we treated the saturated pixels within the images. In the first approach, saturated pixels were not used to compute the DHI_SONA. In the second, saturated pixels were considered by using the maximum radiative value that the camera can measure. Therefore, we have two sets of DHI_SONA in all kinds of sky scenarios. Both sets of DHI results derived from the sky-camera were compared for validation against those acquired by a CMP21 pyranometer mounted on a sun tracker with a shadow ball (Figure 11). In addition, a third dataset, using only images with no saturated pixels, was used in the comparison (grey points in Figure 11).

The SONA and pyranometer DHI show a high degree of correlation in all cases, with R > 0.98 independently of the way that saturated pixels were considered to compute DHI_SONA (Figure 11a,b). In general terms, all fittings produced a very similar slope of around two (

Table 5. Linear fit parameters obtained from the validation of the broadband diffuse irradiance.

Saturated Pixels for DHI Calculation	Sky Conditions	Slope	Offset	R	RMSE	N
Not considered	All	2.10	−20.60	0.98	21.40	287,398
Considered as max Irradiance	All	2.04	−19.61	0.98	19.89	28,7398
Not present	All	1.97	−7.28	1.00	6.84	86,108

). This slope largely accounts for the different spectral range covered by the instruments and their different spectral response. Our results for the slope are as expected, since the SONA sky-camera is a three broadband channel sensor covering the 300−1100 nm spectral range (Figure 5), which represents approximatively half of the quasi-flat spectral response of the pyranometer in the 340−2800 nm range. This fact was also supported by the simple test of calculating the ratio between the pyranometer spectral response and the sum of the spectral responses of the sky-camera channels. A ratio of 2.06 was obtained, when the spectral responses of the devices were convolved with the sun irradiance model ASTM G173-03 and an air mass of 1.5 (not shown).

No substantial differences in the fit parameters were observed regarding the way saturated pixels were used to obtain the DHI_SONA. This is mainly due to the fact that, from a statistical point of view, the images containing saturated pixels do not represent a large proportion of the number of considered images. However, the best fit is obtained using data from the images with no saturated pixels, since the offset and RMSE were significantly reduced with respect to the other two cases (Table 5). These images correspond to a pyranometer DHI of lower than 500 Wm⁻² (grey points in Figure 11).

In spite of that, DHI_SONA presents a notable dependence on the percentage of saturated pixels, especially at higher values of DHI (>400 Wm⁻²) measured by the pyranometer (Figure 11a,b). In fact, an increasing underestimation of DHI_SONA with respect to the pyranometer values was observed as the fraction of saturated pixels increased. As expected, this is especially relevant when the saturated pixels are not used in the computations (Figure 11a,c), since a substantial part of the sky radiance may not have been considered in the retrieval. This underestimation effect is more evident for images with percentages of saturated pixels of larger than 3.5%, for which most of these DHI results fall out of the 95% confidence interval for the predictions (±50 Wm⁻²) (Figure 11c).

Conversely, the dependence on the percentage of saturated pixels was significantly reduced when the DHI_SONA was obtained following our second approach (Figure 11b,d). In this case, the 95% confidence interval was reduced to ± 40 Wm⁻². In spite of the inherent uncertainty of the radiative value of the saturated pixels, they had a significant effect in reducing the DHI_SONA, even for low saturated pixel percentages.

Usually, a large occurrence of saturated pixels is associated with the presence of a high aerosol load that produces a strong scattering of sun light in the circumsolar area or under partially cloudy sky conditions, in which very reflective clouds do not cover the sun. Both atmospheric conditions produced an increase in DHI that was not captured well by the SONA.

The best results of the comparison were obtained when the images were free of saturated pixels, even though the data used in this model did not cover all possible DHI values. It may represent an issue in the attempt to obtain the DHI_SONA in all sky conditions, since the underestimation effect of the saturated pixels appears for higher DHI values. Therefore, for operational purposes, we chose the model that accounted the saturated pixels as the maximum value. It allowed the SONA to determine the DHI in a more reliable way by accounting for the missing spectral region and regarding the saturated pixels, with less dependency on the sky conditions (not shown in this paper).

Figure 12 shows the predictions of the DHI_SONA by applying the spectral response correction factor and the offset provided by the chosen fit model, compared to the pyranometer DHI from 28 March to 27 July 2021 (not used in Figure 11). The comparison results indicate a high degree of agreement between data, with a slope of close to one and the offset relative to the maximum irradiance being well below the pyranometer uncertainty (~5%). The maximum difference between the pyranometer and DHI_SONA is 2% considering all the sky conditions in the validation. In addition, 95% of the images (±2 standard deviation) show DHI values within the ±50 Wm⁻². These results suggest an optimal capability of our methodology to determine the DHI from a well-calibrated sky-camera.

Our validation results are satisfactory, since they improve those reported by [19] that showed lower correlations in comparison to the pyranometer measurements using a sky-camera and a different methodological approach.

5.3. Validation of Broadband Diffuse on an Arbitrarily Oriented Plane

The broadband diffuse radiation on an arbitrary oriented plane DTI_β,ϕ can be determined from sky-cameras. In our case, we used β = 30° and ϕ = 0° to determine DTI_30,0 from SONA since the validation was carried out against the measurements of a CM11 pyranometer set up on a 30° tilted and South-oriented plane, from 4 August to 5 September 2021. As discussed in Section 5.2, the DHI_SONA was corrected by a spectral response range factor to match the pyranometer values, giving saturated pixels the maximum possible measured radiative value.

The CM11 pyranometer also measured the global irradiance on the selected tilted plane GTI_β,ϕ. However, to be able to carry out the comparison we determined the DTI_β,ϕ from the CM11 by subtracting the projection of the DNI, measured by a pyrheliometer on a solar tracker, from the beta tilted plane DNI·cosω following Equations (19) and (20).

$$DT{I}_{\beta ,\varphi }=GT{I}_{\beta ,\varphi }-DNI·\cos \omega$$

(19)

$$\cos \omega =\cos {\theta }_s\cos \beta +\sin {\theta }_s\sin \beta \cos \left(\pi -{\phi }_s-\varphi \right)$$

(20)

The results of this comparison are shown in Figure 13. The SONA retrieval of DTI_β,ϕ offered an extremely high correlation with the independent pyranometer measurements (R = 0.99). The slope deviated 5% from the unity, indicating a slight overestimation of the SONA DTI_β,ϕ with respect to the pyranometer. Moreover, 95% of the data used in the comparison falls within a ± 50 W/m² interval around the fit line. No particular differences were observed regarding the variability of the percentage of the saturated pixels.

After the necessary processing and corrections, the sky-camera may be a useful instrument for measuring the solar radiation in different geometries.

6. Conclusions

In this work, we have presented a methodological approach that is able to provide accurate and calibrated measurements of sky radiance and broadband solar irradiance using HDR images from a sky-camera.

Firstly, we performed a very detailed characterization of the SONA, in terms of electronics and image processing, and a significant improvement of the dynamic range of the images, conferring on them an extremely linear behavior under different illumination scenarios. In addition, we developed a methodology to automatically perform the geometric calibration of the camera. In that way, an optical characterization of the camera was carried out, obtaining the center of the image, the focal length of the fisheye lens and consequently, the observation geometry (in polar coordinates of the sky-dome). Then, the HDR images were calibrated using a radiative transfer model to obtain the sky radiance of each pixel, with an uncertainty smaller than 10%. Finally, we constructed a 1-min temporal resolution calibrated HDR imagery database, which was available from February 2020 and updated once a day.

The resulting sky radiance from these images was compared to the AERONET spectral radiances in the almucantar plane in different sun geometries and sky conditions. The results of this validation presented a very good agreement between our estimations and the AERONET measurements in all sky conditions. This agreement is especially relevant in clear sky cases, showing correlation coefficients higher than 0.96 and an RMSD < 10% for all spectral channels.

In addition, we developed methods to determine the shortwave diffuse radiation over arbitrary geometries and planes, starting from calibrated sky radiance images obtained from our sky-camera. In this sense, we obtained the diffuse radiation over a horizontal plane DHI_SONA. It showed a very good correlation with the independent pyranometer measurements in all sky conditions. The results of this DHI_SONA and pyranometer comparison were used to correct DHI_SONA retrievals due two facts: the limited spectral range of SONA sky-camera that does not cover the whole solar spectrum and how the saturated pixels should be treated in the retrieval of DHI_SONA. When these corrections were applied, the difference between pyranometer and our DHI_SONA retrieval was 2% and, in 95% of the images (±2 standard deviations), the DHI values were within the ±50 Wm⁻².

Finally, we tested our methodology to estimate the diffuse irradiance over an arbitrarily tilted and oriented plane from SONA by comparing it to measurements collected by a pyranometer tilted 30° with a South orientation. Our results showed an extremely high correlation between both estimates (R = 0.99) in all sky conditions, which also indicates a slight overestimation of 5% of the SONA DTI_β,ϕ with respect to the pyranometer.

All results suggest an optimal capability of our methodologies to determine the sky radiance and the diffuse irradiance over arbitrary geometries from a well-calibrated sky-camera. In turn, they represent a significant contribution to the improvement of the capabilities of the sky-cameras to provide accurate and well-calibrated radiometric images of the whole sky, making them suitable for solar and atmospheric research. In addition, the algorithms and retrievals developed in this work may be of great interest for both solar energy producers and energy network operators, allowing them to monitor and estimate solar radiation incidents on an arbitrary plane with a single and affordable instrumentation on a steady setup and in all sky conditions.