Terraced artefacts in Fmask cloud shadow data

[robbibt]

Hi all, we've noticed some strange "terraced" artefacts in the cloud shadow data produced by Fmask in Digital Earth Australia's Landsat Collection 3 datasets:

The example above is from 2020-02-26 offshore of Kurnell, Sydney, however the artefacts are common and affect most timesteps for that location and others (I haven't exhaustively tested this, but I've seen it come up regularly in other locations too). I can provide a reproducible DEA datacube query if it would help track down the issue.

The Fmask data was produced using v0.5.4 of python-fmask:

Does anyone know what is causing these artefacts?

[neilflood]

I have never seen anything like this before.

I just checked the same date and location on the Qld RSC system (path/row 089/084, date 2020-02-26). It was also run using python-fmask 0.5.4, and does not show any artefacts of this type.

I suggest that you take this up with the people in GA. At this point, I have no reason to believe that it is inherent in python-fmask. I have never looked into what GA have done with the code to run for DEA.

[robbibt]

Hey @neilflood, thanks very much for the quick reply! Any chance you could post a screenshot of the equivalent Qld result for this area, or point me to where I could access this? It would help me demonstrate that this might be a DEA implementation issue.

[neilflood]

This is the cloud (pink) and shadow (yellow)

And here is the underlying imagery, just to check I am in the right location

I have never before looked at the DEA cloud masks, but from your snapshot, it does appear that they might have switched off all the buffering. This is quite optimistic, so I am not sure what other parameters they might have changed.

If you wanted to test in isolation from everything else, you could do worse than to download the original image file from USGS and run the original python-fmask code on it directly, as per instructions at http://www.pythonfmask.org/en/latest/#usgs-landsat. This will use the default behaviours, and I would expect it to be roughly the same as the snapshot I show above.

[robbibt]

That's really helpful - yes, DEA's Fmask layers do deliberately have buffering turned off so that users can make their own decisions on buffers for their own applications (we had issues for coastal applications where false positive clouds caused entire coastal strips to be removed from the data after buffering). It could be that those artefacts are just a standard output of Fmask, they're just usually hidden by the default buffering selection.

[neilflood]

I guess it might be, but I suggest you look deeper. I don't recall ever seeing anything like it during development.

[robbibt]

Thanks for your help @neilflood. I've been able to install python-fmask using the instructions in INSTALL.txt - no obvious errors have come up, but when it comes to actually running it on a USGS C1 L1 scene I get the following Cannot find a *_MTL.txt file in specified dir error.

I'm working off the DE Africa Sandbox which is a Python 3.6.9 Linux-based environment - can send through any other environmental details if it would help.

jovyan@jupyter-robbibt:~$ unzip -q python-fmask-pythonfmask-0.5.4.zip 

jovyan@jupyter-robbibt:~$ cd python-fmask-pythonfmask-0.5.4

jovyan@jupyter-robbibt:~/python-fmask-pythonfmask-0.5.4$ python setup.py build
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "_fillminima" sources
building extension "_valueindexes" sources
building data_files sources
build_src: building npy-pkg config files
running build_py
creating build
creating build/lib.linux-x86_64-3.6
creating build/lib.linux-x86_64-3.6/fmask
copying fmask/valueindexes.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/sen2meta.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/fmask.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/fillminima.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/landsatangles.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/landsatTOA.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/__init__.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/zerocheck.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/fmaskerrors.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/saturationcheck.py -> build/lib.linux-x86_64-3.6/fmask
copying fmask/config.py -> build/lib.linux-x86_64-3.6/fmask
creating build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/usgsLandsatStacked.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/sentinel2Stacked.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/usgsLandsatMakeAnglesImage.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/sentinel2makeAnglesImage.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/usgsLandsatSaturationMask.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/usgsLandsatTOA.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
copying fmask/cmdline/__init__.py -> build/lib.linux-x86_64-3.6/fmask/cmdline
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
building '_fillminima' extension
compiling C sources
C compiler: x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC

