Data fusion of complementary data sources using Machine Learning enables higher accuracy Solar Resource Maps

Numerical model data are a convenient input product in case they are available for the full 30-year period over which we generate the resource map, and there are relevant products that cover the whole spatial domain of interest. Unfortunately, requirement R2) implies that many high-quality operational models that are not rerun in hindcast or reanalysis mode upon model update cannot be used here since, as discussed above, this leads to changes in the properties (including possible biases and error distribution) of the models.

Therefore, in this study, we use only two numerical model products as input for the ML models:

• •

  ERA5 [16]: this is the fifth-generation ECMWF reanalysis for global climate and weather, which is available globally over the period 1940–present. ERA5 data are available at approximately 31 km spatial resolution and 1-hour temporal resolution.

• •

  NORA3 [17]: this is a high-resolution numerical mesoscale weather reanalysis dataset from MET Norway. NORA3 data are available at 3 km spatial resolution and 1-hour temporal resolution.

We use several variables from both ERA5 and NORA3, as summarized in Table 1. As presented there, we collect variables that contain information about i) the incoming solar energy both at the top of the atmosphere and at the surface, ii) the atmospheric conditions focusing on cloud properties (cloudiness, height of the clouds, cloud types i.e. convective or stratiform), and iii) the ground conditions, i.e. surface albedo, type of vegetation and precipitation type. Data fields belonging to the category i) are naturally useful for the present task, as they constitute the models best prediction for either the GHI, or incoming solar radiation at the top of the atmosphere. Data fields ii) contain information about atmospheric conditions that influence the transfer of solar radiation such as absorption, scattering and reflection. Finally, data fields iii) contain information about the surface conditions that influences how much solar energy is reflected by the ground into the atmosphere, and hence contributes to indirect reflected irradiance (and to errors and uncertainties in the satellite products). Moreover, this surface state information is also potentially useful in combination with the satellite data presented in the next section, since the accuracy and bias of these satellite data is expected to be dependent on the state of snow and ice coverage.

Although numerical models are taking great care to estimate quantities such as GHI as accurately as possible, it is challenging to accurately model all the physics involved in determining the GHI, so we can expect that there are still biases, uncertainties, and inaccuracies in the GHI estimates the models provide, which will be shown in the results later in this paper, and as previously reported by Perez et al. [8].

We do not include information on the latitude, longitude, or altitude of the stations, nor on the time or date. This is because the number of pyranometers in Norway, and hence the number of locations used for training, is limited in both time and space. Therefore, including this information (which would otherwise be relevant) may lead to ML algorithms partially managing to identify each station individually and, hence, overfitting. The influence of these parameters on relevant quantities is still available to the ML models through their influence on the predictors, so that we expect that not providing these parameters has only a weak influence on the performance of the models.

The model data are nearest-neighbor interpolated to the position of the pyranometers to generate the training dataset. When generating the final product, these data are nearest-neighbor interpolated to the grid points where each GHI prediction is performed.

Satellites are key tools for Earth monitoring [19], including measuring solar irradiance [20, 21]. We use solar irradiance measurements from geostationary satellites as input data products. The advantage of these data is that they are available regularly (every 30 minutes for the product selected here) in the entire satellite coverage area. This makes geostationary data easy to incorporate into ML models that need collocated input data in space and time to perform predictions. The downside is that geostationary satellites are not well suited to measure irradiance further north than 65° north (for the product used over Norway) because of the high satellite viewing angle at these high latitudes. This means that satellite data are not available for the northern half of Norway. We do not consider the use of polar orbiting satellites in this study, due to the challenges associated with collocating their swaths in time and space with other measurements and the irregular time coverage at any given location.

More specifically, we use geostationary measurements of solar irradiance as provided by the EUMETSAT CM SAF SARAH-3 product [22]. The data are openly available from https://navigator.eumetsat.int/product/EO:EUM:DAT:0863, at a time resolution of 30 minutes and a spatial resolution of 0.05°. Using the SARAH-3 dataset ensures that, while the actual satellites being used and the instruments they carry vary within the 30-year period, the data processing methodology is homogeneous over the full time span considered. In particular, best efforts are made by the data producer to avoid sensor drift in time (see the self calibration procedure highlighted in [23, 24]). When either training or prediction are performed at a given hour $H$, we provide as input to the models the three closest SARAH-3 estimates, corresponding to the set of time instants $\{H$ - 30 min, $H$, $H$ + 30 min$\}$. Interpolation to the location of the pyranometers (which is performed for the generation of both the training dataset and the validation dataset), or to each grid point (for the GHI product map), is used to obtain the SARAH-3 data at the set of locations needed.

