Cloud Detection and Prediction with All Sky Cameras

Atmospheric monitoring is a field of special importance for astroparticle physics, especially for Imaging Atmospheric Cherenkov Telescopes (IACTs) as clouds will absorb and scatter the Cherenkov photons of air showers. Conventional tools used for atmospheric monitoring (e.g. LIDAR) are very expensive and monitor only a small part of the sky at once. Therefore, they are not suitable to perform a wide scan of the sky which is necessary to detect clouds in advance. This article gives a short overview about a method that uses an all sky camera with a 180 ◦ field of view to identify the cloud distribution by measuring the absorption of star light. It can be used to assign a sky quality rating to single spots, arbitrary regions or the whole sky at once within a 1 min exposure time. A cloud map can be created from the available data that can be used to determine shape and dimension of clouds and to predict their movement. The resulting data can be used by a scheduling algorithm or the operating crew to point the telescope to a different source before the current source gets covered by clouds. The all sky cameras used so far are located on La Palma at the observatory Roque de los Muchachos close to the telescopes FACT and MAGIC and the planned northern CTA site.


Introduction
A precise cloud forecast would be very valuable for the operation of IACTs.If the time was known when a source will cover up, the telescope could switch to another source in advance and continue taking data.The detection of clouds at night using a camera is difficult as they do not emit light on their own and are as dark as the night sky unless they refract moon light and appear as white objects in the camera image.The number of visible stars in the image is therefore a good approximation of the cloud coverage as their appearance does not depend on ambient light or wind conditions.
Regions with visible stars can be marked as not cloudy but a cloud will not always cover stars completely.A partial absorption of light will not be recognized if stars have only two states: 'visible' or 'not visible'.This effect makes it very difficult to reconstruct the smooth shape of a cloud or detect thin clouds in the first place.

Image Analysis
Initially, a star catalog is used to find the positions in the image where a star is expected to show.The actual detection of a star is done by convolving the image with a Laplacian of Gaussian filter [3].This filter is ideal for star detection in long exposure images as it is rotational invariant (anisotropic) and robust against image noise.Figure 1a shows a plot of the filter response R of stars and their e-mail: jan.adam@tu-dortmund.demagnitude during a starry night.The filter response is correlated to the star magnitude.This allows to use an exponential function: to define an upper and a lower limit for the response of a star (green and red line).The parameters of the green curve were choses such that all stars lie above this curve in a starry night (cf.Fig. 1a).
Stars above this threshold are unaffected by clouds while stars below indicate that a cloud absorbed a portion of their light.
A cloudy night is depicted in Figure 1b and defines the lower threshold for the filter response.The parameter of the red curve were chosen such that the majority of stars lies below this curve in a cloudy night.
Stars below the red curve are not visible at all.Combining both thresholds defines the visibility v of a star as its relative position between upper and lower threshold in logarithmic space.For example: in Figure 1a the star with magnitude ≈ 0 lies above the green curve, thus its visibility is 1.In Figure 1b its visibility is 0 because it lies below the red curve.If it had a kernel response of 10 4 its position between red and green curve would be ≈90 %, thus its visibility would be 0.9.Equation 4 describes this Figure 1: Distribution of the kernel response of stars in a good and a cloudy night.Upper limit is defined such that the majority of stars lie above this line in a starry night and in a cloudy night they should be below the lower limit.Stars with mag > 6.7 are always colored red because once the upper limit crosses the lower limit, the visibility is no longer defined.
In the actual analysis, these stars were removed from the catalog and had no influence.
in a mathematical way.The visibility is with the slope of the upper threshold m u and lower one m l and their y-axis intersection b u and b l , respectively, and the star magnitude mag.
Averaging the visibility of the stars i in a certain area, makes the algorithm more robust (cf.eq.5).This area can be chosen arbitrary and may correspond to the telescope's field of view but choosing a wider area that contains more stars yields more robust results.Bright stars are less likely to be affected by a cloud than faint stars therefore a weighting factor 1 of 2.5 mag was used to differentiate between thin and thick clouds as faint stars will be covered completely by both.
The averaged star visibility A is defined as:

High Zenith Angle Correction
At high zenith angles the star's filter response decreases and even drops below the upper threshold.Therefore, it appears that the horizon is always covered by clouds.The cause of this phenomenon is most likely either vignetting of the fish-eye lens [4] or absorption and scattering by the atmosphere's natural aerosols.Although, the real cause was not identified yet, the light loss has to be corrected for the algorithm to work reliably.At high zenith angles the light has to travel a longer way through the atmosphere and thus more photons can be 1 Brightness factor between consecutive magnitudes 5 √ 100 ≈ 2.5.
absorbed or scattered.If the absorption is very small, the transmission can be described by the Beer-Lambert law: with I 0 the initial intensity, I(X) the remaining intensity after distance X and k the absorption coefficient.X is given in units of air mass and by applying the law of cosine on a spherical earth model, the airmass X can be calculated by: with observer altitude y obs = 2.2 km, earth radius R E = 6378 km and atmosphere height y atm = 9.5 km.
To estimate the atmospheric absorption coefficient k, images of three consecutive nights with good weather were used.The images were taken at an interval of 15 min.These measurements of stars at different zenith angles allow to observe the decreasing filter response.A fit of equation 6 to a star's data points yields a k for every single star.The absorption coefficient should be equal for all stars but is spread with k = 0.35 ± 0.27.No correlation to the star's magnitude was observed.The 75 % quantile of k was used to correct the response of all stars.This leads to an over correction of ≈ 3  4 of all stars but a cloud will still block enough light for a notable effect.Under corrected stars, however, induce fake clouds all the time.Therefore, an over correction is the better choice if one corrects all stars with a constant k.The correction of the kernel response was performed before calculating the star visibility.

