diff --git a/odc/stats/plugins/l34_utils/l4_bare_gradation.py b/odc/stats/plugins/l34_utils/l4_bare_gradation.py
index da675d46..45e08804 100644
--- a/odc/stats/plugins/l34_utils/l4_bare_gradation.py
+++ b/odc/stats/plugins/l34_utils/l4_bare_gradation.py
@@ -7,28 +7,20 @@
 
 def bare_gradation(xx: xr.Dataset, bare_threshold, veg_cover):
 
+    # Address nodata
     bs_pc_50 = expr_eval(
         "where((a!=a), nodata, a)",
         {"a": xx.bs_pc_50.data},
-        name="mark_veg",
+        name="mark_bare_gradation_nodata",
         dtype="float32",
         **{"nodata": NODATA},
     )
 
-    # Map any data > 100 ---> 100
-    bs_pc_50 = expr_eval(
-        "where((a>100)&(a!=nodata), 100, a)",
-        {"a": bs_pc_50},
-        name="mark_veg",
-        dtype="uint8",
-        **{"nodata": NODATA},
-    )
-
     # 60% <= data  --> 15
     bs_mask = expr_eval(
         "where((a>=m)&(a!=nodata), 15, a)",
         {"a": bs_pc_50},
-        name="mark_veg",
+        name="mark_bare",
         dtype="uint8",
         **{"m": bare_threshold[1], "nodata": NODATA},
     )
@@ -37,7 +29,7 @@ def bare_gradation(xx: xr.Dataset, bare_threshold, veg_cover):
     bs_mask = expr_eval(
         "where((a>=m)&(a<n), 12, b)",
         {"a": bs_pc_50, "b": bs_mask},
-        name="mark_veg",
+        name="mark_very_sparse_veg",
         dtype="uint8",
         **{"m": bare_threshold[0], "n": bare_threshold[1]},
     )
@@ -46,7 +38,7 @@ def bare_gradation(xx: xr.Dataset, bare_threshold, veg_cover):
     bs_mask = expr_eval(
         "where(a<m, 10, b)",
         {"a": bs_pc_50, "b": bs_mask},
-        name="mark_veg",
+        name="mark_sparse_veg",
         dtype="uint8",
         **{"m": bare_threshold[0]},
     )
diff --git a/odc/stats/plugins/l34_utils/l4_cultivated.py b/odc/stats/plugins/l34_utils/l4_cultivated.py
index d7e4f46f..8ae8f8c3 100644
--- a/odc/stats/plugins/l34_utils/l4_cultivated.py
+++ b/odc/stats/plugins/l34_utils/l4_cultivated.py
@@ -3,95 +3,106 @@
 NODATA = 255
 
 
-def lc_l4_cultivated(l34, level3, lifeform, veg_cover):
+def lc_l4_cultivated(l34, level3, woody, veg_cover):
+    
+    woody = expr_eval(
+        "where((a!=a), nodata, a)",
+        {"a": woody.data},
+        name="mask_woody_nodata",
+        dtype="float32",
+        **{"nodata": NODATA},
+    )
 
+      # the 4-8 classes can't happen in LC since cultivated class will not be classified if vegetation doesn't exist.
+    # skip these classes in level4
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==10)&(c==1), 9, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l34},
+        "where((a==111), 1, d)",
+        {"a": level3, "d": l34},
         name="mark_cultivated",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==12)&(c==1), 10, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==113), 2, d)",
