Corrected the implementation of water seasonality and water persistence. Also made a range of other mappings at Level 4 based on the latest understanding of Collection 2

Author: tebadi

odc/stats/plugins/l34_utils/l4_natural_aquatic.py

@@ -4,12 +4,24 @@
 
 from odc.stats._algebra import expr_eval
 
+NODATA = 255
 
-def natural_auquatic_veg(l4, veg_cover, water_seasonality):
+def natural_auquatic_veg(l4, veg_cover, water_season):
 
     # mark woody/herbaceous
     # mangroves -> woody
     # everything else -> herbaceous
+
+    water_seasonality = expr_eval(
+        "where((a==a), a, nodata)",
+        {
+            "a": water_season,
+        },
+        name="mark_water_season",
+        dtype="float32",
+         **{"nodata": NODATA},
+    )
+  
     res = expr_eval(
         "where((a==124), 56, a)",
         {
@@ -102,7 +114,7 @@ def natural_auquatic_veg(l4, veg_cover, water_seasonality):
     )
 
     res = expr_eval(
-        "where((a==252)&(b==10), 80, a)",
+        "where((a==252)&(b==10), 79, a)",
         {
             "a": res,
             "b": veg_cover,
@@ -112,7 +124,7 @@ def natural_auquatic_veg(l4, veg_cover, water_seasonality):
     )
 
     res = expr_eval(
-        "where((a==251)&(b==10), 81, a)",
+        "where((a==251)&(b==10), 80, a)",
         {
             "a": res,
             "b": veg_cover,
@@ -271,7 +283,7 @@ def natural_auquatic_veg(l4, veg_cover, water_seasonality):
         name="mark_final",
         dtype="uint8",
     )
-
+    
     res = expr_eval(
         "where((a==251)&(b==16), 92, a)",
         {
@@ -282,4 +294,24 @@ def natural_auquatic_veg(l4, veg_cover, water_seasonality):
         dtype="uint8",
     )
 
+    # There are cases where a tile falls over water. 
+    # In these cases, the PC will have no data so we map back 251-254 to their corresponding classes
+    res = expr_eval(
+        "where((a>=251)&(a<=252), 57, a)",
+        {
+            "a": res,
+        },
+        name="mark_final",
+        dtype="uint8",
+    )
+    res = expr_eval(
+        "where((a>=253)&(a<=254), 58, a)",
+        {
+            "a": res,
+        },
+        name="mark_final",
+        dtype="uint8",
+    )
+    
+
     return res

odc/stats/plugins/l34_utils/l4_veg_cover.py

@@ -7,6 +7,7 @@
 
 def canopyco_veg_con(xx: xr.Dataset, veg_threshold):
 
+    
     # Mask NODATA
     pv_pc_50 = expr_eval(
         "where(a==a, a, nodata)",
@@ -88,7 +89,7 @@ def canopyco_veg_con(xx: xr.Dataset, veg_threshold):
 
     # 65-100 --> 10
     veg_mask = expr_eval(
-        "where((a>=m)&(a<n), 10, b)",
+        "where((a>=m)&(a<=n), 10, b)",
         {
             "a": pv_pc_50,
             "b": veg_mask,

odc/stats/plugins/l34_utils/l4_water.py

@@ -5,43 +5,62 @@
 
 def water_classification(xx, intertidal_mask, water_persistence):
 
+    # Replace nan with nodata
     l4 = expr_eval(
-        "where(((a==223)|(a==221))&(b==1), 101, a)",
-        {"a": xx.level_3_4.data, "b": water_persistence},
+        "where((a==a), a, nodata)", 
+         {"a": xx.level_3_4.data}, 
+         name="mark_water", 
+         dtype="uint8",
+         **{"nodata": NODATA},
+    )
+
+    l4 = expr_eval(
+        "where((a==223)|(a==221), 98, a)", 
+         {"a": l4}, 
+         name="mark_water", 
+         dtype="uint8"
+    )
+    
+    l4 = expr_eval(
+        "where((a==98)&(b!=1), 99, a)",
+        {"a": l4, "b": intertidal_mask},
         name="mark_water",
         dtype="uint8",
     )
-
+    
     l4 = expr_eval(
-        "where(((a==223)|(a==221))&(b==7), 102, a)",
-        {"a": l4, "b": water_persistence},
+        "where((a==98)&(b==1), 100, a)",
+        {"a": l4, "b": intertidal_mask},
         name="mark_water",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where(((a==223)|(a==221))&(b==8), 103, a)",
+        "where((a==99)&(b==1), 101, a)",
         {"a": l4, "b": water_persistence},
         name="mark_water",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where(((a==223)|(a==221))&(b==9), 104, a)",
+        "where((a==99)&(b==7), 102, a)",
         {"a": l4, "b": water_persistence},
         name="mark_water",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where(((a==223)|(a==221))&(b==1), 100, a)",
-        {"a": l4, "b": intertidal_mask},
+        "where((a==99)&(b==8), 103, a)",
+        {"a": l4, "b": water_persistence},
         name="mark_water",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where((a==223)|(a==221), 98, a)", {"a": l4}, name="mark_water", dtype="uint8"
+        "where((a==99)&(b==9), 104, a)",
+        {"a": l4, "b": water_persistence},
+        name="mark_water",
+        dtype="uint8",
     )
 
     return l4

odc/stats/plugins/l34_utils/lc_water_seasonality.py

@@ -11,7 +11,6 @@
 
 from .l34_utils import (
     l4_water_persistence,
-    lc_water_seasonality,
     l4_veg_cover,
     lc_level3,
     l4_cultivated,
@@ -27,7 +26,6 @@
 
 NODATA = 255
 
-
 class StatsLccsLevel4(StatsPluginInterface):
     NAME = "ga_ls_lccs_Level34"
     SHORT_NAME = NAME
@@ -39,7 +37,6 @@ def __init__(
         veg_threshold: Optional[List] = None,
         bare_threshold: Optional[List] = None,
         watper_threshold: Optional[List] = None,
-        water_seasonality_threshold: int = None,
         **kwargs,
     ):
         super().__init__(**kwargs)
@@ -51,9 +48,7 @@ def __init__(
         self.watper_threshold = (
             watper_threshold if watper_threshold is not None else [1, 4, 7, 10]
         )
-        self.water_seasonality_threshold = (
-            water_seasonality_threshold if water_seasonality_threshold else 3
-        )
+        
 
     def fuser(self, xx):
         return xx
@@ -64,19 +59,15 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
         water_persistence = l4_water_persistence.water_persistence(
             xx, self.watper_threshold
         )
-
-        water_seasonality = lc_water_seasonality.water_seasonality(
-            xx, self.water_seasonality_threshold
-        )
-
+   
         intertidal_mask = lc_intertidal_mask.intertidal_mask(xx)
 
         # #TODO WATER (99-104)
         l4 = l4_water.water_classification(xx, intertidal_mask, water_persistence)
 
         # Generate Level3 classes
         level3 = lc_level3.lc_level3(xx)
-
+        
         # Vegetation cover
         veg_cover = l4_veg_cover.canopyco_veg_con(xx, self.veg_threshold)
 
@@ -88,18 +79,19 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
 
         # Apply terrestrial vegetation classes [19-36]
         l4 = l4_natural_veg.lc_l4_natural_veg(l4, level3, lifeform, veg_cover)
-
+ 
         # Bare gradation
         bare_gradation = l4_bare_gradation.bare_gradation(
             xx, self.bare_threshold, veg_cover
         )
-        l4 = l4_natural_aquatic.natural_auquatic_veg(l4, veg_cover, water_seasonality)
+
+        l4 = l4_natural_aquatic.natural_auquatic_veg(l4, veg_cover, xx.water_season)
 
         level4 = l4_surface.lc_l4_surface(l4, level3, bare_gradation)
 
         level3 = level3.astype(np.uint8)
         level4 = level4.astype(np.uint8)
-
+   
         attrs = xx.attrs.copy()
         attrs["nodata"] = NODATA
         dims = xx.level_3_4.dims[1:]
@@ -113,4 +105,4 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
         return leve34
 
 
-register("lc_l3_l4", StatsLccsLevel4)
+register("lc_l3_l4", StatsLccsLevel4)

odc/stats/plugins/l34_utils/utils.py

@@ -17,8 +17,8 @@ def image_groups():
             [
                 [210, 210, 210],
                 [210, 210, 210],
-                [210, 210, 210],
-                [210, 210, 210],
+                [223, 210, 210],
+                [221, 221, 221],
             ]
         ],
         dtype="uint8",
@@ -89,13 +89,25 @@ def image_groups():
             [
                 [1, 3, 2],
                 [4, 5, 6],
-                [9, 2, 11],
+                [2, 2, 11],
                 [10, 11, 12],
             ]
         ],
         dtype="uint8",
     )
 
+    water_season = np.array(
+        [
+            [
+                [1, 2, 1],
+                [2, 1, 2],
+                [1, 1, 2],
+                [2, 2, 1],
+            ]
+        ],
+        dtype="uint8",
+    )
+
     tuples = [
         (np.datetime64("2000-01-01T00"), np.datetime64("2000-01-01")),
     ]
@@ -141,6 +153,11 @@ def image_groups():
             dims=("spec", "y", "x"),
             attrs={"nodata": 255},
         ),
+        "water_season": xr.DataArray(
+            da.from_array(water_season, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
     }
 
     xx = xr.Dataset(data_vars=data_vars, coords=coords)
@@ -149,9 +166,9 @@ def image_groups():
 
 
 def test_l4_classes(image_groups):
-    expected_l3 = [[216, 216, 215], [216, 216, 216], [215, 215, 215], [215, 215, 215]]
+    expected_l3 = [[216, 216, 215], [216, 216, 216], [220, 215, 215], [220, 220, 220]]
 
-    expected_l4 = [[95, 97, 93], [97, 96, 96], [93, 93, 93], [93, 93, 93]]
+    expected_l4 = [[ 95, 97, 93], [97, 96, 96], [100, 93, 93], [101, 101, 101]]
     stats_l4 = StatsLccsLevel4(measurements=["level3", "level4"])
     ds = stats_l4.reduce(image_groups)