Cloud Map
The star visibility is only defined at the location of a star.To visualize clouds, it is necessary to develop a parameter that is homogeneous in space.One could achieve this by  calculating the averaged star visibility for every pixel of the image within some radius r but this is computationally very expensive.
A less expensive algorithm that also uses a filter kernel is visualized in Figure 3. Every star was filled into a 2D histogram based on its x and y position in the image.Each entry then was weighted by the star's visibility and the magnitude factor 2.5 mag (green bins).Then, the histogram was treated as a regular image and a Gaussian filter was applied to smear the discrete bin values (green curve).
However, this 'map' is biased to the density of stars.To normalize this curve, another 2D histogram was created with v = 1 for all stars (red bins and curve).After applying the Gaussian kernel, the normalized difference between both curves can be calculated (black curve in lower plot).This is the cloud map used in all further steps.It visualizes the cloud coverage and differentiates between thick and thin clouds but is not biased by changing light conditions like the original camera image.
This implementation takes only seconds to compute a cloud map for the whole sky whereas a LIDAR measurement is in the order of minutes for a spot-measurement.This map can then be made available for the telescope's operating crew once a minute each night.An example is displayed in Figure 2. The cloud map visualizes the structure of clouds with a very good resolution and even very thin clouds that only block a portion of the star light are visible.

Cloud Tracking and Prediction
Identifying single clouds in the cloud map could be done by using a clustering algorithm.Then their direction could be reconstructed by tracking each cloud separately in consecutive images.Clouds, however, can not only change their shape but also split and merge if they move at different velocities.
As a proof of concept, the movement of all clouds was estimated as a unified object instead of identifying and tracking single clouds.The previous and the current was calculated between both maps u and v to rate their similarity.By shifting the previous map in x and y direction, the best translation was found by minimizing the Euclidean distance.This translation vector should correspond to the current wind speed and direction and can be used to predict the next position of the clouds if the wind is constant to some degree.Wind speed and direction can also be obtained from wind profiles [5] and could be used to verify the results.The process is displayed in Figure 4 for a cloud that passes through the image.
The algorithm manages to extract the direction of the cloud in the first three rows (yellow arrow).The reconstruction fails if the cloud enters or leaves the image (row four) because then no valid match can be achieved between both maps.

Conclusion
Cloud tracking and prediction are promising but need some improvements: Changes of the wind direction could be limited to a few degree within a short period of time such that rapid changes like in row four of Fig. 4 do not occur anymore.
Projecting the cloud map onto the model of the spherical atmosphere will increase the accuracy of the prediction because right now, clouds move faster in the image center due to the fish eye lens.
The usage of wind profiles might be sufficient to calculate wind speed and direction with a good precision such that the cloud map will only be used for detecting clouds and determining their shapes.This also increases the forecast time because calculating the speed of a cloud from the cloud map requires some time to accumulate enough images.By that time the cloud already moved further into the image center and by then it might be too late to switch the telescope to a different source.
The cloud maps provide a detailed visualization of the sky's cloud coverage.In combination with a LIDAR for altitude measurements they can be a valuable tool for future IACTs.

Figure 2 :
Figure 2: All sky camera image (left) and the corresponding cloud map (right) that results from the algorithm explained in section 4 and Figure 3. Dark pixels show cloudy regions.The shape of the cloud was reconstructed very detailed and even thin cloud areas can be identified.Clouds close to the horizon were not reconstructed correctly because the number of available stars was too low.In the lower left, the MAGIC telescopes were removed from the image.

Figure 3 :
Figure 3: Explanation of the Cloud map algorithm.The upper plot is a slice of the 2D histogram discussed in the text.The bin size correlates to the star magnitude and the level of the green color displays the fraction of their visibility.The curves are the result of Gaussian smoothing.The normalized difference between both curves is the cloud map displayed in the lower plot.Here, two clouds are visible as peaks with different intensities.

Figure 4 :
Figure 4: Cloud tracking and prediction.The left column shows cloud maps containing a cloud that moved within approx.40 min from the lower right of the image to the upper left.The yellow arrow shows the translation that yielded the smallest Euclidean distance between this map and the previous one (wind speed and direction).The right column shows how the previous cloud map was translated to get the best matching position.The tracking fails in the final row because the cloud left the image and no good match was achieved.