Add water frequency to A1 module. (#155)

* Added water frequency range to a1 module

* Apply specific data types.

* Change the output type of water seasonality band to int and ensure no data is 255.

Author: tebadi

odc/stats/plugins/lc_veg_class_a1.py

@@ -71,6 +71,7 @@ def __init__(
         saltpan_threshold: Optional[int] = None,
         water_threshold: Optional[float] = None,
         veg_threshold: Optional[int] = None,
+        water_seasonality_threshold: Optional[float] = None,
         **kwargs,
     ):
         super().__init__(**kwargs)
@@ -83,11 +84,14 @@ def __init__(
         )
         self.water_threshold = water_threshold if water_threshold is not None else 0.2
         self.veg_threshold = veg_threshold if veg_threshold is not None else 2
+        self.water_seasonality_threshold = (
+            water_seasonality_threshold if water_seasonality_threshold else 0.25
+        )
         self.output_classes = output_classes
 
     @property
     def measurements(self) -> Tuple[str, ...]:
-        _measurements = ["classes_l3_l4"]
+        _measurements = ["classes_l3_l4", "water_seasonality"]
         return _measurements
 
     def native_transform(self, xx):
@@ -193,6 +197,7 @@ def l3_class(self, xx: xr.Dataset):
         # issues:
         # - nodata information from non-indexed datasets missing
 
+        # Mask nans with NODATA
         l3_mask = expr_eval(
             "where((a!=a)|(b>=nodata), nodata, e)",
             {
@@ -205,17 +210,51 @@ def l3_class(self, xx: xr.Dataset):
             **{"nodata": NODATA},
         )
 
-        return l3_mask
+        # Now add the water frequency
+        # Divide water frequency into following classes:
+        # 0 --> 0
+        # (0,0.25] --> 1
+        # (0.25,1] --> 2
 
-    def reduce(self, xx: xr.Dataset) -> xr.Dataset:
-        l3_mask = self.l3_class(xx)
+        water_seasonality = expr_eval(
+            "where((a > 0) & (a <= wt), 1, a)",
+            {"a": xx["frequency"].data},
+            name="mark_wo_fq",
+            dtype="float32",
+            **{"wt": self.water_seasonality_threshold},
+        )
+
+        water_seasonality = expr_eval(
+            "where((a > wt) & (a <= 1), 2, b)",
+            {"a": xx["frequency"].data, "b": water_seasonality},
+            name="mark_wo_fq",
+            dtype="float32",
+            **{"wt": self.water_seasonality_threshold},
+        )
+
+        water_seasonality = expr_eval(
+            "where((a != a), nodata, a)",
+            {
+                "a": water_seasonality,
+            },
+            name="mark_nodata",
+            dtype="uint8",
+            **{"nodata": NODATA},
+        )
+
+        return l3_mask, water_seasonality
 
+    def reduce(self, xx: xr.Dataset) -> xr.Dataset:
+        l3_mask, water_seasonality = self.l3_class(xx)
         attrs = xx.attrs.copy()
         attrs["nodata"] = int(NODATA)
         data_vars = {
             "classes_l3_l4": xr.DataArray(
                 l3_mask[0], dims=xx["veg_frequency"].dims[1:], attrs=attrs
-            )
+            ),
+            "water_seasonality": xr.DataArray(
+                water_seasonality[0], dims=xx["veg_frequency"].dims[1:], attrs=attrs
+            ),
         }
         coords = dict((dim, xx.coords[dim]) for dim in xx["veg_frequency"].dims[1:])
         return xr.Dataset(data_vars=data_vars, coords=coords, attrs=xx.attrs)

tests/test_landcover_plugin_a1.py

@@ -149,10 +149,55 @@ def test_l3_classes(dataset):
         dtype="uint8",
     )
 
-    res = stats_l3.l3_class(dataset)
+    res, water_seasonality = stats_l3.l3_class(dataset)
     assert (res == expected_res).all()
 
 
+def test_l4_water_seasonality(dataset):
+    stats_l3 = StatsVegClassL1(
+        output_classes={
+            "aquatic_veg": 124,
+            "terrestrial_veg": 110,
+            "water": 221,
+            "intertidal": 223,
+            "surface": 210,
+        },
+        optional_bands=["canopy_cover_class", "elevation"],
+    )
+
+    wo_fq = np.array(
+        [
+            [
+                [0.0, 0.021, 0.152, 255],
+                [0.249, 0.273, 0.252, 0.0375],
+                [0.302, 0, 0.789, 0.078],
+                [0.021, 0.243, 255, 0.255],
+            ]
+        ],
+        dtype="float32",
+    )
+    wo_fq = da.from_array(wo_fq, chunks=(1, -1, -1))
+
+    dataset["frequency"] = xr.DataArray(
+        wo_fq, dims=("spec", "y", "x"), attrs={"nodata": np.nan}
+    )
+
+    expected_water_seasonality = np.array(
+        [
+            [
+                [0, 1, 1, 255],
+                [1, 2, 2, 1],
+                [2, 0, 2, 1],
+                [1, 1, 255, 2],
+            ]
+        ],
+        dtype="float32",
+    )
+
+    res, water_seasonality = stats_l3.l3_class(dataset)
+    assert np.allclose(water_seasonality, expected_water_seasonality)
+
+
 def test_reduce(dataset):
     stats_l3 = StatsVegClassL1(
         output_classes={