Mapping Solar Global Radiation and Beam Radiation in Taiwan

Abstract

Data for solar radiation resources play a pivotal role in assessing the energy yield capability of solar applications. A nationwide database for the typical meteorological year from the 30 weather stations of the Central Weather Bureau (CWB) in Taiwan is used to determine the spatial distribution of global radiation over the terrain of Taiwan. There is no available beam radiation information in daily reports from all CWB stations. Information on the diffuse fraction for all CWB stations is estimated using three available correlation models that account for topographical and geographical effects in Taiwan. The databases for beam radiation are generated using these estimated diffuse fractions. The mappings of global and beam radiation on the Taiwanese mainland are performed with databases from 24 CWB stations using the residual kriging method. There are no mappings of the remote islands, where six CWB stations are located. The databases for global and beam radiation for these six CWB stations are applied to nearby remote islands. The effects of topography and geography on the distributions of global and beam radiation are discussed. The spatial distributions of solar radiation presented are good scientific references for assessing the performances of solar energy systems in Taiwan.

1. Introduction

Imported fossil fuels have long been Taiwan’s main energy source, but public support for renewable energy is growing. Additionally, renewable energy is essential for decarbonizing Taiwan’s heavily fossil-fuel-dependent energy system. The Taiwanese government actively promotes the development of renewable energy to achieve a goal of 20% of energy generated being renewable in 2025. The framework for a Sustainable Energy Policy was announced in 2008 [1]. The Renewable Energy Development Act was promulgated in 2010 and was amended on 1 May 2019. The Greenhouse Gas Emission Reduction and Management Act was enacted in 2015. These two acts encourage energy consumers in Taiwan to produce and use green energy instead of fossil fuel energy.

Taiwan’s key renewable energy resources are solar and wind energy. Two phases of the purchased-based subsidy program for solar water heaters (SWH) were implemented in the periods 1986–1991 and 2000–2017 for promoting the use of solar energy [2]. The cumulative installed SWH area reached a peak of 1882 ${\mathrm{k}\mathrm{m}}^2$ in 2019 [3]. However, more than 98% of SWH systems were installed in the residential sector [2]. To ensure sustainable development, large-scale SWHs will be the key elements for expansions of the market for solar thermal applications, particularly after the termination of the second-phase subsidy program in 2018. The promotion of large-scale SWHs in the commercial and industrial sectors relies upon how the Greenhouse Gas Emission Reduction and Management Act is enforced in this context.

The purchased-based subsidy program for photovoltaics (PV), which is the other key strategy for the application of solar energy, was also implemented by the Taiwanese government in the period 1992–2009. After 2010, the feed-in-tariff (FIT) program for PV replaced the purchased-based subsidy program when the Renewable Energy Development Act came into effect. The cumulative installed PV capacity until November 2023 was 13,622 MW [3], most of which was installed in the residential sector. In order to achieve the goal of 20% of energy generated being renewable by 2025, it is estimated that the installed PV capacity must be at least 20 GW [4], so the cumulative installed PV capacity must be more than a 47% expansion of the present quantity by the end of 2025. However, to achieve the long-term goal of producing and using green energy in Taiwan, large-scale PV systems that are land-based or water-based are key to the promotion of solar PV applications.

Recent years have seen significant advancements in harnessing solar energy. These advancements go beyond improvements in solar panel efficiency and energy storage technologies, encompassing innovations in analyzing and utilizing solar radiation data. Current and reliable solar radiation data are essential for energy assessments in large-scale solar projects in Taiwan. Global solar radiation ($I_{global}$) in the sky comprises beam radiation ($I_{beam}$) and diffuse radiation ($I_{diffuse}$):

$$I_{global}=I_{beam}+I_{diffuse}$$

(1)

Beam radiation data are essential for energy assessments in concentrating solar applications, which are sometimes used in large-scale systems. However, no beam and diffuse radiation information is provided in the daily reports from all weather stations of the Central Weather Bureau (CWB) of Taiwan. Recently, a nationwide database has been established for global solar radiation for the typical meteorological year (TMY), using data from 30 CWB weather stations across the terrain of Taiwan [5]. Of these, twenty-four stations are located on the Taiwanese mainland (Stations 1–24), while the remaining six (Stations 25–30) are situated on various remote islands, as shown in Figure 1 and

Table 1. Information for the 30 CWB weather stations in Taiwan.

Number	Name	Latitude (°N)	Longitude (°E)	Altitude (m)
1	Anbu	25.18	121.53	837.6
2	Zhuzihu	25.16	121.54	607.1
3	Tamsui	25.16	121.45	19.0
4	Keelung	25.13	121.74	26.7
5	Taipei	25.04	121.51	5.3
6	Banqiao	25.00	121.44	9.7
7	Hsinchu	24.83	121.01	26.9
8	Wuqi	24.26	121.52	31.7
9	Taichung	24.14	120.68	84.0
10	Sun Moon Lake	23.88	120.91	1017.5
11	Alishan	23.51	120.81	2413.4
12	Yushan	23.49	120.96	3844.8
13	Chiayi	23.50	120.43	26.9
14	Chigu	23.15	120.09	2.9
15	Yongkang	23.04	120.24	8.1
16	Tainan	22.99	120.20	40.8
17	Kaohsiung	22.57	120.32	2.3
18	Hengchun	22.00	120.75	22.3
19	Dawu	22.36	120.90	8.1
20	Taitung	22.75	121.15	9.0
21	Chenggong	23.10	121.37	33.5
22	Hualien	23.98	121.61	16.1
23	Su-Ao	24.60	121.86	24.9
24	Yilan	24.76	121.76	7.2
25	Pongjiayu	25.63	122.08	101.7
26	Lanyu	22.04	121.56	324.0
27	Penghu	23.56	119.56	10.7
28	Dongjidao	23.26	119.67	43.0
29	Matsu	26.17	119.92	97.8
30	Kinmen	24.41	118.29	47.9

. This development marks a significant step forward in providing the necessary data for accurate energy assessments. Ground measurements for global solar radiation in monotonous, low-lying landscapes are generally reliable within a few tens of kilometers. However, in areas near the sea, urban centers, or rapidly changing landscapes such as mountains, this reliability may be restricted to just a few kilometers. Figure 1 illustrates that mountains and hills cover nearly two-thirds of the Taiwanese mainland. Five mountain ranges, referred to as the Central Mountain Ranges, extend from the north-northeast to the south-southwest, dividing the mainland into eastern and western regions. The land slopes gently towards broad plains and basins adjacent to the Taiwan Strait in the west. The steep Central Mountain Ranges feature over 200 peaks exceeding 3000 m and descends sharply to the Pacific Ocean in the east. A study by Hsieh et al. [5] revealed a significant topographic effect on global radiation at weather stations between the eastern and western parts of Taiwan. Additionally, there is a notable geographic effect on global radiation between the mainland and the islands in the Pacific Ocean and Taiwan Strait. Given the large area of the Taiwanese mainland and the limited number of weather stations, understanding the spatial variability of global radiation data is crucial for accurate solar energy assessments.

