add urban mask for landcover

odc/stats/plugins/_utils.py

@@ -0,0 +1,26 @@
+import dask
+import fiona
+from rasterio import features
+
+
+def rasterize_vector_mask(shape_file, transform, dst_shape, threshold=None):
+    with fiona.open(shape_file) as source_ds:
+        geoms = [s["geometry"] for s in source_ds]
+
+    mask = features.rasterize(
+        geoms,
+        transform=transform,
+        out_shape=dst_shape[1:],
+        all_touched=False,
+        fill=0,
+        default_value=1,
+        dtype="uint8",
+    )
+
+    # if valid area >= threshold
+    # then the whole tile is valid
+    if threshold is not None:
+        if mask.sum() > mask.size * threshold:
+            return dask.array.ones(dst_shape, name=False)
+
+    return dask.array.from_array(mask.reshape(dst_shape), name=False)

odc/stats/plugins/l34_utils/lc_level3.py

@@ -5,7 +5,7 @@
 NODATA = 255
 
 
-def lc_level3(xx: xr.Dataset):
+def lc_level3(xx: xr.Dataset, urban_mask):
 
     # Cultivated pipeline applies a mask which feeds only terrestrial veg (110) to the model
     # Just exclude no data (255 or nan) and apply the cultivated results
@@ -23,10 +23,11 @@ def lc_level3(xx: xr.Dataset):
     # Mask urban results with bare sfc (210)
 
     res = expr_eval(
-        "where(a==_u, b, a)",
+        "where((a==_u)&(c>0), b, a)",
         {
             "a": res,
-            "b": xx.urban_classes.data,
+            "b": xx.artificial_surface.data,
+            "c": urban_mask,
         },
         name="mark_urban",
         dtype="float32",

odc/stats/plugins/lc_level34.py

@@ -6,8 +6,10 @@
 
 import numpy as np
 import xarray as xr
+import os
 
 from ._registry import StatsPluginInterface, register
+from ._utils import rasterize_vector_mask
 
 from .l34_utils import (
     l4_water_persistence,
@@ -36,13 +38,21 @@ class StatsLccsLevel4(StatsPluginInterface):
 
     def __init__(
         self,
+        urban_mask: str = None,
+        mask_threshold: Optional[float] = None,
         veg_threshold: Optional[List] = None,
         bare_threshold: Optional[List] = None,
         watper_threshold: Optional[List] = None,
         water_seasonality_threshold: int = None,
         **kwargs,
     ):
         super().__init__(**kwargs)
+        if urban_mask is None:
+            raise ValueError("Missing urban mask shapefile")
+        if not os.path.exists(urban_mask):
+            raise FileNotFoundError(f"{urban_mask} not found")
+        self.urban_mask = urban_mask
+        self.mask_threshold = mask_threshold
 
         self.veg_threshold = (
             veg_threshold if veg_threshold is not None else [1, 4, 15, 40, 65, 100]
@@ -58,6 +68,7 @@ def __init__(
     def fuser(self, xx):
         return xx
 
+    # pylint: disable=too-many-locals
     def reduce(self, xx: xr.Dataset) -> xr.Dataset:
 
         # Water persistence
@@ -75,7 +86,13 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
         l4 = l4_water.water_classification(xx, intertidal_mask, water_persistence)
 
         # Generate Level3 classes
-        level3 = lc_level3.lc_level3(xx)
+        urban_mask = rasterize_vector_mask(
+            self.urban_mask,
+            xx.geobox.transform,
+            xx.artificial_surface.shape,
+            threshold=self.mask_threshold,
+        )
+        level3 = lc_level3.lc_level3(xx, urban_mask)
 
         # Vegetation cover
         veg_cover = l4_veg_cover.canopyco_veg_con(xx, self.veg_threshold)

odc/stats/plugins/mangroves.py

@@ -9,10 +9,9 @@
 import dask
 import os
 from odc.algo import keep_good_only, erase_bad
-import fiona
-from rasterio import features
 
 from ._registry import StatsPluginInterface, register
+from ._utils import rasterize_vector_mask
 
 NODATA = 255
 
@@ -45,22 +44,6 @@ def measurements(self) -> Tuple[str, ...]:
         _measurements = ["canopy_cover_class"]
         return _measurements
 
-    def rasterize_mangroves_extent(self, shape_file, transform, dst_shape):
-        with fiona.open(shape_file) as source_ds:
-            geoms = [s["geometry"] for s in source_ds]
-
-        mangrove_extent = features.rasterize(
-            geoms,
-            transform=transform,
-            out_shape=dst_shape[1:],
-            all_touched=False,
-            fill=0,
-            default_value=1,
-            dtype="uint8",
-        )
-
-        return dask.array.from_array(mangrove_extent.reshape(dst_shape), name=False)
-
     def fuser(self, xx):
         """
         no fuse required for mangroves since group by none
@@ -73,7 +56,7 @@ def reduce(self, xx: xr.Dataset) -> xr.Dataset:
         mangroves computation here
         it is not a 'reduce' though
         """
-        extent_mask = self.rasterize_mangroves_extent(
+        extent_mask = rasterize_vector_mask(
             self.mangroves_extent, xx.geobox.transform, xx.pv_pc_10.shape
         )
         good_data = extent_mask == 1

tests/test_lc_l34.py

@@ -3,13 +3,76 @@
 import pandas as pd
 import xarray as xr
 import dask.array as da
+import json
+import tempfile
+import os
+import fiona
+from fiona.crs import CRS
+from datacube.utils.geometry import GeoBox
+from affine import Affine
+
 
 import pytest
 
 
 NODATA = 255
 
 
+@pytest.fixture(scope="module")
+def urban_shape():
+    data = """
+    {
+   "type":"FeatureCollection",
+   "features":[
+      {
+         "geometry":{
+            "type":"Polygon",
+            "coordinates":[
+               [
+                  [
+                     0,
+                     0
+                  ],
+                  [
+                     0,
+                     100
+                  ],
+                  [
+                     100,
+                     100
+                  ],
+                  [
+                     100,
+                     0
+                  ],
+                  [
+                     0,
+                     0
+                  ]
+               ]
+            ]
+         },
+         "type":"Feature"
+      }
+   ]
+}
+    """
+    data = json.loads(data)["features"][0]
+    tmpdir = tempfile.mkdtemp()
+    filename = os.path.join(tmpdir, "test.json")
+    with fiona.open(
+        filename,
+        "w",
+        driver="GeoJSON",
+        crs=CRS.from_epsg(3577),
+        schema={
+            "geometry": "Polygon",
+        },
+    ) as dst:
+        dst.write(data)
+    return filename
+
+
 @pytest.fixture(scope="module")
 def image_groups():
     l34 = np.array(
@@ -100,18 +163,22 @@ def image_groups():
         (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]),
-    }
+
+    affine = Affine.translation(10, 0) * Affine.scale(
+        (20 - 10) / l34.shape[2], (5 - 0) / l34.shape[1]
+    )
+    geobox = GeoBox(
+        crs="epsg:3577", affine=affine, width=l34.shape[2], height=l34.shape[1]
+    )
+    coords = geobox.xr_coords()
 
     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(
+        "artificial_surface": xr.DataArray(
             da.from_array(urban, chunks=(1, -1, -1)),
             dims=("spec", "y", "x"),
             attrs={"nodata": 255},
@@ -148,11 +215,13 @@ def image_groups():
     return xx
 
 
-def test_l4_classes(image_groups):
+def test_l4_classes(image_groups, urban_shape):
     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(measurements=["level3", "level4"])
+    stats_l4 = StatsLccsLevel4(
+        measurements=["level3", "level4"], urban_mask=urban_shape, mask_threshold=0.3
+    )
     ds = stats_l4.reduce(image_groups)
 
     assert (ds.level3.compute() == expected_l3).all()