+        {"a": level3, "b": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==13)&(c==1), 11, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==114), 3, d)",
+        {"a": level3, "b": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==15)&(c==1), 12, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==10)&(c==113), 9, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
+    
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==16)&(c==1), 13, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==12)&(c==113), 10, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==10)&(c==2), 14, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==13)&(c==113), 11, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==12)&(c==2), 15, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==15)&(c==113), 12, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==13)&(c==2), 16, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==16)&(c==113), 13, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==15)&(c==2), 17, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==10)&(c==114), 14, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==16)&(c==2), 18, d)",
-        {"a": level3, "b": veg_cover, "c": lifeform, "d": l4},
+        "where((a==111)&(b==12)&(c==114), 15, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==1), 2, d)",
-        {"a": level3, "b": lifeform, "d": l4},
+        "where((a==111)&(b==13)&(c==114), 16, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==111)&(b==2), 3, d)",
-        {"a": level3, "b": lifeform, "d": l4},
+        "where((a==111)&(b==15)&(c==114), 17, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
 
-    # the 4-8 classes can't happen in LC since cultivated class will not be classified if vegetation doesn't exist.
-    # skip these classes in level4
-
     l4 = expr_eval(
-        "where((d==110)&(a==111), 1, d)",
-        {"a": level3, "d": l4},
+        "where((a==111)&(b==16)&(c==114), 18, d)",
+        {"a": level3, "b": veg_cover, "c": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
 
+    
+
+  
+
     return l4
diff --git a/odc/stats/plugins/l34_utils/l4_natural_veg.py b/odc/stats/plugins/l34_utils/l4_natural_veg.py
index e21c4edc..4b40a5ca 100644
--- a/odc/stats/plugins/l34_utils/l4_natural_veg.py
+++ b/odc/stats/plugins/l34_utils/l4_natural_veg.py
@@ -3,128 +3,140 @@
 NODATA = 255
 
 
-def lc_l4_natural_veg(l4, l3, lifeform, veg_cover):
+def lc_l4_natural_veg(l4, l3, woody, veg_cover):
 
+    woody = expr_eval(
+        "where((a!=a), nodata, a)",
+        {"a": woody.data},
+        name="mask_woody_nodata",
+        dtype="float32",
+        **{"nodata": NODATA},
+    )
+    
     l4 = expr_eval(
-        "where((a==110)&(b==nodata), nodata, a)",
+        "where((b==nodata), nodata, a)",
         {"a": l4, "b": l3},
         name="mark_cultivated",
         dtype="uint8",
         **{"nodata": NODATA},
     )
+    
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==10)&(b==1), 27, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112), 19, d)",
+        {"a": l3, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+    
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==12)&(b==1), 28, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(b==113), 20, d)",
+        {"a": l3, "b": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==13)&(b==1), 29, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(b==114), 21, d)",
+        {"a": l3, "b": woody, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==15)&(b==1), 30, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==10), 22, d)",
+        {"a": l3, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==16)&(b==1), 31, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==12), 23, d)",
+        {"a": l3, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==10)&(b==2), 32, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==13), 24, d)",
+        {"a": l3, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==12)&(b==2), 33, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==15), 25, d)",
+        {"a": l3, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==13)&(b==2), 34, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==16), 26, d)",
+        {"a": l3, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==15)&(b==2), 35, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==10)&(b==113), 27, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==16)&(b==2), 36, d)",
-        {"a": l3, "b": lifeform, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==12)&(b==113), 28, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(b==1), 20, d)",
-        {"a": l3, "b": lifeform, "d": l4},
+        "where((a==112)&(c==13)&(b==113), 29, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
-
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(b==2), 21, d)",
-        {"a": l3, "b": lifeform, "d": l4},
+        "where((a==112)&(c==15)&(b==113), 30, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
 
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==10), 22, d)",
-        {"a": l3, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==16)&(b==113), 31, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==12), 23, d)",
-        {"a": l3, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==10)&(b==114), 32, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==13), 24, d)",
-        {"a": l3, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==12)&(b==114), 33, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==15), 25, d)",
-        {"a": l3, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==13)&(b==114), 34, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((d==110)&(a==112)&(c==16), 26, d)",
-        {"a": l3, "c": veg_cover, "d": l4},
+        "where((a==112)&(c==15)&(b==114), 35, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
+
     l4 = expr_eval(
-        "where((d==110)&(a==112), 19, d)",
-        {"a": l3, "d": l4},
+        "where((a==112)&(c==16)&(b==114), 36, d)",
+        {"a": l3, "b": woody, "c": veg_cover, "d": l4},
         name="mark_cultivated",
         dtype="uint8",
     )
 