creating build/temp.linux-x86_64-3.6/src
compile options: '-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -I/env/lib/python3.6/site-packages/numpy/core/include -I/env/include -I/usr/include/python3.6m -c'
x86_64-linux-gnu-gcc: src/fillminima.c
x86_64-linux-gnu-gcc -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.6/src/fillminima.o -o build/lib.linux-x86_64-3.6/fmask/_fillminima.cpython-36m-x86_64-linux-gnu.so
building '_valueindexes' extension
compiling C sources
C compiler: x86_64-linux-gnu-gcc -pthread -DNDEBUG -g -fwrapv -O2 -Wall -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 -fPIC

compile options: '-DNPY_NO_DEPRECATED_API=NPY_1_7_API_VERSION -I/env/lib/python3.6/site-packages/numpy/core/include -I/env/include -I/usr/include/python3.6m -c'
x86_64-linux-gnu-gcc: src/valueindexes.c
x86_64-linux-gnu-gcc -pthread -shared -Wl,-O1 -Wl,-Bsymbolic-functions -Wl,-Bsymbolic-functions -Wl,-z,relro -Wl,-Bsymbolic-functions -Wl,-z,relro -g -fstack-protector-strong -Wformat -Werror=format-security -Wdate-time -D_FORTIFY_SOURCE=2 build/temp.linux-x86_64-3.6/src/valueindexes.o -o build/lib.linux-x86_64-3.6/fmask/_valueindexes.cpython-36m-x86_64-linux-gnu.so
running build_scripts
creating build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_sentinel2Stacked.py -> build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_usgsLandsatStacked.py -> build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_sentinel2makeAnglesImage.py -> build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_usgsLandsatSaturationMask.py -> build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_usgsLandsatTOA.py -> build/scripts.linux-x86_64-3.6
copying and adjusting bin/fmask_usgsLandsatMakeAnglesImage.py -> build/scripts.linux-x86_64-3.6
changing mode of build/scripts.linux-x86_64-3.6/fmask_sentinel2Stacked.py from 644 to 755
changing mode of build/scripts.linux-x86_64-3.6/fmask_usgsLandsatStacked.py from 644 to 755
changing mode of build/scripts.linux-x86_64-3.6/fmask_sentinel2makeAnglesImage.py from 644 to 755
changing mode of build/scripts.linux-x86_64-3.6/fmask_usgsLandsatSaturationMask.py from 644 to 755
changing mode of build/scripts.linux-x86_64-3.6/fmask_usgsLandsatTOA.py from 644 to 755
changing mode of build/scripts.linux-x86_64-3.6/fmask_usgsLandsatMakeAnglesImage.py from 644 to 755

