Types of heat generation in Ukraine in 2016 and their cost
Январь 31st, 2016
The calculation of the direct normal irradiance is performed in different steps. In a first step cloud information is derived from the satellite measurements using an enhanced version of the Heliosat method [1]. As a measure of cloudiness a dimensionless clear sky index k* is derived. Using a clear sky model the global irradiance is immediately related to the clear sky index, while our new beam fraction model allows for calculating the DNI. Fig. 1 represents the necessary process steps, which are described in this section.
Fig. 1: Significant process steps of the Heliosat method. 
The model HELIOSAT2 arises as a natural evolution of the model called HELIOSAT developed by Cano and Diabate during the 80s. The progress introduced in this model corresponds to modifications on the model of clear sky used by means of the introduction of some descriptive parameter of the atmospheric quality. In both models estimation of a dimensionless parameter (n), called cloud index, is proposed on a basis of using the brilliance detected by the satellite.
The cloud index is estimated for every pixel of the image according to so called apparent or
instantaneous albedo (p ), albedo of the clouds (Pn) and ground albedo (pg) as shown in the following expression:
n = (p‘ — Pg ){p’n — Pg )1 (1)
The evaluation of the different albedos necessary for the calculation of the cloud index is obtained from the radiance (L), the geographical parameters that define the pixel in question (latitude, length and altitude) and the information about the atmospheric local quality contained in the coefficient of atmospheric turbidity of Linke (TL).
To conclude, the model HELIOSAT2 calculates the global horizontal irradiance by means of a relation between the so called clear sky index Kc (global / global of clear sky), and the cloud index (n). The relation used in this work is a modification of the original proposal in which the local median of the cloud index is used [5]. So that the global irradiance estimated in every pixel can be calculated as:
where GCs is the global irradiance of clear sky calculated by the ESRA [6] algorithm. The final daily irradiation is calculated by integration of the different images available in a day.
This activity focuses on collection of high quality ground measurements which can be used as a reference in the benchmarking exercises. The measurements should be conducted with high accuracy, high frequency and traceable maintenance of the equipment. Data has been collected from the Baseline surface radiation network (BSRN), International Daylight Measurement programme (IDMP), the meteomedia network, the World Radiation Data Center (WRDC) and the Global Atmospheric Watch (GAW) programme. In addition further measurements were collected from scientific institutions, providing they fulfil the quality criteria above.
A common quality control procedure has been defined for all broadband time series data. The parameters for the quality assessment have been deducted from the Baseline Surface Radiation Network Operation Manual [3] and operational experience of the partners involved. All data is flagged if:
• It is within the physical possible limits. The applied physical limits are shown in table 1.
• It is lower than it would be in a dry and clean atmosphere. We apply the e. g. the Bird clear sky model (Model C in [4]) with aerosols and water vapour set to zero.
• The components of direct, diffuse and global parameters are coherent (if all are available). The conditions and limits are described in table 2.
Only values passing all tests will be used in the benchmarking exercises.
Table 1: Physical limits for quality control, with: G0: irradiance available at the top of atmosphere; GV0: illuminance available at the top of atmosphere; 9z: solar zenith angle

Parameter 
condition 
Limits 

GGlobal 
Oz < 75°, Effuse + GDirect cos Oz 
V О 
1.0±8% 
GDiffuse + GDirect cos Oz 

GGlobal 
93° > dz > 75°, Gf + GDirect 
W cosOz > 50—j — m 
1.0±15% 
GDiffuse + GDirect cos Oz 

GDiffuse GGlobal 
Oz < 75°, GOMaI > 50W m 
< 1.05 

