Day-Ahead Hourly Solar Irradiance Forecasting Based on Multi-Attributed Spatio-Temporal Graph Convolutional Network

Solar irradiance forecasting is fundamental and essential for commercializing solar energy generation by overcoming output variability. Accurate forecasting depends on historical solar irradiance data, correlations between various meteorological variables (e.g., wind speed, humidity, and cloudiness), and influences between the weather contexts of spatially adjacent regions. However, existing studies have been limited to spatiotemporal analysis of a few variables, which have clear correlations with solar irradiance (e.g., sunshine duration), and do not attempt to establish atmospheric contextual information from a variety of meteorological variables. Therefore, this study proposes a novel solar irradiance forecasting model that represents atmospheric parameters observed from multiple stations as an attributed dynamic network and analyzes temporal changes in the network by extending existing spatio-temporal graph convolutional network (ST-GCN) models. By comparing the proposed model with existing models, we also investigated the contributions of (i) the spatial adjacency of the stations, (ii) temporal changes in the meteorological variables, and (iii) the variety of variables to the forecasting performance. We evaluated the performance of the proposed and existing models by predicting the hourly solar irradiance at observation stations in the Korean Peninsula. The experimental results showed that the three features are synergistic and have correlations that are difficult to establish using single-aspect analysis.


Introduction
Extensive growth in the global population has led to an increase in the use of fossil fuels and greenhouse gas emissions, leading to worsening environmental pollution and global warming problems [1]. In 2015, the United States and China pledged to achieve 100% reliance on renewable energy to tackle climate change [2]. In addition, the European Union has decided to reduce greenhouse gas emissions and transition to renewable energy entirely by 2050 [3]. Among renewable energy sources, solar (photovoltaic) energy is estimated to meet a quarter of the electricity demand by 2050 [4]. However, because various factors, such as solar position, time, geographical location, and meteorological conditions, affect solar power generation, the efficiency of solar power plants is highly volatile [5,6]. Volatility also causes problems such as output instability of solar power plants and overloads of power grids, which should be addressed for commercializing solar energy [7][8][9][10]. Therefore, methods for accurately forecasting solar energy production in a specific region have become essential [11,12]. Various solar irradiance prediction models, from statistical to neural network-empowered, have been proposed for providing a scientific basis for managing solar power generation and power grid overloads [13][14][15][16].
Conventional solar irradiance forecasting models can be classified as physical, empirical, and statistical models. The physical approach represents meteorological conditions in a and analyzing the influences than grids and CNN models, respectively. Jiao et al. [43] composed a graph with observation stations and their adjacencies. Adjacency was defined based on Pearson correlation coefficients (PCC) between the historical solar irradiance at each station and the minimum threshold. They then predicted future solar irradiance by analyzing the spatial influences of the GNN layers and applying LSTM layers to the spatial features extracted by the GNN layers. Khodayar and Wang [44] applied a similar approach to wind-speed forecasting. They defined spatial adjacency using the mutual information between historical wind speeds and directions. In contrast to Jiao et al. [43], Khodayar and Wang first extracted temporal features from each station using LSTM and then used the GNN for analyzing the spatial correlations between the temporal features. Dong et al. [45] used dilated CNN layers instead of RNN layers for analyzing the temporal characteristics of wind power production. Because CNN layers are more effective in parallelizing its computation than recurrent models, this approach could reduce the time consumption for both model training and forecasting services. However, CNN layers have limitations in explicitly considering the time-sequential characteristics of the meteorological data. Muthukumar et al. [46] predicted PM 2.5 concentration by combining graph convolutional network (GCN) and ConvLSTM and constructed a graph using the distance between fine dust sensors and interpolated data from unobserved locations using the GCN layer. The output of the GCN for spatial interpolation was converted into an image and input into ConvLSTM. Spatiotemporal models performed better than models without fusing spatial and temporal information [47,48]. However, these models cannot utilize a variety of meteorological data (e.g., humidity, temperature, air pressure, and cloudiness) observed with prediction targets. This makes the conventional models far from understanding the weather conditions of observation stations and their spatiotemporal correlations. Cheng et al. [49] used GNN to analyze correlations between atmospheric variables, and Huang et al. [50] extracted the temporal features of each variable and analyzed feature correlations using MLP. However, these studies omitted the spatial influences between observation stations.
The weather conditions of spatially adjacent observation stations influence each other; for example, clouds move with wind. The influences occur with non-uniform time lags, and weather conditions have temporal patterns. Sunrise and sunset create daily patterns, and yearly patterns are correlated with the regional climate. Prediction targets and a few meteorological variables related to the targets (e.g., wind speed and direction) are insufficient in providing contextual information on the weather in a region. Thus, analyzing spatiotemporal correlations between various meteorological variables with an end-toend network will improve the performance of weather forecasting models. In addition, as discussed above, models that are not based on atmospheric knowledge have limited model interpretability. Although several existing studies have attempted to combine multiple features, they did not closely examine the effects of combining the three features on weather forecasting with a case study of solar irradiance. Therefore, we first developed a novel solar irradiance forecasting model that considers (i) temporal patterns of meteorological variables, (ii) spatial influences between observation stations, and (iii) correlations among a variety of meteorological variables. The weather data were represented as a graph, with the observation stations as nodes, the spatial adjacency of the stations as edges, and meteorological variables as attributes. Thus, the graph exhibits static structures and dynamic attributes. Then, we extended the attribute-augmented spatiotemporal GCN (AST-GCN) model [51], which considers both static and dynamic attributes, to analyze spatiotemporal correlations between multiple meteorological variables. The AST-GCN model consists of graph convolution layers and recurrent layers, which extract the spatial and temporal features of the dynamic networks, respectively. However, the graph convolution layers also analyze the temporal changes in the dynamic attributes using a fixed-length window. We call the proposed model for multivariate spatiotemporal analysis of meteorological data "Multi-attributed Spatio-Temporal GCN (MST-GCN)." To examine the effects of the feature combination, we compared the performance of the proposed model with baseline models, which are based on each part of the three features, by adjusting the prediction sequence lengths, seasons, weather conditions, etc. Based on this comparison, we attempted to validate the following research questions: • RQ1. Weather conditions of spatially adjacent observation stations influence each other, and the influence is significant in predicting solar irradiance. • RQ2. Temporal changes in historical weather data are effective in solar irradiance forecasting. • RQ3. Meteorological variables observed at a station have correlations with future solar irradiance of the station.
The performance of the proposed and existing models demonstrated the contribution of each feature to the aspects of weather forecasting. The performance comparison between the models showed that the spatial, temporal, and multivariate features complemented each other and were synergistic. The main contributions of this study can be summarized as follows: • We propose MST-GCN, which allows for spatiotemporal analysis of dynamic multiattributed networks to conduct day-ahead hourly solar irradiance forecasting for multiple stations. Our proposed model consists of GCN layers for spatial features, GRU layers for temporal features, and multi-attribute fusion modules for multivariate features to fuse the three features of meteorological data. • We demonstrated the superiority of MST-GCN in terms of forecasting performance and stability over the baseline models, including T-GCN (spatiotemporal), GRU (temporal), GCN (spatial), and MLP (multivariate) with intensive experiments. Furthermore, we verified the above research questions, RQ1, RQ2, and RQ3, by comparing T-GCN with GRU, T-GCN with GCN, and MST-GCN with T-GCN, respectively.
The remainder of this paper is organized as follows: Section 2 describes the acquisition and pre-processing of the meteorological data used in this study. Section 3 presents the proposed solar irradiance forecasting model. Section 4 details the experimental procedures and results used to evaluate the performance and practicality of the proposed model. Section 5 presents the concluding remarks and discusses the limitations and future directions of this study.