Geostationary satellite-based products, such as SARAH-3, also have several challenges associated with the method deriving surface irradiance. Indeed, this method is quite complex: typically, satellite products first characterize the clouds, then assume a clear sky irradiance which involves both climatology data and modeling, and finally they use these data to estimate the surface irradiance [23]. This involves several steps, and complex calibration and geophysical transfer functions that are based partly on heuristics. As a consequence, even satellite products are not a true, direct observation of the GHI. In particular, it is well known that these products, despite continuous improvements and the development of novel algorithms for the SARAH-3 product version, still have some level of difficulty distinguishing low-altitude clouds from snow or ice cover [25, 26, 27, 28, 29], and that they also struggle in foggy conditions [23]. Such challenging conditions, specifically the presence of snow and ice covers, are particularly common and relevant in Norway, given the local weather and climate conditions. As a consequence, we expect it to be possible to improve upon SARAH-3, even though it is already a high-quality product for solar irradiance, by applying ML-based post-processing. Moreover, effects due to the high viewing angle at the latitude of continental Norway, as well as the extent of the satellite pixels footprint, can also contribute to added errors and uncertainties in the SARAH-3 data.

In addition to the predictor data from both models and satellite observations described above, we also use in-situ data from pyranometers as a ground truth to train and validate our models. Surface measurements of hourly global horizontal irradiance were collected at 106 locations in Norway between 2016 and 2020. This time period was selected due to good data coverage (see [30]). Earlier than 2016 much less measuring stations were available in Norway. Therefore, this is the 5-year period that is used for training and validation of our ML models. Most of these data are openly available through the Frost API of MET Norway at https://frost.met.no. These have also been supplemented with additional data from four locations collected in the course of this project.

Since the pyranometer data had not undergone prior quality control, visual inspection and quality flagging routines have been implemented to only use high-quality station and measurement data to train and validate our ML models [30]. More specifically, stations with insufficient data coverage were excluded first, leaving 69 stations that met the criteria of at least 80% overall coverage and 60% average coverage for any month of the year in 2016–2020. Then, the data were inspected using both automatic and visual techniques such as time series analysis and comparisons with clear-sky GHI data to identify spurious trends and other significant shortcomings in any station data. As a result, 22 stations were rejected due to issues such as shadowing effects or data inconsistencies, leaving 47 stations in the final dataset. The data from the 47 remaining stations have been flagged for exceedances of theoretical limits, temporal variability, and data distributions [30]. Only valid data have been kept in the final training and validation dataset and the hourly GHI data have been collected in a netCDF file. The dataset is freely available at https://zenodo.org/records/8082726 [31], together with visualization and quality control plots, scripts, and files that include flags and uncleaned data, as well as the original raw data of the 106 stations.

When training and validation of our ML models are performed, we split our data between two disjoint sets of stations used for each purpose, so that we can detect (and avoid) any overfitting. The split between the stations used for training and validation is illustrated in Fig. 2. We use 8 stations for validation purposes. The split between stations used for training versus validation is performed so as to represent the general range of conditions over Norway as well as possible in both subsets, based on the general topography and climate conditions observed over Norway.

As visible in Fig. 2, there are in situ pyranometer stations in the southern and northern parts of continental Norway. Moreover, the different geographic and climatic areas of continental Norway are well covered: there are stations both on the coast, in the plains of southern Norway, and in the mountain ranges that run along the country. Although there is an area between approximately 64° N and 67° N, where unfortunately no station data of sufficient quality is present, there is nothing particular about the geography and climatic conditions in this area relative to those immediately north and south of it. Therefore, we expect that the models trained in this work are still able to perform well in this area, since they have seen conditions similar to those encountered there in their training dataset. This, combined with the fact that simple and robust models are used, as will be described in the next section, mitigates the lack of in situ ground truth data in this specific area. We hope that studies in a few years will be able to remediate this limitation, as more pyranometers are now being deployed over Norway; however, at present, we can only use the data available to us. Since, as will be described in more detail in the next sections, the models rely on pyranometer data only for their training and validation, we are still able to generate GHI data over the corresponding area. Naturally, we expect that having more and better-distributed stations in the future will enable to develop better, larger, and more expressive correction models before overfitting happens, which will allow even better and more accurate ML-corrected maps to be developed.

