Projecting Cloud Shadows from Clouds
Can we project cloud masks to efficiently detect their shadows on Earth’s surface?
This blog post is the result of activities within the AgriDataValue project, where the aim is to use satellite imagery in the scope of various use cases to understand and model specific phenomena happening on the ground. Removing cloud shadows from the processing steps improves the quality of the input data and thus the downstream models and results.
Clouds have been a point of fascination for many people throughout history, and a point of frustration for many scientists, including those working in the field of Remote Sensing. They pose a significant problem due to obstructing the view of the Earth’s surface from satellites and other sensing platforms. Since such data is often useless, this leads to reduced spatial and temporal resolution in imagery and hinders accurate data collection and interpretation.
To make matters worse, clouds cast large shadows on the surface of the Earth, and these shadows present a unique challenge, as clouds change position and shape with varying cloud cover and sun angles. The resulting shadows distort the appearance of objects and affect the accuracy of data analysis by creating inconsistencies in temporal observations.
So, although the beauty of clouds cannot be denied, sometimes their beauty stands in the way of science, and this is where we draw the line, roll up our sleeves, and start looking for solutions.
Detecting Clouds
Today, cloud detection is mostly a solved problem. Solutions are plentiful, varying in complexity from band math value masking to machine learning and computer vision models. The choice of the right tool depends on the specific problem at hand, taking into account the source and the scale of the data, as well as the required accuracy for your application.
A few years back we designed a pixel-based classifier for detecting clouds on Sentinel-2 imagery, called s2cloudless. We’ve started using it in the majority of our applications, and we’ve even produced it (and continue to do so) for the full Sentinel-2 dataset, so you can use it too! According to the the last Cloud Mask Intercomparison eXercise (CMIX), our classifier was among the top best performing models.
Tackling Cloud Shadow Detection
We asked ourselves if we could solve the problem of detecting cloud shadows with the same approach as mentioned above. After trying it out, we ran into some problems. It turns out that while pixel-based classifiers are effective at distinguishing between different spectral signatures associated with clouds and clear sky, they may struggle to differentiate between shadows cast by clouds and other features on the Earth’s surface.
The main issue is that areas covered by cloud shadows still exhibit spectral characteristics similar to those of non-shadowed areas (bare soil, water bodies, dense vegetation, …). Additionally, the appearance of shadows can vary depending on factors like terrain elevation, surface orientation, and time of day, further complicating their detection using pixel-based classifiers. These problems lead to misclassifications, where shadows are incorrectly detected as non-shadowed areas or vice versa.
To address these challenges, more advanced techniques are often employed, such as taking into account the temporal stability of an area, utilizing dedicated object-detection methods, or incorporating ancillary data sources like terrain models or sun angle information. These approaches allow for analyses beyond simple spectral signatures, improving the accuracy of shadow detection and reducing false positives in remote sensing applications.
Detecting Cloud Shadows via Cloud Projection
We have decided to explore the possibilities of cloud shadow detection via projecting cloud masks. The main inspiration for this blog post was an article with a similar approach to projection called A Newly Developed Algorithm for Cloud Shadow Detection — TIP Method. If you really need a good and reliable cloud shadow detector, we would suggest looking into it, as this blog post is only for demonstration purposes.
The process can be broken down into the following steps:
- Derivation of projection equations
- Generation of cloud masks
- Identification of shadowed areas
- Projection onto the surface of the Earth
- Post-processing
There is also a Jupyter Notebook available here so you can follow along on your own, or if you’re just interested in the technical parts of the tool presented here. Let’s take a deeper look at each of the steps.
Derivation of Projection Equations
This is the point where we start to get a bit more technical. The math we’ll be dealing with is not too hard, it’s high-school-level trigonometry. The problem is that it’s in 3D, which might break our brains, but that’s a good thing!
…Wait, where are you going? Don’t run away!
Did you accidentally scroll below and see the scary equations? It’s OK, come back, don’t be afraid… You got this! We got this!
Baby steps… Let’s start with a simple schematic of how a satellite takes an image of a cloud. Below we see a satellite which is oriented at a specific zenith and azimuth angle. The zenith angle θᵥ defines the angle between the vertical line and the view direction of the satellite (0° means looking down), while the azimuth angle 𝜑ᵥ is the angle between the north vector and the view direction on the horizontal plane (0° means looking north).
Due to the way how satellites see the surface of the Earth, a cloud at the height h above point P would seem to be located at point C in the reconstructed 2D image.
But all we would see in this case is darkness if it weren’t for the Sun. The Sun illuminates the surface and the objects so that we can see them, but the Sun also casts shadows of these objects. Let’s say that the Sun is oriented at some solar zenith (θₛ) and solar azimuth (φₛ) angles which are different from the ones of the satellite.
The cloud at the same height h above point P would cast a shadow down to the Earth’s surface at point S, as shown in the schematic below.
Now we have all the variables necessary to derive the projection equations as a function of the cloud height h. What we actually need to do is move the cloud mask at point C to the cloud shadow at point S, and this is easiest to do through point P, like so:
Clouds cover certain pixels in the image, and since the whole image area is defined with a bounding box in a specific coordinate system, it is easy to extract the global coordinates of a point (x, y).
The hard part now is breaking down the CP and the PS vectors. Looking at the main schematic above and with the help of defining angles α and β, we can define and simplify the remaining pieces of the puzzle:
Finally, we can put it all back together to obtain the formula for the projected point S.
And that’s it! See, it wasn’t that hard. I told you we got this :) Let’s continue with the next section.
Generation of Cloud Mask
There is not much to add here. As mentioned above, cloud masking is more or less a solved problem, and for the purpose of this task, we have generated cloud masks using our s2cloudless model.
Let’s focus on an example area and use it throughout the process. The image below contains clouds and cloud shadows, and we’ve already applied the process of cloud masking.
Identification of Shadowed Areas
We can use the obtained cloud masks to project them onto the Earth’s surface, but we don’t know when to stop projecting, because we don’t know the height of the clouds. One idea for how to know when to stop the projection is to get a mask of cloud shadow candidates and project to the point when the two are best aligned.
However, constructing a mask of cloud shadow candidates is not a simple task, since we stumble upon the already mentioned issue of shadowed areas having spectral similarities with non-shadow areas. No need to be a perfectionist here. Let’s just define a luminosity image and look at the darkest pixels. The luminosity can be defined as the root-mean-square (RMS) of the band values.
Since images can vary in brightness both at different locations and at different times, “dark pixels” is a relative term, so we determine the threshold as the 15% point on the range defined by the cloudless pixels in the image. To remove the outliers, we only consider the values between the 1st and the 99th percentile of the luminosity values of cloudless pixels.
The approach seems to work well and is independent of the type of ground cover, at least as long as the image is large enough that the shadow and non-shadow areas are well-represented. If you look closely, though, you will notice some non-cloud-shadow areas in the dark-pixel mask. The most obvious ones are water bodies and topographic shadows. Ideally, we would want to exclude these regions as much as possible, but it’s fine if the end result is not perfect, as long as the general area of the cloud shadows is represented in the constructed mask.
Detecting water is easy enough — we simply use the Normalized-Difference Water Index (NDWI). Index values greater than ~0.4 usually correspond to water bodies. Topographic shadows — on the other hand — are not so simple to obtain. Assuming you have the Digital Elevation Model (DEM) available at hand, you can write functions for calculating the slope, aspect, and then also hill shadows.
There is a really nice tutorial available here, describing how to obtain the hill shadow mask from DEM values. Here is the extracted code in Python, for your convenience.
Lastly, we can also reduce some noise by removing groups of connected pixels which are too small. Since clouds are large, we also expect the cloud shadows to be large, so we can drop all connected regions which are smaller than, e.g., 50 pixels.
Now we can correct the dark-pixel mask above by excluding small groups of pixels and other areas classified as water bodies and topographic shadows.
Projection onto the Surface of the Earth
Let’s stop for a second, catch our breath, and take stock of what we’ve done so far:
Time to finally run our projection!
Below you can see how the projected cloud masks move across the image. For evaluation, we use the intersection-over-union (IoU) and plot it as a function of cloud height. The optimal cloud height is obtained when the IoU function reaches the maximum value, which corresponds to the point when the two masks are best aligned. By using your favorite optimizer method, you can determine the optimal cloud height for a given image with relative ease.
We’re almost there! The projected mask already looks really good, but we’re still noticing a few issues. Let’s see how to handle them in the next few sections.
Post-processing: Part 1— Buffering the projected mask
As you may have seen above, the cloud masks don’t perfectly match the cloud shadow candidates. Due to the projection and the complicated shape of the 3D clouds, we cannot expect that the top-down satellite view of the shape of the cloud will match the shape of its shadow.
The second reason is that we are using a single, average estimate of the cloud height, but the clouds generally populate a range of heights in a single layer. If, for example, the clouds were present at layers of significantly different cloud heights, our algorithm would fail to take this into account.
To somewhat remedy this, we can apply a mask buffer to improve the match between the two masks, and we can then re-calculate the intersection-over-union. This second part can also be solved with an optimizer in a similar way as above.
Post-processing Part 2— Ghost Clouds
We don’t need to call Ghostbusters for this one, it has more to do with how the clouds cast shadows. For example, it is possible that a cloud is not present in the image, but still casts a visible shadow. Such cloud shadows are not detected with our approach, because we don’t have the cloud mask to project in the first place.
Here is a simple solution — by using the information about the x and y projection offsets, we can define the region where the potential shadows of the hypothetical clouds outside of the image would fall (the green region in the image below). Then we simply take the intersection of this area with the cloud shadow candidates to fill out the problematic edge of our mask.
We’ll stop here with the post-processing and say we’re satisfied with the results. Now it’s time to take a look at the results in other areas and compare them with ground-truth labels.
Validation
For properly validating the method, we need reliable cloud shadow labels. The labels dataset which we decided to use is the CESBIO dataset, developed by Louis Baetens, under the direction of Olivier Hagolle at CESBIO/CNES. The dataset was generated by an active learning method and provides reference cloud mask data for 38 Sentinel-2 scenes. The dataset can be downloaded from here.
We obtained the Sentinel-2 L1C band data for the same areas where the labels are available and ran our method of cloud shadow detection. We also compared the performance of our method with the cloud shadow layer in the Scene Classification map available from the Sentinel-2 L2A processing.
The metrics we have chosen to compare are Matthew’s correlation coefficient (MCC), precision, and recall. MCC is used due to the fact that it takes into account true and false positives and negatives, and is generally regarded as a balanced measure which can be used even if the classes are of very different sizes, which is true in our case. We additionally focused on precision and recall, because precision tells us something about false positives (wrong classifications of shadows) in the sample, while recall contains information on false negatives (not detecting shadows) which weren’t picked up. Below are the results for all of the 38 scenes.
We see that for some of the labeled scenes the L2A data is missing or is not available (data from 2015). In general, it seems that our model prefers fewer false negatives (higher recall) on account of having more false positives (lower precision). We believe this makes sense because false negatives are more costly when you are constructing a mask to filter out data.
The two methods don’t seem to be two equivalent points on the same ROC curve, since, according to the MCC metric, our classification model outperforms the Scene Classification (SC) from Sentinel-2 L2A in 15 out of 29 scenes where the SC data is available. In 6 scenes the SC method seems to be superior, while in the case of the remaining 8 the two methods seem comparable.
Let’s take a closer look at a few examples (excluding the cases where the SC layer is not available).
Example 1 – SC outperforms our method (largest differences)
A few examples where the SC method outperforms our model. The 2nd and 3rd row show a case with clouds at multiple altitudes, so our method fails. The last two examples show a good positional alignment, but our model fails to cover the cloud shadows as well as the SC layer does.
Example 2 — Our method outperforms the SC (best examples)
These are cases where our method outperforms the SC model, sorted by the difference in the MCC metric. You can see how well our detector performs by looking at the visual coverage of the clouds, or by comparing the masks.
Example 3 — Our method outperforms the SC (worst examples)
These are the bottom few cases where our model still outperforms the SC results, however, the differences between the two are getting smaller in terms of the MCC metric.
Caveats of Our Shadow Detector
Some of the caveats became obvious already as we were constructing the detector, so let’s write down the most problematic ones here.
Chronologically speaking, the first problem of our detector is that it depends on the quality of the cloud masks we are using in the projection. Any biases or problems in detecting clouds are propagated to the cloud shadow detector, so the choice of the cloud detector plays an important role.
Secondly, a little less obvious issuse is that in our projection we used spatially averaged values of the view- and Sun angles to make the computation easier and the resulting masks more consistent. While this is completely acceptable for the Sun angles (due to the Sun being far away), doing this for the satellite view angles is a bit more troublesome, since the rate of change across an image can be significant, so such an approach is perhaps more suitable for smaller scenes and not the 100 km × 100 km Sentinel-2 scenes we used here. Therefore we would suggest running the process on smaller patches, with a side length between 5 to 20 km.
The final issue we will mention here (but not the final in general) is the issue with scenes containing clouds at various altitudes. In our approach, we use a single, average cloud height value, but we have seen in the validation section that this assumption is not always correct and alternative approaches might yield better results. One such example is shown below, where a layer of high-altitude clouds is seen on top of a layer of low-altitude clouds, with the additional complexity that high-altitude clouds are cast on top of the low-altitude clouds, confusing the cloud detector as well.
Conclusion
In wrapping up our exploration of cloud shadow detection, we’ve taken quite the journey — from geeky trigonometry to a practical application. However, our goal was not to construct the new best cloud shadow detector, but rather to explore the data, have fun, and see if we learn something new.
What we’ve discovered is both promising and challenging. We’ve come face-to-face with the reality that the nature of cloud shadows is quite frustrating and that to really make a successful classification, one would need a multi-step detector looking at various aspects of tackling the problem at hand, potentially also taking temporal information and spatial context into account.
On the other hand, we have shown that a lot can be done even with a simple and straightforward approach with no involvement of machine learning (if you assume the availability of cloud masks, that is). Such approaches might be quicker to plug into your pipelines, and — as we have seen — have the potential to produce results which are good enough, or often even better than what’s available out there.
Special thanks to Tom Aucler and Urh Primožič for their initial contributions to this work and laying the groundwork for our exploration of cloud shadow detection.
The project has received funding from European Union’s Horizon 2020 “Research and Innovation Programme” under the Grant Agreement №101086461.