Data Acquisition and Preprocessing
This section describes the procedures for acquiring meteorological data used to evaluate the proposed model and validate the research questions.

Resource Data
There are two methods for measuring solar irradiance. The first method uses a pyrometer, and the other indirectly estimates solar irradiance by analyzing satellite images. Although the pyrometer can accurately measure the amount of insolation per hour, it has disadvantages in terms of the high cost of the measurement system and the limited measurable range [30,52,53]. Satellite image analysis has advantages in observing solar irradiance over a wide area. However, this method also has difficulties in real-time estimation, owing to the characteristics of satellite imaging [54,55]. In addition, because satellite images are taken over clouds, cloudiness can cause images to be much less accurate than images from ground observations. Owing to this inherent limitation, satellite image analysis showed a relatively lower accuracy than the pyrometer. Uncertainties in the observation data also affect the performance of solar irradiance forecasting models using measurements. Most of the existing forecasting models based on satellite images have lower and less stable accuracies than pyrometer-based models [56,57]. Thus, we acquired meteorological data collected by automated surface observing systems (ASOS), which are based on a pyrometer and more accurate than satellite image analysis, to reduce uncertainties caused by input variables.
The ASOS Programme is a joint effort of the National Weather Service (NWS), Federal Aviation Administration (FAA), and Department of Defense (DOD). The ASOS serves as the nation's primary weather-observing surface network. This system was designed to support weather forecasting and aviation operations. Simultaneously, the ASOS supports the needs of meteorological, hydrological, and climatological research communities [58,59]. The ASOS conducts time-synchronized ground observations at every participating observatory for obtaining time-sequential data for atmospheric conditions. In addition, the system automatically measured meteorological variables using synoptic meteorological observation equipment. The observational data were accessible through a public repository (https://data.kma.go.kr/cmmn/main.do (accessed on 15 August 2022)). Among the ASOS observation stations, we selected 42 stations measuring solar irradiance since 2017 and located in the Korean Peninsula, as shown in Figure 1.

Meteorological Variables
Solar irradiance (S r ) is closely related to the geographical factors of observatories, the date and time of observations, and other meteorological variables (e.g., cloudiness and precipitation). From a geographical perspective, solar irradiance varies with latitude and longitude. As discussed in Section 1, the spatial adjacency between observatories indicates that weather contexts can influence each other. In addition, time-sequential analysis can establish daily and seasonal patterns of solar irradiance. Therefore, utilizing these spatial and temporal features can improve the performance of solar irradiance forecasting models [60,61]. From an atmospheric standpoint, weather contextual information, which can be inferred from meteorological variables, is significant for predicting solar irradiance [30]. For example, if it is a cloudy day, the sun is blocked and the amount of insolation reaching the ground decreases. Therefore, to consider weather contexts with high variability, we gathered meteorological parameters correlated with solar irradiance, and spatiotemporal parameters. Among the variables in the ASOS data, we selected 17 variables in the three categories as input parameters for the proposed model, as listed in Table 1.

Data Preprocessing
The ASOS data have a significant number of missing values, and interpolating the omitted observations can cause uncertainties and affect the performance of the forecasting models. Thus, for fair evaluation and validation, we removed the variables and adjusted the observation period for avoiding missing values. However, a few values are significantly correlated with solar irradiance and are not difficult to reliably substitute for omitted values. For precipitation, we checked records from the Korea Meteorological Administration for regions where observation stations with missing precipitation values were located. If there was no precipitation when missing values occurred, we replaced them with zero. We examined sunrise and sunset times in cases of missing sunshine duration and solar irradiance. Because insolation cannot exist between sunset and sunrise (e.g., 21:00 KST to 05:00 KST), we replaced the missing sunshine duration and solar irradiance values in the period with zero. As a result, we gathered hourly observation data for four years (from 1 January 2017 to 31 December 2020), including the 17 meteorological variables observed at the 42 observatories. The first three years of data were used to train the proposed and baseline models, and the remaining year was used for model evaluation.

Methods
We propose a novel solar irradiation forecasting model that considers (i) spatial features, (ii) temporal features, and (iii) correlations between meteorological variables. First, we represented the ASOS data as undirected networks with multiple dynamic attributes. We then modified and extended the existing spatiotemporal GCN models [51,62] to analyze the spatial and temporal correlations between these dynamic attributes.

