How to increase mask buffer and then offset it using r.mapcalc in GRASS?

Question

I'm in the process of learning on how to remove clouds and cloud-shadows from Landsat images using GRASS. So far I was able to remove a lot of clouds from those images using i.landsat.acca, but it seems like they still require a lot of refinements:

I followed the GRASS manual for i.landsat.acca to obtain the image above: http://grass.osgeo.org/manuals/html70_user/i.landsat.acca.html.

Basically, the process is to run different Landsat band images for the same scene with i.landsat.toar, and then processed with i.landsat.acca to get the cloud areas. Then use the cloud area as a mask using r.mapcalc:

r.mapcalc "MASK = if(isnull(clouds))"

Problem#1: the cloud edges is still visible. I need to increase the cloud-mask buffer to catch those renegade cloud edges.

Problem#2: the cloud-shadows is still visible. If I include the -s flag to include a category for cloud-shadow when running i.landsat.acca, a large amounts of shadow from the terrain will consequently be removed as well... So to catch the cloud-shadows, I believe I need to make a new mask for them - by offsetting the current cloud mask until they cover the cloud-shadow areas.

To achieve both objectives, I believe I need to modify and play with the r.mapcalc command, but I can't figure out how to write it (due to my bad Math and my-even-worse programming skills)..

So, my ultimate question:

Question#1: How do I edit r.map.calc to increase the mask area (let's say by a buffer of 10 units)? Question#2: How do we offset a mask in r.map.calc (in this case a little to the south-west direction)?

If I know the answer to both of these: I can combine both cloud-mask and cloud-shadow-mask, to get a relatively cloud-free image. Then it's a matter of overlaying it to images of the same scene of different dates.

I appreciate any help that I can get. Many thanks!

Answer 1

Thanks to @markusN and @whuber for pointing out some of my glaring oversight on my approach of dealing with this problem.

At first I thought what I needed was a mathematical r.mapcalc formula that could magically offset a mask around - which is the wrong thought-process. What I really need is to first apply specific raster commands to the cloud raster, and then include the processed raster in the r.mapcalc as a mask.

So, the answer to my problems:

Problem#1: the cloud edges is still visible. I need to increase the cloud-mask buffer to catch those renegade cloud edges.

Solution: Use the cloud-cover generated from i.landsat.acca and apply r.buffer or r.grow to it. For example:

r.buffer input=clouds output=clouds_buffered distances=30 unit=meters

Problem#2: the cloud-shadows is still visible. If I include the -s flag to include a category for cloud-shadow when running i.landsat.acca, a large amounts of shadow from the terrain will consequently be removed as well... So to catch the cloud-shadows, I believe I need to make a new mask for them - by offsetting the current cloud mask until they cover the cloud-shadow areas.

Solution: Duplicate the cloud-cover and offset it to its shadow using r.region:

r.region map=shadow n=n+100 e=e+100 w=w+100 s=s+100

or just by using 2-directions with + and - sign:

r.region map=shadow n=n-100 e=e+100

The final mask:

We can finally combine both cloud and cloud-shadow cover using r.patch:

r.patch input=clouds,shadow output=totalcloudcover

and then create a mask from that with r.mapcalc:

r.mapcalc "MASK = if(isnull(totalcloudcover))"

Final thoughts:

As @whuber pointed out, it's not as simple as reprojecting clouds around to mask its own shadow, because of the range of terrain, cloud-elevation and illumination angle effects.. The mask resulting from lateral movements overestimates all the areas with potential cloud shadows, and could include some underlying usable data - and should be further differentiated. At the very least, this technique is a start.