Hsiao et al. [6] estimated downward solar irradiance (DSI) over Taiwan using 2006–2007 data from a geostationary multifunctional transport satellite (MTSAT). The DSI estimates utilized the MTSAT’s image data and 40-meter-resolution digital terrain data for Taiwan. These calculations had to account for two factors: (1) atmospheric effects from suspended aerosols and clouds and (2) the blocking and diffusion effects of terrain gradients at each pixel-point of the MTSAT data. These lead to the result that the accuracies for the satellite measurements were less than the accuracies for the ground measurements. Indeed, a comparison of the hourly integrated DSI results with the available CWB’s hourly global data by Hsiao et al. [6] showed that the correlation coefficients were at best slightly greater than 0.9 for the low-altitude weather stations but dropped to about 0.8 for the high-altitude weather stations. In addition, the DSI data covered two years (2006 and 2007) only, which was insufficient to represent long-term typical solar conditions in Taiwan.

Beam radiation data are vital for conducting energy assessments in concentrating solar applications. The direct measurement of beam radiation is conducted using a pyrheliometer, which relies on a collimated detector to capture solar radiation at a normal incidence within a small area around the sun, and it necessitates a sun-tracking device in the operation. In contrast, measuring diffuse radiation is less complex and involves two pyranometers: one measures global radiation, while the other, equipped with a shadow band, measures diffuse radiation. With these measurements, beam radiation can be calculated using Equation (1).

Unfortunately, daily reports from CWB weather stations do not provide diffuse radiation data. While global maps of solar diffuse radiation exist, they do not offer the detailed resolution required for small regions. This limitation has led to the development of correlation models to estimate diffuse radiation where direct ground measurements are unavailable. By utilizing meteorological data in Taiwan and applying regression techniques to develop correlation models, these models can estimate diffuse radiation for Taiwan, enabling the calculation of beam radiation through Equation (1) using the measured global radiation from all CWB weather stations.

Most correlation models are designed to estimate the solar diffuse fraction (d), as defined below, rather than solar diffuse radiation by using predictors such as global radiation and other meteorological factors:

$$d={\displaystyle \raisebox{1ex}{$I_{diffuse}$}\!\left/ \!\raisebox{-1ex}{$I_{global}$}\right.}$$

(2)

It is widely acknowledged that existing correlation models cannot be universally applied to all geographical regions and climatic conditions [7,8,9,10,11,12]. Chang et al. [13] conducted a series of monitoring studies for solar diffuse radiation on the eastern and western sides of the Taiwanese mainland and on a remote island (the Penghu archipelago), revealing significant seasonal variations in diffuse fractions at these locations due to topographical and geographical influences. These findings are consistent with the variations in global radiation observed at the same sites, as reported by Hsieh et al. [5]. Recently, Chang’s research team [11,12] developed three modified Boland–Ridley–Lauret-type correlation models [14]. These models utilized in situ data from a TMY for global and diffuse radiation, collected separately from the western and eastern regions of the Taiwanese mainland and a remote island in the Penghu archipelago. This approach addresses the impacts of topography and geography across Taiwan’s terrain. These three correlation models for the diffuse fraction are, thus, suitable for estimating beam radiation throughout Taiwan.

Kuo et al. [15] developed a Liu–Jordan-type correlation model [16] using the data for global and diffuse radiation for two years (2011 and 2012), which was measured in the western part of the Taiwanese mainland. They then generated datasets for the hourly diffuse fraction for 20 CWB stations on the Taiwanese mainland (except for Stations 3, 6, 14, and 16 in Figure 1) by individually inputting each hourly data point for the global radiation that was reported by these weather stations into their developed correlation model. Maps of the diffuse fraction for each month in a year were then created using a simple regression (spline interpolation) technique in terms of two geographical parameters: latitude and altitude. The spline interpolation method estimates the value at a specific point using data for neighboring stations to calculate the monthly diffuse fraction on the Taiwanese mainland. However, as reported in the recent study of Lin et al. [11], the multiple-predictor model has a better long-term performance than the Liu–Joran-type (single predictor) model that was developed by Kuo et al. [15] A simple interpolation method may be used to generate a climate-variable map using ground-based measurements, provided that the number of the measurement stations is large enough to ensure that the grid size for the measurement stations on the map provides a sufficient spatial resolution for the studied climate variable. Park and Park [17] compared four different spatial interpolation methods, including spline, distance weighting, simple kriging, and ordinary kriging [18], to construct the distribution for the average daily solar global radiation in South Korea. The mean area density per station for the study was 1282 ${\mathrm{k}\mathrm{m}}^2$ (a total area of 100,032 ${\mathrm{k}\mathrm{m}}^2$ divided by 78 stations), which is slightly less than the 1492 ${\mathrm{k}\mathrm{m}}^2$/station for Taiwan (see

Table 2. Summary of spatial distributions for previous studies using the kriging family.

Country or District	Area/Mean Area Density per Station (km2)	Grid Size (km × km)	(Number of Stations for Model Development)/(Number of Stations for Validation)	Data Sampled	Type of Kriging
Saudi Arabia [19]	2,150,000/52,439	55 × 33	40/1 *	Since 1971 (incomplete information)	Ordinary kriging [18]
Turkey [20]	783,356/4927	0.5 × 0.5	124/35	Mean values in 1968–2004	Universal kriging [21]
Andalusia, Spain [22,23]	87,000/482	1 × 1	112/54	Mean values in 2003–2006	Ordinary kriging and residual kriging [24]
South Korea [25]	100,032/1235	0.03 × 0.03	80/1 *	Mean values in 2001–2012	Residual kriging
Greece [26]	131,957/3384	No information	39/39 **	TMY (1985–1999)	Empirical Bayesian kriging [27]
Taiwan (present study)	35,808/1492	1 × 1	23/1 *	TMY (2004–2018)	Residual kriging

*: cross validation, a leave-one-out estimate of each observation; **: a comparison of each measured data point with the estimated value at each location using the regression model that is developed using all data locations.

). The comparison shows that ordinary kriging is the most suitable spatial interpolation method. In addition, the effect of topography on the diffuse fraction is mainly due to the high-rise mountain ranges that run from north-northeast to south-southwest on the Taiwanese mainland, as observed by the study of our research team [12,13], so the geographical parameter of longitude must be considered in the mapping formulation of the diffuse fraction on the Taiwanese mainland.

Two spatial interpolation methods are found in the literature to determine spatially continuous databases using measurements from discrete stations. The first method involves deterministic interpolation techniques, such as spline functions or weighted averages, which calculate unknown values using mathematical functions based on the similarity to known points but do not provide reliability assessment errors for the predictions. The second method includes stochastic approaches, like kriging methods, which combine analytical and statistical techniques to predict unknown values based on the spatial correlation among known data points. Kriging methods are particularly popular for applications in solar radiation interpolation (see Table 2). Kriging methods have demonstrated significant advantages over deterministic interpolation techniques in estimating the spatial distribution of global solar radiation [17].