jovyan@jupyter-robbibt:~/python-fmask-pythonfmask-0.5.4$ python setup.py install
running install
running build
running config_cc
unifing config_cc, config, build_clib, build_ext, build commands --compiler options
running config_fc
unifing config_fc, config, build_clib, build_ext, build commands --fcompiler options
running build_src
build_src
building extension "_fillminima" sources
building extension "_valueindexes" sources
building data_files sources
build_src: building npy-pkg config files
running build_py
running build_ext
customize UnixCCompiler
customize UnixCCompiler using build_ext
running build_scripts
running install_lib
creating /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/valueindexes.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/sen2meta.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/fmask.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/fillminima.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/landsatangles.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/_fillminima.cpython-36m-x86_64-linux-gnu.so -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/landsatTOA.py -> /env/lib/python3.6/site-packages/fmask
creating /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/usgsLandsatStacked.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/sentinel2Stacked.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/usgsLandsatMakeAnglesImage.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/sentinel2makeAnglesImage.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/usgsLandsatSaturationMask.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/usgsLandsatTOA.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/cmdline/__init__.py -> /env/lib/python3.6/site-packages/fmask/cmdline
copying build/lib.linux-x86_64-3.6/fmask/_valueindexes.cpython-36m-x86_64-linux-gnu.so -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/__init__.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/zerocheck.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/fmaskerrors.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/saturationcheck.py -> /env/lib/python3.6/site-packages/fmask
copying build/lib.linux-x86_64-3.6/fmask/config.py -> /env/lib/python3.6/site-packages/fmask
byte-compiling /env/lib/python3.6/site-packages/fmask/valueindexes.py to valueindexes.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/sen2meta.py to sen2meta.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/fmask.py to fmask.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/fillminima.py to fillminima.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/landsatangles.py to landsatangles.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/landsatTOA.py to landsatTOA.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatStacked.py to usgsLandsatStacked.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/sentinel2Stacked.py to sentinel2Stacked.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatMakeAnglesImage.py to usgsLandsatMakeAnglesImage.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/sentinel2makeAnglesImage.py to sentinel2makeAnglesImage.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatSaturationMask.py to usgsLandsatSaturationMask.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatTOA.py to usgsLandsatTOA.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/cmdline/__init__.py to __init__.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/__init__.py to __init__.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/zerocheck.py to zerocheck.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/fmaskerrors.py to fmaskerrors.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/saturationcheck.py to saturationcheck.cpython-36.pyc
byte-compiling /env/lib/python3.6/site-packages/fmask/config.py to config.cpython-36.pyc
running install_scripts
copying build/scripts.linux-x86_64-3.6/fmask_sentinel2Stacked.py -> /env/bin
copying build/scripts.linux-x86_64-3.6/fmask_usgsLandsatStacked.py -> /env/bin
copying build/scripts.linux-x86_64-3.6/fmask_sentinel2makeAnglesImage.py -> /env/bin
copying build/scripts.linux-x86_64-3.6/fmask_usgsLandsatSaturationMask.py -> /env/bin
copying build/scripts.linux-x86_64-3.6/fmask_usgsLandsatTOA.py -> /env/bin
copying build/scripts.linux-x86_64-3.6/fmask_usgsLandsatMakeAnglesImage.py -> /env/bin
changing mode of /env/bin/fmask_sentinel2Stacked.py to 755
changing mode of /env/bin/fmask_usgsLandsatStacked.py to 755
changing mode of /env/bin/fmask_sentinel2makeAnglesImage.py to 755
changing mode of /env/bin/fmask_usgsLandsatSaturationMask.py to 755
changing mode of /env/bin/fmask_usgsLandsatTOA.py to 755
changing mode of /env/bin/fmask_usgsLandsatMakeAnglesImage.py to 755
running install_data
copying LICENSE.txt -> /env/lib/python3.6/site-packages/
running install_egg_info
Writing /env/lib/python3.6/site-packages/python_fmask-0.5.4.egg-info
running install_clib
customize UnixCCompiler

jovyan@jupyter-robbibt:~/python-fmask-pythonfmask-0.5.4$ fmask_usgsLandsatStacked.py -o cloud.img --scenedir LC08_L1TP_089084_20200226_20200313_01_T1

Traceback (most recent call last):
  File "/env/bin/fmask_usgsLandsatStacked.py", line 22, in <module>
    usgsLandsatStacked.mainRoutine()
  File "/env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatStacked.py", line 203, in mainRoutine
    makeStacksAndAngles(cmdargs)
  File "/env/lib/python3.6/site-packages/fmask/cmdline/usgsLandsatStacked.py", line 112, in makeStacksAndAngles
    raise fmaskerrors.FmaskFileError("Cannot find a *_MTL.txt file in specified dir")
fmask.fmaskerrors.FmaskFileError: Cannot find a *_MTL.txt file in specified dir
jovyan@jupyter-robbibt:~/python-fmask-pythonfmask-0.5.4$

[neilflood]

Well, the most obvious question is: Is there a *_MTL.txt file in that directory? If not, then you have a problem.

[robbibt]

My apologies - I misread the docs and thought fmask_usgsLandsatStacked.py accepted a tarred (not untarred) scene - I have it running now I've untarred it. Here's the same scene with no cloud shadow buffering (cloud_0.zip): looks like the artefact is produced here too

