Evaluating the spatio-temporal performance of sky-imager-based solar irradiance analysis and forecasts

Clouds are the dominant source of small-scale variability in surface solar radiation and uncertainty in its prediction. However, the increasing share of solar energy in the worldwide electric power supply increases the need for accurate solar radiation forecasts. In this work, we present results of a very short term global horizontal irradiance (GHI) forecast experiment based on hemispheric sky images. A 2-month data set with images from one sky imager and high-resolution GHI measurements from 99 pyranometers distributed over 10 km by 12 km is used for validation. We developed a multi-step model and processed GHI forecasts up to 25 min with an update interval of 15 s. A cloud type classification is used to separate the time series into different cloud scenarios. Overall, the sky-imager-based forecasts do not outperform the reference persistence forecasts. Nevertheless, we find that analysis and forecast performance depends strongly on the predominant cloud conditions. Especially convective type clouds lead to high temporal and spatial GHI variability. For cumulus cloud conditions, the analysis error is found to be lower than that introduced by a single pyranometer if it is used representatively for the whole area in distances from the camera larger than 1–2 km. Moreover, forecast skill is much higher for these conditions compared to overcast or clear sky situations causing low GHI variability, which is easier to predict by persistence. In order to generalize the cloud-induced forecast error, we identify a variability threshold indicating conditions with positive forecast skill.


Introduction
As a result of worldwide growing photovoltaic electricity production, the energy sector is facing new challenges. One major issue is solar variability (Sayeef et al., 2012), which, on short timescales, is mainly caused by changes in cloud cover. With an increased share of solar power in the electricity grid, balancing power production and consumption is getting more and more challenging for power plant and grid operators. Consequently, flexibility options such as demandside management, backup capacities, inverter control, storage and strengthening of the grid are in the focus of research. In order to control and manage the flexibility options, the expected solar power production is an important piece of information. Although the large variety of cloud characteristics (opacity, motion, height, spatial distribution) makes these cloud-induced fluctuations difficult to predict, solar irradiance forecasting techniques have been successfully developed (a comprehensive overview is given in Inman et al., 2013, andHeinemann, 2012). The spectrum comprises numerical weather models (NWPs) (Perez et al., 2013), satellite-based forecasts using cloud motion vectors (Kühnert et al., 2013;Lorenz et al., 2004;Hammer et al., 1999), statistical methods based on machine learning (Wolff et al., 2013) and time series analysis (Reikard, 2009) predominantly developed for intra-day and day-ahead forecasts. For very short term forecasts with time horizons of up to 30 min, both NWP and satellite image-based models lack spatial and temporal resolution regarding cloud-induced small-scale variability (Inman et al., 2013).
With the aim to fill this gap of local high-resolution and very short term forecasts, recent research has made use of ground-based (whole/all/total) sky imagers. Sky imagers have been used for years for monitoring cloud cover characteristics (Pfister et al., 2003;Long et al., 2006;Cazorla et al., 2008) and aerosol properties Cazorla, 2010).
By analysing distribution, movement and optical properties of clouds, the incoming solar irradiance can be forecasted. Most of the cameras that are used are equipped with fisheye lenses capturing the whole visual sky. Evidently, the possible spatial coverage of irradiance analysis and consequently the temporal forecast horizon are variable and depend on daytime and, accordingly, sun position, cloud distribution (type and altitude) and cloud motion (speed and direction). Two types of forecast experiments have been reported in recent work. Point forecasts only predict the occlusion of the sun with clouds and therefore can only process forecasts for the location of the camera (e.g. West et al., 2014). Area forecasts, on the other hand, incorporate cloud base height in order to calculate the location of clouds and shadows on the ground (e.g. Yang et al., 2014).
This work presents results of short-term area forecast experiments. We developed and applied a multi-step model on a large data set of sky images and processed forecasts up to 25 min ahead for 99 locations distributed in the surrounding area. Global horizontal irradiance (GHI) measurements at these locations are used for evaluating the forecast performance.
The sky images, the ceilometer-based cloud base height measurements and the pyranometer data used in this study were collected during the High Definition of Clouds and Precipitation for advancing Climate Prediction (HD(CP) 2 ) Observational Prototype Experiment (HOPE) in spring 2013. The data set provides both a high spatial density of solar radiation measurements and the necessary temporal resolution of pyranometer and ceilometer data as well as sky images.
This work focuses on the investigation of the performance under different cloud conditions. A cloud classification scheme is used to categorize the data set into seven different cloud conditions. This differentiation is helpful in the comparison with reference models, such as persistence, which are almost perfect for short forecast horizons and low variability in cloud cover. The spatial distribution of pyranometers is used to identify differences in performance for locations distant from the camera. This analysis is helpful for investigating the usefulness of sky-imager-based irradiance field analysis instead of using several expensive pyranometers. As our cloud detection scheme only provides two binary states (sky/cloud) and no cloud transmissivity infor-mation, the irradiance retrieval for GHI has weaknesses in the case of deviations from this simplification. Therefore, we also evaluate binary forecasts in order to identify the contribution to the overall forecast error caused by the irradiance retrieval.

