mask_scl_dilation: refactor to create the full mask first, for consistency

[jdries]

mask_scl_dilation often has issues on tile edges, because clouds in neighbouring products are not visible. This was solved by the large overlap in Sentinel-2 products, but returns as an issue when loading with large tile sizes.

I propose to first build the entire mask, then apply it using the regular masking procedure. This has most chances of ensuring consistency.

[JeroenVerstraelen]

@jdries When we construct the entire mask this will also have to be done in tiles. How is this different from the current implementation? It seems like the same issues will arise.

[JeroenVerstraelen]

I did several tests where with a tilesize of 2048 and the new approach solves the issues at the tile borders:

Current New

In the top right of the first image you can see that the mask is clearly cut off by the tile borders. In the second image this border is no longer present and the dilated clouds of the neighbouring tiles have removed much of that mask.

[jdries]

@JeroenVerstraelen what's the setting to enable this? We should also make an assessment of performance impact. Can be done by looking at cpu use, but also spark-UI inspection.

[jdries]

I get an issue on the output edge still, this is my code for the mask:

scl = connection.load_collection(
            S2_collection,
            bands=["SCL"],
            spatial_extent=bbox,
            temporal_extent=[start, end],
            max_cloud_cover=95
        ).resample_spatial(projection=target_crs, resolution=10.0)
        mask = scl.process(
            "to_scl_dilation_mask",
            data=scl,
            kernel1_size=17, kernel2_size=77,
            mask1_values=[2, 4, 5, 6, 7],
            mask2_values=[3, 8, 9, 10, 11],
            erosion_kernel_size=3)
        bands = bands.mask(mask)

[JeroenVerstraelen]

@jdries is it possible to get the full code with bbox and temporal extent? I used what I think is your process graph from the CDSE logs but the result appears to be correct.

This is the code I'm currently using:

bbox = {'west': 3740000.0, 'south': 3020000.0, 'east': 3740150.0, 'north': 3020150.0, 'crs': 'EPSG:3035'}
bands = ["B08", "SCL"]
start = "2021-03-01"
end = "2021-10-31"
S2_collection = "TERRASCOPE_S2_TOC_V2"
connection = openeo.connect("https://openeo-dev.vito.be/openeo/1.0").authenticate_oidc()
target_crs = "EPSG:3035"

s2 = connection.load_collection("TERRASCOPE_S2_TOC_V2", spatial_extent = bbox, temporal_extent = [start, end],
    bands = bands).resample_spatial(projection = target_crs, resolution = 10.0)

scl = connection.load_collection(
            S2_collection,
            bands=["SCL"],
            spatial_extent=bbox,
            temporal_extent=[start, end],
            max_cloud_cover=95
        ).resample_spatial(projection=target_crs, resolution=10.0)
mask = scl.process(
    "to_scl_dilation_mask",
    data=scl,
    kernel1_size=17, kernel2_size=77,
    mask1_values=[2, 4, 5, 6, 7],
    mask2_values=[3, 8, 9, 10, 11],
    erosion_kernel_size=3)
masked = s2.mask(mask)
masked.download("mask_scl_dilation.nc")

[JeroenVerstraelen]

When using a large area I get these results: mask_scl_dilation (current):

to_scl_dilation_mask (new):

Overlayed:

So it looks like:

• The current process is more aggressive than the new process
• Everything that has been masked by the new process was also masked by the current one, they align perfectly.
• The new process masked away some of the south and east edges, while the current process didn't.

[JeroenVerstraelen]

For a far larger extent they appear to be almost identical:

current:

new:

overlay (current in red):

There is only a slight difference at the edge where the current process (red) masked away more than the new one (black):

[jdries]

@JeroenVerstraelen the problem is also visible in your code, but note that QGis does not show these nodata pixels at the edge. When I enlarge the bbox of the mask, the problem goes away. I think we need to do something like this: https://github.com/Open-EO/openeo-python-driver/blob/779b690183932f5f4acb84925167bb6306f8f44f/openeo_driver/dry_run.py#L617 To ensure that load_collection is aware that a pixel buffer is needed.

[jdries]

added a commit to get fix going. The buffer is normally applied here: https://github.com/Open-EO/openeo-geotrellis-extensions/blob/7e8521ce39b6f82f7526cbdce4b49f48842cb410/openeo-geotrellis/src/main/scala/org/openeo/geotrellis/layers/FileLayerProvider.scala#L621

[jdries]

works better now, I only still see a one pixel difference on the mask boundary when working on a very small (15x15 pixel) area.