Retrieval of Solar Shortwave Irradiance from All-Sky Camera Images

Abstract

The present work proposes a new model based on a convolutional neural network (CNN) to retrieve solar shortwave (SW) irradiance via the estimation of the cloud modification factor (CMF) from daytime sky images captured by all-sky cameras; this model is named CNN-CMF. To this end, a total of 237,669 sky images paired with SW irradiance measurements obtained by using pyranometers were selected at the following three sites: Valladolid and Izaña, Spain, and Lindenberg, Germany. This dataset was randomly split into training and testing sets, with the latter excluded from the training model in order to validate it using the same locations. Subsequently, the test dataset was compared with the corresponding SW irradiance measurements obtained by the pyranometers in scatter density plots. The linear fit shows a high determination coefficient ($R^2$) of $0.99$. Statistical analyses based on the mean bias error (MBE) values and the standard deviation (SD) of the SW irradiance differences yield results close to $-2\%$ and $9\%$, respectively. The MBE indicates a slight underestimation of the CNN-CMF model compared to the measurement values. After its validation, model performance was evaluated at the Antarctic station of Marambio (Argentina), a location not used in the training process. A similar comparison between the model-predicted SW irradiance and pyranometer measurements yielded $R^2=0.95$, with an MBE of around $2\%$ and an SD of approximately $26\%$. Although the precision provided by the SD at the Marambio station is lower, the MBE shows that the model’s accuracy is similar to previous results but with a slight overestimation of the SW irradiance. Finally, the determination coefficient improved to $0.99$, and the MBE and SD are about $3\%$ and $11\%$, respectively, when the CNN-CMF model is used to estimate daily SW irradiation values.

1. Introduction

Total solar shortwave (SW) irradiance is the radiative power emitted by the Sun per unit surface integrated over the solar wavelength spectrum [280–3500 nm]. The SW irradiance that reaches the Earth’s surface is usually defined over a horizontal plane at ground level and is called global horizontal irradiance (GHI). Then, the incoming GHI is modulated by the relative position of the Sun and the observer, with the main factor being the solar zenith angle (SZA). SW irradiance is attenuated by different interactions with atmospheric components, such as absorption or scattering processes involving gas molecules, aerosols, and cloud water droplets [1,2]. The Earth also emits irradiance but in the thermal infrared range [4–40 $\mathsf{\mu }$m] (longwave radiation; LW). Hence, the Earth’s radiative budget, corresponding to the difference between SW and LW net irradiances, determines the temperature of the planet and its atmosphere since the equilibrium of this budget represents the fundamental driver of the climate system. Therefore, incoming radiation must be balanced with outgoing radiation to achieve this equilibrium [2,3].

The monitoring of GHI in polar regions is crucial since it is well known that these areas play a critical role in the global energy balance through the maintenance of the north–south temperature gradients that drive Earth’s oceanic and atmospheric circulations [4]. These areas are highly sensitive to changes in climate forcing, and warming of more than twice the global average has been reported in the northern regions; in contrast, the southern latitudes have experienced slower warming or even cooling in some areas (see [5] and references therein). The scientific understanding in these regions is particularly challenging due to the complexity of the underlying processes of changes in their climates [5], highlighting the need for further research on this topic. Given recent developments, solar energy technology has emerged as one of the leading renewable energy sources currently available. Knowledge of solar energy enabled by GHI reaching the Earth’s surface is essential for research on this topic as input for short-term irradiance forecasting and for other applications in the field of photovoltaics [6].

GHI measurements are typically obtained using pyranometers; however, in locations where there is a lack of available data, GHI data series are usually reconstructed from models. Radiative transfer models (RTMs) are frequently used to estimate GHI under cloud-free conditions, an estimate that can also be obtained by using other kinds of empirical models. For these cloud-free conditions, the results retrieved by RTMs are generally very accurate. However, difficulties appear when cloud effects are considered. One way to take cloud effects into account is by using the cloud modification factor (CMF), which quantifies cloud transmittance as the ratio between the GHI that reaches the surface under cloudy and cloud-free conditions [7]. GHI can be accurately estimated as the product of the CMF and the simulated cloud-free GHI if the CMF is known. Some works have used auxiliary measurements, such as sunshine duration measured by pyrheliometers, as a proxy for clouds to calculate the CMF and reconstruct long-term GHI series. Thanks to these long-term series of measured and reconstructed GHI, global phenomena, such as global dimming and brightening [8,9,10,11,12,13], have been proven to occur in the 20th century.

Sky images from all-sky cameras can be used as a proxy for cloud effects, as these images capture the position and brightness of clouds across the entire sky. Furthermore, other works have shown that these sky images contain information about other atmospheric components, such as aerosols [14,15]. Recently, different image treatment techniques have emerged, such as convolutional neural networks (CNNs), which are able to extract different characteristics from images, such as, for instance, cloud cover (CC) from all-sky camera images [16], among others. This framework motivates the present work, the main objective of which is to develop a model that applies these new CNN techniques to all-sky camera images to estimate the CMF and reconstruct the GHI values. In addition, the accuracy and precision of the proposed model are studied by comparing the computed GHI and CMF predictions with real measurements obtained at a remote location.