GDiffuse GGlobal 
93°> Oz > 75°, GOMaI > 50W m 
< 1.10 
The maps give good insight into regional differences, however due to the coarse grid resolution (5 arc minutes) they smooth out local effects. Some of the integrated systems incorporate specific algorithms, such as an application of highresolution Digital Elevation Model (DEM), which may increase the difference between results for a given site. Thus, for a set of 37 randomly selected points we analyse the values and resulting differences that a user obtains for a particular site when consulting directly each data source/system. In areas with higher agreement of all 6 databases (standard deviation lower than 4%) we selected 15 points, and in areas with higher disagreement (standard deviation between 5% and 11%) another 22 points.
For each point we calculated the yearly sum of global irradiation for (1) southoriented surface inclined at 34 degrees (close to the typical optimum angle for PV systems in Europe), and (2) for 2axis suntracking systems.
While the map of yearly sum of global horizontal irradiation (Fig. 1) shows the average of 6 databases, the map of standard deviation (Fig. 2) indicates how much the datasets differ in various regions. Higher disagreement between the databases can be found not only in higher mountains (the Alps, Pyrenees, Carpathians, etc.), but also along some coastal zones and in flat regions, e. g. around the Baltic and North Seas. Higher standard deviation in Bulgaria and Romania, in South Scandinavia, and Po plain relates to uneven distribution of the input ground data and questionable quality of the older ground measurements that are used in ESRA and PVGIS. However, cumulative distribution of the map values (Fig. 4 right) shows that standard deviation does not exceed 6% in 78% of the study region, and only in rare (extreme) cases is it higher than 12%.
The analysis on 37 selected points shows some cases with higher differences (some of them are not present in the coarseresolution maps) with average standard deviation 5.6%. Mean bias difference (see yellow box in the Fig. 3) indicates that satellite databases (SatelLight and HelioClim2) show generally higher values compared to the products relying on ground measurements (ESRA, Meteonorm and PVGIS). This difference may be partially explained by the time period of data covered. While SatelLight and HelioClim2 represent last 10 years, ESRA and PVGIS are based on one decade of data (1980s), and Meteonorm represents two decades (1980s and 1990s). The NASA SSE product partially spans both periods and also generally gives results that are intermediate between the two groups.
Fig. 1. Yearly sum of global irradiation on horizontal surface — average of 6 databases: Meteonorm v.6, ESRA, PVGIS, NASA SSE v.6, SatelLight and HelioClim2 [kWh/m2]. Results for 37 randomly selected sites are shown in Fig. 3. 
Fig. 2. Yearly sum of global horizontal irradiation — standard deviation of the values from 6 databases relative to the overall average shown in Fig. 1 [%]. Results for 37 randomly selected sites are shown in Fig. 3. 
Part of the contribution from higher values obtained from HelioClim2 database can be probably explained by Mean Bias Deviation calculated by comparison with data from 11 highquality BSRN stations ([9] see Tab. 1).
For the same 37 points, global irradiation has been calculated for a southoriented surface inclined at an angle of 34 degrees. Results from five systems are available (all except HelioClim2). From linear fit of the values on the scattergram (Fig. 4 left) one can see that the standard deviation of
estimations for an inclined surface increases in average by 21% compared to the horizontal irradiation (correlation coefficient 0.98). The differences between values for 2axis tracking surface are calculated only for 3 systems: Meteonorm, NASA SSE/RETScreen and PVGIS. Here the standard deviation increases by 36% compared to the estimates for horizontal surface (correlation coefficient 0.97). Both simulations show a good agreement in the implementation of the underlying algorithms in the compared software.
Fig. 3. Yearly sum of global horizontal irradiation — differences of the values from 6 databases relative to the overall average. First 15 points represent areas with higher agreement between databases; the other 22 points are randomly selected in areas where the difference between the databases is higher. 
Fig. 4. Estimates of yearly sum of global irradiation. Left: Scattergram of 37 points showing relation between relative standard deviation calculated for horizontal surface and those calculated for 34° inclined and 2axis tracking system [%]; Right: Cumulative distribution of relative standard deviation [%] for global horizontal irradiation summarised from the map and the standard deviation estimated for inclined, and 2axis tracking surface by linear fit of points from the left. 
By integration of findings from the map and point analyses (Figs. 2 and 4 left) we can estimate uncertainty for any point on the map for yearly values for horizontal, inclined and 2axis tracking surfaces. Thus the cumulative distribution of user’s uncertainty can be calculated from the map data for the whole study region. Fig. 4 (right) shows that for 90% of the study region the
uncertainty of estimates of yearly global irradiation (expressed by standard deviation) is lower than 7% for a horizontal surface, 8.3% for a southfacing surface inclined at 34°, and 10% for a 2axis tracking surface.
Differences in the estimates using several solar radiation databases relate to a number of factors.
Analysis of interannual variability for 70 meteorological stations over Europe having data for at least 10 years shows that the standard deviation from the longterm average ranges typically between 3% and 6%. To capture this variability, a database should include data for at least 5 to 10 years depending on the climate region.
There is an inherent difference between in situ (ground) and satellite observations, and the methods how these data are processed. Databases relying on the interpolation of ground observations (Meteonorm, ESRA, and PVGIS Europe) are sensitive to the quality of measurements and density of the measuring stations (which is not satisfactory in many regions), and they typically represent only statistical values. The satellitederived databases (e. g. HelioClim, NASA SSE, and Satellight) are more affected by higher uncertainty of the cloud cover assessment when the ground is covered by snow and ice and for low sun angles, but they offer time series with high time resolution (e. g. hourly data) and provide spatiallycontinuous coverage.
Quality and the spatial detail of spatial database is determined by input data used in the models, mainly parameters describing the optical state of the atmosphere (such as Linke atmospheric turbidity, ozone, water vapour, aerosol optical depth), and Digital Elevation Models. The community of data developers lacks high quality aerosol data. Highresolution DEM is presently considered only in Meteonorm and PVGIS. Databases with coarser spatial resolution (e. g. NASA SSE) provide good regional estimates, however for studies at local level they may show deviations as they ignore local climate and terrain features.
The crosscomparison approach presented in this study assumes that each database has equal weight. Uncertainty can be reduced by weighting each contribution that is based on indicators of quality, number of datayears included, spatial resolution and incorporation of terrain effects. Including another two European satellitederived databases, those owned by German Aerospace Center (DLR) and University of Oldenburg (both in Germany), and application of weighting will provide more complete picture of the uncertainty in the solar resource assessment.
This study provides a first insight into spatial distribution of uncertainty of solar radiation estimates by relative cross comparison of six data sources. In this stage only the yearly sum of global irradiation is considered, and all databases are assumed to give an equal contribution to the overall average. The map of standard deviation from the average indicates combined effect of differences between the databases, and in this study it is used as an indicator of the user’s uncertainty.
Differences at the regional level indicate that within 90% of the study region the uncertainty of yearly global irradiation estimates expressed by standard deviation does not exceed 7% for horizontal surface, 8.3% for surface inclined at 34°, and 10% for 2axis tracking surface. A user has to expect higher differences in the outputs from the studied databases in complex climate conditions of mountains, some coastal zones and in areas where solar radiation modelling cannot rely on sufficient density and quality of input data.
[1] J. Remund, S. Kunz, C. Schilter, Handbook of METEONORM Version 6.0, Part II: theory. Meteotest, 2007, http://www. meteonorm. com/.
[2] J. Greif, K. Scharmer, Eds., European Solar Radiation Atlas (ESRA), 4th edition. Scientific advisors: R. Dogniaux, and J. Page; Authors: L. Wald, M. Albuisson, G. Czeplak, B. Bourges, R. Aguiar, H. Lund, A. Joukoff, U. Terzenbach, H.G. Beyer, and E. P. Borisenko, Paris: Presses de l’Ecole des Mines de Paris, 2000.
[3] D. Dumortier, M. Fontoynont, D. Heinemann, A. Hammer, J. A. Olseth, A. Skartveit, P. Ineichen, C. Reise. SatelLight, a www server which provides high quality daylight and solar radiation data for Western and Central Europe. CIE 24th Session, Warsaw, pp. 277281, 1999. http://www. satellight. com/core. htm.
[4] Surface meteorology and Solar Energy (SSE) Release 6.0 Methodology, draft version 1.005, 2008, http://eosweb. larc. nasa. gov/sse/. Accessible also through RETScreen software http://www. retscreen. net/.
[5] C. Rigollier, M. Lefevre, L. Wald, The method Heliosat2 for deriving shortwave solar radiation from satellite images, Solar Energy, 77 (2004) 159169, http://www. helioclim. org/radiation/. Accessible also through SoDa portal http://www. sodais. com/.
[6] M. Shri, T. Huld, T. Cebecauer, E. D. Dunlop, Geographic Aspects of Photovoltaics in Europe: Contribution of the PVGIS Web Site. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 1 (2008), in press, http://re. irc. ec. europa. eu/pvgis/.
[7] J. Remund, Quality of Meteonorm Version 6.0. Proceedings of 10th World Renewable Energy Conference, 1925th July 1008, Glasgow. (2008).
[8] H.G. Beyer, G. Czeplak, U. Terzenbach, L. Wald, Assessment of the method used to construct clearness index maps for the new European Solar Radiation Atlas (ESRA), Solar Energy, 61 (1997) 389397.
[9] P. Inneichen, Irradiance products validation and comparison against 13 ground stations, Working document of the MESoR project, 2008, http://www. mesor. net/.
The Danish King Christian VIII sent the first Galathea Expedition to the Far East in 18451847. The second Galathea Expedition took place 100 years later from 19501952, devoted primarily to ethnographic research. The present paper describes data collected on the Galathea III Expedition which was undertaken from August 2006 through April 2007. Figure 1 shows the Royal Danish Navy (RDN) vessel Vsdderen which was converted to a research platform for the expedition.
Figure 1: The RDN vessel Vsdderen served as the research vessel. On the right the SolData Instruments optics table can be seen with the Palmer Peninsula in the background. 
Data was integrated over 10minute intervals and recorded with a Cambridge Scientific CR10X data logger along with date, time and position information. Data was uploaded via a satellite link to the internet weekly. Thanks to excellent support from the Danish Department of Fisheries a complete data set is available for every single day of the 8 month expedition. This high level of data integrity is partly because the data logger was battery powered with a continuous trickle charge and could be independent of external power during for up to several days if necessary. The complete data set has been placed on the internet, and it is available in Excel format at the home page www. soldata. dk. Click on the link “Galathea data…”.
The results of this proposed model are compared with three other previously reported models which are believed to be applicable anywhere in the world and yield superior results for cloudy
conditions.
( і) Rietveld model [2]
Rietveld examined several published values of a and b, and noted that a is related linearly and b hyperbolically to the mean value of S / S0 such that
(9)
Substituting a and b in Eq (1), one
(10)
(ii) Ogelman et al model [2, 4]
Ogelman et al proposed the use of a correlation which relates the global solar radiation to S / S0 in a quadratic form as


