Added masking to binning implementation

Author: Hneuschmidt

tests/test_binning.py

@@ -5,7 +5,8 @@
 import time
 from efast import eoprofiling
 
-from efast.binning import BBox, SynProduct, bin_to_grid, create_geogrid, get_reflectance_filename, bin_to_grid_numpy
+from efast.binning import FILL_VALUE, BBox, SynProduct, bin_to_grid, create_geogrid, get_reflectance_filename, bin_to_grid_numpy
+from efast.constants import S3L2SYNClassificationAerosolFlags, S3L2SYNCloudFlags
 
 SYN_PRODUCT_PATH = (
     Path(__file__).parent.parent
@@ -31,22 +32,37 @@ def test_create_geogrid():
     lat, lon = create_geogrid(np.array([], dtype=np.float32))
     return True
 
+CLOUD_FLAGS_COMBINED = (
+    S3L2SYNCloudFlags.CLOUD
+    | S3L2SYNCloudFlags.CLOUD_AMBIGUOUS
+    | S3L2SYNCloudFlags.CLOUD_MARGIN
+)
 
 def test_bin_to_grid_dev():
     reflectance_bands = [get_reflectance_filename(i) for i in (4, 6, 8, 17)]
-    band_names = [*reflectance_bands, "CLOUD_flags"]
+    band_names = [*reflectance_bands, "CLOUD_flags", "SYN_flags"]
     prod = SynProduct(SYN_PRODUCT_PATH, band_names=band_names)
     ds = prod.read_bands()
     bbox = BBox.from_wkt(REGION)
     grid = create_geogrid(bbox, num_rows=66792)
 
-    #binned = bin_to_grid_numpy(ds, reflectance_bands, grid, super_sampling=2)
+    cloud_no_land_mask = ((ds["CLOUD_flags"] & CLOUD_FLAGS_COMBINED.value) != 0) & (ds["SYN_flags"] & S3L2SYNClassificationAerosolFlags.SYN_land.value > 0)
+
+    for reflectance_band in reflectance_bands:
+         band_values = ds[reflectance_band].data
+         band_values[cloud_no_land_mask] = -10000
+         ds[reflectance_band] = (["lat", "lon"], band_values)
+
+
     binned = bin_to_grid_numpy(ds, reflectance_bands, grid, super_sampling=2)
     pixel_size = grid.lat[1] - grid.lat[0]
-    # TODO why last lat?
+
+    # TODO
+    binned[binned < 0] = np.nan
+
     transform = Affine.translation(grid.lon[0], grid.lat[-1]) * Affine.scale(pixel_size, pixel_size)
 
-    with rasterio.open(TEST_DATA_OUTPUT_PATH / "binned.tif", "w", driver="GTiff", height=binned.shape[1], width=binned.shape[2], count=len(reflectance_bands), dtype=binned.dtype, crs="+proj=latlong", transform=transform) as ds:
+    with rasterio.open(TEST_DATA_OUTPUT_PATH / "binned.tif", "w", driver="GTiff", height=binned.shape[1], width=binned.shape[2], count=len(reflectance_bands), dtype=binned.dtype, crs="+proj=latlong", transform=transform, nodata=np.nan) as ds:
         ds.write(binned[:, ::-1, :])
 
     return True