Detecting cloud shadows

Question

I am trying to detect cloud shadows in Landsat images. I came up with producing good shadows and cloud masks. Now I'd like to implement my cloud mask adding cloud shadows.

Do You have a good idea to do that?

OPTION 1: I thought combining the two masks I have to be a good idea, but I would need an oriented buffer to do that, given that I know the solar azimuth. I found a How to create an oriented buffer using arcpy? with a super good explanation, but I have no idea how to translate the R prototype code provided there into Python (which I use), as I have too few programming skills.

Does anybody know how to deal with this (maybe @whuber himself who provided the code)?

OPTION 2: As I know sun azimuth and angle (from .MTL file of Landsat), I thought to be a good idea to create a fake DEM with the cloud mask inside. Clouds would have their elevations as value in the DEM. Then I would cast the shadows against this "fake DEM" to obtain my cloud shadow mask (or to "play with" it and my shadow mask). Anyway, I don't know a good way to estimate cloud elevations, and, having tried with 5000, 8000, and 10000 meters values with very poor results, I am not totally convinced this is the best method.

Answer 1

As you've probably found, removing clouds from Landsat imagery is not a trivial problem. There is a decent amount of scholarly research about methods for implementing masking and correction for clouds and their shadows. So due to the complexity of this issue, you probably won't find any one size fits all solutions that work perfectly out of the box.

That being said, I've found a paper that may be helpful: Huang, C., Thomas, N., Goward, S. N., Masek, J. G., Zhu, Z., Townshend, J. R., & Vogelmann, J. E. (2010). Automated masking of cloud and cloud shadow for forest change analysis using Landsat images. International Journal of Remote Sensing, 31(20), 5449-5464. A PDF copy is available here. The authors in this paper estimate cloud height based on the thermal infrared data of the clouds, as described in this paragraph:

Once the cloud pixels are identified, their shadows are detected geometrically and spectrally. Specifically, the projected location of each detected cloudy pixel is calculated according to its location, cloud height and solar illumination geometry. Cloud height is estimated by applying a normal lapse rate of 6.4 C km-1 (Christopherson 2002, p. 67, Smithson et al. 2008) to the temperature difference between a cloudy pixel and nearby surface air.

Answer 2

You may check the implementation in GRASS GIS's i.landsat.acca which performs Landsat TM/ETM+ Automatic Cloud Cover Assessment (ACCA). The source code is in C but may give you some ideas.