(iv) Page model [8]
Page quoted a correlation which was believed to be suitable anywhere.
H c
The three models listed above were applied to the sunshine data in China. The calculated and measured values of average daily global radiation on the horizontal surface were compared using the following parameters.
In order to compute the deviation and compare the results of above models, this paper makes use of following parameters.










N
MBE = [I (Hestimated — Hmeasured )] / N
i=1
RMSE: Root Mean Square Error MABE: Mean Absolute Bias Error
MBE: Mean Bias Error
Hestimated : the estimated value of the global solar radiation Hmeasured: the measured value of the global radiation
He: the arithmetic mean value of the m estimated values of the global solar radiation
Hm : the arithmetic mean value of the m measured values of the global solar radiation N: the number of the data
The NWP uncertainties of the rainfall forecast is based on nonstationary, nonlinear and dynamic effects as stated by Todini [18]. The solar radiation forecast, also a function of the cloud cover, is expected to be subjected to the same underlying effects. Time series models which using the DWT are capable to interpret non stationary effects as stated in [27] and time series models using the autoregression (AR) models are capable to interpret dynamic effects [28]. Applying the AR model to the Haar and the biorthogonal DWT family, shows that the highest order biorthogonal wavelet transform of the Matlab™ tool set (bior6.8) :, has the best performance for the DWTAR MOS of the solar radiation correction as shown in [29]. Thus, it can be considered, that the DWT, which uses the highest order biorthogonal wavelet corrects best the dynamic and non stationary characteristics of statistic MOS model. Improving the correction of the non linear effects, the AR model is substituted by an ANN in the present article.
3.1 Statistical correction of the NWP model output
The idea of the proposed statistical correction method is based on the time series correction MOS as presented in [15] for local correction. The Kalman Filter is substituted by the ANNDWT (section 3.5). This method is utilized to estimate the correction of the solar radiation for the forecasted day, based on the recognized pattern of the NWP residual obtained from twelve previous days. The residual {sA} between the measured {H} and the forecasted daily mean solar radiation {HA} is estimated to correct the NWP output with this estimations {sA}. The forecasted {HA} and the measured {H} solar radiation time series are in a first step both decomposed in its subseries {HAs} and {Hs} by the DWT. For mx time scales one obtains s = 1 … (mx + 1) subsignals of residuals {sas} of the ARPS forecasts by the equation (2). For s = 1 … mx, the vectors {sas} are the subsignals which represent the details and for s = (mx + 1), the {sas} is the subsignal which represents the approximation of the ARPS residuals.
{sas} = {Has} — {Hs} (2)
The independent {sAS} subsignal vectors can also be obtained by the direct decomposition of the residual vector {sA}. Each of these time series is utilized for the training of its corresponding ANN, minimizing the squared error (eqn. 3) in a supervised learning process [30] to estimate the predictand sAs i for the day i, based on the characteristic pattern of the predictors, the sAS — values of (i1) to (ik) previous days. The symbol E stands for the energy of the error [30].
E _ 2 ^^ij (s As, i (sAs, (i1) ■■■ sAs, (ik)) — sAs, i) ^ тІП (3)
This function has a support length of thirteen sampled discrete values.
The forecast output of the NWP model is corrected with the sum of the (mx+1) estimated residual values ё а5,і.
In this work, we have only considered explanatory variables that can be derived from DEM. This makes the results more general and useful, since these variables are readily available elsewhere. Particularly, in this study we have evaluated the effect of the elevation and shadows cast as external variables in the kriging procedure. To deal with the seasonality of the solar radiation variability, the study is carried out for each month independently. The topographic characteristics were derived based on a 1 km spatial resolution DEM and the IDRISI GIS software [6] was used to carry out the regression analysis and the whole residual kriging procedure.
Regarding the elevation, it is clear that this variable is related to the atmospheric attenuation of the incoming solar radiation at the earth surface. The higher the elevation, the lower the atmospheric layer thickness and, therefore, the lower the atmospheric attenuation. Nevertheless, the importance of this explanatory variable depends on the elevation gradient. In the present study, this gradient is relative small (minimum elevation is 4 m and maximum 1212 m).
The other external explanatory variable that has been explored in this study is related to the shadow cast caused by the topography. The slope and aspect effect on the measured values can not be evaluated in this study since the stations are located on horizontal surfaces. The skyview factor influence, nevertheless, can be evaluated. Particularly, to assess this effect, a skyview factor can be used. This skyview factor can be defined as the ratio between radiation actually received by a planar surface to that received from the entire hemispheric radiating environment (without any obstruction). It represents the shadows cast by topographic features. Sky view factor can be derived from DEM using several available procedures. In the present study the skyview factor was computed using the GRASS GIS [6] environment in the following way: given a location, the maximum elevation angle of sky obstruction were computed for 36 directions around this location. Then, a continuous curve of angle values was computed using linear interpolation. The area under this continuous curve, normalized to the entire hemisphere is considered the skyview factor. Note that the value of the skyview factor varies from 1, when the whole sky is obscured, to 0, when no obstructions are present. In actual cases, skyview factor seldom are higher than 0.20.
In this study, a modified version of this skyview factor was used. Particularly, the skyview factor was computed just using the topographic information in the south direction. That is, the angles were only computed for the half of the sky view, from solar azimuth angles between 90° (eastern ) to 90° (western), instead of for the complete circumference (solar azimuth angles from 180° to 180°). We will call this new skyview factor semiskyview factor hereinafter. As for the complete circumference, the semiskyview factor may take values between 0 and 1. Actual values for all the stations used in this study ranged between 0 and 0.195. This new explanatory variable reveled to be much more effective in explaining the solar radiation variability than a complete skyview factor. The rationale behind this result is that, as it is known, in the northern hemisphere, the amount of radiation received from solar azimuth angles between below 90° and above 90° is scarce. Therefore, the topographic information associated with these angles is not relevant in explaining the solar radiation variability and just increases the noise of the explanatory variable, without adding new information. This new semiskyview factor was used as independent variable in the residual kriging procedure.
The Heliosat method is a technique of determining the global radiation at the ground by using data from a geostationary satellite.
In the Heliosat method it is assumed that the intensity of the visible solar irradiance which is scattered back to the satellite from the Earth and the atmosphere, behaves proportional to the atmospheric reflection. Also the isotropy of the atmospheric reflection is suppositional in order to disable the influence of the atmosphere’s heterogeneity on the reflection properties.
Due to the dominating dependence of the reflection on the cloudiness, it is feasible to derive an important quantity characterizing the degree of cloudiness existing within a solid angle from the radiation measured by the satellite. Via this quantity it is possible to deduce the transmission properties of the atmosphere and then determine the global irradiance.
Figure 1C shows a typical case of a clear day with the snow covering the Sierra Nevada Mountains. In these particular images, the method Heliosat2 can lead to subestimations, because the pixels are considered to be covered by clouds instead of by snow. In this work, a preprocessing method is applied to perform a detection of possible snow covers. Consisting in a comparison between the mountain pixels versus the surrounding ones, where a very low albedo in the surrounding if compared with the mountain is considered to belong to a clear day with snow covers. The method offers a success of near 85% in the detection of snow covers and it will be used. Error reduction due to this detection is discussed in the Results section.
1.2. Horizon calculation
Computational cost of the horizon calculation can be a problem when dealing with large areas; that issue has been addressed by developing an algorithm that reduces drastically the time spent in this process, without loosing much information about the actual horizon. The used DTM file contains around two million values of altitude, and the horizon for a specific point should be calculated by using all that information, but in this work most of the points are skipped and only around 70000 are utilized. The skipping step is defined as a function of the distance to the studied point, and a random contribution is also added. Finally, the horizon is only calculated for those selected points which altitude is higher than the studied point one, reducing the time consumption. Figure 3 shows the random and the selected points in the calculation of the horizon.
A comparison of the calculated horizon following the described algorithm and the real horizon (calculated by using the whole set of DTM points) appears in Figure 4, where it can be noted the little influence of the selection in the final calculation. Finally, for each day, the horizon will be overlaid with the sun path, obtaining the actual sunrise and sunset hours from the solar altitudes that will be introduced to the model.
Fig. 4 . Differences between real and used horizons. 
Daily irradiation values were calculated with and without inclusion of local effects; afterwards, the error was calculated in terms of Mean Bias Error (MBE) and Root Mean Square Error (RMSE), as a percentage of the mean measured value. Table 2 shows the obtained errors for the different stations and the three considered methods: A. Heliosat2 with no horizon and using a mean altitude for the pixel, B. Heliosat2 considering the calculated horizon and the altitude of the point and C. same as B but including the detection of clear days with snow covers.
Table 2. Obtained errors for the 14 stations and three different methods