The first model kind we use is a simple linear regression model, taking as input all the parameters presented within Table 1. LR models are well known and can be seen as the simplest form of machine learning models [32, 33]. Despite their simplicity, LR models are known to be effective in many real-world applications [34], in particular when several estimators with different and separate error sources (i.e., estimators with non-fully correlated errors) are combined into a Best Linear Unbiased Estimator (BLUE) [35, 34]. This is naturally the case in the present study, as the ERA5 and NORA3 numerical models have quite different setups, and the satellite data provide a separate estimate that is derived from observations.

In the present study, our LR model is simply obtained by applying the $LinearRegression$ function from the sklearn software package [36] to our set of predictors using the ordinary least square optimizer, while selecting pyranometer data as target. While we have investigated the use of more sophisticated optimizers, e.g. ridge regression (see e.g. [37, 38]; ridge regression works by imposing a penalty on the regression coefficient magnitudes using a L2 norm), we found during testing that the plain ordinary least-square optimizer is providing the best results on our task, and it is hence used in the following.

The second kind of models that we use is a fully connected neural network (FCNN). Neural networks have been known for a long time (for example, their property of universal approximators was already discussed in the 1980s [39, 40]). However, their massive popularity increase happened arguably only recently, concomitantly with the drastic rise in dataset size and available computational power that we have witnessed since the early 2010s [41]. Our motivation for using a simple FCNN is that these models are effective and flexible non-linear function approximators that are supported by a wide range of software packages and easy to train. FCNNs do not require specific prior knowledge about the problem being considered. While large scale models are attracting most attention over the last couple of years (e.g., large image analysis networks [42, 43], the large language models (LLMs) revolution [44, 45], diffusion models [46, 47], large reinforcement learning models [48], or large data-driven global weather models [49, 50]), even simple FCNNs of modest size and complexity are known to be effective nonlinear function approximators that can help on tasks such as error correction [51, 52, 53] or data fusion [54], given that enough training data are available to perform a representative training.

As simple neural network theory, including the theory describing FCNNs, is now well known across many fields and covered by a very large number of books (e.g., [41]) and reviews (e.g., [55]), we will not reproduce boilerplate materials that discuss simple neural network theory here. Instead, in the following, we consider that the reader has a basic knowledge of the general technicalities of NNs and FCNNs, and we refer the reader without such a background to, for example, [41]. For all practical matters in the present study, FCNN models can be thought of as simple nonlinear function approximators, that can be effectively tuned through gradient descent on the training set to represent the mapping from the combination of model and satellite predictors to the target pyranometer GHI measurements.

We expect FCNN models to work similarly to the LR model described above, except that they have more freedom in designing the mapping function from the predictors to the model output, which can be expected a priori to translate into better accuracy. Indeed, the FCNN can include non-linearity into its mapping, which is not possible for the LR model. Given the relatively moderate amount of data available in the present study, and the fact that both numerical models and satellite observations already provide estimators of the GHI target that have been tuned and optimized by the geoscience community over several decades, we expect moderate accuracy improvements to be obtained with the LR and FCNN models relative to the set of GHI predictors. The improvements obtained by the FCNN model will be the result of both the BLUE effect (similar to the LR model), and possibly the compensation of minor remaining non-linear biases and distorted couplings in the physics of the models and geophysical transfer function of the satellite system, if these exist.

From a practical implementation point of view, neural network methods are known to be potentially sensitive to class imbalance [56]. Therefore, since over half the raw training data correspond to night or very low light conditions, we exclude samples for which the GHI is lower than $P_{min}$=50 W/m² (which is typically around 5% of the mid-day summer value) from our training and validation dataset. This means that the FCNN that we train can only be used when the incoming solar radiation is equal or higher to $P_{min}$. In other cases, the LR model will be used.