Ordinary kriging can produce reliable estimates of solar radiation for uniform terrains with consistent climatic conditions. However, its accuracy diminishes as topographical complexity increases or when the Earth’s surface is heterogeneous, such as in areas where the land meets the sea. An early study of Becker and Boyd [28] showed that there is no significant increase in solar global radiation with an altitude up to an elevation of 305 m (1000 ft) above sea level. However, there is an approximately quadratic percentage increase in solar global radiation with an altitude between the elevations of 610 m (2000 ft) and 3050 m (10,000 ft) above sea level in both the summer and winter seasons. The complex topography of the Taiwanese mainland creates difficulties for the kriging procedure. Discontinuity in elevation and shadowing due to topographical features creates significant local gradients in solar radiation for which the ordinary kriging interpolation procedure cannot properly account. Many different kriging processes (named the kriging family hereafter) have been developed for different applications in spatial interpolation, as summarized in Table 2. There are several approaches to incorporating external variables into the kriging process to account for topographic effects. One approach involves treating external variables, such as altitude, as covariables during interpolation [29]. However, this method becomes complicated when multiple external variables are involved [30]. Alternatively, external variables can be addressed in a preliminary step before the interpolation. Such a technique is known as residual kriging [24]. Alsamamra et al. [22] found that residual kriging was particularly effective for mapping complex terrains, which is why it was chosen for this study.

The solar radiation climate at a specific location reflects the levels and trends of global, diffuse, and beam radiation components over an extended period, typically at least ten years. Due to significant annual variation in solar radiation data, a typical year of solar radiation is used to simplify simulations and data management, rather than using multiple years of data. This study uses data from a TMY, a concept introduced by Hall et al. [31]. A TMY is constructed from twelve typical meteorological months (TMMs) selected from different years and compiled into a single representative year. This compilation provides hourly values of meteorological and solar radiation parameters, encapsulating the climatic specifics of the location for the period it covers. The benefit of a TMY over other methods, such as averaging parameter values (see Table 2), is that it retains original, unaveraged data.

This study aims to generate the spatial distributions of global and beam radiation over the entire terrain of Taiwan, which are not available so far but are necessary references for assessing the performances of solar energy systems in Taiwan. Maps for global and beam radiation on the Taiwanese mainland are developed using the residual kriging method. The dataset for global radiation is excerpted directly from the TMY database (2004–2018) for the 30 CWB stations established by Hsieh et al. [5]. The information on the diffuse fraction for the 30 CWB stations is calculated using the three available correlation models [11,12], which account for the topographical and graphical factors in Taiwan, together with the TMY database for global radiation [5]. The dataset for beam radiation can be then generated through Equations (1) and (2) for this study.

2. Database

Two databases for global radiation and beam radiation from the 30 CWB weather stations were established prior to mapping.

2.1. Global Radiation Database

The monthly global radiation (MGR) database for the TMY from 2004 to 2018, derived from 30 CWB weather stations, is based on recent research by Hsieh et al. [5] and is listed in