+
+
     return l4
diff --git a/odc/stats/plugins/l34_utils/l4_surface.py b/odc/stats/plugins/l34_utils/l4_surface.py
index 6ac1cb68..cb96e0cb 100644
--- a/odc/stats/plugins/l34_utils/l4_surface.py
+++ b/odc/stats/plugins/l34_utils/l4_surface.py
@@ -4,33 +4,34 @@
 def lc_l4_surface(l4, level3, bare_gradation):
 
     l4 = expr_eval(
-        "where((a==210)&(b==10)&(c==216), 95, a)",
-        {"a": l4, "b": bare_gradation, "c": level3},
+        "where((c==215), 93, a)",
+        {"a": l4, "c": level3},
         name="mark_surface",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((a==210)&(b==12)&(c==216), 96, a)",
-        {"a": l4, "b": bare_gradation, "c": level3},
+        "where((c==216), 94, a)",
+        {"a": l4, "c": level3},
         name="mark_surface",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((a==210)&(b==15)&(c==216), 97, a)",
+        "where((b==10)&(c==216), 95, a)",
         {"a": l4, "b": bare_gradation, "c": level3},
         name="mark_surface",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((a==210)&(c==215), 93, a)",
-        {"a": l4, "c": level3},
+        "where((b==12)&(c==216), 96, a)",
+        {"a": l4, "b": bare_gradation, "c": level3},
         name="mark_surface",
         dtype="uint8",
     )
     l4 = expr_eval(
-        "where((a==210)&(c==216), 94, a)",
-        {"a": l4, "c": level3},
+        "where((b==15)&(c==216), 97, a)",
+        {"a": l4, "b": bare_gradation, "c": level3},
         name="mark_surface",
         dtype="uint8",
     )
+    
     return l4
diff --git a/odc/stats/plugins/l34_utils/l4_veg_cover.py b/odc/stats/plugins/l34_utils/l4_veg_cover.py
index 2346ed76..60bc8201 100644
--- a/odc/stats/plugins/l34_utils/l4_veg_cover.py
+++ b/odc/stats/plugins/l34_utils/l4_veg_cover.py
@@ -16,17 +16,6 @@ def canopyco_veg_con(xx: xr.Dataset, veg_threshold):
         **{"nodata": NODATA},
     )
 