Experimental setup and database
The data sets used in this short-term forecast experiment have been collected during the HOPE measurement campaign in 2013. For this work, data from a network of 99 irradiance sensors, 1 ceilometer and 1 sky imager were used (Fig. 1). The following subsections give a short description of the data sets used. Here, measurements from 1 April to 31 May were used. The measurement site is located in Jülich, Germany. The area is rather flat and surrounded by two large open-cast lignite mines (Fig. 1).

Sky imager
A sky imager developed at the GEOMAR Helmholtz Centre for Ocean Research (Kalisch and Macke, 2008) was used for continuous sky observations. The imager was part of the LACROS supersite within the HOPE measurement campaign; see Madhavan et al. (2015) for the location and details. The Canon digital CCD camera equipped with a Raynox fisheye lens realized a field of view of 183 • . The hemispheric sky images with 2592 pixel × 1744 pixel resolution were sampled at a rate of 15 s.

Irradiance sensor network
A irradiance measurement network with 99 pyranometer stations was set up around Jülich, Germany, on an area of 10 km×12 km. Each station was equipped with an EKO ML-020VM photodiode pyranometer. The 10 bit data-logging system was synchronized with the GPS time. The irradiance was measured with 10 Hz resolution and was averaged to 1 Hz. Maintenance and checks for cleaning and tilt were performed on a weekly basis. Based on the maintenance protocol, data are provided with quality flags ("good", "okay, but sometimes spurious", "bad or ignore completely" or "missing or no observations"). Madhavan et al. (2015) give a description of the pyranometer network within the HOPE campaign, details of the hardware, quality flags and an investigation of measurement uncertainties.

Additional data
Further information about cloud base height and sun position for ray tracing and following cloud shadow mapping is needed for processing of sky images for solar irradiance area forecasts. Clear sky irradiance information is necessary for reference cloud-free sky conditions and irradiance retrieval. Information about cloud base height is retrieved from a Jenoptik CHM15k-x ceilometer that was located next to the sky imager. Ceilometers are recognized by the WMO as the most accurate, reliable and efficient means of measuring cloud base height from the ground when compared with alternative equipment (World Meteorological Organization, 2008). One measurement was done every 20 s. As a ceilometer provides only point measurements, the median of the last 30 measurements was used in order to smooth the signal. Although multi-layer cloud height information is available, only lower-level cloud height was used, because the sky imager algorithm that was used does not yet support multilayer clouds. Clear sky irradiance is estimated with the clear sky model of Dumortier (Fontoynont et al., 1998) and turbidity values according to Bourges (1992) and Dumortier (1998). The model is also used in internal operational services for photovoltaic plant energy yield monitoring and evaluated continuously at more than 100 sites in Germany within our group. For a recent independent validation of the model, see Ineichen (2013).
The solar zenith and azimuth angle are calculated with the solar geometry2 (SG2) algorithm (Blanc and Wald, 2011).

Methods
In order to determine and predict the surface patterns of global horizontal irradiance distribution from sky images, several preprocessing steps on the image have to be done. This section subdivides the processing chain into image analysis, irradiance analysis and irradiance forecasting. Figure 2 gives an overview of the workflow, which is described in  more detail in the following sections. Furthermore, a cloud classification scheme is introduced.