The present work is structured in the following way: Section 2 shows the instrumentation used and the sites where they are installed. The dataset classification and the model to retrieve the GHI, with the training process and validation of the test set, are described in Section 3. Section 4 contains the analysis of the results at the Antarctic station of Marambio. At the end, the main conclusions are summarized in Section 5.

2. Instrumentation and Sites

2.1. Instruments

Two kinds of instruments were used in this work to collect the used data: all-sky cameras and pyranometers. The all-sky camera model used in this work is the OMEA-3C (Alcor System, Lyon, France). This camera model is composed of a SONY IMX178 CMOS sensor, which works with a Bayer RGGB mosaic with an infrared filter, coupled to a fisheye lens with a 180° × 180° field of view covering all the sky. It also incorporates external sensors to measure humidity and temperature, an internal heating system to avoid condensation, and a BK7 glass dome on top [17,18]. The sensor takes pictures of 6.44 megapixels at a 14-bit resolution with a size of 3096 × 2080 pixels.

This camera model was available in the four stations described in Section 2.2. These cameras were controlled and configured by the GOA-OMEA Capture-1.0 software, which was developed by the Group of Atmospheric Optics of the University of Valladolid (GOA-UVa). These all-sky cameras were configured to take (every 5 min during the daytime) a sequence of consecutive raw pictures with different exposure times. The images of each sequence are combined to obtain a high dynamic range (HDR) image [19]. These HDR pictures are tone mapped and converted into images of size 2000 × 2000 pixels at an 8-bit resolution [20].

Regarding GHI, this was recorded every 1 min at the mentioned stations using four different types of pyranometers: CMP-11, CM-21, and CM-22 from Kipp & Zonen (Delf, The Netherlands) and MS-802F from EKO Instruments (Tokyo, Japan). The CMP-11 model has a response time of less than 5 s in a spectral range of 285–2800 nm; the achievable uncertainty (confidence level of $95\%$) for the daily totals of this instrument is expected to be lower than $2\%$ [21]. CM-21 has an uncertainty (with a $95\%$ confidence level) of $2\%$ for daily totals, with a time response of 5 s in a spectral range of 305–2800 nm [22]. CM-22, with a spectral range of 200–3600 nm, has a time response of 5 s at a $95\%$ confidence level and an achievable uncertainty of $2\%$ for daily totals [23]. Finally, MS-802F is a high-precision instrument with a spectral range of 285–3000 nm, with less than 5 s of time response (confidence level of $95\%$) and an expected uncertainty of lower than $\pm 1\%$ for daily totals [24,25].

2.2. Sites

The all-sky cameras and pyranometers that were used, as described in the previous section, were installed at four different locations: Valladolid (Spain), Izaña (Spain), Lindenberg (Germany), and Marambio (Argentinian Antarctic Station). All these all-sky cameras are managed by the GOA-UVa, and the pyranometers are directly managed at each station.

The Valladolid station is located on the rooftop of the Science Faculty of the University of Valladolid (41.66°N, 4.71°W, 705 m asl), which is a medium-sized city located in north-central Spain. It is an urban city with approximately 400,000 inhabitants in all of the metropolitan area, surrounded by rural areas. Its climate is classified as clean continental, with hot summers, cold winters, and occasional Saharan dust episodes [26,27,28]. The GOA-UVa manages the instrumentation platform installed at this station, which is composed of different instruments, such as sunphotometers, a ceilometer, all-sky cameras, and various radiometers [14,16,29,30]. The OMEA-3C camera used in this work has been operating in Valladolid since 5 October 2021, and the pyranometer is a CM-21 model installed on a SOLYS-2 solar tracker (Kipp & Zonen).