-    # Map any data > 100 ---> 100
-    pv_pc_50 = expr_eval(
-        "where((a>100) & (a!=nodata), 100, a)",
-        {
-            "a": pv_pc_50,
-        },
-        name="mark_veg",
-        dtype="uint8",
-        **{"nodata": NODATA},
-    )
-
     # data < 1 ---> 0
     veg_mask = expr_eval(
         "where(a<m, 0, a)",
diff --git a/odc/stats/plugins/l34_utils/l4_water.py b/odc/stats/plugins/l34_utils/l4_water.py
index c1d32144..478dd4be 100644
--- a/odc/stats/plugins/l34_utils/l4_water.py
+++ b/odc/stats/plugins/l34_utils/l4_water.py
@@ -3,7 +3,8 @@
 NODATA = 255
 
 
-def water_classification(xx, intertidal_mask, water_persistence):
+def water_classification(xx, water_persistence):
+
 
     # Replace nan with nodata
     l4 = expr_eval(
@@ -13,6 +14,14 @@ def water_classification(xx, intertidal_mask, water_persistence):
         dtype="uint8",
         **{"nodata": NODATA},
     )
+    
+    intertidal_mask = expr_eval(
+        "where(a==_u, 1, 0)",
+        {"a": l4},
+        name="mask_intertidal",
+        dtype="uint8",
+        **{"_u": 223},
+    )
 
     l4 = expr_eval(
         "where((a==223)|(a==221), 98, a)", {"a": l4}, name="mark_water", dtype="uint8"
diff --git a/odc/stats/plugins/l34_utils/lc_intertidal_mask.py b/odc/stats/plugins/l34_utils/lc_intertidal_mask.py
deleted file mode 100644
index ba44033b..00000000
--- a/odc/stats/plugins/l34_utils/lc_intertidal_mask.py
+++ /dev/null
@@ -1,18 +0,0 @@
-import xarray as xr
-
-from odc.stats._algebra import expr_eval
-
-NODATA = 255
-
-
-def intertidal_mask(xx: xr.Dataset):
-
-    res = expr_eval(
-        "where(a==_u, 1, 0)",
-        {"a": xx.level_3_4.data},
-        name="mask_intertidal",
-        dtype="uint8",
-        **{"_u": 223},
-    )
-
-    return res
diff --git a/odc/stats/plugins/l34_utils/lc_level3.py b/odc/stats/plugins/l34_utils/lc_level3.py
index f8071a50..c9cd646c 100644
--- a/odc/stats/plugins/l34_utils/lc_level3.py
+++ b/odc/stats/plugins/l34_utils/lc_level3.py
@@ -51,4 +51,12 @@ def lc_level3(xx: xr.Dataset):
         dtype="uint8",
     )
 
+    # Combine woody and herbaceous aquatic vegetation
+    res = expr_eval(
+        "where((a==124)|(a==125), 124, b)",
+        {"a": xx.level_3_4.data, "b": res},
+        name="mark_aquatic_veg",
+        dtype="uint8",
+    )
+
     return res
diff --git a/odc/stats/plugins/l34_utils/lc_lifeform.py b/odc/stats/plugins/l34_utils/lc_lifeform.py
deleted file mode 100644
index 53378be1..00000000
--- a/odc/stats/plugins/l34_utils/lc_lifeform.py
+++ /dev/null
@@ -1,33 +0,0 @@
-from odc.stats._algebra import expr_eval
-import xarray as xr
-
-NODATA = 255
-
-
-def lifeform(xx: xr.Dataset):
-
-    # 113 ----> 1 woody
-    # 114 ----> 2 herbaceous
-
-    lifeform_mask = expr_eval(
-        "where((a!=a)|(a>=nodata), nodata, a)",
-        {"a": xx.woody.data},
-        name="mark_lifeform",
-        dtype="float32",
-        **{"nodata": NODATA},
-    )
-
-    lifeform_mask = expr_eval(
-        "where(a==113, 1, a)",
-        {"a": lifeform_mask},
-        name="mark_lifeform",
-        dtype="uint8",
-    )
-    lifeform_mask = expr_eval(
-        "where(a==114, 2, a)",
-        {"a": lifeform_mask},
-        name="mark_lifeform",
-        dtype="uint8",
-    )
-
-    return lifeform_mask
diff --git a/odc/stats/plugins/lc_level34.py b/odc/stats/plugins/lc_level34.py
index 69a51b5d..b6b30ea4 100644
--- a/odc/stats/plugins/lc_level34.py
+++ b/odc/stats/plugins/lc_level34.py
@@ -19,8 +19,6 @@
     l4_surface,
     l4_bare_gradation,
     l4_water,
-    lc_lifeform,
-    lc_intertidal_mask,
 )
 
 
@@ -60,10 +58,8 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
             xx, self.watper_threshold
         )
 
-        intertidal_mask = lc_intertidal_mask.intertidal_mask(xx)
-
         # #TODO WATER (99-104)
-        l4 = l4_water.water_classification(xx, intertidal_mask, water_persistence)
+        l4 = l4_water.water_classification(xx, water_persistence)
 
         # Generate Level3 classes
         level3 = lc_level3.lc_level3(xx)
@@ -71,14 +67,11 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
         # Vegetation cover
         veg_cover = l4_veg_cover.canopyco_veg_con(xx, self.veg_threshold)
 
