Started work on s3 pre-processing, some tests

Author: Hneuschmidt

efast/openeo/preprocessing/s3.py

@@ -0,0 +1,88 @@
+from datetime import datetime
+import openeo
+from openeo import processes, UDF
+
+from efast.constants import S3L2SYNClassificationAerosolFlags, S3L2SYNCloudFlags
+from efast.openeo.preprocessing.general import _extract_bit_mask
+
+from . import extract_mask
+
+DATE_FORMAT = "%Y-%m-%d"
+
+
+def extract_clear_land_mask(s3_cube: openeo.DataCube):
+    CLOUD_FLAGS_COMBINED = (
+        S3L2SYNCloudFlags.CLOUD
+        | S3L2SYNCloudFlags.CLOUD_AMBIGUOUS
+        | S3L2SYNCloudFlags.CLOUD_MARGIN
+    )
+    cloud_mask = (s3_cube.band("CLOUD_flags") & CLOUD_FLAGS_COMBINED.value) != 0
+    land_mask = (
+        s3_cube.band("SYN_flags") & S3L2SYNClassificationAerosolFlags.SYN_land.value
+    ) != 0
+    return ((cloud_mask == 0) & land_mask) != 0
+
+
+# TODO add land mask, broken
+# TODO debug as a batch job, to access logs
+def _extract_clear_land_mask(s3_cube: openeo.DataCube) -> openeo.DataCube:
+    mask = _extract_bit_mask(
+        s3_cube,
+        {
+            "CLOUD_flags": (
+                S3L2SYNCloudFlags.CLOUD
+                | S3L2SYNCloudFlags.CLOUD_AMBIGUOUS
+                | S3L2SYNCloudFlags.CLOUD_MARGIN
+            ).value
+        },
+        invert=True,
+    )  # | extract_bit_mask(
+    #    s3_cube,
+    #    {
+    #        "SYN_flags": S3L2SYNClassificationAerosolFlags.SYN_land.value,
+    #    },
+    # )
+    return mask
+
+
+def composite(cube):
+    # What we want: sliding window computation
+    # distance_to_cloud_score = compute_distance_to_cloud_score(distance_to_cloud)  # add normalization
+    # merged = distance_to_cloud_score.merge_cubes(cube)
+    # merged.apply_neighbourhood(process=produce_weighted_composite)
+    pass
+
+
+# TODO rename D to something more useful
+def compute_distance_to_cloud_score(
+    distance_to_cloud: processes.ProcessBuilder, D: float
+) -> processes.ProcessBuilder:
+    score = (distance_to_cloud - 1) / D
+    return processes.clip(score, 0.0, 1.0)
+
+
+def compute_time_weighted_cube(
+    unweighted_cube: processes.ProcessBuilder, target_date: datetime
+) -> processes.ProcessBuilder:
+    udf_compute_temporal_score = UDF.from_file(
+        "efast/openeo/udf/compute_temporal_score",
+        context={"target_date": target_date.strftime(DATE_FORMAT)},
+    )
+    return unweighted_cube.apply_dimension(
+        process=udf_compute_temporal_score, dimension="t"
+    )
+
+
+def produce_weighted_composite(
+    s3_cube: processes.ProcessBuilder,
+    distance_to_cloud_score: processes.ProcessBuilder,
+    target_date: datetime,
+) -> processes.ProcessBuilder:
+    """
+    :param distance_to_cloud_score:
+    """
+    # input1: all relevant Sentinel-3 images to produce a single output composite
+    # input2: distance to cloud scores for all Sentinel-3 observations to avoid
+    # recomputing them for each composite
+
+    pass

tests/test_preprocessing_openeo.py

@@ -18,10 +18,12 @@
 
 from openeo.udf import execute_local_udf
 
+from efast.constants import S3L2SYNClassificationAerosolFlags
 from efast.openeo import preprocessing
 # TODO this should live somewhere else
 from efast.openeo.preprocessing import connect
 from efast import s2_processing, s3_processing
+from efast.openeo.preprocessing.s3 import extract_clear_land_mask
 
 
 TEST_DATA_ROOT = Path(__file__).parent.parent / "test_data"
@@ -40,8 +42,9 @@
 
 SKIP_LOCAL = True
 DOWNLOAD_INTERMEDIATE_RESULTS = True
-SMALL_AREA = True
+SMALL_AREA = False
 
+VISUAL_OUTPUT_PATH = Path(__file__).parent.parent / "visual_test_results"
 
 @contextmanager
 def create_temp_dir_and_copy_files(
@@ -316,3 +319,38 @@ def apply_datacube(cube: XarrayDataCube, context: dict) -> XarrayDataCube:
     #return XarrayDataCube(xr.DataArray(array, dims=["t", "x", "y", "bands"]))
     return XarrayDataCube(xr.DataArray(array, dims=["bands", "x", "y"]))
 """)
+
+def test_extract_clear_land_mask_s3():
+    out_path = VISUAL_OUTPUT_PATH / "extract_clear_land_mask"
+    out_path.mkdir(exist_ok=True, parents=True)
+    conn = connect()
+    conn = conn.authenticate_oidc()
+    bounds = {
+        "west": 399960.0,
+        "south": 1590240.0,
+        "east": 509760.0,
+        "north": 1700040.0,
+        "crs": 32628,
+    }
+    if SMALL_AREA:
+        dist = 3600
+        bounds["east"] = bounds["west"] + dist
+        bounds["north"] = bounds["south"] + dist
+
+    bands = [
+        "Syn_Oa17_reflectance",
+        "CLOUD_flags",
+        "SYN_flags",
+    ]
+    test_area = preprocessing.TestArea(
+        bbox=bounds, s3_bands=bands, temporal_extent=(TEST_DATE_DASH, TEST_DATE_DASH)
+    )
+    cube = test_area.get_s3_cube(conn)
+
+
+    print("Downloading input")
+    #cube.download(out_path / "input.nc")
+    mask = extract_clear_land_mask(cube)
+    print("Downloading mask")
+    mask.download(out_path / "mask.nc")
+    print("Done")