Added udf to compute temporal score and combine with distance score for

weighted composites.

The computation still expects a target date to be provided explicitly,
I am not sure if that is possible with the moving window approach
that I was targeting.

Author: Hneuschmidt

efast/openeo/udf/__init__.py

@@ -0,0 +1,52 @@
+import numpy as np
+import xarray as xr
+from openeo.udf import XarrayDataCube
+from openeo.udf.debug import inspect
+from datetime import datetime
+
+DATE_FORMAT = "%Y-%m-%d"
+CONTEXT = {
+    "target_date": "2018-06-18",
+    "sigma_doy": 10,
+    }
+
+def apply_datacube(cube: XarrayDataCube, context: dict=CONTEXT) -> XarrayDataCube:
+    distance_score = cube.get_array()
+    inspect(distance_score["t"], level="INFO", message="inspecting time dimension")
+    target_date = datetime.strptime(get_from_context(context, "target_date"), DATE_FORMAT)
+
+    # TODO see if this is the correct order
+    time_differences = distance_score.indexes["t"] - target_date
+    sigma_doy = get_from_context(context, "sigma_doy")
+
+    doy_score = np.exp(
+        -1
+        / 2
+        * (np.array(time_differences.days) ** 2 / sigma_doy**2)
+    )
+    # TODO this division is likely redundant
+    #doy_score = doy_score / np.max(doy_score)
+
+    inspect(distance_score.dims, message="distance score dims")
+    inspect(doy_score.shape, message="doy score dims")
+    score = distance_score * doy_score.reshape((-1, 1, 1, 1))
+    score = score / (np.sum(score, axis=0) + 1e-5)
+    inspect(score.shape)
+
+    score = xr.DataArray(
+            score,
+            #dims=["t", "x", "y", "bands"],
+            dims=["t", "bands", "y", "x"],
+            coords={
+                "t": distance_score["t"],
+                "x": distance_score["x"],
+                "y": distance_score["y"],
+                "bands": ["weight"]
+            }
+        )
+
+    return XarrayDataCube(score)
+
+
+def get_from_context(context, key):
+    return context.get(key, CONTEXT[key])

efast/openeo/udf/compute_s3_composite_weights.py

@@ -0,0 +1,70 @@
+from openeo import UDF
+from openeo.udf import XarrayDataCube
+from openeo.udf.run_code import execute_local_udf
+import xarray as xr
+import pandas as pd
+import numpy as np
+
+from efast.openeo.udf.compute_s3_composite_weights import apply_datacube
+
+
+DATE_FORMAT = "%Y-%m-%d"
+
+
+def create_dummy_data():
+    time = pd.date_range("2018-06-01", periods=30, freq="D")
+    x = np.arange(7)
+    y = np.arange(6)
+    bands = ["red", "green", "blue"]
+    data = np.random.rand(len(time), len(x), len(y), len(bands))
+
+    ds = xr.DataArray(
+            data,
+            dims=["t", "x", "y", "bands"],
+            coords={
+                "t": time,
+                "x": x,
+                "y": y,
+                "bands": bands,
+                }
+            )
+    return ds
+
+
+def create_dummy_dtc():
+    time = pd.date_range("2018-06-01", periods=30, freq="D")
+    x = np.arange(7)
+    y = np.arange(6)
+    bands = ["distance_to_cloud"]
+    #data = np.random.rand(len(time), len(x), len(y), len(bands))
+    data = np.random.rand(len(time), len(bands), len(y), len(x))
+
+    da = xr.DataArray(
+            data,
+            #dims=["t", "x", "y", "bands"],
+            dims=["t", "bands", "y", "x"],
+            coords={
+                "t": time,
+                "x": x,
+                "y": y,
+                "bands": bands,
+                }
+            )
+    return da
+
+
+def test_time_difference_score():
+    data = create_dummy_dtc()
+    udf = UDF.from_file("efast/openeo/udf/compute_s3_composite_weights.py")
+    target_date = "2018-06-18"
+    context = {
+        "target_date": target_date,
+        "sigma_doy": 10,
+        }
+    apply_datacube(XarrayDataCube(data), context)
+    res = execute_local_udf(udf, data, fmt="netcdf")
+
+    return True
+
+
+