-        # Define life form
-        lifeform = lc_lifeform.lifeform(xx)
-
         # Apply cultivated Level-4 classes (1-18)
-        l4 = l4_cultivated.lc_l4_cultivated(l4, level3, lifeform, veg_cover)
+        l4 = l4_cultivated.lc_l4_cultivated(l4, level3, xx.woody, veg_cover)
 
         # Apply terrestrial vegetation classes [19-36]
-        l4 = l4_natural_veg.lc_l4_natural_veg(l4, level3, lifeform, veg_cover)
+        l4 = l4_natural_veg.lc_l4_natural_veg(l4, level3, xx.woody, veg_cover)
 
         # Bare gradation
         bare_gradation = l4_bare_gradation.bare_gradation(
diff --git a/tests/test_lc_l4_ctv.py b/tests/test_lc_l4_ctv.py
index 8592c9dc..6f0035ab 100644
--- a/tests/test_lc_l4_ctv.py
+++ b/tests/test_lc_l4_ctv.py
@@ -7,7 +7,6 @@
     l4_cultivated,
     lc_level3,
     l4_veg_cover,
-    lc_lifeform,
 )
 
 import pandas as pd
@@ -131,10 +130,8 @@ def test_ctv_classes_woody():
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
 
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
-
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
 
     assert (l4_ctv.compute() == expected_cultivated_classes).all()
 
@@ -211,10 +208,9 @@ def test_ctv_classes_herbaceous():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
     assert (l4_ctv.compute() == expected_cultivated_classes).all()
 
 
@@ -290,10 +286,9 @@ def test_ctv_classes_woody_herbaceous():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
-    veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
+    veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
 
     assert (l4_ctv.compute() == expected_cultivated_classes).all()
 
@@ -370,8 +365,8 @@ def test_ctv_classes_no_vegcover():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
+   
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
-
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
+    
     assert (l4_ctv.compute() == expected_cultivated_classes).all()
diff --git a/tests/test_lc_l4_natural_surface.py b/tests/test_lc_l4_natural_surface.py
index 81306a18..819155cf 100644
--- a/tests/test_lc_l4_natural_surface.py
+++ b/tests/test_lc_l4_natural_surface.py
@@ -15,7 +15,6 @@
     l4_natural_aquatic,
     l4_surface,
     l4_bare_gradation,
-    lc_lifeform,
 )
 
 import pandas as pd
@@ -191,14 +190,13 @@ def test_ns():
     )
 
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
 
     veg_threshold = [1, 4, 15, 40, 65, 100]
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, veg_threshold)
 
     # Apply cultivated to match the code in Level4 processing
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
-    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
+    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, xx.woody, veg_cover)
 
     l4_ctv_ntv_nav = l4_natural_aquatic.natural_auquatic_veg(
         l4_ctv_ntv, veg_cover, xx.water_season
diff --git a/tests/test_lc_l4_nav.py b/tests/test_lc_l4_nav.py
index 87022208..f91ac6fb 100644
--- a/tests/test_lc_l4_nav.py
+++ b/tests/test_lc_l4_nav.py
@@ -13,7 +13,6 @@
     l4_veg_cover,
     l4_natural_veg,
     l4_natural_aquatic,
-    lc_lifeform,
 )
 
 import pandas as pd