Meteorological Networks
The proposed model conducts solar irradiance forecasting by analyzing (i) spatial correlations between ASOS stations, (ii) historical patterns of meteorological variables, and (iii) correlations of solar irradiance with the variables. These three viewpoints will enable the proposed model to establish weather contexts at each ASOS station and to predict future weather by understanding the spatiotemporal influences between the stations. First, we represented the spatial correlations as an undirected network and historical meteorological variables observed at each ASOS station as the dynamic node attributes of the network.
Most of the existing studies defined correlations between meteorological observation sites by using mutual information [49,63] and the PCC [43,45]. However, the influence between the two observation sites will have inconsistent time lags according to distances, landforms, weather contexts, and so on. Moreover, mutual information and correlation coefficients are not proper metrics for detecting the influence of dynamic time lags. Therefore, we defined the correlations between ASOS stations using geographical distances. Then, the relative correlations among the adjacent stations can be learned by the GCN layers in the proposed model. For the same reason, although existing studies [43,45,49,63] defined the adjacency between stations by using minimum thresholds for correlations, we searched for N-nearest neighborhoods of the stations according to the distances. We call the network representing dynamic weather data 'meteorological network', and it can be defined as follows: Definition 1 (Meteorological Network). The geographical adjacency of the ASOS stations is described as N = (V, E ), where V = {v 1 , · · · , v V } is the set of stations and V is the number of stations. E e i,j is the set of edges, where v j is one of the N-nearest neighbors of v i . In addition, each node had meteorological variables (listed in Table 1) collected by the corresponding ASOS stations as dynamic attributes. Thus, the structure of the static meteorological network can be represented as an adjacency matrix A ∈ R V×V . Then, the dynamic attributes can be represented as a sequence of matrices, X = X 1 , · · · , X T , where T denotes the number of time points and X t ∈ R K×V refers to node attributes at time t when K is the number of meteorological variables. Figure 1b presents an example of a meteorological network when the number of neighborhoods (N) is two. The proposed approach assigns at least N candidate stations that can be correlated with the target station, assuming that we do not know the degree of correlation at this moment. Figure 1c shows the case with the minimum threshold (θ R ) for the PCC ( Figure A1). This approach allows stations to have a flexible-size neighborhood, but there can be isolated stations; models cannot learn the spatial correlations of those stations. Similarly, when N is too large or θ R is too small, we can miss the spatial correlations between stations. In the opposite case, the model is confused by overabundant information. In Section 4.4.2, we discuss the advantages and disadvantages of the two approaches by evaluating the proposed model based on N and θ R .
The node attributes at time point t (X t ) consist of solar irradiance, which is the forecasting target, and other meteorological variables related to solar irradiance. The solar irradiance is defined as follows: Definition 2 (Solar Irradiance). Solar irradiance at a station at time t is viewed as one of the node attributes of a meteorological network. Thus, X t,S r ∈ R V is a row vector of X t that represents the solar irradiance of all stations at time t. In addition, the ith component of X t,S r (S r (t, i)) corresponds to the solar irradiance degree at the ith station at time t.
The proposed model predicts future solar irradiance by analyzing previous solar irradiance and meteorological variables. The spatiotemporal correlations of meteorological variables with solar irradiance will enable the proposed model to understand weather contexts that can affect solar irradiance. Although we acquired the 16 variables listed in Table 1 from the ASOS data, the last column of the table says that not all the variables have explicit correlations with solar irradiance. Similar to the adjacency of observation stations, the PCC will be insufficient for establishing the spatiotemporal correlations of solar irradiance with meteorological variables. However, omitting highly correlated variables can hinder the proposed model from recognizing weather contexts, and appending extraneous variables can confuse the model. Thus, the meteorological variables can be defined based on the threshold for their PCC with solar irradiance (θ V ) as follows: Definition 3 (Meteorological Variables). The remaining node attributes are multiple variables that correlate with solar irradiance and reflect the weather context. When K = {k 1 , · · · , k K } is the set of all available meteorological variables and r(·, ·) indicates the PCC between two variables, the node attributes can be formulated as K * = {k j |k j ∈ K, r(k j , S r ) ≥ θ V }. In addition, similar to solar irradiance (S r ), when k j (t, i) refers to the value of k j at the ith station at time t, X t,k j indicates a vector representing the values of k j at time t at every observation station. By concatenating X t,k j , ∀k j ∈ K * , we can compose an attribute matrix at time t, X t .
In Section 4.4.1, we evaluate the proposed method to compose a set of meteorological variables by adjusting θ V . Therefore, solar irradiance forecasting, which models temporal and spatial dependencies between solar irradiance and meteorological variables, can be defined as learning a mapping function f based on the meteorological network N = A, X that consists of the static adjacency matrix A and the dynamic node attributes X . When L p and L o are the prediction and observation sequence lengths, respectively, the forecasting procedure can be formulated as: (1)