We have tested several different FCNN architectures, optimizers, loss functions, and metaparameters in the early phases of the training process. We found that a simple FCNN with 4 layers of 60 neurons each, using the rectified linear unit (ReLu) activation function, the mean absolute error loss, and the Adagrad optimizer, gives the best validation loss over the present task. This results in a small network with slightly more than 13k trainable parameters in total. We use non-trainable linear scaling normalization layers to scale the input and output to the FCNN, so that the input and output of the hidden neuron layers, as well the neuron activation values inside the FCNN, are typically scaled to the range [-1,1]. The network and its training and prediction routines are implemented in Tensorflow/Keras [57, 58]. We find through trial and error that an initial learning rate of 2e-4 (which is then tuned by the Adagrad optimizer as learning progresses) results in satisfactory training within around 2000 epochs. These parameters are typical of what is expected for a small and simple FCNN as implemented here. We have also tried introducing more complexity in the FCNN setup in exploration phases of the training process, such as inserting dropout layers, adding a regularization term to the loss function, inserting batch normalization layers, and using more and larger layers, but this did not provide meaningful improvement, and sometimes even degraded (due to overfitting), the validation loss. This is typical of training processes that correspond to a relatively simple task, for which a moderate amount of data are available, so that adding too much additional complexity in the FCNN does not help improving the validation loss (and can even lead to overfitting). This is consistent with our experience on previous similar projects where moderate amounts of in-situ data were available, and models were already providing a reasonably good prediction of the phenomena considered, so that the FCNN only needs to add small changes on top of the predictors already provided [53].

The main operational goal of this work is to generate a ML-optimized GHI map that covers continental Norway over a 30-year period. For this, we need to be able to produce an estimate of the GHI even in cases where a specific input product is unavailable (because it is not available over the whole domain, as is the case with the SARAH-3 data, or because there were some operational production issues for a specific time and place). As a consequence, we train several flavors of both the LR and FCNN models, using different sets of input data. These are summarized in Table 2. We then produce GHI estimates using each of the models that can be run (i.e., for which the predictors required are available) at each point in space and time, and we select the best (i.e., the one with the highest model_index following the validation presented below) that can be successfully produced as the corresponding best estimate for the GHI. Hence, if all data sources are available at a given location and time, the best available model (i.e. the FCNN model using the predictors from all the data sources) is used to estimate the corresponding solar irradiance. If some data are not available for a given point in space or time to apply the best ML model, a simpler model is used instead for the corresponding time and location estimate.

Note that more models could be built and trained, using different subsets of predictors in addition to what is presented in Table 2. For example, one could train both LR and FCNN models that use only the ERA5 or only the NORA3 predictors (effectively creating 4 more models, LR_ERA5, FCNN_ERA5, LR_NORA3, FCNN_NORA3). Due to how the present project developed, and owing to the more modest gains that can be obtained from applying model correction to a single category of predictors at a time (see Appendix C), this is not implemented in the final GHI product that we provide. Moreover, since ERA5 and NORA3 are reanalysis, both these model data are available with good quality for the vast majority of samples in time and space, so that it is seldom that the model_index 3, or a fortiori 2, cannot be produced. Therefore, in our context, leveraging additional LR and FCNN models based on NORA3 or ERA5 alone would provide a limited improvement only (this is documented in more detail in Appendix C). Although such models could also be used in the future during the generation of the GHI map, this is not the case here.

Note that the pyranometer data used in the present study are only covering continental Norway, over land. Therefore, sea areas and other land locations, such as Svalbard, are not covered by the final GHI map product, since the ML models are not trained and validated over these areas. This is enforced by masking out the areas where the ML models are not applicable in the final GHI product. As mentioned previously, the SARAH-3 data are not available further than 65° north, since this data source is based on geostationary observations that cannot be made close to the poles. Therefore, only model indexes 3 and lower are obtained in the northernmost part of Norway.

The full input FCNN model (model_index 6) relies on all data sources, and is the most sensitive to the input source data quality. More generally, FCNN models, being (non-linear) neural networks, are potentially the models that are the most sensitive to unusual or previously unseen (during training) input data. Therefore, we prefer not to risk extrapolating the FCNN models outside of the strict range of inputs for which they were trained. By contrast, the LR models, being simple linear models, are more robust to inputs that are slightly unusual. As a consequence, we take extra precautions with regard to the FCNN models. More specifically, we check that i) the incoming solar radiation is greater than 50 W/m² (since we set a lowest value cutoff of 50 W/m² to filter the training data to avoid class imbalance due to the night, as highlighter higher up), ii) all the input variables are valid within the range seen during training, before accepting that FCNN models can be used. This may not be the case if, e.g., strong precipitations happen that temporarily reduce the amount of solar radiation, or lead to extreme values on some of the predictors to the FCNN models – in which case, the LR models will be used instead. This is also not the case in some specific locations, such as glaciers in the Norwegian mountains.