Table 3. TMY data for monthly and annual global radiation (${\displaystyle \raisebox{1ex}{$\mathrm{M}\mathrm{J}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.}$) for the 30 CWB weather stations (excerpted from [5]).

Station Number/Name	January	February	March	April	May	June	July	August	September	October	November	December	Sum
1/Anbu	179.62	208.55	282.76	349.30	354.37	374.06	494.62	440.12	330.78	248.30	186.55	179.65	3268.7
2/Zhuzihu	164.20	179.73	334.12	218.92	335.81	416.09	342.42	429.70	341.63	353.64	197.38	208.13	3421.8
3/Tamsui	190.44	213.37	308.74	321.57	435.22	441.43	580.10	490.31	469.13	309.13	247.17	200.91	4208.1
4/Keelung	141.76	178.63	252.12	313.84	397.27	445.67	630.05	565.90	354.30	250.23	173.73	137.00	3840.5
5/Taipei	186.92	209.80	249.61	371.37	438.10	421.69	547.94	429.09	368.22	273.89	248.54	202.06	3947.2
6/Banqiao	214.77	224.03	257.07	318.59	451.30	434.64	603.76	509.89	419.07	304.32	231.09	189.91	3968.5
7/Hsinchu	231.26	224.59	276.61	353.47	422.42	460.72	570.72	506.62	456.43	405.70	270.44	212.41	4391.1
8/Wuqi	263.62	303.24	374.62	429.72	507.76	486.03	560.38	516.30	471.91	453.63	297.14	281.21	4945.6
9/Taichung	355.97	332.78	416.85	417.28	484.45	487.79	573.82	496.54	486.37	483.21	376.07	326.58	5237.7
10/Sun Moon Lake	311.97	310.18	342.28	318.25	359.59	390.65	457.56	400.69	367.36	388.29	329.75	276.01	4252.6
11/Alishan	359.17	331.43	387.33	379.50	372.37	376.03	408.91	330.54	326.02	374.32	359.04	348.19	4352.8
12/Yushan	400.85	385.51	456.67	412.03	454.42	471.16	547.76	468.82	411.77	524.81	376.57	394.10	5304.5
13/Chiayi	351.67	332.32	434.81	474.24	566.83	483.25	550.59	480.96	508.74	469.35	320.19	324.11	5297.1
14/Chigu	380.60	391.62	453.95	439.48	527.92	489.00	579.36	483.06	533.14	459.80	348.86	325.76	5412.6
15/Yongkang	336.12	351.04	383.23	445.50	539.81	470.92	511.31	464.71	440.55	461.17	341.67	315.96	5062.1
16/Tainan	366.81	369.90	482.78	498.13	550.66	537.61	572.01	481.67	457.07	464.28	348.68	350.34	5479.9
17/Kaohsiung	342.94	363.43	463.19	484.57	548.18	518.64	546.12	488.83	447.62	465.65	329.77	327.45	5326.4
18/Hengchun	335.57	414.00	447.52	486.01	518.02	505.44	523.87	460.64	456.74	449.28	375.64	323.97	5296.7
19/Dawa	276.22	319.26	345.00	411.85	499.42	531.17	551.18	516.73	462.86	434.28	325.72	251.87	4925.6
20/Taitung	289.65	297.98	338.25	430.25	560.46	588.17	678.33	585.66	484.83	451.15	330.13	302.22	5337.1
21/Chengong	232.33	259.12	296.97	328.47	402.78	561.51	661.58	552.57	488.42	416.07	312.95	256.76	4769.5
22/Hualien	201.38	220.22	286.63	303.06	420.74	499.39	623.50	552.67	438.94	342.35	255.14	221.22	4365.2
23/Su-Ao	166.25	189.63	252.48	332.22	385.58	505.51	581.22	582.96	399.17	273.58	191.91	170.51	4031.0
24/Yilan	171.19	212.05	282.50	337.25	402.22	458.92	630.85	559.38	435.54	277.87	189.70	156.94	4114.4
25/Pongjiayu	162.32	186.79	290.64	377.78	477.04	533.40	748.67	663.66	516.22	399.85	230.50	186.63	4773.5
26/Lanyu	219.97	240.67	307.95	364.64	419.61	431.12	541.76	423.75	371.95	368.67	251.38	229.60	4171.1
27/Penghu	218.22	235.57	343.08	424.50	498.01	527.09	611.54	529.94	478.05	410.44	298.95	256.43	4811.8
28/Dongjidao	283.88	293.42	380.38	470.03	575.94	532.77	621.03	557.39	510.96	455.64	307.64	261.34	5250.4
29/Matsu	195.47	247.70	289.66	358.42	419.64	407.83	601.57	518.26	379.56	339.09	215.17	194.03	4166.4
30/Kinmen	287.77	258.95	324.43	384.80	442.75	456.21	622.04	545.34	466.73	417.45	311.93	281.43	4799.8

. Out of these 30 stations, 24 (Stations 1–24, see Figure 1) are situated on the Taiwanese mainland, and their data are utilized to create spatial maps of solar radiation resources. The remaining six stations (Stations 25–30) are positioned on various remote islands, as illustrated in Figure 1.

2.2. Beam Radiation Database

Three correlation models developed specifically for the western part [11] and the eastern part [12] of the Taiwanese mainland and a remote island [12] (hereafter named the MW, ME, and MR correlation models, respectively) were used to generate the diffuse fraction data over the entire terrain of Taiwan. These models are all piecewise linear multiple-predictor functions in terms of predictors including the hourly and daily sky clearness indices ${(k}_t$ and $K_T$), global radiation persistence ($\psi$), solar altitude angle ($\alpha$, in radians), and apparent solar time (AST, in hours), as outlined below.

The MW correlation model:

$$\left\{\begin{array}{c}d=1.0 0\le k_t<0.2412\\ d=\min \left(1.0, 1.2207-0.6179k_t-0.1993\psi \right) 0.2412\le k_t<0.4091\\ d=\min (1.0, 1.6215-1.5379k_t+0.1486\alpha +0.006AST-0.3186K_T-0.4051\psi \\ 0.4091\le k_t<0.7222\\ d=0.9718-0.3805K_T-0.7169\psi 0.7222 \le k_t\le 1\end{array}\right.$$

(3)

The ME correlation model:

$$\left\{\begin{array}{c}d=1.0 0\le k_t<0.2876\\ d=\min \left(1.0, 1.3527-0.8494k_t-0.1469K_T-0.1304\psi \right) 0.2876\le k_t<0.4366\\ d=\min (1.0, 1.8375-1.7883k_t-0.4991K_T+ 0.4366\le k_t<0.7548\\ 0.1691\alpha -0.2669\psi ) \\ d=-0.2284+1.3803k_t-0.4153K_T-0.5520\psi 0.7548\le k_t\le 1\end{array}\right.$$

(4)

The MR correlation model:

$$\left\{\begin{array}{c}d=1.0 0\le k_t<0.2247\\ d=\min \left(1.0, 1.1588-0.4836k_t-0.0557K_T\right) 0.2247\le k_t<0.3863\\ d=\min (1.0, 1.7653-1.6041k_t-0.4417K_T+ 0.3863\le k_t\le 0.7641\\ 0.1012\alpha -0.2239\psi ) \\ d=0.9886-0.2061\alpha -0.7357\psi 0.7641\le k_t\le 1.0\end{array}\right.$$

(5)

with

$$k_t={\displaystyle \raisebox{1ex}{$I_{global}$}\!\left/ \!\raisebox{-1ex}{$I_0$}\right.}$$

(6)

$$K_T={\displaystyle \raisebox{1ex}{$\sum _{i=1}^{24}{I}_{global, i}$}\!\left/ \!\raisebox{-1ex}{${\sum }_{i=1}^{24}{I}_{0,i}$}\right.}$$

(7)

$$\psi =\left\{ \begin{array}{c}(k_{t-1}+k_{t+1})/2 \\ { k}_{t+1} for t=sunrise \\ k_{i-1} for t=sunset \end{array}\right.$$

(8)

where $I_0$ is the hourly extraterrestrial global horizontal radiation between the hour angles ${\omega }_1$ and ${\omega }_2$, defined as follows [32]:

$$I_0={\displaystyle \frac{12\times 3600}{\mathsf{\pi }}}{G}_{sc}\left(1+0.033 cos{\displaystyle \frac{360 n}{365}}\right)\left[cos\varphi cos\delta \left(sin{\omega }_2-sin{\omega }_1\right)+{\displaystyle \frac{\pi ({\omega }_2-{\omega }_1)}{180}}sin\varphi sin\delta \right]$$

(9)

where $G_{sc}$ denotes the solar constant, $n$ is the day of the year, $\varphi$ is the latitude, $\delta$ is the solar declination angle, and $\omega$ is the hour angle (in degrees). The performance differences between the MW and ME correlation models are mainly due to cloud coverage, influenced by the interaction of the Central Mountain Ranges with the Northeast Monsoon during the monsoon season (November–March). On the other hand, the performance differences between the MR and MW (or ME) models are attributed to atmospheric constituents, such as the humidity over remote islands, which is carried by moist sea air and higher than that over the inland areas of the Taiwanese mainland. Therefore, the MW and ME models are applicable to the leeward (western) and windward (eastern) sides of the Central Mountain Ranges, respectively, while the MR model is applicable to all remote islands in Taiwan’s terrain. Further details on these aspects can be found in the study of [12].

It was shown that, for better estimations of the diffuse fraction over the entire terrain of Taiwan, three d correlation models were needed. For CWB Stations 19–24 in the narrow eastern coastal plains (see Figure 1), Equation (4) was used to determine the d values. For the remaining CWB stations on the Taiwanese mainland, CWB Stations 1–18 (see Figure 1), Equation (3) was applied to determine the d values. Equation (5) was applied to all CWB stations on the remote islands to determine the d values, including CWB Stations 25–30, as shown in Figure 1.

The hourly diffuse fraction (d) data for a TMY (2004–2018) for each CWB station were determined using one of Equations (3)–(5), depending on its location in Taiwan, together with the available TMY hourly global radiation data in a publicly accessible repository, http://gitub.com/p4706115/Taiwan TMY (accessed on 20 March 2024), which were generated by a previous study by this research group in [5]. The hourly data for d and $I_{global}$ were used to determine the hourly diffuse radiation ($I_{diffuse})$ using Equation (2), which was next used to determine the hourly beam radiation ($I_{beam}$) through Equation (1). The daily, monthly, and annual beam radiation were, respectively, calculated for a day, a month, and a year based on the $I_{beam}$ data.

Table 4. TMY data for monthly and annual beam radiation (${\displaystyle \raisebox{1ex}{$\mathrm{M}\mathrm{J}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.}$) for 30 CWB weather stations.

Station Number/Name	January	February	March	April	May	June	July	August	September	October	November	December	Sum
1/Anbu	64.61	75.61	91.82	118.06	103.40	89.79	176.76	162.40	105.61	86.04	54.67	61.35	1190.1
2/Zhuzihu	43.04	40.40	137.67	15.07	94.66	110.23	26.91	124.94	120.35	72.32	21.60	87.75	849.9
3/Tamsui	55.07	81.59	95.78	79.60	171.67	83.96	196.30	141.24	195.48	110.06	68.45	77.29	1356.5
4/Keelung	43.53	70.96	108.77	93.87	116.81	151.96	278.76	265.72	131.87	80.23	54.03	28.54	1425.0
5/Taipei	59.48	78.51	61.32	111.43	121.69	90.29	185.31	110.11	148.22	89.07	81.02	72.70	1209.2
6/Banqiao	85.88	87.91	51.86	105.01	141.20	143.24	289.14	185.54	182.38	100.34	78.53	75.77	1526.8
7/Hsinchu	78.89	60.58	73.26	105.63	135.63	149.22	200.38	177.58	174.42	150.89	91.56	53.98	1452.0
8/Wuqi	67.76	85.64	74.51	133.50	134.63	115.50	171.98	206.87	238.78	173.25	105.05	106.57	1614.0
9/Taichung	191.73	152.09	178.57	115.09	161.36	146.76	207.16	133.16	194.48	214.44	187.33	168.94	2051.1
10/Sun Moon Lake	155.51	138.96	122.18	59.57	100.13	77.98	114.05	82.72	80.23	125.67	136.63	131.09	1324.7
11/Alishan	208.65	143.56	166.19	128.10	130.46	84.65	117.63	67.74	94.79	163.09	176.84	184.80	1666.5
12/Yushan	261.68	228.94	231.90	170.56	138.37	140.31	229.04	183.52	164.18	309.55	214.71	256.11	2528.9
13/Chiayi	161.71	118.53	156.54	180.09	226.90	185.47	181.43	187.81	236.46	216.46	114.42	127.97	2093.8
14/Chigu	205.97	163.06	157.14	117.33	227.56	203.20	203.81	149.35	220.88	156.12	145.35	124.04	2073.8
15/Yongkang	151.62	150.64	104.92	127.28	210.09	123.37	156.14	203.73	151.08	206.23	139.29	147.58	1872.0
16/Tainan	150.59	152.68	191.10	202.89	252.17	265.16	204.15	186.28	158.56	218.32	142.08	154.58	2278.6
17/Kaohsiung	131.03	133.16	176.20	168.66	198.21	192.38	227.37	180.43	146.02	182.92	115.97	118.23	1970.6
18/Hengchun	121.87	185.86	143.32	171.55	206.42	186.07	176.93	144.26	148.74	176.45	161.74	132.75	1956.0
19/Dawa	105.38	136.19	87.90	132.41	147.67	183.54	197.70	274.24	183.54	189.38	129.34	81.53	1848.8
20/Taitung	119.57	115.72	124.62	214.74	283.12	276.15	368.11	320.64	272.97	238.04	143.77	126.33	2603.8
21/Chengong	74.97	90.42	52.87	97.85	139.78	250.65	349.27	281.63	280.68	215.39	139.45	101.71	2074.7
22/Hualien	73.16	70.74	105.11	81.47	103.20	199.70	262.10	229.48	164.06	128.03	98.89	93.79	1609.7
23/Su-Ao	59.09	67.89	69.18	108.11	103.56	182.11	267.18	298.37	134.43	93.78	69.79	68.97	1522.5
24/Yilan	53.30	86.03	76.26	92.36	140.26	121.90	297.94	275.27	191.14	114.14	63.58	45.64	1557.8
25/Pongjiayu	59.43	79.60	129.66	151.95	183.04	210.98	463.17	310.33	226.38	182.12	71.29	77.73	2145.7
26/Lanyu	56.80	75.39	72.53	95.13	117.95	95.10	214.58	127.39	134.64	119.75	68.59	72.50	1250.4
27/Penghu	55.80	78.96	78.65	112.47	155.53	171.29	244.56	207.58	214.30	128.76	99.92	66.76	1614.6
28/Dongjidao	87.96	102.50	102.03	127.75	206.42	227.03	253.87	204.23	201.38	183.23	115.03	78.96	1890.4
29/Matsu	83.67	90.69	57.17	89.80	91.76	104.64	206.25	159.21	111.79	96.85	57.80	54.62	1204.2
30/Kinmen	92.87	89.40	86.17	143.89	165.70	124.34	244.68	178.44	175.39	138.87	107.58	101.36	1648.7

lists the TMY database for the monthly and annual beam radiation for the 30 CWB weather stations.

3. Spatial Interpolation Methodology

3.1. Interpolation

Kriging is a group of least-square linear regression techniques used to estimate values for a variable at unsampled locations by analyzing the spatial patterns of observed data. The ordinary kriging method [18] utilizes two key pieces of information: the variation in the attribute and the distance between data points. This method assumes that the data are a continuous realization of a random function, $z\left(x\right)$, in two-dimensional space, and that it is second-order stationary, meaning that the mean and variance are dependent on the lag distance, $h$. However, most datasets do not meet this criterion, as variance increases with domain size. Instead, a less strict intrinsic hypothesis is applied, which requires that the variance of the first-order increment, $z(x+h)-z\left(x\right)$, be finite and the second-order be stationary. This variance of the increment is represented by a semivariogram, $\gamma \left(h\right),$ defined as follows:

$$\mathsf{\gamma }\left(\mathrm{h}\right)={\displaystyle \frac{1}{2}}\mathsf{{\rm E}}\left\{{[z\left(x+h\right)-z\left(h\right)]}^2\right\}$$

(10)

where $\mathsf{{\rm E}}$ denotes the expectation. The exponential model, which has been commonly used for studies with the kriging family [22,23,25,31], and a linear model, as expressed, respectively, in the two following equations, are used in the study:

$$\gamma \left(h\right)=C_0+C_1\left[1-exp\left({\displaystyle \frac{-h}{r}}\right)\right]$$

(11)

$$\gamma \left(h\right)=bh$$

(12)

where $C_0, C_1, r,$ and b are, respectively, the nugget, sill, range, and slope.

3.2. Residual Kriging

Residual kriging is a mixed technique, which isolates the first-order constituent trend from the other stochastic model using a multiple linear regression operation:

$$z\left(x\right)=z^{\ast }\left(x\right)+r\left(x\right)=\sum {\alpha }_i{a}_i\left(x\right)+r\left(x\right)$$

(13)

where x is located within the domain to which the observations, $a_i\left(x\right)$, pertain. This study uses three independent variables, latitude ($x_1$), longitude ($x_2),$ and altitude ($x_3$), so Equation (13) is rewritten as follows:

$$z\left(x\right)={\alpha }_0+{\alpha }_1{x}_1\left(x\right)+{\alpha }_2{x}_2\left(x\right)+{\alpha }_3{x}_3\left(x\right)+r\left(x\right)$$

(14)

Another nonlinear regression operation that was devised by Wu and Li [33] is written as follows:

$$\begin{gathered} z^{\ast }\left(x\right)={\beta }_0+{\beta }_1{x}_1\left(x\right)+{\beta }_2{x}_2\left(x\right)+{\beta }_3{x}_1^2\left(x\right)+{\beta }_4{x}_2^2\left(x\right)+{\beta }_5{x}_1\left(x\right)x_2\left(x\right) \\ +{\beta }_6{x}_3\left(x\right)+{\beta }_7{x}_1\left(x\right)x_3\left(x\right)+{\beta }_8{x}_2\left(x\right)x_3\left(x\right) \end{gathered}$$

(15)

These two regression operators are used for this study. The coefficients ${\alpha }_i$ and ${\beta }_i$ are fitted using a standard least-square procedure. The spatial field $z^{\ast }\left(x\right)$ is determined for each location, x, where ${\alpha }_i\left(x\right)$ or ${\beta }_i\left(x\right)$ are known using Equation (14) or Equation (15), respectively, so a residual error, r(x), is determined as the difference between the observation, z(x), and the regression estimation, $z^{\ast }\left(x\right)$, which is given by

$$r\left(x\right)=z\left(x\right)-z^{\ast }\left(x\right)$$

(16)

Ordinary kriging is used to obtain the final residuals, $\widehat{r}\left(x\right)$, representing the corrections to apply to the regression model. The final estimates $\widehat{z}\left(x\right)$ are obtained as follows:

$$\widehat{z}\left(x\right)=z^{\ast }\left(x\right)+\widehat{r}\left(x\right)$$

(17)

3.3. Cross Validation

Cross validation is the leave-one-out estimate for each observation and is a technique for assessing the performance of a statistical analysis method. It removes one data location and performs the regression estimate using the data for the rest of the locations. This estimate is compared with the data for the excluded location. The process iterates until all data locations are used. Five different measures of fit for the monthly global radiation were used to determine how well an interpolated result, ${MGR}_{est}\left(x_i\right),$ represents the observed data, ${MGR}_{obs}\left(x_i\right),$ at a location, $x_i$, including the mean error (ME), the mean percentage error (MPE), the mean absolute error (MAE), the mean absolute percentage error (MAPE), and the root-mean-square error (RMSE), which are, respectively, defined as follows:

$$\mathrm{M}\mathrm{E}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{∆}_{bias, i}$$

(18)

$$\mathrm{M}\mathrm{P}\mathrm{E}={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left\{1-{\displaystyle \frac{{MGR}_{est}\left(x_i\right)}{{MGR}_{obs}\left(x_i\right)}}\right\}\times 100\%$$

(19)

$$\mathrm{M}\mathrm{A}\mathrm{E}={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left|{∆}_{bias,i}\right|$$

(20)

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left|1-{\displaystyle \frac{{MGR}_{est}\left(x\right)}{{MGR}_{obs}\left(x_i\right)}}\right|\times 100\%$$

(21)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\left[{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left({∆}_{bias,i}\right)}^2\right]}^{1/2}$$

(22)

where

$${∆}_{bias,i}={MGR}_{obs,i}-{MGR}_{est,i}$$

(23)

There are 24 CWB weather stations within the Taiwanese mainland (see Figure 1), so n has a value of 24 for the five measures.

3.4. Determination of Semivariogram Model, Regression Formula, and Grid Size

There were two semivariogram models employed for this study: the exponential, Equation (11), and the linear, Equation (12), models. A multiple linear regression formula, Equation (14), and a multiple nonlinear regression formula, Equation (15), were considered in the study. To determine the best combination of semivariogram model and regression formula, four cases in

Table 5. Test cases using four different combinations of the semivariogram model and the multiple regression formula.

Case	Semivariogram Model		Multiple Regression Formula
	Exponential, Equation (11)	Linear, Equation (12)	Linear, Equation (14)	Nonlinear, Equation (15)
1		√	√
2		√		√
3	√		√
4	√			√

were tested using the residual kriging method with a grid size of 1 km × 1 km. Validation error statistics for the four combination cases using the MPE, MAPE, ME, MAE, and RMSE are listed in

Table 6. Comparison of monthly error statistics for the four combination cases in Table 5 using the residual kriging method.

Month	January	February	March	April	May	June	July	August	September	October	November	December	Mean **
MPE-1 * (%)	0.51	0.38	0.70	1.16	0.79	0.51	1.64	0.95	0.96	0.43	0.22	0.58	0.74
MAPE-1 (%)	5.73	4.57	6.96	8.24	6.70	5.24	9.64	7.72	8.75	3.57	1.99	6.06	6.27
ME-1$(\times {10}^5 \mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	−4.90	−3.30	3.77	1.08	0.144	−1.10	−1.90	−0.035	−0.130	−11.0	−23.0	−1.10	−3.50
MAE-1$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.469	0.455	0.780	0.945	0.926	0.788	1.550	1.144	1.212	0.433	0.173	0.472	0.780
RMSE-1$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.623	0.568	0.917	1.146	1.234	1.040	1.980	1.436	1.335	0.647	0.263	0.566	0.980
MPE-2 (%)	0.36	0.28	0.61	0.98	0.31	0.10	0.69	0.24	0.25	0.05	0.32	0.54	0.39
MAPE-2 (%)	4.55	3.95	6.14	7.75	3.91	2.59	5.63	3.53	3.09	2.24	3.92	5.95	4.44
ME-2$(\times {10}^5 \mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	−4.70	−5.50	0.674	~ 0	1.06	8.88	1.10	1.45	−1.60	1.09	−2.00	0.470	0.074
MAE-2$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.392	0.384	0.674	0.895	0.562	0.410	0.925	0.561	0.447	0.260	0.380	0.464	0.530
RMSE-2$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.508	0.471	0.858	1.072	0.812	0.505	1.270	0.761	0.634	0.320	0.519	0.571	0.692
MPE-3 (%)	0.41	0.04	0.11	0.06	0.73	0.42	0.89	0.70	0.80	0.62	0.73	0.62	0.51
MAPE-3 (%)	4.74	0.59	0.42	0.43	6.34	4.45	5.22	5.45	6.85	5.00	5.82	6.32	4.30
ME-3$(\times {10}^5 \mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	−9.40	−3.51	81.74	−3.90	1.04	0.585	6.59	−3.90	−9.22	−9.70	−7.12	−1.21	64.2
MAE-3$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.393	0.0547	0.0464	0.0476	0.877	0.672	0.840	0.792	0.940	0.613	0.511	0.492	0.523
RMSE-3$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.528	0.0721	0.0571	0.0622	1.147	0.882	1.081	1.072	1.112	0.861	0.668	0.593	0.678
MPE-4 (%)	0.28	0.20	0.05	0.05	0.21	0.09	0.25	0.21	0.02	0.06	0.02	0.30	0.14
MAPE-4 (%)	3.35	2.69	0.36	0.43	2.66	2.47	1.89	3.20	0.09	2.29	0.20	3.16	1.90
ME-4$(\times {10}^5 \mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	−9.10	−9.50	11.07	−6.30	4.22	1.06	20.1	2.95	23.3	0.587	1.06	−1.30	3.99
MAE-4$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.288	0.256	0.0397	0.0485	0.376	0.392	0.312	0.508	0.0123	0.266	0.0187	0.246	0.230
RMSE-4$(\mathrm{M}\mathrm{J}/{\mathrm{m}}^2$)	0.373	0.308	0.0567	0.0641	0.525	0.485	0.428	0.680	0.0174	0.328	0.0238	0.306	0.300

* The digit behind the segment line denotes the case number. ** Mean monthly value.

. Note that the values for the ME are increased by an order of ${10}^5$ in Table 6. Kriging is an unbiased estimator that assumes the expected value of the estimate to be equal to the observed value, so all the ME values for the four cases were close to zero, as shown in Table 6, and were not used for comparison. A comparison of the other four statistical errors shows that Case 4 (a combination of the exponential semivariogram model, Equation (11), and the multiple nonlinear regression function, Equation (15), see Table 5) achieved the best performance and was used for this study.

A grid independence test was made with the cross validation for the Tainan weather station (Station 16 in Table 1 and Figure 1) using three grid sizes of $0.1 \mathrm{k}\mathrm{m}\times 0.1 \mathrm{k}\mathrm{m}$, 1 km × 1 km, and 10 km × 10 km, and the results are listed in

Table 7. Grid independence test: statistical (MPE, %) errors between observed and estimated values for Tainan CWB weather station.

Grid Size (km by km)	January	February	March	April	May	June	July	August	September	October	November	December
10 × 10	−10.25	3.92	−5.92	−5.14	−0.74	−10.01	−2.14	−0.03	5.16	−0.13	−2.63	−5.13
1 × 1	−10.09	3.86	−5.24	−5.11	−0.75	−10.16	−2.16	−0.04	5.10	−0.10	−2.61	−5.10
0.1 × 0.1	−10.08	3.86	−5.20	−5.11	−0.74	−10.15	−2.15	−0.04	5.09	−0.10	−2.61	−5.09

. The statistical (MPE) error for each TMM for the Tainan weather station approaches its asymptotic value for a grid size of 1 km × 1 km for the residual kriging computation, so a grid size of 1 km × 1 km was used for this study.

4. Results and Discussion

The predicted global radiation was validated for each CWB weather station on the Taiwanese mainland using a cross-validation procedure and TMM data. For the sake of conciseness, only results for the statistical indicator of the MPE, Equation (19), are shown in Figure 2. Figure 2 presents the MPEs of 12 TMMs (from January to December) for 24 CWB weather stations located on the Taiwanese mainland. Except for four stations, Anbu (Station 1 in Table 1 and Figure 1), Zhuzihu (Station 2), Alishan (Station 11), and Yushan (Station 12) which will be elaborated on later, the MPEs for the cross-validation for each weather station fall within ($-$20%, 25%) for all months of each TMY. It is noted that the month here is a TMM which is the representative month between 2004 and 2018 for each CWB weather station. A given TMM at different CWB weather stations is usually selected from different years in between 2004 and 2018. It certainly leads to a more scattering spatial distribution for global radiation among the 24 CWB stations on the Taiwanese mainland for a given TMM as compared to a case under the condition of selecting a specified month of the same year for all 24 CWB stations, which would have a smoother spatial continuity for the global radiation among them. A similar situation was observed for the development of a global-radiation map using a TMY (1985–1999) for Greece [26]. The study showed relatively large deviations and greater scattering between the predicted and observed annual global radiation data for 39 locations with an $R^2$ value of 0.87 (see Figure 2 of [26]).

The Alishan (Station 11) and Yushan (Station 12) weather stations are the two outlying stations in the Central Mountain Ranges with an elevation of more than 2000 m (see Table 1). Mountainous areas feature complicated topography in terms of slope and aspect as the sun rays irradiate the ground [25,34], which may blockade some beam and diffuse radiation from the sky dome or receive more reflected radiation from the surroundings to the ground [35]. Therefore, the data quality for the measured global radiation at these two outlying stations is not as good as that for other stations on a monotonous landscape at a low altitude. In particular, Yushan station (altitude of 3844.8 m, Table 1) is the highest station on the Taiwanese mainland. Its $z^{\ast }$ variables are extrapolated from the lower $x_3$ (latitude) values for nearby stations using Equation (15), so there are greater scatterings for its MPE than that for other stations. An extreme outlier for the MPE (around $-$100%) occurs at this station in February, as shown in Figure 2.

Stations 1 and 2 (Anbu and Zhuzihu, respectively) are located in a volcanic mountainous area of the northern Taiwanese mainland. The horizontal distance between the Anbu (at ${25.18}^{° }\mathrm{N},{121.53}^{° }\mathrm{E}$) and Zhuzihu (at ${25.16}^{° }\mathrm{N},{121.54}^{° }\mathrm{E}$) stations is short, but their elevations are significantly different, at 837.6 m for Anbu and 607.1 m for Zhuzihu. The annual TMY global radiation for both stations differs by 200 ${\displaystyle \raisebox{1ex}{$\mathrm{M}\mathrm{J}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{m}}^2$}\right.}$ [5]. Variability in elevation due to topographic features creates a high local gradient in the solar radiation [22]. The reliability of the kriging estimate decreases for this case because the simulation using the exponential semivariogram model is affected by sudden local variations in radiation with a very short horizontal distance, which looks mathematically like a jump condition. The higher deviations in the MPE for Stations 1 and 2 in Figure 2 are attributed to this factor. However, except for these four weather stations, the remaining CWB weather stations on the Taiwanese mainland are mostly located in a monotonous landscape with a low terrain (see Table 1), where the solar energy applications in Taiwan are primarily promoted to, and their MPEs for the cross-validation are within ($-$20%, 25%).

The spatial distributions of monthly global radiation for the Taiwanese mainland were created using residual kriging on basis of the TMY global radiation data for the 24 CWB weather stations in Table 3, and the results are shown in Figure 3. The annual average daily quantity was determined by summing the 12 values for the monthly global radiation for each grid and dividing it by 365 days. Figure 4a shows the spatial distribution of the annual average daily global radiation. The spatial distributions for monthly beam radiation for the Taiwanese mainland were also created using the TMY data for the beam radiation for the 24 CWB weather stations in Table 4, and the results are shown in Figure 5. The spatial distribution for the annual average daily beam radiation is shown in Figure 4b.

Seasonal trends for the monthly variations in the global radiation (Figure 3) and beam radiation (Figure 5) are very similar. The Northeast Monsoon from November to March dominates the weather in the northeastern Taiwanese mainland. During the monsoon season, moist ocean air is blown southwest-bound and is lifted upward by the Central Mountain Ranges. Therefore, cloudy skies and orographic rain are frequent features in the windward side of the mountains (i.e., the eastern coastland) during the monsoon season. This results in low levels of monthly global and beam radiation in the northern Taiwanese mainland and the eastern coastland during November–March, as, respectively, shown in Figure 3 and Figure 5. Orographic rain sometimes lasts until April on the eastern coastland and results in relatively low monthly beam radiation, as shown in Figure 5d. The Northeast Monsoon has a smaller influence on the monthly global and beam radiation on the leeward side of the Central Mountain Ranges (i.e., the western part of the Taiwanese mainland), as shown in Figure 3d–j and Figure 5e–j, respectively. The relative lower levels of monthly global and beam radiation observed in the regions of mountains and hills in comparison to the regions of plains/basins for the western part of the Taiwanese mainland after the monsoon season are attributed to the altitude factor [28]. It shows significant topographic effects on both global and beam radiation between the western and eastern parts of the Taiwanese mainland during the Northeast Monsoon season. For the other months in a year, the regions to the east (sunrise) side of the Central Mountain Ranges experience higher global radiation and beam radiation than regions to the west (sunset) side, as, respectively, shown in Figure 3e–j and Figure 5e–j. The spatial distributions of the annual results for global radiation and beam radiation are shown in Figure 4. It reveals that the eastern Taiwanese mainland has more solar resources than the western part. It is agreed that solar resources become more abundant as latitude decreases toward the equator. The results in Figure 4 exhibit this trend.

Extraterrestrial radiation on a horizontal plane is a function of latitude and is independent of longitude, as shown in Equation (9). This value theoretically defines the upper limit of global radiation that is incident on the ground. In order to investigate the geographic effect on the solar energy resources in Taiwan, three weather stations located on remote islands, including the Kinmen (Station 30) and Penghu (Station 27) stations situated in the Taiwan Strait as well as the Lanyu (Station 26) station situated in the west Pacific, were selected for comparing their annual global and beam radiation with three corresponding weather stations each located on the Taiwanese mainland and with closer latitude to one of these three remote-island stations, as summarized in

Table 8. Comparison between the annual global and beam radiation observed for the weather station on a remote island and for their closest-latitude weather station located on the Taiwanese mainland: at Kinmen station versus Taichung station.

Station Number	$\mathbf{Global} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference ***	$\mathbf{Beam} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference
30 * (24.41° N, 118.29° E)	4799.8	$-$8.36%	1648.7	$-$19.62%
9 (24.14° N, 120.68° E)	5237.7		2051.1

* located in the Taiwan Strait (see Figure 1). *** on the basis of the data for the station on the Taiwanese mainland.

,

Table 9. Comparison between the annual global and beam radiation observed for the weather station on a remote island and for their closest-latitude weather station located on the Taiwanese mainland: at Penghu station versus Chiayi station.

Station Number	$\mathbf{Global} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference	$\mathbf{Beam} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference
27 * (23.56° N, 119.56° E)	4811.8	$-$9.16%	1614.6	$-$22.88%
9 (23.50° N, 120.43° E)	5297.1		2093.8

* located in the Taiwan Strait (see Figure 1).

and

Table 10. Comparison between the annual global and beam radiation observed for the weather station on a remote island and for their closest-latitude weather station located on the Taiwanese mainland: at Lanyu station versus Hengchun station.

Station Number	$\mathbf{Global} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference	$\mathbf{Beam} \mathbf{Radiation} \left(\raisebox{1ex}{$\mathbf{M}\mathbf{J}$}\!\left/ \!\raisebox{-1ex}{${\mathbf{m}}^2$}\right.\right)$	Relative Difference
26 ** (22.04° N, 121.56° E)	4171.1	$-$21.25%	1250.4	36.07%
18 (22.00° N, 120.75° E)	5296.7		1956.0

** located in the western Pacific (see Figure 1).

. The comparison results in Table 8, Table 9 and Table 10 show remarkable effects on solar radiation, particularly on beam radiation. This is attributed to the moist ocean (or sea) air over the sky dome of the remote islands. Water vapor can absorb some solar radiation (spanned in the wavelength range of 0.3 to 2.8 $\mathsf{\mu }$m) when sun rays pass through the atmosphere and thus reduces the incident solar (beam radiation) intensity on ground. The excited water molecules in the atmosphere have to return to the ground (stable) state associated with the process of emitting diffuse radiation outward by moving half-upward (i.e., escaping to space) and half-downward (i.e., being incident on ground). This explains why the high humidity in the sky dome of a remote island has a higher impact on the reduction in beam radiation than that of global radiation, which is the sum of beam and diffuse radiation. The results in Table 8, Table 9 and Table 10 also show that the geographical effect on solar radiation is more apparent for the remote islands in the Pacific Ocean than the remote islands in the Taiwan Strait.

5. Conclusions

A nationwide TMY (2004–2018) database for global radiation from the 30 CWB weather stations, of which twenty-four are located on the Taiwanese mainland and six are located on six different remote islands, was used to determine the spatial distribution for global radiation over the terrain of Taiwan. There is a lack of beam (and diffuse) radiation information in daily reports from all CWB weather stations. Therefore, information on the diffuse fraction for all CWB weather stations was estimated using three available correlation models that were developed for Taiwan and account for the topographic and geographical effects on it on basis of the TMY data. A database for beam radiation for the 30 CWB weather stations was then generated using the estimated diffuse fractions together with the TMY database for global radiation.

Mapping of the global radiation or beam radiation on the Taiwanese mainland used the databases from the 24 CWB weather stations on it and the residual kriging method. Seasonal patterns for monthly variations in global radiation and beam radiation due to the interaction between the Central Mountain Ranges and the Northeast Monsoon are similar. The Central Mountain Ranges have a significant topographic effect on the spatial distributions of global radiation and beam radiation for the western and eastern Taiwanese mainland. The geographic effect is more significant for beam radiation than for global radiation because there is constant moist ocean (or sea) air in the sky dome over the remote islands.

In terms of the databases for global radiation (Table 3) and beam radiation (Table 4) for the six CWB weather stations on remote islands, no mapping actions were performed. The databases in Table 3 and Table 4 for Station 25 (Pongjiayu) were applied to all the remote islands off the northeast coast of the Taiwanese mainland, and the databases in Table 3 and Table 4 for Station 26 (Lanyu) were applied to the two inhabited islands: Lanyu and Green (at at 22.67° N, 121.48° E) Islands, off the southeast coast of the Taiwanese mainland. There are two CWB stations in the Penghu archipelago: Stations 27 (Penghu) and 28 (Dongjidoa). The databases in Table 3 and Table 4 for Station 27 and for Station 28 were, respectively, applied to the nearby remote islands located beyond and beneath the Tropic of Cancer (latitude: 23.5° N). The databases in Table 3 and Table 4 for Station 29 (Matsu) and Station 30 (Kinmen) were, respectively, applied to all islands in the Matsu archipelago and to those in the Kinmen archipelago, as shown in Figure 1.

As solar energy continues to play a pivotal role in the transition to sustainable energy sources, the spatial distributions of credible, long-term databases for global radiation and beam radiation presented are good scientific references for performing assessments of solar energy systems in Taiwan.