@@ -169,12 +168,11 @@ def test_ntv_classes_woody_herbaceous():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
     # Apply cultivated to match the code in Level4 processing
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
-    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
+    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, xx.woody, veg_cover)
 
     l4_ctv_ntv_nav = l4_natural_aquatic.natural_auquatic_veg(
         l4_ctv_ntv, veg_cover, xx.water_season
@@ -280,12 +278,11 @@ def test_ntv_herbaceous_seasonal_water_veg_cover():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
     # Apply cultivated to match the code in Level4 processing
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
-    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
+    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, xx.woody, veg_cover)
 
     l4_ctv_ntv_nav = l4_natural_aquatic.natural_auquatic_veg(
         l4_ctv_ntv, veg_cover, xx.water_season
@@ -391,12 +388,11 @@ def test_ntv_woody_seasonal_water_veg_cover():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
     # Apply cultivated to match the code in Level4 processing
-    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, lifeform, veg_cover)
-    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, lifeform, veg_cover)
+    l4_ctv = l4_cultivated.lc_l4_cultivated(xx.level_3_4, level3, xx.woody, veg_cover)
+    l4_ctv_ntv = l4_natural_veg.lc_l4_natural_veg(l4_ctv, level3, xx.woody, veg_cover)
 
     l4_ctv_ntv_nav = l4_natural_aquatic.natural_auquatic_veg(
         l4_ctv_ntv, veg_cover, xx.water_season
diff --git a/tests/test_lc_l4_ntv.py b/tests/test_lc_l4_ntv.py
index 9cea5d55..0d89005e 100644
--- a/tests/test_lc_l4_ntv.py
+++ b/tests/test_lc_l4_ntv.py
@@ -11,7 +11,6 @@
     lc_level3,
     l4_veg_cover,
     l4_natural_veg,
-    lc_lifeform,
 )
 
 import pandas as pd
@@ -134,9 +133,9 @@ def test_ntv_classes_herbaceous():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
+
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
-    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, xx.woody, veg_cover)
     assert (l4_ntv.compute() == expected_natural_terrestrial_veg_classes).all()
 
 
@@ -212,10 +211,9 @@ def test_ntv_classes_woody():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
 
-    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, xx.woody, veg_cover)
     assert (l4_ntv.compute() == expected_natural_terrestrial_veg_classes).all()
 
 
@@ -291,9 +289,9 @@ def test_ntv_classes_no_veg():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
+
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
-    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, xx.woody, veg_cover)
     assert (l4_ntv.compute() == expected_natural_terrestrial_veg_classes).all()
 
 
@@ -369,7 +367,7 @@ def test_ntv_classes_no_lifeform():
 
     stats_l4 = StatsLccsLevel4()
     level3 = lc_level3.lc_level3(xx)
-    lifeform = lc_lifeform.lifeform(xx)
+
     veg_cover = l4_veg_cover.canopyco_veg_con(xx, stats_l4.veg_threshold)
-    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, lifeform, veg_cover)
+    l4_ntv = l4_natural_veg.lc_l4_natural_veg(xx.level_3_4, level3, xx.woody, veg_cover)
     assert (l4_ntv.compute() == expected_natural_terrestrial_veg_classes).all()
diff --git a/tests/test_lc_l4_water.py b/tests/test_lc_l4_water.py
index 3a70db88..b612ee50 100644
--- a/tests/test_lc_l4_water.py
+++ b/tests/test_lc_l4_water.py
@@ -10,7 +10,6 @@
 from odc.stats.plugins.l34_utils import (
     l4_water_persistence,
     l4_water,
-    lc_intertidal_mask,
 )
 
 import pandas as pd
@@ -163,7 +162,6 @@ def test_water_classes():
     )
 
     stats_l4 = StatsLccsLevel4()
-    intertidal_mask = lc_intertidal_mask.intertidal_mask(xx)
 
     # Water persistence
     water_persistence = l4_water_persistence.water_persistence(
@@ -171,7 +169,7 @@ def test_water_classes():
     )
 
     l4_water_classes = l4_water.water_classification(
-        xx, intertidal_mask, water_persistence
+        xx, water_persistence
     )
 
     assert (l4_water_classes.compute() == expected_water_classes).all()
@@ -272,7 +270,6 @@ def test_water_intertidal():
     )
 
     stats_l4 = StatsLccsLevel4()