As a consequence, GHI predictions with model_index 6 are only available south of 65° north (i.e., where geostationary satellite data are available), when all data sources are of high quality and when the GHI is higher or equal to 50 W/m². In all other cases, predictions of lower model_index are provided (in which case, we provide the output corresponding to the highest model_index model that could be produced given available data). Finally, due to challenges with some of the data sources over the period 2012-05 to 2012-08, only ERA5 data (corresponding to model_index 1) are used on the corresponding period North of 65° N.

The key metrics (computed on the validation stations) comparing the most important models are summarized in Table 3. The different columns indicate the values obtained from observations (i.e., the ground truth pyranometer data), GHI estimate fields available in numerical models (NORA3 and ERA5) and in satellite data (SARAH-3), as well as from ML models. ML models with the _NOSAT postfix (i.e., LR_NOSAT and FCNN_NOSAT) indicate models flavors that do not use the SARAH-3 data as predictors, and these are trained and evaluated over all of Norway. ML models without this postfix (i.e., LR and FCNN) indicate model flavors that use SARAH-3 data as predictors, and these models are trained and evaluated using only the data available south of 65° north, where SARAH-3 data are available. This is reflected in the number of samples used to compute each statistics, which is also reported in Table 3.

The results presented in Table 3 are in good agreement with what can be expected based on the characteristics of each product or model. The two numerical models (NORA3 and ERA5) have the highest level of both mean absolute error (MAE) and root mean square error (RMSE) relative to observations. NORA3 has slightly higher error levels compared to ERA5, but has also significantly higher and more realistic levels of variability in the predicted data (indicated by the "std. of data" column). This corresponds well to the fact that NORA3 is a higher resolution model compared to ERA5. Said otherwise, ERA5 has a tendency to spatially smooth out physical processes and the predictions that it provides, which produces slightly better error metrics at the cost of reduced intrinsic variability in its output. The SARAH-3 satellite data have better MAE, RMSE, and variability levels compared to the data from numerical models. Naturally, the SARAH-3 product still has non-negligible levels of MAE and RMSE, as well as slightly lower intrinsic variability, compared to the actual observations. This is a natural consequence of the finite footprint of the satellite data acquisition pixels, as well as of the post-processing that is necessary to derive an actual GHI observation from the physical signals acquired by the satellite, and a variety of error sources induced by high angle of incidence at the latitude considered. As noted earlier, snow and ice cover also add to the uncertainty of the satellite-derived irradiance.

Regarding the ML models, the FCNN models have slightly better performance compared with the LR models, both when satellite data are used and when they are not. Both the LR and FCNN models using all data inputs (including SARAH-3) have lower levels of MAE and RMSE compared to the SARAH-3 data. In the case of the LR model, this comes at the cost of a non-negligible reduction in the intrinsic variability of the predicted GHI, while the FCNN model has an intrinsic variability that is closer to that of SARAH-3 and the pyranometer observations. This can naturally be understood by the "purely averaging" effect obtained applying a linear regression model, compared to the non-linear properties of neural network based techniques. ML models that do not use SARAH-3 satellite data as an input (i.e., LR_NOSAT and FCNN_NOSAT) have, unsurprisingly, higher levels of MAE, RMSE, and slightly lower levels of intrinsic variability. This is understandable, as the SARAH-3 data input was the one that brought the highest accuracy GHI estimates, with the most realistic level of data variability. However, these _NOSAT models still significantly outperform the NORA3 and ERA5 GHI estimates by themselves, which indicates that they manage to extract meaningful information from the numerical model predictors provided, and to use this to improve on the GHI predictions.

The data presented in Table 3 can be extended and shown in Taylor plots, as shown in Fig. 3. Taylor plots show the amount of normalized intrinsic variability (radial distance) and the correlation coefficient (angle) between model and ground truth data. These are computed between the normalized ground truth pyranometer data (where normalization is performed individually for each station so that the normalized intrinsic variability reported for each station is matched to 1), and the normalized GHI prediction obtained at each individual station with each model (there, also, the same normalization is used for each station as was used for the normalized pyranometer data). This enables effective visualization of both the accuracy and the amount of intrinsic variability in each model output for each station (small dots), as well as the average performance of each model across each set of stations (thick crosses). The Taylor plots are shown for both the training (left panel) and the validation (right panel) subsets of stations.