Cloud detection
To identify clouds, we apply a binary classification (cloud or sky) of each image pixel (Fig. 2). As a consequence, we do not account for variations in cloud optical thickness (from thin semi-transparent to thick opaque). Here, we use the concept of the red-blue ratio (RBR), first developed by Scripps Institution of Oceanography (Johnson et al., 1989(Johnson et al., , 1991Shields et al., 1998). RBR is the ratio between the red and blue channel of the image. The RBR indicates whether the scattered light comes from a cloud (value close to 1) or blue sky (value 1). Based on an empirically determined threshold of RBR = 0.82, each pixel is classified as cloudy or non-cloudy. Cloud detection based on RBR has been used in several sky-imager-based forecast applications (e.g. Chow et al., 2011;Yang et al., 2014;Urquhart et al., 2015). The RBR is not homogeneously distributed over the whole field of view for the same sky conditions. RBR has an angular dependency (Pfister et al., 2003) and the area close to the sun (circumsolar region) is affected by the bright sun (RBR ≈ 1). Consequently, misclassifications are likely when one single global threshold is applied to the image. Optically dense clouds, which appear quite dark in the centre of their base (West et al., 2014), are another source of errors. Here, the RBR is very low and clouds can be misclassified as sky.
To overcome these disadvantages, we correct the RBR with a set of clear sky images similar to Chow et al. (2011) and Shields et al. (2009). Here, the clear sky library (CSL) contains RBR images from one clear sky day (4 May 2013) of the measurement period. The database serves as a reference for clear sky conditions (see Fig. 3 for an example). The reference image (Fig. 3c) is selected by calculating the angular distance of the currently sun position from the references and choosing the closest one.
A modified RBR (R mod ) is given for each pixel at the image position i, j by the following equation: It first accounts for the difficult circumsolar area. We defined the grade of saturation (S ∈ [0, 1]) as the average pixel intensity in the disc up to an angular distance of 5 • to the centre of the sunspot. A value of S = 1 would correspond to a completely saturated sun area (each pixel's intensity I = 255). Weighted by the grade of saturation S multiplied with a constant factor a, we subtract the clear sky RBR (R CSL ) from the original RBR (R orig ). Moreover, a correction based on the pixel intensity I 1 multiplied with a second constant factor b and clear sky RBR R CSL is applied, which increases RBR in the case of dark clouds and decreases RBR in the case of bright clouds (Fig. 3).
The values for a = 0.1, b = 0.0018 and R thres = 0.82 that discriminates clouds and sky were determined empirically on a test data set of 40 images with different sky conditions. They were adjusted with the aim of achieving good results for all possible sky conditions, including in particular thick and dark clouds, semi-transparent cirrus clouds and clear sky.
Note that the used CSL introduces errors on days where solar zenith and azimuth angles deviate from the reference day. Moreover, days with different atmospheric conditions (aerosol load, scattered light) from those of the reference day will lead to errors not quantified in the RBR corrections (Ghonima et al., 2012).
The proposed approach aims to reduce the mentioned misclassifications in the circumsolar area and in the case of thick and dark clouds.

Camera calibration and image undistortion
In order to project an image pixel from a fisheye lens image in geometric coordinates, two types of parameters are needed. First, intrinsic parameters describe the geometric distortion introduced by the optics used to project 2-D image pixel points onto a unit sphere. Next, extrinsic parameters describe the transformation from the unit sphere in the real world. This can be expressed with a rotation matrix accounting for orientation errors.
The intrinsic parameters are determined by a calibration of the fisheye lens following Scaramuzza (2014). The method detects straight known lines on photographs of a checkerboard and retrieves the distortion (Scaramuzza et al., 2006). Assuming a radial symmetrical distortion, a fifth-degree polynomial function with coefficients k in Eq. (2) is fitted on the detected data points. It assigns each pixel's distance r from the centre of the image to the corresponding incidence angle θ .
Extrinsic parameters are estimated by a visual comparison of the reprojected sun position (azimuth and zenith angle) to image coordinates and their visual appearance in the image. In this case, we assume a perfectly horizontally mounted camera and define a rotation matrix which rotates the top of the image to geographic north. Equation (2) and the rotation matrix are used for undistorting the image.

Image masking
Static artificial objects in the field of view are masked out. Furthermore, the field of view had been limited to an incidence angle of 80 • in order to reduce perspective errors at large incidence angles.

Cloud mapping
Determination of the three-dimensional position of a cloudy pixel with incidence angle θ i,j needs the clouds' base height, h, as a further input. The geometric distance of a single pixel d i,j from the position of the camera is calculated with The clouds' position is then calculated using the measured cloud base height from the ceilometer and the pixels' incidence θ i,j and azimuth angles φ i, j retrieved from the camera calibration (Sect. 3.1.2).

Shadow mapping
With the information about current sun position (azimuth angle φ sun and zenith angle θ sun ) and cloud base height h, a sun ray tracing is applied to map the cloud layer as a shadow layer on the ground. Eq. (4) gives the basic formula for calculating the horizontal distance d of a cloud's shadow on the ground from the camera.
A topographically flat surface was assumed. The application of a more realistic topography could lead to better results, but considering the almost flat surface at the measurement site, the introduced error will be small related to other error sources.

Gridding
In order to analyse the cloud shadow field at the location of the pyranometer stations, image pixels are mapped on a regular grid of 20 km×20 km with a resolution of 20 m. One has to consider that, depending on cloud base height, the raw image pixel resolution is higher than the final grid resolution in the image centre and lower in the outer region. In the former case the central pixel is used, while nearest-neighbour interpolation is used for interpolating in regions where the image resolution is below the grid resolution. Afterwards, a Gaussian filter with σ = 3 is applied on the gridded binary data to smooth cloud edges. The effect of smoothed cloud edges is illustrated in the time series of the forecast example in Fig. 6.

Cloud classification
We apply a cloud classification algorithm in order to classify each image instance in different cloud condition categories. This data separation is used for evaluating forecast performance under different cloud conditions. A review of existing cloud detection and classification methodologies is given by Tapakis and Charalambides (2013). Here, we modified the cloud classification scheme introduced by Heinle et al. (2010). The modified classification algorithm uses "support vector classification" (SVC) as it outperforms k-nearest neighbours (k-NN) in our application. We also extended the number of features to 16 image-based features and trained on a data set of 600 images manually classified into seven categories. The seven categories are meteorologically justified according to Heinle et al. (2010) -cumulus (Cu); Three of the additional features include image texture properties derived from the grey-level co-occurrence matrix (GLCM) and defined by (Haralick et al., 1973). The angular second-moment (ASM) feature is a further measure for homogeneity, correlation is a measure of grey-tone linear dependencies and dissimilarity is a measure that defines the variation of grey level pairs in an image (Gebejes and Huertas, 2013). Furthermore, as possible informative features, the ratio of the number of saturated pixels (all channels have intensities of 255) to all non-masked pixels, the average pixel intensity, and the average RBR value are used as input.
The classification was validated by 10-fold crossvalidation with an accuracy of 92 %. Note that only the dominant cloud type according to the classification model is determined.

Irradiance analysis
The transformation from surface shadow fields to irradiance fields is based on past records of clear sky indices measured at each pyranometer station. The clear sky index k * is the ratio of measured global horizontal irradiance GHI meas and a clear sky reference value GHI clear (Eq. 5).
A typical histogram of measured k * has two peaks for overcast and clear sky conditions. Here, this information is used for the irradiance retrieval for the two states: shadow and no shadow (see Fig. 4).
We calculate the histogram for each station for the previous 30 min to account for changing atmospheric conditions. The method takes the global peak below for k * < 0.5 for shadow state and k * > 0.9 for no shadow. We decided to use 100 bins for 0.2 ≤ k * ≤ 1.4. If no peaks can be determined (in the case of homogeneous irradiance conditions in the previous 30 min), default values of k * hist = 0.4 and k * hist = 1.0, respectively, have been assigned for the two states. See Sect. 2.3 for the clear sky irradiance model used. The corresponding GHI can then be calculated with The spatial smoothing (introduced in Sect. 3.1.6) of the shadow field leads to smoothed cloud shadow edges. This could be regarded as more realistic for transitions from nonshaded to completely shaded conditions. Despite the adaptation to the situation of the previous 30 min, irradiance levels for clear sky (diffuse + direct irradiance) and cloudy sky (diffuse irradiance only) may deviate from measured values (see Fig. 6). A more realistic retrieval can benefit from the estimation of direct and diffuse components, e.g. by considering additional image features with machine learning (Schmidt et al., 2015).

Cloud motion
The fundamental information needed for cloud forecasts is cloud movement and cloud transformation. As the transformation (development and dissolution) of clouds is a very complex task, our algorithm does not yet account for that. As a consequence, predicted cloud scenes are the result of a translation of the current analysed cloud scene. Cloud movement is determined by applying the optical flow algorithm available in OpenCV (Open source Computer Vision Library, http://opencv.org). Optical flow calculations have been used in other sky imager applications by West et al. (2014) and Wood-Bradley et al. (2012). The first step is to determine good features to track in the image (Shi and Tomasi, 1994). These objects -mostly found on strong gradients like cloud edges -serve as input for the Lucas-Kanade tracking algorithm (Lucas and Kanade, 1981;Bouguet, 2001). The algorithm yields cloud motion vectors (CMVs). In this study, new features are determined every 2 min as old features change too much or move out of the visible image. The algorithm is applied to the original greyscale image, where artificial objects are masked out.
Similar to the cloud shadow projection, each single CMV is transformed to the underlying metric grid by projecting the image coordinates of the vectors initial and terminal point. (Sect. 3.1.6). Figure 5 shows an example of the transformation from the circular fisheye image to the grid. This scene illustrates the rectification of the CMVs, which is important for quality control and averaging to a global CMV.
To increase the CMV quality, we first mask out the circumsolar area in the feature detection step, as its brightness disturbs the algorithm. Next, we apply a quality control. Initial CMVs are flagged as invalid if their speed is lower than 0.2 m s −1 to avoid tracking artificial objects in the image. If clouds are moving at a speed below that threshold and all CMVs are flagged as invalid, a persistent cloud mask is assumed. If sudden changes in direction and speed (changes in cloud speed > 2 m s −1 ) of follow-up vectors occur, which can happen if brightness in the image changes rapidly, these vectors are removed. The final CMVs are then averaged to one global vector which determines the principal movement of the cloud scene for the forecast. Due to recalculation of CMV positions every 2 min, discontinuities in the global CMV may occur. As these are unlikely to represent true changes in cloud motion, which is rather inert, the last four global vectors are also averaged in time. Furthermore, each change in the average CMV will affect the forecasted cloud distribution and the irradiance forecast. An approach that uses the uncertainty in cloud motion for an estimation of uncertainty irradiance forecasts is being developed.

Solar irradiance prediction
Irradiance forecasts are calculated for each pyranometer station with a horizon of a maximum of 1500 s and a resolution of 1 s. A forecast run is computed for each image (every 15 s). This is done by advecting the "frozen" cloud field with the global CMV (Sect. 3.4.1) and calculating the surface shadow maps (Sect. 3.1.6) and irradiance maps (Sect. 3.3). We considered the varying sun position in the 25 min forecast horizon by computing its position for each forecast step. Afterwards, the irradiance forecast at each pyranometer station is retrieved.
As an example, Fig. 6 illustrates a forecast run for a pyranometer in the north of the sensor arrangement. The thick coloured line represents the forecast path along the opposite direction of the global CMV, indicating a mean cloud motion from a southern direction. Here, cloud speed is low enough for processing a full forecast up to 25 min ahead for this location. The binary pattern of the forecast is a result of the measured GHI in the previous 30 min (Sect. 3.3). Although the binary pattern is represented in the forecast time series, slight smoothing at the cloud edges is pronounced as well.

Concept of evaluation
In order to evaluate the forecast data set we focused on two main aspects: 1. How accurate is the sky-imager-based analysis during different cloud conditions and with respect to distance from the camera?
2. How accurate are sky-imager-based forecasts in different cloud conditions especially compared to persistence?
For answering the first question, we analyse mean bias error (MBE) and root mean square error (RMSE) spatial distribution (see Sect. 3.5.2) for each cloud class. By sorting the stations by distance from the camera position, we were able to compare the analysis (forecast lead time = 0) error to the error introduced if a single pyranometer at the location of the camera was representative of the whole area. The second question is investigated by evaluating the forecast performance depending on the forecast lead time. As a reference forecast we use persistence. Persistence forecasts account for changing sun angles but assume no change in cloudiness described by a constant clear sky index k * : We keep the raw resolution of 1 s for the persistence definition. As a consequence, persistence forecasts have no initial error, but it increases with time. To evaluate performance in different cloud conditions the data set is separated in the seven analysed classes. Forecast error and skill are then calculated for each of the classes (for definition of error metrics see Sect. 3.5.2).
During the processing chain, several assumptions and simplifications are made which contribute to final analysis and forecast errors. One error source is the irradiance retrieval (Sect. 3.3) based on binary cloud maps processed before. Particularly, cloud irradiance enhancements due to reflections at cloud edges, irradiance reductions due to semitransparent clouds, and changes in diffuse irradiance levels due to a changing cloud distribution cannot be accurately addressed with the proposed methods. Therefore, we evaluate the ability of the forecast to distinguish between the two states (sunny and cloudy) by introducing a threshold of k * = 0.7. The time series in Fig. 6 illustrate the error introduced by GHI values deviating from the average.

Data selection
To analyse the performance of our forecasting system, we had to attend to data availability and quality. The total number of processed forecast runs is 138 912, corresponding to the number of available images processed for sun elevations greater than 10 • . The number of forecasts used for the evaluation is reduced by non-available measurements or forecasts. We decided to use only measurements which were flagged by the data provider as "good". As stations were maintained once a week and quality flags were given for the whole week, data gaps most of the time cover a whole week (Madhavan et al., 2015). As a consequence, a reduced subset of 50 stations with at least 70 % of the maximal possible number of measurements available was used when comparing performance for different stations (Sect. 4.2). Forecast availability for each location is limited by several factors. The size of the underlying grid, the field of view of the camera (we masked out the area beyond 80 • lens angle of incidence), current cloud base height, cloud speed and direction, and the sun position lead to a varying maximum forecast horizon. Figure 7 illustrates the data availability for the evaluation depending on the forecast horizon as well as the spatial distribution for a forecast horizon of 10 min.

Error metrics
For measuring the accuracy and performance of the forecast system we used mean bias error (MBE, Eq. 8), root mean square error (RMSE, Eq. 9), forecast skill (FS, Eq. 10) and accuracy (ACC, Eq. 11) in this analysis.
MBE is the average deviation of the forecast or analysis y from the measurement x: where subscript i refers to a single forecast or analysis y or measurement x. By definition, RMSE is given by FS is given by A positive FS means that the sky-imager-based forecast outperforms persistence (Eq. 7). ACC is used for measuring the ratio of the number of correctly predicted states (sunny and cloudy) in all instances (Metz, 1978): where TS represents "true sunny", TC "true cloudy", FS "false sunny" and FC "false cloudy". For example, a forecast is true sunny if measured and predicted k * are > 0.7. A forecast is false sunny if measured k * > 0.7 and predicted k * ≤ 0.7. These error metrics are calculated for each station and forecast horizon separately. Table 1 shows the results of the cloud classification. This table gives an overview of the predominant cloud conditions Table 1. Results of cloud type classification. The first two rows give the frequency of occurrence and its share of the total number of analysed images. k * and V represent the corresponding k * statistics (average, variability) of station 33, which had the highest number of valid measurements during the period. Average and standard deviation are given for the measured cloud base height (CBH), the analysed hemispheric cloud coverage (CC) and cloud motion speed (Cspd). 2200 ± 1500 1600 ± 800 2700 ± 1400 1300 ± 700 3400 ± 2500 1300 ± 800 NaN CC (%) 55 ± 32 95 ± 11 59 ± 31 99 ± 3 64 ± 34 100 ± 1.0 5 ± 8 Cspd (m s −1 ) 9.6 ± 6.0 11.6 ± 5.5 10.8 ± 6.6 7.8 ± 5.2 7.1 ± 6.5 7.0 ± 5.7 NaN and their characteristics mainly affecting the surface solar irradiance and its variability in space and time. The GHI statistics are calculated for a single station which had the highest availability. Variability V is defined according to Marquez and Coimbra (2012),

Cloud type distribution
with the number of images in each class N and t set to 5 min. As expected, the convective cloud type classes Cu, Ac/Cc and Sc have the highest variability. Sc, in contrast to Cu and Ac/Cc, has a high cloud coverage and therefore causes a lower average clear sky index. St/As causes low variability close to that of clear sky. The non-intuitive variability for scenes classified as clear sky can be traced back to scenes not fully clear but predominantly clear (not shown here). Cb/Ns and Ci/Cs also cause low variability compared to the first three classes. Cu, Sc, and Ac/Cc occurred about 37 % of the time, while low-variability classes except for clear sky occurred 53 % of the time; 10 % were clear. No big differences can be seen in the cloud motion statistics for all non-clear situations. An average cloud speed of 10 m s −1 has the effect that a cloud will move across the domain in about 33 min from east to west or from north to south. This number illustrates one aspect of the limits to the forecast horizon. For further evaluation purposes we group the convective type clouds Cu, Sc, Ac/Cc together into a new category, "heterogeneous" clouds, while the cloud types St/As, Ci/Cs and the clear sky situations build the category "homogeneous" clouds, as they cause rather low variability in surface solar irradiance.

Irradiance analysis accuracy
Irradiance analysis is evaluated depending on the distance of the stations from the camera and according to the different cloud classes.
The spatial distribution of the mean bias error MBE of the GHI analysis (forecast lead time t = 0) is shown in Fig. 8 for Cu and clear sky situations. Here, the MBE is given for each of the stations of the subset introduced in Sect. 3.5.1. The MBE distribution for Cu shows a negative MBE of about −80 W m −2 for stations close to the camera, increasing with distance to positive values around 70 W m −2 . A similar overestimation for stations close to the camera can also be found for Ac/Cc, Sc and Ci/Cs (not shown here). This is probably explained by the fact that the correction of RBR (Sect. 3.1.1) is too strong in the circumsolar region (affecting these locations) in the presence of the mentioned clouds. As a result, clouds in the circumsolar region are possibly too often misclassified as clear sky and surface irradiance is overestimated in the area around the camera. This phenomenon is not found for St and Ns/Cb situations dominated by (dark) overcast sky not affected by the correction. Moreover, the clear sky MBE distribution in Fig. 8 shows that the correction performs on average well in clear sky situations as no significant MBE for stations surrounding the camera is present.
An increasing tendency in MBE with distance from the camera is also found for the aforementioned types Ac/Cc, Sc, and Ci/Cs, while it is not present during clear sky or overcast stratus clouds (only clear sky is shown here). A similar tendency is identified for RMSE in Fig. 9 showing the cumulus conditions again. However, even if there is a large (absolute) MBE for stations close to the camera, no increased RMSE is present. This makes it clear that the main contribution to RMSE is the standard deviation of the analysis error and not the MBE.
Several possible explanations for these results can be identified. First, the perspective error increases with distance from the centre of the image. As a result, convective clouds with vertical extent (mostly cumulus), which are interpreted as horizontally flat in our scheme, are projected incorrectly if they are seen from their side near the edge of the field of view. Among other things, this leads to an underestimation of gaps in the cloud layer contributing to RMSE and a positive MBE. Furthermore, uncertainties in cloud base height lead to higher errors in the shadow mapping the more distant the clouds are (can be derived from Eq. 4). Moreover, cloud base height was measured at the position of the camera. Therefore, its representativeness for locations more distant is reduced depending on the cloud situation. This displacement of shadow patterns contributes mainly to RMSE. As the temporal and spatial resolution of 1 Hz and 20 m, respectively, is quite high, double penalties in the case of small cumulus or broken cloud layers are likely (Gilleland et al., 2009) and increase RMSE even more. Furthermore, the pixel resolution is reduced for larger lens incidence angles. This leads to a reduced spatial resolution for locations distant from the camera, which affects the accuracy of the camera-based irradiance analysis.
Moreover, Fig. 9 shows the RMSE that is introduced if a single pyranometer is used representatively for the whole area. It is assumed that the pyranometer closest to the camera is the reference sensor, and the RMSEs of its measure- ments compared to the remaining pyranometers are calculated. As expected, the error increases very fast with distance as the cross-correlation between the sensor pairs is reduced especially in conditions with high GHI variability. It can be stated that the "break-even" distance, where the sky-imagerbased irradiance analysis outperforms a single sensor spatial extrapolation for these highly variable cloud conditions, is found at a distance between 1 and 2 km from the camera. For other convective cloud types a distance of 2-3 km for Sc and Ac/Cc and 6 km for Ns/Cb is found. In the case of St/As and Ci/Cs clouds and in clear sky conditions, the analysis error is always larger due to the high sensor pair correlation in these less variable situations. Figure 10 shows the RMSE of the sky imager forecast and its corresponding persistence forecasts depending on the forecast horizon for the different cloud conditions. Here, the average RMSE of all evaluated pyranometer stations is shown. As expected, the overall forecast error is higher in situations with more variability in cloud cover and therefore in surface solar irradiance. forecast horizon depending on cloud base height, cloud speed and direction, sun position, and the location of each individual pyranometer. The limits of the underlying domain are a fixed constraint. A detailed analysis (not shown here) revealed that forecast runs with large forecast horizons have a lower RMSE value, probably caused by lower cloud speed resulting in less GHI variability. This is maximally expressed in the cumulus cloud class. From Fig. 10 it can also be seen that the difference between sky imager forecast RMSE and persistence RMSE is much more pronounced for the stratiform cloud types St/As and Ci/Cs. Figure 11 underlines this result as it shows the forecast skill FS for the categories of "homogeneous" and "heterogeneous" clouds defined in Sect. 4.1. While the skyimager-based forecasts are able to outperform persistence under heterogeneous conditions for at least a few stations after about 10 min, the forecast skill under homogeneous conditions is much worse.

Forecast performance
In order to determine the influence of the irradiance retrieval based on the binary cloud/sky decision on the forecast error, binary forecasts with the accuracy metric (Eq. 11) are also evaluated. It is expected that the forecast performance is higher if only the two main states, sun and shadow, are considered, as GHI forecast errors are introduced into our algorithm during conditions in which the measured GHI distribution deviates strongly from our simplified binary model from Sect. 3.3.
The evaluated accuracy for both sky-imager-based forecasts and persistence is shown in Fig. 12. Obviously, the accuracy for stratus (St) and nimbostratus/cumulonimbus (Ns/Cb) clouds is very high for both forecasts, indicating that irradiance is constantly lower than k * = 0.7 and that this state is predicted accurately. The clear sky case can also be predicted with an accuracy of more than 90 %. Forecasts in times of semi-transparent cirrus/cirrostratus (Ci/Cs) clouds still have a low skill, indicating that misclassifications in  cloud detection are preceding the irradiance retrieval. Moreover, from our experience we know that the RBR threshold used for cloud detection is not able to distinguish well between thin cirrus clouds and blue sky. Stratocumulus (Sc) also achieves high accuracy larger than 80 % for the whole horizon. In contrast to RMSE, forecast accuracy can outperform persistence from a forecast lead time of 3-4 min on. This indicates that GHI forecast errors for Sc conditions, which are dominated by high cloud coverage, can be attributed a considerable amount to irradiance retrieval errors. For Cu and Ac/Cc, only low improvements can be stated compared to RMSE for continuous forecast verification. As a consequence, other error sources like spatial mismatch dominate the error in this case. Furthermore, this result is of interest for applications focusing on binary events, which is the case for concentrated solar power (CSP), dealing mainly with variations in direct normal irradiance.
In the previous section, we identified a spatial "breakeven" distance for GHI analysis for different cloud classes. Such a "break-even" point can also be identified for an increased temporal GHI variability. Figure 13 displays the RMSE (based on clear sky index k * ) of a 10 min forecast depending on the prevailing variability (Eq. 12) for 10 min k * increments. Here, no distinction in cloud classes is made. RMSE and variability are calculated for short moving time windows of 25 min each. The time step between two time windows is 1 min, resulting in an overlapping database. The lines in Fig. 13 represent the average values of each bin. With that definition, persistence forecast errors fall on the diagonal line of the plot. This analysis summarizes the previous investigations of forecast errors under specific cloud conditions. In situations of low GHI variability there is only low forecast skill, which is increased with increased GHI variability. As these situations are much less frequent (see dashdotted line in Fig. 13), this skill is not visible in the average error statistics. Therefore, the strength of deterministic sky-imager-based forecasts for changing cloud conditions is made visible here. It can be stated that there is a specific value of 0.3 to 0.4 k * variability in the given case where a skyimager-based forecast can have skill against persistence. For other locations in the area covered, the results are similar, with a slightly different "break-even" value (not shown here).

Conclusions
A short-term GHI forecast experiment based on hemispheric images of the visible sky was conducted on a large data set of spatially distributed pyranometers. A processing chain comprising cloud detection, cloud motion, cloud and shadow mapping, and irradiance retrieval was proposed and applied to sky images retrieved during April and May 2013. The results show that the forecast performance and the benefit of sky-imager-based forecasts vary a lot depending on the given cloud conditions. A cloud classification scheme was used to determine seven different cloud conditions in order to evaluate the performance in more detail. Even though the overall forecast performance is quite low compared to persistence, one has to point out that the skill increases in heterogeneous cloud conditions, leading to increased variability in surface solar irradiance.
The evaluation of the GHI analysis shows the potential of sky imagers for areal irradiance monitoring. The study shows that the sky imager retrieval for distances of more than 1-2 km from the camera under cumulus cloud conditions outperforms a single pyranometer representing the spatial irradiance distribution. This value is increased for stratocumulus and altocumulus/cirrocumulus to 2-3 km and for nimbostratus/cumulonimbus to 6 km. As setting up a pyranometer network with a comparable density or resolution to a sky imager is more expensive than a camera, a camera-based areal irradiance monitoring can be beneficial. The impact of irradiance retrieval on forecast errors is shown by comparing standard GHI forecast errors to a binary forecast evaluation. This indicates potential for improvements by enhancing the irradiance retrieval. We also see potential to improve the model in the handling of multi-layer clouds (accurate cloud base height and cloud motion) and in a better cloud detection (assigning transmissivity to each cloud pixel instead of simple binary states). With these improvements in image processing and forecast methods, forecast error will be reduced continuously in the future. Skyimager-based analysis and forecast methods can then contribute to site monitoring and short-term forecasting, especially in highly variable cloud conditions.