-    intertidal_mask = lc_intertidal_mask.intertidal_mask(xx)
 
     # Water persistence
     water_persistence = l4_water_persistence.water_persistence(
@@ -280,7 +277,7 @@ def test_water_intertidal():
     )
 
     l4_water_classes = l4_water.water_classification(
-        xx, intertidal_mask, water_persistence
+        xx, water_persistence
     )
 
     assert (l4_water_classes.compute() == expected_water_classes).all()
diff --git a/tests/test_lc_water_seasonality.py b/tests/test_lc_water_seasonality.py
new file mode 100644
index 00000000..397fe761
--- /dev/null
+++ b/tests/test_lc_water_seasonality.py
@@ -0,0 +1,159 @@
+from odc.stats.plugins.lc_level34 import StatsLccsLevel4
+import numpy as np
+import pandas as pd
+import xarray as xr
+import dask.array as da
+
+import pytest
+
+
+NODATA = 255
+
+
+@pytest.fixture(scope="module")
+def image_groups():
+    l34 = np.array(
+        [
+            [
+                [210, 210, 210],
+                [210, 210, 210],
+                [210, 210, 210],
+                [210, 210, 210],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    urban = np.array(
+        [
+            [
+                [216, 216, 215],
+                [216, 216, 216],
+                [215, 215, 215],
+                [215, 215, 215],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    woody = np.array(
+        [
+            [
+                [113, 113, 113],
+                [113, 113, 255],
+                [114, 114, 114],
+                [114, 114, 255],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    pv_pc_50 = np.array(
+        [
+            [
+                [1, 64, 65],
+                [66, 40, 41],
+                [3, 61, 78],
+                [4, 23, 42],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    bs_pc_50 = np.array(
+        [
+            [
+                [1, 64, NODATA],
+                [66, 40, 41],
+                [1, 40, 66],
+                [NODATA, 1, 42],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    cultivated = np.array(
+        [
+            [
+                [255, 255, 255],
+                [255, 255, 255],
+                [255, 255, 255],
+                [255, 255, 255],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    water_frequency = np.array(
+        [
+            [
+                [1, 3, 2],
+                [4, 5, 6],
+                [9, 2, 11],
+                [10, 11, 12],
+            ]
+        ],
+        dtype="uint8",
+    )
+
+    tuples = [
+        (np.datetime64("2000-01-01T00"), np.datetime64("2000-01-01")),
+    ]
+    index = pd.MultiIndex.from_tuples(tuples, names=["time", "solar_day"])
+    coords = {
+        "x": np.linspace(10, 20, l34.shape[2]),
+        "y": np.linspace(0, 5, l34.shape[1]),
+        "spec": index,
+    }
+
+    data_vars = {
+        "classes_l3_l4": xr.DataArray(
+            da.from_array(l34, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "urban_classes": xr.DataArray(
+            da.from_array(urban, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "cultivated_class": xr.DataArray(
+            da.from_array(cultivated, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "woody_cover": xr.DataArray(
+            da.from_array(woody, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "pv_pc_50": xr.DataArray(
+            da.from_array(pv_pc_50, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "bs_pc_50": xr.DataArray(
+            da.from_array(bs_pc_50, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+        "water_frequency": xr.DataArray(
+            da.from_array(water_frequency, chunks=(1, -1, -1)),
+            dims=("spec", "y", "x"),
+            attrs={"nodata": 255},
+        ),
+    }
+
+    xx = xr.Dataset(data_vars=data_vars, coords=coords)
+    return xx
+
+
+def test_l4_classes(image_groups):
+    expected_l3 = [[216, 216, 215], [216, 216, 216], [215, 215, 215], [215, 215, 215]]
+
+    expected_l4 = [[95, 97, 93], [97, 96, 96], [93, 93, 93], [93, 93, 93]]
+    stats_l4 = StatsLccsLevel4()
+    ds = stats_l4.reduce(image_groups)
+
+    assert (ds.level3.compute() == expected_l3).all()
+    assert (ds.level4.compute() == expected_l4).all()