On 4 February 2022, an OMEA-3C camera from GOA-UVa was installed at The Izaña Atmospheric Research Center (28.30°N, 16.49°W, 2400 m asl), which is managed by the Meteorological State Agency of Spain (AEMet). This meteorological observatory is a high-mountain station above a quasi-permanent strong temperature inversion layer that prevents the arrival of local anthropogenic influence with stable and pristine conditions, located near the Teide peak in Tenerife (Canary Islands, Spain) [24,31,32,33,34]. This observatory has been enrolled in the Baseline Surface Radiation Network (BSRN) since 2009, implemented by the World Climate Research Program (WCRP) under the support of the World Meteorological Organization (WMO). The GHI is one of the basic BSRN measurements of the BSRN program at Izaña. The BSRN GHI measurements at Izaña used in this work (available at https://doi.pangaea.de/10.1594/PANGAEA.968676; URL access on 25 August 2024) were recorded by an MS-802F pyranometer, which has been working since 11 November 2016 at this station [24].

MOL-RAO (Meteorologisches Observatorium Lindenberg—Richard-Aßmann-Obser-vatorium; 52.21°N, 14.12°E, and 122 m asl) is located in the village of Lindenberg (Tauche) in the north-east German region of Brandenburg (green flat land and partly woody). Lindenberg is located about 50 km southeast of the outskirts of the city of Berlin, far from the coast, with predominantly continental influence. The observatory belongs to the Department of Research and Development of the DWD (Deutscher Wetterdienst—German Weather Service). It is a historical station (since 1905) for meteorological and atmospheric measurements; today, it is a supersite for meteorological measurements and measurements of atmospheric parameters and the physical processes of the atmosphere, from the boundary layer to the stratosphere. MOL-RAO hosts the GRUAN leading center, which is a WMO testbed for atmospheric measurements and has been a BSRN station since 1994. GHI measurements have been obtained from a CM-22 pyranometer at this station. Two OMEA-3C all-sky cameras from GOA-UVa were installed at this station. One camera has been operating since 29 July 2020 until 12 February 2021, when it was substituted by another OMEA-3C, which is still working.

Marambio is an Argentinian Antarctic station located in the ice-free Marambio island (64.24°S, 56.62°W, and 200 m asl) at the north-east tip of the Antarctic Peninsula [35]. In January 2018, the GOA-UVa installed an OMEA-3C all-sky camera at this site, the images of which have been used since August 2020, when the GOA-OMEA Capture-1.0 software started to be used to record sky images. This all-sky camera is installed on the rooftop of the Scientific Laboratory, which is located in the middle of the station between the airfield and station buildings. This laboratory belongs to the National Meteorological Service of Argentina (SMN) and is shared with the Argentinian National Direction of the Antarctic (DNA) [35]. The station is also equipped with instrumentation focused on atmospheric studies using different techniques, such as a CMP-11 pyranometer, the GHI records of which were used for this work.

3. Method

A neural network model to retrieve GHI via CMF estimation from sky images is described in this section. The datasets used for training, validation, and testing are detailed as well. Data from Valladolid, Izaña, and Lindenberg were used for model training and validation, whereas the Marambio data were only used to test model performance in a different location than those used for training.

3.1. Dataset

We obtained the dataset retrieved from the pyranometers described in Section 2, containing measurements of GHI every minute. Additionally, we have daytime images captured by the aforementioned all-sky cameras every 5 min. The picture datasets contain all-sky HDR images taken using the different cameras operating at the different stations. These images were captured every 5 min, combining a sequence of various raw images captured consecutively at different exposure times, usually doubling the time exposure; this helped obtain the final HDR image. All the available raw images of a sequence were normalized to the same lowest exposure time and averaged (discarding saturated pixels). The averaged image was tone mapped by using a root square function, and they were white balance-corrected using the same fixed values for all images. After that, the signal of each pixel was normalized by dividing each pixel signal by the 99.5 percentile of all the pixel signals. This image was multiplied by 255 and rounded to be converted to 8-bit when it was finally demosaiced. This process is explained in detail in [20]. The final HDR images are 2000 × 2000 pixels at an 8-bit resolution. In Figure 1, some examples of these sky HDR images are shown for different sky conditions and for all cameras and stations. These datasets have been filtered and classified, as explained in the next sections.

3.1.1. Initial Dataset and Data Filtering

In order to classify the data for model training and validation, the first step consisted of matching the GHI measurements with the all-sky camera images. To conduct this matching, we used two different datasets: the data file with the pyranometer measurements and a list of the available sky images, both for each station. These two datasets contained a column with the timestamps when the measurements were taken, so we selected the GHI and SZA values where the timestamp coincides for both datasets. For this purpose, the selected periods of the coincident measurements of the GHI and sky images are the following: from 1 January 2022 to 23 August 2023, from 4 February 2022 to 31 January 2023, and from 30 July 2020 to 12 September 2023 for Valladolid, Izaña, and Lindenberg, respectively. The number of images in the initial dataset corresponding to these periods at each station is shown in the second column of

Table 1. Number of sky images classified to train and test the model for different sites.

Site	Initial Dataset	Filtered Dataset	Train Dataset	Validation Dataset	Test Dataset
Valladolid	82,865	70,962	25,654	7096	38,212
Izaña	46,509	41,327	6104	4133	31,090
Lindenberg	146,575	125,380	54,679	12,539	58,072
Total	275,949	237,669	86,527	23,768	127,374

.

All the data were filtered to avoid very low GHI values because these measurements represent noisy results close to 0; this introduces high uncertainty in the CMF value since it will be the ratio between two values, with both being close to zero, especially for high SZA values. For that, the solar zenith angle (SZA) was limited to be below 85° due to the very low GHI measurements between 85° and 90°, even under cloud-free conditions, and the GHI was above 5 Wm⁻². The sky images were captured quasi-instantaneously, while the pyranometer values were recorded over a longer integration time (1 min). As a result, inconsistencies may arise between the GHI reaching the camera and the pyranometer, especially when the sky conditions are highly variable. For example, the Sun may be obstructed by a cloud in an instantaneous image but remain unobstructed for most of the integration time during the pyranometer recordings. Therefore, the matched camera–pyranometer data pairs were discarded if the difference between the recorded GHI and the corresponding GHI measurements 1 min before and after were higher than $30\%$. In summary, it is remarkable that this filtering process avoids possible issues regarding GHI measurements (as outliers or noise values) and sky images, as overexposure could occur. After this data filtering, the CMF calculations and the image classifications follow in the next section.

3.1.2. CMF Calculation and Data Classification

Generally, a tone map [36,37] is applied to all-sky camera images to display all elements present in the sky with sufficient brightness and clarity. Unfortunately, it removes the proportionality between the pixel signals of two images and the incoming SW radiation to both images. On the other hand, sky images do provide information about the quantity and type of clouds present, their position in the sky, and their brightness, as well as whether they obstruct the Sun or not; in other words, the effect of clouds can be inferred from the images. These facts suggest that it may not be very accurate to directly adjust the sky images to GHI. Instead, they were adjusted to an intrinsic value of the clouds, their transmittance, which is commonly parameterized by the cloud modification factor. In the case of GHI, the CMF_GHI was calculated using Equation (1):

$${\mathrm{CMF}}_{\mathit{GHI}}={\displaystyle \frac{{\mathrm{GHI}}_{\mathit{meas}}}{{\mathrm{GHI}}_{\mathit{cf}}}}$$

(1)

The sub-indices “meas” and “cf” refer to the measured GHI under cloudy conditions (GHI_meas) and the same under cloud-free conditions (GHI_cf), respectively. GHI, under cloud-free conditions, can be estimated with a radiative transfer model; hence, if CMF_GHI is known, the real GHI can be calculated by multiplying that value by those simulated under cloud-free conditions. Therefore, this work looks to derive CMF_GHI directly from the sky images since GHI can be obtained from these values.

Then, once the datasets were matched and filtered, the measured GHI values were converted into CMF_GHI by dividing them by the cloud-free GHI simulations [27]. These cloud-free GHI values were obtained using the libRadtran-2.0.5 package. LibRadtran is a library of radiative transfer routines and programs, the main program of which is the radiative transfer tool uvspec, which allows for calculating irradiance after specifying the sky conditions in the input files [38,39]. These GHI values were simulated using the radiative transfer equation solver twostr (two-stream) [40]. The input values were the following: the extraterrestrial spectrum from Kurucz [41]; the date and coordinates for each station; the monthly average ozone value for each location, calculated using the daily average of ozone total column of GIOVANNI provided by OMI-DOAS (Ozone Monitoring System Differential Optical Absorption Spectroscopy; https://giovanni.gsfc.nasa.gov/giovanni; URL access on 25 August 2024); the spectral surface albedo, obtained from the climatology of MODIS bi-directional reflectance distribution function (BRDF) and albedo [42]; the total column water vapor, obtained from AERONET’s daily climatology [43]; and the aerosol optical depth by using the Angstrom exponent and turbidity obtained from AERONET’s daily climatology [43]; the single scattering albedo and the asymmetry parameter were set to a constant value given by AERONET’s monthly climatology. In addition, the spectral resolution was set by using 2 nm steps from 280 to 800 nm, by 5 nm steps until 1600 nm, and by 10 nm steps from there to 2800 nm. The SBDART model [44] was used for gas absorption parametrization. Finally, the irradiance values were interpolated every 1 nm from 280 nm to 2800 nm, and the obtained spectrum was integrated to obtain GHI in a similar way to [45]. According to the results obtained in [45], the comparison of GHI measurements and simulations under clear conditions shows an MBE value of 0.2% and an SD of 6.3%. This agreement between measurements and simulations allows us to use libRadtran’s GHI simulations for cloud-free conditions with high confidence.

After calculating the corresponding CMF_GHI values, the dataset was filtered again to remove data pairs with CMF_GHI values above 1.3 because they are unrealistic and could indicate some problems with the CMF_GHI measurements, such as snow covering the pyranometers, among others. The CMF_GHI values above 1 are due to the enhancement effect, where the radiation measured is greater than the cloud-free value [46,47]. The enhancement cases, which present CMF_GHI values above 1.3 are very unusual, being between 1.3% [47] and 3% [46] of the total of the enhancement effect cases. For this reason, it has been decided that a threshold of 1.3 for CMF_GHI should be used so as to not misunderstand the results with outliers, as the values above this threshold are very unusual according to the literature. Then, the sky images of the final datasets were classified in different folders, grouping them according to 0.01 intervals of their corresponding CMF_GHI for each site, having a total of 130 different classes of CMF_GHI. As a result of the filtering step and grouping according to CMF_GHI, a total of 237,669 images were classified: 70,962 corresponding to Valladolid, 41,327 to Izaña, and 125,380 to Lindenberg, as shown in the third column of Table 1. These images were divided into three sets: training, validation, and test. The model, training, and validation procedures are described in Section 3.2, and the test of its performance is in Section 3.3.

The Marambio measurements were reserved for the application of the model as a case study and were excluded from the training process. This case is detailed in Section 4, where the results are presented. In Marambio, the selected period of data ranges from 16 November 2021 to 31 October 2023. The data before this initial date were discarded due to technical issues detected in the images and the pyranometer. After all, a total of 59,444 all-sky images were collected for application in the model at this place. These images were matched and filtered following the same procedure as at the other stations, obtaining a total of 57,852 camera–pyranometer data pairs from this station.

3.2. Model

In order to retrieve the GHI (through CMF estimation) from the mentioned daytime sky images, it was proposed that a model should be based on a convolutional neural network. The proposed model, named CNN-CMF, is an adaptation of an existing CNN model developed to retrieve cloud cover (CNN-CC; [16]). The CNN-CMF model architecture is exactly the same as that of CNN-CC (more details can be found in [16]), but it uses 130 CMF_GHI classes in the training instead of 9 cloud cover classes. This CNN-CMF model is designed to provide the CMF_GHI as a float number between 0.01 and 1.30 using a sky image as input.

The sky images from Valladolid, Izaña, and Lindenberg were divided into three datasets, as described in the previous section. The training and validation sets are part of the model training process, where the validation set corresponds to the sample of data used to provide an unbiased evaluation of the model fit on the training dataset after each epoch while tuning the model parameters. The test set is the sample of data used to provide an unbiased evaluation of the final model, which is detailed in the next section. The classification of images into these sets was carried out for each station as follows. First, the number of available images in each CMF_GHI class was calculated, and after that, the 25 percentile $\left(P_{25}\right)$ of these 130 numbers was calculated and considered as the maximum number of images in the training dataset for each CMF_GHI class. Then, for each CMF_GHI class, (1) $10\%$ of the images were randomly added to the test set; (2) another $10\%$ of the images were randomly selected for the validation set; (3) a total of $P_{25}$ images were randomly added into the training set; (4) the rest of available images were added to the test set. This warrants at least 10% of the images of each class for testing and validation. After that, the datasets of the three locations were joined.

As a result, a total of 86,527 pictures belonging to 130 classes formed the training set. From this set, 25,654 images correspond to Valladolid, 6104 to Izaña, and 54,769 to Lindenberg. The validation set is formed from 23,768 sky images distributed into the same 130 classes, of which 7096 are from Valladolid, 4,133 are from Izaña, and 12,539 are from Lindenberg. This gives us a total of 127,374 pictures for the test set, constituting 38,212 from Valladolid, 31,090 from Izaña, and 58,072 from Lindenberg. This data classification is summarized in the last three columns of Table 1.

Once the datasets were formed, the model was trained using the mean square error (mse) as the loss function, considering that the model must penalize more those cases where the predicted CMF_GHI values are farther from the labeled reference than the cases where they are closer even when they are not the same. The training also had to be carried out using batches, with a batch size of 32, as each one is formed by 32 images. The optimizer is the Adam algorithm, a similar method to the one used for the training of the CNN-CC model explained in [16]. The initial learning rate was the same as in the reference model, which is dynamically reduced by dividing it by 10 when the loss metric of the validation set did not improve over 5 consecutive epochs. The training stopped when the loss value of the validation set did not improve after 10 consecutive epochs, trying to avoid overfitting. The CNN-CMF model was trained, reaching a total of 34 epochs before stopping, where epoch 24 corresponds to the chosen version of the model.

3.3. Test Set

Once the CNN-CMF model had been trained, its performance was evaluated at the three corresponding stations. For this purpose, the model was run using the images of the reserved test dataset, providing an estimation of the CMF_GHI values, which were transformed into GHI values by multiplying their corresponding cloud-free GHI simulations. These predicted GHI values were compared against the measured GHI values with the pyranometers, as shown in the density scatter plots of Figure 2 for the three locations and all of them together. These graphs also present the linear fit line and some statistics of the comparison, such as the standard deviation (SD) and the mean bias error (MBE) of the GHI differences between CNN-CMF and the pyranometer values.

In general, the density scatter plots of Figure 2 show a high correlation between predicted and measured GHI values, with a determination coefficient ($R^2$) value of $0.99$ for the four cases. The results reveal that high values of GHI are more frequent in Izaña, while low values are more common in Lindenberg. This discrepancy is primarily due to the high frequency of cloud-free days in Izaña and quite cloudy days in Lindenberg, as well as lower Sun elevation in Lindenberg. Valladolid exhibits an intermediate behavior between the other two stations. Additionally, in the panel that includes data from all three stations, a higher frequency of low values of GHI is observed due to the greater number of available data from Lindenberg (58,072 out of a total of 127,374). Regarding the linear fit, Valladolid shows the slope value closest to 1 (1.00), followed by Izaña (0.96) and Lindenberg (0.95), which also shows the y-intercept closest to zero (−3.80 Wm⁻²). The CNN-CMF model in Valladolid also exhibits the lowest MBE (−1%) and SD (7%) values. However, it slightly underestimates the measured GHI in the three stations with MBE values between −1% (Valladolid) and −3% (Izaña). This better fitting in Valladolid may be attributed, at least in part, to the mentioned availability of a similar number of data for all GHI intervals. This is also noticed when comparing the absolute results, as the absolute value of the MBE is notably larger at Izaña (where the available data present bigger GHI values) than the others, but the relative results are more similar (although it is also larger). Finally, for the total case of the three stations together, the MBE is about −2% and the SD close to 9%, indicating a slight underestimation of the CNN-CMF model compared to the measurements.

4. Results

In order to study the performance of the developed CNN-CMF model in a location not used in training, this model was applied to the sky images recorded in the mentioned station of Marambio to estimate the GHI. A total of 59,444 images recorded since 16 November 2021 at this site were used for this task.

4.1. SW Irradiance

Firstly, for a qualitative comparison, the intra-daily evolution of the predicted GHI values was plotted together with the measured ones for four different days in Figure 3. It must be noted that the predicted values are shown every 5 min (camera capture frequency), whereas the measured ones are every 1 min. This means that data filtering has not yet been performed for this comparison. In general, the CNN-CMF predicted and measured values are in concordance with all the days shown. The first day (21 December 2021) corresponds to a day with almost cloud-free conditions. In fact, both data series match with the cloud-free simulations, except between 13:00 and 14:00 UTC, when some clouds appeared. This cloudy period presented enhancement effects (SW higher than under cloud-free conditions; [46,48,49,50]) but also subperiods with a significant attenuation of GHI by clouds; both the predicted and measured values are capable of reproducing these behaviors. Something similar happens for panel (b), where both the predicted and measured values detect similar enhancement effect conditions between 13:00 and 14:00 UTC but for a cloudy day in this case. Moreover, during this day, the predicted values fit with the measured ones, especially in the afternoon, when the GHI variation with time is not so abrupt. A similar behavior can be appreciated in panel (c). This is explained by the fact that, as mentioned above, rapid and abrupt changes in incoming GHI can provide differences between the values predicted from images (instantaneous) and the values measured by the pyranometer (higher integration time). Finally, panel (d) presents a day with the presence of cloud-free and cloudy periods with different levels of GHI attenuation. In this case (4 February 2022), the predicted values overestimate the measurements between 17:00 and 19:00 under low CMF_GHI conditions.

Once the daily evolution of the CNN-CMF predicted GHI values was observed, even for abrupt temporal changes in GHI values, the dataset was filtered, as described in Section 3.1, for a quantitative comparison against the pyranometer measurements. As a result, Figure 4 shows the density scatter plot between the CNN-CMF predicted GHI values and the measured ones. The density of points indicates a high occurrence of low GHI values related to cloudy conditions, similar to those observed at Lindenberg in Figure 2, but also with the prevalence of high SZA values in this area. The predicted values correlate with the measurements, with an $R^2$ value of about $0.95$. Some data pairs in Figure 4 present high dispersion, which provides a higher value of SD of about 26%, indicating the lower precision of the CNN-CMF model than that obtained in the previous comparison with Valladolid, Izaña, and Lindenberg. In contrast, the obtained MBE is about 2%, similar to the absolute value obtained in Figure 2d when the results for Valladolid, Izaña, and Lindenberg are combined. This points out that the CNN-CMF model is accurate, but it slightly overestimates the GHI measured at Marambio.

GHI is modulated by clouds and also by the SZA due to the importance of the inclination of the incoming solar radiation. In addition, the sky images usually present a radial distortion caused by fisheye lens projection, which can mean that the Sun looks different in the image depending on the SZA. Therefore, in order to observe the effect of SZA in the performance of the CNN-CMF model, the MBE and SD of the GHI differences ($\Delta$GHI) between the model predictions and the measurements have been calculated for SZA ±2.5° bins. The absolute and relative values of these MBE and SD data have been plotted as points and error bars, respectively, in Figure 5 as a function of SZA. The 25th and 75th percentiles are shown as boxes in the figure, and the median values are depicted as straight lines. These results have also been calculated by taking into account cloud cover, splitting the data into the following categories: cloud-free (0–1 oktas in red), partly cloudy (2–6 oktas in green), and totally covered (7–8 oktas in orange). The total results for all data are depicted in blue. The absolute MBE for all data generally ranges between 0 and 25 Wm⁻², pointing out an overestimation of 20 Wm⁻², except for SZA $>{70}^{\circ }$, which ranges between −10 and 0 Wm⁻², showing a bit of an underestimation of the predicted values. In this plot, the SD decreases as the SZA increases, ranging around ±100 Wm⁻² for lower-angle values, around ±50 Wm⁻² for medium SZA values, and around ±25 Wm⁻² for higher angles. The variation in the relative MBE with SZA does not show any clear pattern, with the MBE ranging between 0 and 5% for SZA below 72.5° and between $-10\%$ and $0\%$ for the rest. In contrast, the SD increases with SZA, ranging between 20% (for the shortest) to 30% for the largest SZA values. These results reveal a slight overestimation of the CNN-CMF model to the measurements for SZA values below 72.5° and with a precision of about 20%; however, for higher SZA values, this model underestimates the measurements with a precision close to 30%.

As can be noticed from Figure 5, the presence of clouds affects the results for both absolute and relative values. For cloud-free sky conditions, the MBE differences are closer to 0 for SZA’s between 45° and 60°. In the range from 65° to 70°, the results are quite similar for all-sky conditions, which are also around 0. For larger angles of between 75° and 80°, the cloud-free conditions present higher results, especially for the relative values. The absolute MBE is around −10 Wm⁻² and between −10% and −15% for relative values. For 40° of SZA, the absolute MBE of the cloud-free condition is larger than for the rest of the sky conditions (near −10 Wm⁻²); however, for the relative values, the results are also closer to the rest of the sky conditions (around −2%). In general, the cloud-free cases present more negative values than the rest of the cases. It is also noticeable that totally covered sky conditions generally show positive relative MBE values, except at 80° where all results are negative. The fully covered and partly cloudy conditions are quite similar to the overall results for MBE, both in relative and absolute terms, except at 45°, where the partly cloudy results are closer to the cloud-free values, and the fully covered values are slightly larger than the others; additionally, at 75°, where the fully covered results are higher than the others, and the partly cloudy values fall between the total and cloud-free values. The SD error bars are shorter, indicating higher precision, under cloud-free conditions. The absolute SD shows similar results for fully covered conditions compared to the overall results, while partly cloudy cases present the largest values. Conversely, for relative values, the SD results show the opposite pattern, with the partly cloudy error bars generally more similar to the overall error bars and the bars for fully covered cases being larger. The median values generally exhibit similar behavior to the MBE values, showing comparable differences between the overall results and the different sky conditions and similar values to the corresponding MBE results. In absolute terms, the median shows more significant differences compared to the mean values for partly cloudy conditions with lower SZA values, as well as for cloud-free conditions, and the overall results at 45°. For the relative results, the largest differences from the mean values for each angle bin occur under fully covered conditions; for the other sky conditions, these differences increase at higher SZA values. Despite these differences, both the absolute and relative results for the median fall within very similar ranges to the MBE results. The boxes representing the 25th and 75th percentiles exhibit similar behavior to the error bars but with lower values. The relative MBE and SD were calculated and are shown in more detail in

Table 2. Mean bias error (MBE) and standard deviation (SD) of the shortwave global horizontal irradiance (GHI) differences between the predicted CNN-CMF model and measured data (in relative values) for different solar zenith angle (SZA) bins.

SZA [°]	40	45	50	55	60	65	70	75	80
oktas	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]	MBE ± SD [%]
0	−1.1 ± 4.1	0.4 ± 6.4	1.3 ± 6.7	−0.5 ± 6.4	−0.7 ± 7.7	−1.9 ± 10.1	−1.8 ± 13.0	−6.5 ± 15.1	−12.8 ± 22.8
1	−2.8 ± 7.5	0.4 ± 11.1	−1.3 ± 9.8	−0.1 ± 12.1	−0.2 ± 14.3	−1.9 ± 15.7	−3.6 ± 20.4	−11.8 ± 23.8	−18.2 ± 26.8
2	−0.9 ± 11.3	−2.1 ± 17.2	−0.8 ± 13.9	−1.0 ± 17.3	−2.0 ± 19.4	−1.9 ± 22.2	−2.5 ± 22.8	−6.1 ± 27.9	−8.6 ± 29.4
3	−0.0 ± 16.8	−0.4 ± 18.2	6.3 ± 18.2	2.9 ± 15.7	−0.0 ± 20.1	2.7 ± 20.1	3.5 ± 21.8	−4.5 ± 29.6	−10.5 ± 31.7
4	−3.2 ± 26.3	−2.4 ± 19.3	0.0 ± 21.7	1.8 ± 23.8	3.9 ± 21.0	4.1 ± 22.9	−4.6 ± 29.8	−3.7 ± 25.6	−9.7 ± 28.6
5	−1.0 ± 23.6	2.8 ± 23.4	3.8 ± 18.8	2.7 ± 24.1	3.9 ± 26.2	5.4 ± 25.1	−0.5 ± 26.8	−1.3 ± 27.0	−9.4 ± 32.7
6	3.0 ± 22.6	1.3 ± 25.9	6.3 ± 21.1	2.6 ± 23.4	3.7 ± 25.2	−1.3 ± 24.2	−1.5 ± 28.6	−5.9 ± 27.0	−10.6 ± 32.9
7	5.4 ± 24.1	8.9 ± 22.3	6.7 ± 22.7	2.0 ± 24.9	2.0 ± 27.0	−3.9 ± 28.7	−0.9 ± 28.5	−3.7 ± 27.8	−11.7 ± 31.4
8	−0.8 ± 29.1	4.7 ± 29.1	4.1 ± 28.4	4.1 ± 28.5	5.1 ± 26.4	3.8 ± 27.4	3.0 ± 31.8	3.8 ± 29.5	−5.2 ± 37.4
All	−0.0 ± 27.1	3.5 ± 21.8	3.3 ± 19.8	2.7 ± 21.8	3.3 ± 22.1	1.7 ± 23.4	0.9 ± 26.1	−0.7 ± 27.3	−8.5 ± 33.0

, where these results are presented for the individual cloud cover values. In good agreement with Figure 5, the results generally present more negative values for less cloudy conditions, except for larger SZAs (from 75° to 80°). The MBE ranges from −10% to 10%, showing accurate results, except 80° (ranging from −5% to −20%). The SD results indicate that the precision is higher for lower cloud cover and SZA values. The SD is generally larger for SZAs above 65°, ranging from 20% to 30% (and in a few cases, for 80°, above 30%), except for 0 oktas between 70° and 75° that present values of 13% to 15%, respectively. For the rest of the solar zenith angles, the SDs present larger values (around 30%) for eight oktas, ranging from 20% to 30% and between four and seven oktas. For 2–3 oktas, the results range from 15% to 20% and between 10% to 15% percent for 1 oktas. The lowest values correspond to 0 oktas, which varies from 5% to 10%.

4.2. Daily SW Irradiation

In Section 4.1, it was shown that the predicted GHI values with the CNN-CMF model accurately match the 1 min measurements, although with a precision of around 26%. However, in many cases, it is not necessary to know the minutely measurement of irradiance in detail, but rather, it suffices to know the daily value of SW global horizontal irradiation (GHI_d), as is the case in many studies on the temporal trends of SW radiation. By capturing sky images throughout the day, the CNN-CMF model is able to estimate the corresponding GHI, and by temporally integrating these predictions, the estimated value of GHI_d is obtained. This integration process could help to reduce the uncertainty of the inferred GHI_d value. In order to estimate the accuracy and precision of the GHI_d values predicted by CNN-CMF, the GHI_d values were calculated at Marambio as the daily mean of the GHI values multiplied by daytime length. This was carried out for the GHI predicted by CNN-CMF and also for the measurements taken by the pyranometer. In this case, day length was obtained by limiting the SZA to 85° in order to obtain the daily values with the available GHI data. The data have not been filtered by the variation between consecutive data, but the GHI_d values were removed if the amount of GHI data that day was less than the expected 70% of the GHI data.

The predicted GHI_d values have been represented as a function of the measured reference values in Figure 6, which shows a density scatter plot for the 482 available days. As a result, the GHI_d from the CNN-CMF model correlates with the measurements with a high determination coefficient $R^2$ of $0.99$. In addition, the linear fitting is close to the line 1:1, with the slope and the y-intercept equal to $1.00$ and 0.35 MJm⁻², respectively. The MBE value shows an accuracy of about 3% (0.38 MJm⁻²), with a precision given by a SD value of about 11% (1.32 MJm⁻²). These results indicate that, on average, the CNN-CMF model slightly overestimates GHI_d (by 3%) more than GHI (by 2%). This can be explained by the fact that CNN-CMF tends to overestimate GHI more for low SZA values (see Figure 5). Since the GHI values for the lower SZA values are typically higher than for larger ones, they contribute more to the final GHI_d, resulting in an overestimation. The precision of the CNN-CMF model is better for these values, as the SD reduces from 26% for GHI to 11% for GHI_d.

5. Conclusions

In this work, a new method was developed to retrieve GHI via CMF estimation from all-sky camera images. The images were captured and paired with GHI measured using pyranometers to train the proposed model at three different stations: Valladolid, Izaña (Spain), and Lindenberg (Germany). The model, named CNN-CMF, is based on a convolutional neural network, which was validated using a test dataset excluded from the training, that uses the same stations. The obtained results were compared against the GHI measurements from the pyranometers, showing a high correlation with a determination coefficient of $0.99$, which indicates that the model accuracy (MBE) is close to $-2\%$ and the precision (SD) is about $9\%$.

However, these results might be partially biased by the comparison with images similar to but not the same as the training dataset. In order to avoid any kind of autocorrelation between the image training and test sets, the sky images retrieved at another location, the Marambio station, were used to study the CNN-CMF model’s performance. The model can predict GHI with a slight overestimation of around $2\%$ with respect to pyranometer data. This overestimation is higher for lower SZA values and shorter for larger SZA values. The accuracy and precision of the model changed to values close to $3\%$ and $11\%$, respectively, when daily values are analyzed.

Retrieving radiometric measurements using pyranometers still continues to be the best way to monitor solar irradiance and detect possible changes and tendencies. However, in sites where these measurements are not available, this work proves that if sky images can be obtained, they provide a good proxy to estimate the GHI and the GHI_d. We encourage researchers to use the proposed CNN-CMF model. We also urge the scientific community to develop new models or improve existing ones for their benefit. We recommend using our CNN-CMF model for the reconstruction of a GHI time series, as well as also to comparing the model with measurements at other sites and improving the quantification of the uncertainty of this model. Finally, for future scope, we would like to use a long-term series of sky images to reconstruct the SW radiation values and check if this model is able to predict trends as they were observed in the past during global brightening or global dimming.