Multi-Attributed Spatio-Temporal Graph Convolutional Network
The proposed model aims to discover the spatio-temporal correlations of solar irradiance with multiple meteorological variables. The existing spatio-temporal GCN models [51,62,[64][65][66] have barely paid attention to dynamic changes in node attributes. Thus, we propose a novel spatio-temporal GCN model that can consider multiple dynamic node attributes by extending the AST-GCN [51], which considers deals with both static and dynamic attributes. The proposed model mainly consists of GCN layers and GRU layers. The GCN layers focused on extracting spatial features from snapshots of the meteorological network at each time point. Then, the GRU layers then analyze temporal changes in the spatial features for predicting to predict future solar irradiance. We call the proposed model 'MST-GCN (Multi-attributed Spatio-Temporal Graph Convolutional Network)', and the structure of the model are illustrated in Figure 2.

Multi-Attribute Fusion
This study represents multiple meteorological variables observed at each station as attributes of corresponding nodes to infer micro-and macro-weather conditions and their spatiotemporal correlations. Thus, the adjacency matrix A is static and represents only the geographical adjacency of the stations. The node attributes X = X 1 , · · · , X T are dynamic and expressed as a sequence of attribute matrices at each time point. The simplest approach for feeding N = A, X into spatiotemporal GCN models is using N t = A, X t as inputs to GCN layers and learn temporal changes in feature vectors for N t−L o to N t on the GRU layers. However, this approach overlooks the spatial influences of weather conditions on observation stations, which are not immediate. To analyze the spatial influences with agnostic and unfixed time lags, we let the GCN layers observe multiple time points using a fixed-length sliding window. In addition, the GRU layers compress the input features at multiple time points into fixed-length vectors, which causes information loss. This approach can reduce the risk of information loss by sharing the burden of temporal analysis with the GCN layers. When the length of the sliding window is l, the input network of the GCN layer at time t can be formulated as: where X t t−l+1 ∈ R lK×V indicates the concatenation of attribute matrices within the window. Although Zhu et al. [51] suggested this approach, they did not evaluate its effectiveness in analyzing multiple dynamic attributes. Section 4.2 focuses on verifying the effectiveness of enabling the GCN layers to conduct a spatiotemporal analysis by comparing it with T-GCN [62]. We heuristically set the window size l as six. The proposed model aims for day-ahead hourly forecasting, and the 42 ASOS stations used in this study are densely located in the southern part of the Korean Peninsula. Thus, we assumed that l = 6 is a sufficient time period for establishing the inter-station spatial influences required for predicting the weather tomorrow.

Spatial Dependency Modeling
Discovering the spatial influences between the weather contexts of observation stations is significant for predicting future weather contexts and forecasting solar irradiance. Graph convolutional network (GCN) models, which are the generalization of convolutional neural network (CNN) models to graph-structured data, have been shown to be effective for analyzing the propagation of node features between adjacent nodes. The proposed model employs the spectral graph convolution method proposed by Kipf and Welling [67], which improves the computational complexity of the existing spectral GCN models. This method updates node features to smooth (and denoise) the features of neighboring nodes by conducting convolution operations in the spectral domains. Thus, convolution filters extract spatial features between nodes by analyzing them and their first-order neighborhoods. By stacking multiple GCN layers, we can obtain the representation of each node while considering the influence of adjacent nodes. This study initially sets the node features as meteorological variables, including solar irradiance, from t − l + 1 to t. Then, the GCN layers generate representations of weather contexts at each station at time t based on spatial influences between the weather contexts of adjacent stations during the time period [t − l + 1, t]. The GCN layer in the proposed model can be formulated as: is the feature matrix generated by the nth GCN layer at time t. H (0) t represents the initial node feature at time t. θ (n) refers to the convolution filter of the nth layer. Finally, σ(·) denotes the activation function used for nonlinear modeling. Therefore, each GCN layer linearly transforms the feature matrix (H to the GCN layers, it cannot allow the model to determine the temporal order of the variables. However, the model can learn the spatial correlations among variables at each time point in the windows. The GRU layers then discover more distinctive spatiotemporal correlations from temporal changes in the sequence of feature As discussed in the previous section, the meteorological network had 42 nodes (stations), and the out-degrees of the nodes were at least N. The edge density of the meteorological network is far higher than that of conventional networked data such as bibliographic networks. Because aggregating information of a few hops will reach across the network (or its sub-networks), simply stacking a few GCN layers can cause an over-smoothing problem. Thus, we heuristically set the number of GCN layers to 2. In addition, we used a rectified linear unit (ReLU) as the activation function.

Temporal Dependency Modeling
The node representations extracted by the GCN layers reflect the spatiotemporal correlations between the meteorological variables. However, extending the window size (l) as long as the observation period (L o ) provides overabundant information that exceeds the learning capabilities of the GCN layers. Thus, time-varying representations with shortterm spatiotemporal features are fed into GRU layers to establish temporal dependencies between meteorological contexts in the long term. The GRU layers can be regarded as the compositions of the reset and update gates. When r t and u t are the reset gate and update gate at time t, respectively, r t is used to control the amount of information of the previous time points (i.e., the output at the previous time point, h t−1 ) that will be remembered or forgotten. Likewise, u t controls how much past information should be reserved. r t combines the input feature vector on t with h t−1 to compose the current cell state, c t . The final output on t (h t ) is then derived by combining c t and h t−1 using u t . This procedure can be formulated as follows: where σ(·) denotes a sigmoid function. φ(·) denotes the hyperbolic tangent function. * refers to the Hadamard product.
is a spatial feature matrix extracted from N t using the last GCN layer. θ u , θ r , and θ c are learnable weight matrices. We stacked two GRU layers above the GCN layers, and one fully connected layer was used to predict the future solar irradiance from h t with linear activation.
The goal of solar irradiance forecasting is to make the prediction result approximate the actual weather conditions as closely as possible. Thus, the objective of the proposed model was to minimize the prediction error. The error was measured by the L2 loss, and the objective function can be formulated as: where X τ,S r andX τ,S r are the actual and predicted solar irradiances at time τ, respectively, and λ is a hyperparameter that controls the regularization rate.

Evaluation
This section presents the experimental procedures and results for evaluating the prediction performance of the proposed model and validating the research questions underlying the proposed approaches. First, we compared the performance of the proposed model with baseline models, including both conventional regression models (e.g., HA, ARIMA, VAR, and SVR) and neural network models (e.g., MLP, GCN, GRU, and T-GCN). We also examined the performance of the proposed and existing models in terms of longterm predictions. This experiment demonstrates the practicality of the proposed model and shows whether the models understand the dynamic changes in weather contexts. Subsequently, we examined the stability of the forecasting models by comparing their performance variations according to cloudiness and months. Finally, the proposed model has several hyperparameters that determine the meteorological variables and neighboring stations that were used for forecasting. We evaluated the sensitivity of the proposed model by assessing its performance according to the hyperparameters.

Experimental Settings
This section describes the experimental settings, including the datasets, accuracy metrics, hyperparameter settings, and the comparison groups. We acquired meteorological observation data from 42 ASOS stations for four years (1 January 2017 to 31 December 2020), as described in Section 2.1. We collected hourly solar irradiance data during the observation period and the 16 meteorological variables correlated with solar irradiance, as listed in Table 1. The observation data of the first three years (1 January 2017 to 31 December 2019) were used as the training dataset. We then evaluated the proposed model using the remaining years (1 January 2020 to 31 December 2020). For the experiments, every meteorological variable was normalized to [0, 1]. We predicted the solar irradiance of the next L p ∈ {1, 2, 3, 4, 5, 6, 12, 24} h at each ASOS station by analyzing the solar irradiance and meteorological variables during the previous L o ∈ {12, 24} h.
In this study, we used six accuracy metrics to evaluate the performance of solar irradiance forecasting models: root mean square error (E 2 ), mean absolute error (E 1 ), normalized mean square error (NE 2 ), accuracy (A), R-squared (R 2 ), and variance (σ). These metrics can be measured as follows: where Y t andŶ t indicate the observed and predicted solar irradiances at time t, respectively.  [67] is a graph neural network, which extracts spatial features from networked data by transforming and aggregating feature vectors of nodes and their neighborhoods. We applied a two-layered GCN model to predict X t+1,S r , · · · , X t+L p ,S r by analyzing A and X t−L o +1,S r , · · · , X t,S r without considering the temporal order. • GRU [72] is a RNN model that improves the long-term dependency problem of the conventional RNN and reduces parameters of LSTM. We used a three-layered GRU model with 64 hidden units to predict S r (t + 1, i), · · · , S r (t + L p , i) at each station v i by analyzing only temporal changes in S r (t − L o + 1, i), · · · , S r (t, i) . • T-GCN [62] is a spatio-temporal graph neural network, which is a combination of GCN and GRU, in order to analyze networks with static structures and dynamic attributes.
This model uses a GCN layer to extract spatial features (H (N) τ ) from A and X τ,S r on each time point τ ∈ [t − L o + 1, t] and a GRU layer to predict X t+1,S r , · · · , X t+L p ,S r by analyzing temporal changes from H The proposed model was implemented using TensorFlow in Python. We used a hyperbolic tangent function as the activation function of the output layer, ReLu function for the hidden layers, and Adam optimizer [73]. We conducted a grid search for the proposed model's hyperparameters: number of epochs: 500 to 3000 with a step size of +500, learning rate: 0.0001 to 0.01 with a step size of ×10, batch size of 32 to 512 with a step size of ×2, number of GRU hidden units: 8 to 128 with a step size of ×2, and number of neighboring observation sites: 1 to 9 with a step size of +1. The proposed model had the best accuracy for the number of epochs = 1000, learning rate = 0.001, batch size = 128, number of GRU hidden units = 64, and number of neighboring observation sites = 2. The implementation of the proposed model is available in the GitHub repository (https://github.com/higd963/MST-GCN (accessed on 15 August 2022)).

Effectiveness of the Proposed Model
This section evaluates the proposed model by comparing its accuracy with that of the baseline models. The models were trained to predict the solar irradiance at time t by analyzing the previous solar irradiance from t − 1 to t − 12. We also compared the accuracy of the models on multivariate analysis with univariate analysis to demonstrate that the proposed model is more effective for analyzing correlations between multiple meteorological variables than existing models. Additionally, HA and ARIMA cannot deal with multivariate features, and the proposed model and VAR are not designed for univariate analysis. Thus, we could not assess model accuracy for those cases and remained as blanks on Table 2, which lists the experimental results.
The proposed method outperformed the existing models in every evaluation metric. In addition, the existing models exhibited a significant performance decrement in the multivariate analysis compared to the univariate analysis. This result is unexpected because T-GCN [62] and the proposed model do not have significant differences in their model structures. When we modified T-GCN to consider multiple meteorological variables, the major difference between the two models is that the proposed model observes samples from time t − m to t with an m + 1-length window on time t, but T-GCN only considers samples at t. Thus, the graph convolution layers in the proposed model can extract spatiotemporal features from meteorological data. Otherwise, the graph convolution layers and recurrent layers in the T-GCN focus only on spatial and temporal features, respectively. Thus, we can assume that spatiotemporal analysis is more effective for discovering correlations between meteorological variables than aggregating spatial features using recurrent layers (RQ1 and RQ2).
The neural network models with temporal features (e.g., T-GCN and GRU) outperformed the other models in univariate analysis. However, in the multivariate case, GRU exhibited a worse performance than GCN. Although the recurrent layers could be effective for discovering daily patterns of sunshine, stacking the recurrent layers was not sufficient to establish and utilize the correlations between meteorological variables. This point was also shown in that T-GCN underperformed GRU in the univariate case, which was the opposite in the multivariate case. All existing models exhibited significantly worse performance on multivariate analysis than on univariate analysis. This result might be caused by limitations in the learning capabilities of the models, the same as with the GRU.
The deep learning-empowered models significantly outperformed the conventional regression models in both the univariate and multivariate cases, excluding SVR. SVR had the best performance in univariate analysis considering E 2 and NE 2 . On the other metrics, SVR showed the second-best performance, with a slight difference. This is an interesting result in that SVR can comprehend temporal changes in solar irradiance as much as GRU and T-GCN. However, it exhibited a catastrophic performance decrease in the multivariate case. Because the models conducted one-hour-ahead prediction by analyzing previous solar irradiance during 12 h, SVR could be sufficient for applying daily patterns to the meteorological contexts of each day. Meteorological variables have complicated spatiotemporal correlations, requiring forecasting models with highly expressive capabilities, which conventional models cannot support.
As there have been studies on hourly day-ahead forecasting of solar irradiance [74], we examined the practicality of the proposed model by comparing the accuracy of the proposed and existing models according to various prediction and observation periods. Using two observation periods, L o = 12 and 24, we adjusted the prediction periods (L p ) in [1,6] and {1, 3, 6, 12, 24}, respectively. In the above experiment, the conventional regression models significantly underperformed the deep learning-empowered models, and none of the existing models properly conducted multivariate analysis. Thus, the following experiments employed only deep neural networks (e.g., MLP, GCN, GRU, and T-GCN) trained for univariate analysis as baseline methods. Table 3 and Figure A3 show the results of this experiment.
The proposed model exhibited the highest accuracy for most cases and metrics. The performance improvement was more noticeable in the long-term prediction than in the short-term prediction because the proposed model showed consistently high accuracy according to L p in both L o = 12 and = 24 cases. However, in the short-term prediction, GRU performed similarly to the proposed model, even higher for L o = 24 and L p = {1, 3}. As discussed in the previous experiment, the GRU can be effective for discovering periodic patterns of solar irradiance. In short-term prediction (far shorter than a day), observing the discovered patterns would be more significant than recognizing weather contexts by understanding spatiotemporal correlations between the meteorological variables (RQ2). This point would be the opposite in long-term prediction.
Similarly, in both L o = 12 and = 24 cases, T-GCN performed better than GRU in long-term prediction and worse than GRU in short-term prediction. In the L o = 12 case, reversal occurred at L p = 3. However, in the L o = 24 case, the reversal occurred at L p = 24, which is much later than in the L o = 12 case. There are two reasons for this result. First, the GRU can learn daily solar irradiance patterns more deeply at L o = 24 than at L o = 12. In addition, the T-GCN would not be able to recognize weather context information as accurately as the proposed model because it is a univariate model and cannot analyze correlations between the various meteorological variables. Furthermore, a comparison of the GCN with GRU showed a similar result to the above comparisons only in the L o = 12 case from L p = 4. However, GRU had a higher σ than GCN in most cases, excluding when L o = 12 and L p = 5. In addition, the GRU significantly outperformed the GCN in all cases and metrics in the L o = 24 case. The GCN and GRU focus only on spatial correlations and temporal changes in solar irradiance, respectively. The results of the three comparisons indicate that weather contexts discovered from the spatial features are effective for long-term prediction (RQ1). Simultaneously, the spatial features could not show their effectiveness without combining with the temporal features (RQ2), and the combination of spatiotemporal analysis and multivariate analysis could enhance model capabilities for understanding weather contexts (RQ3). Additionally, the performance decrement of the GRU would not come from the long-term dependency problem because the proposed model and T-GCN also learn temporal changes in spatial features extracted by GCN layers using GRU layers.
MLP significantly underperformed the other models. Although MLP exhibited consistent performance according to changes in L p , its prediction results are difficult to be used for forecasting systems, considering that the observed solar irradiance values were normalized into [0, 1]. Although MLP analyzes correlations between the meteorological variables, the multivariate analysis without spatial and temporal features could not effectively recognize weather contexts.
From these experimental results, we can discover that (i) spatial correlations between observation sites are essential for consistent forecasting performance on both long-term and short-term prediction (RQ1), (ii) in short-term prediction, periodic patterns are more effective than the other features (RQ2), (iii) spatial correlations show their worth when used with the periodic patterns (RQ1 and RQ2), and (iv) correlations between multivariate variables could not show high accuracy solely but exhibited its effectiveness when used with the others (RQ3). Table 2. A performance comparison of the proposed model with the baseline methods. The upper, middle, and lower parts of this table present the accuracy of the conventional regression models, the deep learning-empowered models, and the proposed model, respectively. Furthermore, the left and right sides exhibit the model accuracy of the univariate and multivariate analysis, respectively. -indicates cases that we cannot assess the models' accuracy due to their characteristics. Additionally, * and ** denote models with the first and second best performance on each case and metric, respectively. The hyper-parameters were set as L o = 12, L p = 1, N = 4, and θ V = 0.00.

Univariate Analysis
Multivariate Analysis

Stability of the Proposed Model
This section presents the performance stability of the proposed model by comparing its accuracy fluctuation according to weather conditions with those of the baseline models (e.g., GCN, GRU, and T-GCN). Solar irradiance is affected by various weather factors, such as cloudiness, and seasons are correlated with the annual patterns of solar irradiance and weather. Therefore, we first examined the forecasting models' performance at every cloudiness level as a representative factor affecting the solar irradiance. The monthly performance of the models was then evaluated for determining the seasonal influence on solar irradiance and the forecasting models.

Performance Variation according to Cloudiness
Solar irradiance showed relatively consistent patterns on clear days, and sunny days were more frequent than cloudy days. Therefore, we must evaluate whether the proposed model can achieve high accuracy regardless of cloudiness for examining the practicality of the model. We classified cloudiness into 10 degrees, and our data samples were segmented according to the degree of cloudiness. Subsequently, we evaluated the performance of the proposed and existing deep-learning-empowered models within each segment of the dataset. Table 4 and Figure A4 list the experimental results. Table 4. Performance of the proposed and existing models according to cloud cover. The first row indicates the degree of cloud cover from 0 to 10. The second row presents the distribution of cloudiness levels in the experimental dataset. In the third to thirteenth columns, * and ** indicate cloudiness levels that each forecasting model had the best and second-best performances, respectively. Additionally, in the last two columns, * and ** refer to forecasting models that showed the best and second-best performances on average and the lowest and second-lowest standard deviation, respectively. The hyper-parameters were set as L o = 24, L p = 24, N = 4, and θ V = 0.00. The proposed model exhibited the highest accuracy for all cloudiness levels. A performance decrement on cloudy days was commonly observed in all models. However, the decrease in the proposed model was not as severe as that of T-GCN and GRU. Although on a few metrics, the GCN had a similar or lower standard deviation compared to the proposed model, there was a significant difference between the accuracies of the two models. As discussed, the solar irradiance on clear days follows periodic patterns (e.g., daily and yearly). Cloudy days were far less frequent than clear days, as shown in Figure A2. Thus, the high performance of the proposed model supports that the model could predict cloudiness levels the next day by overcoming data imbalance and not merely applying periodic patterns to observations.

Cloud Cover
In the previous experiment, the T-GCN outperformed the GRU for long-term prediction, whereas the opposite was true for short-term prediction. Similar relationships were observed in this study. Although T-GCN outperformed GRU on clear and slightly cloudy days, GRU performed better than T-GCN on extremely cloudy days (CC ≥ 8) in most metrics. Even in the CC = 10 case, the performance of the T-GCN was similar to or worse than that of the GCN, which focused only on spatial features. For the previous results, we assumed that the combination of spatial and temporal features was effective for understanding meteorological contexts in the long term (RQ1 and RQ2). However, in this experiment, the same combination hindered the recognition of the overcast days. Compared to the proposed model, we can assume that considering solar irradiance in adjacent areas could not provide sufficiently deep contexts to the models, which is different from analyzing spatiotemporal correlations between various meteorological variables (RQ3).
Because T-GCN and the proposed model have almost the same architecture, this result would not be achieved from model capabilities for handling the imbalanced distribution of cloudiness. Additionally, the forecasting models commonly showed a rapid performance decrement from CC = 7, although they exhibited a relatively stable performance on 2 ≤ CC ≤ 6. Partly cloudy skies may not significantly influence the solar irradiance.

Performance Variation according to Months
The weather on the Korean Peninsula, which is our experimental subject, has four distinct seasons. According to seasonal changes, the weather in each month might have distinctive patterns. We can examine whether the yearly patterns affect the solar irradiance prediction by assessing the forecasting monthly model performance. In addition, the monthly performance can establish the model that can learn yearly patterns or overcome seasonal differences. As in the previous experiment, we segmented our observation samples into months, and the proposed and existing forecasting models were evaluated for each month. Table 5 and Figure A5 list the experimental results.
The proposed model outperformed existing models in most months and metrics. T-GCN and GRU exhibit lower E 1 values than the proposed model in January and October. The proposed model had a thin lead for the other metrics; we can also see this result from October to December. Considering the previous experiments, T-GCN and GRU exhibited significant performance decrement on cloudy days (Section 4.3.1) and long-term predictions (Section 4.2). Thus, this result could be due to the weather in the Korean Peninsula during winter. (The climate of the Korean Peninsula can be explained as a humid continental climate with humid summers and dry winters. The temperature differences between the hottest part of the summer and the coldest part of the winter are extreme. In addition, most precipitation falls during the summer monsoon period between June and September.) Periodic meteorological patterns were more effective for solar irradiance prediction during winter than the other features. We can also see this point from the result for July, which occurs in the middle of the 'humid summer' on the Korean Peninsula. In July, the proposed model had a relatively thin lead for E 1 but significantly outperformed the other temporal models for normalized accuracy metrics (e.g., A, R 2 , and σ). The similar E 1 values could be caused by the fact that rainy days have low solar irradiance. Thus, although T-GCN and GRU could infer rainy days by analyzing the solar irradiance of adjacent samples, these models failed to predict the solar irradiance on rainy days. Both spatiotemporal and temporal analyses of solar irradiance are insufficient in forecasting nonperiodic and complicated meteorological phenomena. In July, spatiotemporal analysis (T-GCN) underperformed cases in which spatial and temporal features were solely used (GCN and GRU). As discussed for the experiment on cloudiness, considering spatial correlations of solar irradiance only could not make the models understand the meteorological contexts sufficiently deep. However, the proposed model improved this point by analyzing the correlations between the meteorological variables (RQ3).
The T-GCN, GRU, and proposed model exhibited similar tendencies. These models exhibited high normalized accuracy metrics (e.g., NE 2 , A, R 2 , and σ) from February to April and October to December, and low accuracy from January and July to September. However, contrary to these models, the performance of the GCN was worse than its average in spring (March to April). Because the spring climate of Korea is dry and has clear characteristics (https://web.kma.go.kr/eng/biz/climate_01.jsp (accessed on 15 August 2022)), this differ-ence is not due to irregular meteorological phenomena, such as precipitation. This result indicates that temporal features are effective for discovering yearly climate patterns. Table 5. Performance of the proposed and existing models according to months. In the third to fourteenth columns, * and ** indicate months when each forecasting model had the best and secondbest performances, respectively. Additionally, in the last two columns, * and ** denote forecasting models that showed the best and second-best performances on average and the lowest and secondlowest standard deviation, respectively. The hyper-parameters were set as L o = 24, L p = 24, N = 4, and θ V = 0.00.

Parameter Sensitivity Analysis
We assumed that meteorological parameters observed in spatially adjacent areas could influence each other's future meteorological parameters. For example, the wind speed and direction are affected by the atmospheric pressures of adjacent areas. Therefore, we conducted a temporal analysis of meteorological variables in adjacent areas using the spatiotemporal GCN model. However, there are problems in determining (i) spatially adjacent areas and (ii) correlated meteorological parameters. This section evaluates the effectiveness of the proposed methods for defining spatial adjacency and composing a set of input variables. In addition, we assessed the sensitivity of the proposed model to changes in these two factors.

Meteorological Variable Compositions
The proposed model significantly outperformed the T-GCN [62] by analyzing the spatial correlations between various meteorological variables, and solar irradiance. The previous experiments used all the variables that we collected (listed in Table 1). However, variables that are not correlated with the solar irradiance can hinder the forecasting performance of the proposed model. In addition, if the proposed model can exhibit similar or better performance with fewer input variables than with every variable, the practicality of the model will be improved. Table 1 presents the Pearson correlation coefficients of solar irradiance with other meteorological variables. We composed the input variable group using a threshold for the correlation coefficient based on the assumption that variables with higher correlations contribute more to forecasting performance. Subsequently, we validate this assumption by adjusting the threshold, as shown in Table 6. Table 6. Performance of the proposed model according to the variable composition. The first row presents the thresholds of Pearson correlation coefficients of the forecasting target (i.e., solar irradiance) with the other meteorological variables (θ V ). The second row shows the number of variables chosen based on the thresholds. Additionally, * and ** indicate that the proposed model had the best and second-best performances with the corresponding θ V , respectively. The remaining hyper-parameters were set as L o = 12, L p = 1, and N = 4. We assume that not all meteorological variables contribute to the forecasting performance of the proposed model. Variables that are less correlated with solar irradiance provide unnecessary and overabundant information for the forecasting model. In addition, if we choose variables that are too strict (i.e., small θ V ), the model cannot obtain enough information for understanding weather conditions. Therefore, we expect the model performance to exhibit a convex shape for θ V . Nevertheless, unexpectedly, the proposed model had the best performance when we used all 15 variables (θ V = 0.00), with a significant gap. After removing year (YOY) and day (DOY) from 15 (θ V = 0.04), the model performance showed a sharp decrease. Subsequently, when we removed precipitation (P t ) and local pressure (P L ) from the remaining 13 variables (θ V = 0.08), the model exhibited the second highest performance. In addition, this case's performance was similar to that of the θ V = 0.50 case, which used only sunshine duration (S) and solar irradiance itself (S r ). It is not easy to expect the proposed model to understand meteorological phenomena and their spatiotemporal correlations from the two variables. Accordingly, year and day could be key factors for discovering correlations between the variables, despite their low PCC. Therefore, we concluded that the meteorological variables might have non-linear correlations with solar irradiance, and PCC was not sufficient to reflect these correlations. Likewise, the meteorological variables, excluding sunshine duration, air temperature, and relative humidity, have very low PCC with solar irradiance, as shown in Table 1. However, it is difficult to empirically determine the optimal composition of the meteorological variables. Future research will focus on feature selection methods for meteorological variables by combining statistical correlations and domain knowledge in meteorology studies.

The Number of Neighboring Stations
By comparing the ASOS station locations (Figure 1) with the correlation between their historical solar irradiance ( Figure A1), we can see that the stations have higher correlations with the closer stations, but the clusters of correlated stations have various sizes. In addition, every observation station had a high PCC (>0. 75), and this result indicates that solar irradiance at the ASOS stations has similar long-term tendencies. However, we assume that it is difficult to reflect short-term (hourly or daily) changes and differences in the solar irradiance. Thus, this study employed fixed-size distance-based nearest neighborhoods, which is different from existing studies [43,63] that used flexible-size correlation-based neighborhoods. To examine the advantages and disadvantages of both methods, we assessed the performance fluctuations of the proposed model by adjusting the meteorological networks using both static and dynamic neighborhood sizes, as shown in Table 7. Table 7. Performance of the proposed model according to the edge density of meteorological networks. The edge density was adjusted using the number of neighborhood stations (N) and the threshold for the Pearson correlation coefficient between historical solar irradiance of stations (θ R ). # Nodes and # Edges indicate the number of nodes and edges, respectively, in the meteorological networks constructed according to N and θ R . In addition, * and ** denote that the proposed model had the best and second-best performances with the corresponding N and θ R , respectively. The remaining hyper-parameters were set as L o = 12, L p = 1, and θ V = 0.00.

# Neighborhoods (N)
Corr. When we fixed the number of neighborhoods (N), the proposed model exhibited the best performance for N = 2. Then, the performance worsened according to the increment in neighborhood sizes, and suddenly, the model exhibited the second-highest performance at N = 6. After N = 6, the model performance deteriorated again with the neighborhood extension. The proposed model assumes that the weather conditions of neighboring stations influence each other. This result indicates that closer stations do not always have a greater influence on future weather conditions. Because observation stations have different geographical features and cannot be located at equal distances, we should search for the optimal number of neighborhoods according to the compositions of the observation stations. In addition, the number of influential stations is not fixed. Thus, we assessed the performance of the proposed model in cases in which we defined the adjacency of the stations according to correlations of their solar irradiance history.
As shown in Figure A1, the observation stations exhibit high PCC values. This could mean that the solar irradiance at most stations exhibited similar tendencies. However, to handle the output instability of solar power plants, regional differences in solar irradiance should be accurately forecast. Although the threshold-based approach could not outperform the nearest neighborhoods, the θ R = 0.95 case had the second-highest performance among all cases. However, the θ R = 0.93 case underperformed at N = 6, and the θ R = 0.94 case exhibited a similar performance to the N = 9, which is the worst. Considering the values of θ R in these cases, the performance of the proposed model is sensitive to θ R . Furthermore, if PCC can reflect the meteorological influence between observation points, model performance will show consistent tendencies according to θ R . When θ R is lower than a certain value, overabundant information can be provided to the model. In addition, with too high θ R , the model cannot obtain sufficient information to analyze the spatial correlations of meteorological features. However, unexpectedly, the performance of the proposed model had an irregular tendency according to θ R , and among the three cases (θ R = 0.93, 0.94, and 0.95), θ R = 0.94 exhibited a significantly worse performance than the others. This result contradicts the observations above.
Both the distance-based and correlation-based approaches exhibited irregular tendencies. In addition, although the distance-based approach outperformed the correlation-based approach, the difference was not significant. In conclusion, neither approach was sufficient in reflecting the spatial correlations and meteorological influences between the observation areas. Future research should focus on developing measurements of spatial correlations.

Conclusions
This study aims to conduct day-ahead hourly forecasting of solar irradiance by analyzing the spatio-temporal correlations of solar irradiance with multiple meteorological variables. We also evaluated the effectiveness of (i) spatial analysis, (ii) temporal analysis, and (iii) multivariate analysis for solar irradiance forecasting and validated the underlying research questions presented in Section 1. We collected solar irradiance and other meteorological variables (Table 1) from 42 ASOS stations on the Korean Peninsula ( Figure 1). For spatiotemporal analysis of the variables, we modeled the ASOS observation data as a dynamic attribute network, which has the stations as nodes, variables as attributes, and spatial adjacency between the stations as edges. We then developed a novel solar irradiance forecasting model that analyzes the dynamic attributed network and predicts hourly solar irradiance at each station by modifying the AST-GCN model [51].
We evaluated the effectiveness of the proposed model by comparing its prediction accuracy with those of existing deep learning-empowered models and conventional regression models. Subsequently, to validate the practicality of the proposed model, we examined its accuracy according to the prediction sequence lengths (from hour-ahead to day-ahead prediction), cloudiness, months, variable compositions, and edge density of the network. The proposed model outperformed the existing models, especially in terms of long-term prediction. Contrarily, most of the existing studies have been limited in intra-day prediction (1 to 6 h ahead) [18]. The comparison between the performances of the proposed and existing models indicates that the spatial, temporal, and multivariate features of atmospheric data are synergistic for predicting solar irradiance. Although a few previous studies [43,46,50] attempted to combine temporal and spatial features or used both temporal and multivariate features, there have barely been either forecasting models integrating the three factors or validations for synergistic effects among the factors. In addition, the proposed model exhibited higher and more stable performances on most cloudiness levels and months than the existing models. The proposed model exhibited performance decrements on overcast days and summer as with the existing ones. However, more or less, solar irradiance forecasting models are difficult to avoid this problem caused by frequent thunderstorms in the summer of the Korean Peninsula [75,76]. The experiment for variable compositions showed that the correlation coefficients were insufficient in reflecting spatiotemporal correlations between meteorological variables. Likewise, both geographic distances and correlation coefficients were insufficient in establishing spatial influences between the atmospheric contexts of the observatories. In future research, we will address the following limitations: • Prediction sequence length: We evaluated the forecasting performance of the proposed and existing models on multiple prediction sequence lengths (from an hour-ahead to a day-ahead prediction). However, predicting solar irradiance with longer time intervals (e.g., a week or a month) will be helpful for the practical usage of solar power. We assume that the long-term dependency problem caused by adopting GRU layers hindered the long-term prediction performance of the proposed model. In further research, we will improve this problem by applying the attention mechanism to consider relative importance of time points, adjacent stations, and meteorological variables. • Low accuracy on high cloud cover: The proposed model showed performance decrement on cloudy days, although the decrement was not as significant as the existing models. This problem might come from difficulties in predicting solar irradiance on cloudy days but also due to forecasting cloudiness. Wind speeds and directions at high altitudes are closely correlated with cloudiness [77], and future research will attempt to consider these variables in addition to those observed at the ground observatories. • Multi-modal analysis: Atmospheric observation data are collected through various devices (e.g., sensors, radars, cameras, etc.) deployed on ground stations, satellites, observation balloons, aircraft, etc. Despite the variety of observation data, this study has focused on sensor data from ground observatories. Combining the multi-modal and multi-aspect observations will enable forecasting models to discover more accurate information for atmospheric contexts. For example, the ground observatories were not located with a uniform gap, and geographical characteristics in the gaps were also not homogeneous. Thus, covering the gaps by incorporating geographical features [38,78], land usages [79], and satellite data [80] will be effective for analyzing spatial correlations between atmospheric contexts of the observatories.

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