The general pattern for models location in the Taylor plots in Fig. 3 is very similar between the training and validation datasets, with variability for individual stations compared to the average values over each set of stations which comes from the intrinsic variability between stations, as can be expected given the wide range of conditions encountered over Norway (e.g. coastal areas, plains, mountains, and both temperate and sub-polar climate regimes). The very similar localization of the results across the Taylor plots for both the training and validation cases confirms that there is no noticeable overfitting. Indeed, this implies that there is no significant performance degradation observed by applying the models on previously unseen data compared to the training data. Moreover, the results shown in Fig. 3 confirm the general impression provided by the scalar metrics in Table 3. More importantly, the Taylor plots indicate a clear clustering of individual station data (small dots) around individual model cluster centers (thick crosses), confirming that the different predictors and models (different colors) have significantly different accuracy. The GHI estimates obtained directly from numerical models have the lowest coefficients of correlation (though these are already quite high, at over 0.85). ERA5 has higher correlation with pyranometer in-situ observation data than NORA3, though this comes at the cost of too low variability in the data, indicating increased levels of smoothing as pointed out previously. The LR model without satellite data is able to obtain slightly higher correlation than ERA5, but has typically the same low variability. The FCNN model without satellite data provides a bit higher correlation, and variability that is closer to the real data. The SARAH-3 data have higher correlation and better variability than all data sources and all models that do not rely on satellite observations. The LR and FCNN models with satellite information are the best GHI estimates available, with the FCNN being slightly better in both correlation and variability values.

All these results confirm that the LR model flavor is able to provide better GHI estimates than any of its input predictor taken in isolation, which is as expected from the BLUE theory highlighted earlier. In addition, the FCNN model, which also includes non-linearity, performs slightly better than the LR model, as could be expected. It is notoriously difficult to interpret FCNNs to try to understand their inner workings, and even the interpretation of simple models such as LR can be misleading [59]. However, we can still attempt to inspect how each model weights the different predictors available, and try to compare how the LR and FCNN models operate. For this, we present a measure of predictors importance in Fig. 4. For the LR model, the predictor importance (predictor_importance_LR) is defined as:

Computing a predictor importance is more challenging in the case of the FCNN. Since the overall performance of the FCNN is relatively similar (though slightly better) to the one of the LR, we can hypothesize that the FCNN is relatively close to a linear regression model, with some additional non-linearity added on top of the linear LR behavior. Therefore, we expect that the FCNN is relatively close to its linearized approximation at each point in the prediction space. As a consequence, we leverage the local linear sensitivity of the FCNN relative to each predictor, averaged over the whole dataset, to compute the following predictor_importance_FCNN estimate:

The results for the two best performing models, i.e. the LR and FCNN models using all predictors, are presented in Fig. 4. We present both 1) the importance of each predictor for the LR and FCNN models taken individually (top panel), and 2) the aggregated importance for groups of closely related predictors for the same two models (bottom panel). The reason for presenting both the individual and aggregated predictors importance is that when several predictors are closely related (for example, the surface_solar_radiation_downwards_era5 and the surface_net_solar_radiation_era5 predictors are very similar to each other, and differ only due to the effect of relatively minor contributions), the models may have some degree of spurious sensitivity to these predictors. In particular, the models may attribute importances with canceling effects to these closely related predictors. Hence, the true importance of these predictors taken as a whole may look different from their individual (canceling) importances, as visible in Fig. 4.

The general results presented in Fig. 4 confirm our previous impressions about the behavior of the models. Typically, more importance is given to the SARAH-3 predictors, which are also the highest quality ones. This is consistent with what would be expected from the BLUE theory. In addition, some importance is given to both the NORA3 and ERA5 predictors. Interestingly, non-negligible importance is given to the surface state predictors. This may indicate (though this is speculative as it is challenging to interpret the weights and sensitivities of a NN or even a LR model) that this information is used to correct, e.g., biases and errors in SARAH-3 measurements due to snow and ice cover. Generally, Fig. 4 confirms our previous interpretation that the FCNN model is relatively close to the LR model, with some additional non-linear modulation taking place that explains for the slightly better performance. However, we also want to stress that one should not over-interpret such predictor importance map, as similar techniques used in image analysis and that have inspired our present approach, such as Grad-CAM [60] and saliency maps [61], are notorious for being difficult to interpret into details and are known to be potentially misleading.