Error reduction is detected for all the stations, and the importance of horizon effects appears to be larger than the effect of the snow covers. In fact, the consideration of snow covers is almost negligible in the case of the stations 13 and 14, as they are located in the plane area, outside the mountains, where snow is not usually present, but it can decrease up to 2% in the upper stations.
In addition, the station 9 is the one with a major improvement in the estimates, this is ought to the fact that is placed in a gully, having the larger horizon obstructions compared with the other stations. The mean RMSE considering the whole set of stations is reduced from 16% to 11%, verifying the good behaviour of the proposed methodology.
Once the performance of the method has been analyzed and in view of the satisfactory results provided in the global radiation estimates for the 14 stations located on a complex topography area, we generated a map with daily radiation values for the zone. Figure 5 shows a map created using this methodology.
Fig. 5 . Mapping of daily solar irradiation over a complex topography area. 
Results obtained using the proposed method showed a good agreement with the measured ones and it is interesting to note that this procedure can be applied under all kind of sky conditions. A proposed method to calculate the horizon of a point was tested with satisfactory results.
The RMSE was diminished for the whole set of stations with a mean reduction of 4.7%. MBE was also improved, with almost no over or underestimation for most of the stations. It has been observed that the improvement in the estimates is small in plane areas with almost no horizon obstructions, but becomes significant in the locations inside the mountain. The effect of the snow was also studied and it was pointed out that for the stations in the upper area it could lead to a reduction of the RMSE up to 2%, being smaller or nil in the lower stations. Finally, the method was employed over a wide area allowing the generation of an irradiation map on a complex terrain. If no additional information would be considered, the HELIOSAT 2 model would give a fixed value for the whole pixel, taking into account a mean altitude of the pixel and no horizon. On the other hand, using this methodology, the topographic variability is included in the model and a map of irradiation can be made in an easy way for a complex topography site. In future works, this method is intended to be extended to the estimation of different components of the radiation (i. e. Direct Radiation, Photosynthetically Active Radiation PAR); and also to be used together with Artificial Neural Networks.
This work has been financed by the Ministerio de Ciencia y Tecnologia of Spain (projects ENE/2004 0786C0301 and ENE200767849C0202). We would like to thank CIEMAT for the use of the satellite images.
[1] Dubayah and Van Katwijk, 1992. The topographic distribution of annual incoming solar radiation in the grand River basin. Geophys Res. Leter, 19, 22312234.
[2] Fu and Rich, 2000. A geometric solar radiation model and its applications in the agriculture and forestry. Proceedings of Second International Conference on Geospatial Information in Agriculture and Forestry, Lake Buena Vista.
[3] Diabate, L., Demarcq, H., MichaudRegas, N. and Wald, L., 1988. Estimating incident solar radiation at the surface from images of the Earth transmitted by geostationary satellites: the Heliosat Project. International Journal of Solar Energy 5, 261278.
[4] Dribssa, E., Cogliani, E., Lavagno, E. and Petrarca, S., 1999. A modification of the Heliosat method to improve its performance. Solar Energy 65, 369377.
[5] Zarzalejo, L. F., 2006. Irradiancia solar global horaria a partir de imagenes de satelite. Serie: Coleccion Documentos CIEMAT. Editorial CIEMAT, Madrid (Espana).
[6] ESRA, 2000. The European solar radiation atlas. Vol. 1: Fundamentals and maps. Editado por: Scharmer, K. and Reif, J. Les Presses de l’Ecole des Mines, Paris (France).