fmask_usgsLandsatStacked.py -o cloud_0.img --scenedir LC08_L1TP_089084_20200226_20200313_01_T1 --shadowbufferdistance 0

[neilflood]

No worries.

Hmm... that is very interesting, thank you. To answer your original question, no, I don't have an immediate explanation for it. I have a few hypotheses. Given that you have seen this often before, does it generally happen over ocean? I confess I have not paid much attention to the cloud shadow over water, mainly because shadows are very difficult to delineate over water (because everything is dark).

Not sure.

I should ask at this point, is this causing serious problems, or is it just a curiosity? I don't have an immediate solution for it, except to say that buffering of shadows is there for a reason (separate to this), so you might like to turn it on again. Similarly with buffering of the clouds themselves.

[robbibt]

I originally thought this was mainly an ocean issue, but looking into the data I can see that it appears to occur just as frequently inland too. You can browse results for Landsat 5/7/8 here - seems to affect all sensors equally:

https://terria-cube.terria.io/#share=s-lhEIUo45kVkYMdUInhojlVvFVu6

I'm currently trying to evaluate how much of an issue this is for us. For a lot of applications (e.g. inland) buffering is a great solution, but for coastal it's much more problematic as you get problems like this (unbuffered on left, buffered on right):

I've previously used morphological opening operations to first erode then dilate cloud masks to remove long and narrow false positives like bright sandy beaches, then applied the dilation after that. However, if the cloud shadow data is affected by artefacts like this, it's likely that I am inadvertently wiping out the narrow terraced shadow pixels too.

[neilflood]

Thanks @robbibt

If this is also happening over land, then that seems serious. I will do some digging and see if I can understand what is happening. Right now, I don't understand it at all, so this may take a little while. I will get back to you when I know more.

Aside from that, I agree about the problems with bright areas like beaches. This has always been widespread, and I have no perfect solution.

[neilflood]

Hi @robbibt

I have understood what is happening. The problem arises in trying to deduce the 3D structure of the cloud shapes. If you are interested, there are some comments about this problem in the code, from when I first implemented this. You can see my original attempt at a rigorous approach, which turned out to be infeasible and so I commented out, and then the simpler replacement which works reasonably well. https://github.com/ubarsc/python-fmask/blob/063f24cd109883f2b168091f82b0317c8518bc33/fmask/fmask.py#L913-L940

The idea is to use the thermal band to construct the 3D shape, and then calculate the shape of the shadow (to be used later in matching against possible shadows on the ground). If the upper surface of the cloud includes steep gradients on the sun-facing side, and if the sun angle is low enough, then the modelled sun is effectively shining through gaps in the modelled upper surface of the cloud. The terracing you see is effectively the contours of the upper surface of the cloud being sliced into layers by the steep drop interacting with the low sun angle.

It is further complicated by how close to nadir the cloud is, as the sensor only sees the shape of the cloud projected along line of sight, not projected vertically, so cloud shapes are more difficult to infer with off-nadir views.

This all explains why it only occurs in some situations. Mainly, it would be clouds with large variation in height of their upper surface, when the sun elevation is quite low (i.e. worse in winter, and/or further from the equator).

A true fix would be something like the sort of thing which I abandoned in my original code, years ago, or some other thing like ray-tracing of opaque wire-frame shapes. All of which would be expensive in either memory or compute time.

I should add that this problem should not occur at all when used with Sentinel-2, because there is no thermal band, and so we do not even attempt to construct the 3D shape. We just assume the cloud is a flat surface (which is wrong, but the best we can do).

My strongest suggestion is that you turn on the buffering for the shadows, as it blurs these gaps nicely (although this was just fortuitous). If you really don't want to do that, then I would suggest that you use a binary closing operation on the shadows (i.e. dilation followed by erosion), which will have much the same effect within each shadow, without greatly expanding the shadow edges. The original paper also used some buffering of shadows, too (see the 2nd last paragraph of section 3.2 of Zhu et al, 2012), and the algorithm is really not intended to be used without it.