We generate the final product map by applying the models described in Table 2 over continental Norway with a time resolution of 1 hour and a spatial resolution of 0.03° (approximately 3.3 km), covering the 30-year period 1991–2020 (from the start of January 1st, 1991, to the end of December 31st, 2020). In order to do so, we interpolate (using the nearest neighbor method) the input data used by the ML models to the grid on which the final GHI map is produced. The choice of a 0.03° resolution for the final product is motivated by i) the size on disk of the final map product, and ii) the resolution of the input data. Regarding i), the final map is about 45 GB is size at the current resolution. This is a manageable amount of data that allows cheap hosting, downloading a local copy of the data if necessary, fast access and slicing, and is generally convenient to work with. Regarding ii), ERA5 has a resolution of 31 km, NORA3 3 km, and SARAH-3 0.05°. Therefore, the input data with highest resolution, NORA3, typically has a resolution slightly better than 0.03°. Consequently, it would not make sense to generate a GHI map using our methodology while “artificially” interpolating to a resolution much higher than that of the finest data input. Naturally, it is still an interesting question to generate even higher resolution maps, but this would rather be the topic of a project focusing on generation of super-resolution data, which may require additional and different methodology compared to what we use here.

Links to the product files are provided in Appendix A. If all data sources are available at a given location and time, the best available model (i.e., model_index 6) is used to estimate the corresponding solar irradiance. If some data are not available for a given point in space or time to apply the best model, a simpler model (the one with highest model_index that can be applied to the case) is used instead for the corresponding time and location estimate. For each location and time, the best (i.e., a higher "model index") data output that is available is used in the final product, and these data are collected together with the corresponding model_index. The data files of the product are hosted as netCDF-CF files on the THREDDS server of MET Norway: https://thredds.met.no/thredds/catalog/sunpoint/ML-Optimized-Maps/catalog.html.

More technical details are provided in Appendix B. In particular, we want to stress that the data are provided following 2 different data-chunking layouts, and that choosing the right layout is critical for data retrieval performance.

A typical map of the predicted GHI (left panel) and the model_index used at each point (right panel) for a specific point in time is presented in Fig. 5. We see that the FCNN model with SARAH-3 input is used on most of the southern parts of Norway, as expected. Concomitantly, the FCNN model without SARAH-3 input is used on the northern part of Norway. There are a few locations where either the GHI is too low to use the FCNN models, or some predictors are taking unusual values due to specific weather conditions at the time of the snapshot that is presented, and other methods, in particular the LR models output, are used over the corresponding areas. This is especially visible in the area around 64N, where heavy cloud coverage is present which locally reduces the GHI for the specific date and time shown. A map of the averaged accumulated yearly GHI (kWh/m²) for the whole period 1991-2020 is presented in Fig. 6. Values above 1000 kWh/m² are found along the coast of southeastern Norway and in the high mountain region of southern Norway. Lower annual values are seen in the west coast and in central Norway (800-850 kWh/m²) dropping to 700-750 kWh/m² in northern Norway. It is also noted that higher annual values are encountered in the mountainous areas in northern Norway. The values obtained in southern Norway are in line with what is reported by other products, e.g. the GSA shown in Fig. 1. However, with our product, data are available over all of Norway (also where there are no data in Fig. 1 due to the unavailability of geostationary satellite data). There is significant spatial variability in the mean accumulated yearly GHI over Norway. The general variations observed correspond well with physics-based intuition and expectations, with regions in Northern Norway receiving less accumulated solar energy because of the higher latitude and the associated climate. Mountainous areas in Southern Norway are receiving the highest amount of accumulated solar energy in late spring/early summer, likely due to altitude effects (high clear-sky radiation because of few aerosols and low water vapor), high insolation due to low cloud-cover and high albedo values from snow (reflection). Similarly, the southeastern coastal area is characterized by many days with clear skies or little cloud cover in summer, giving rise to high annual average surface irradiance.

Finally, we compare our final GHI ML product with in-situ pyranometer data and other products at the validation stations. This is illustrated in Fig. 8. We see, similarly to what has been presented in Fig. 3, that the GHI ML product compares well with observations. The ML product is able, despite the inherent variability associated with individual stations and the use of a combination of the different model_index predictions at each location, to reduce biases in average compared with other data sources taken individually.