Fix ensemble bug, add clear count

data/raw/intertidal_development_polygons.geojson

@@ -17,6 +17,7 @@
 { "type": "Feature", "properties": { "id": "westernport" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 145.511810576983464, -38.349447692009377 ], [ 145.511810576983464, -38.329670358499861 ], [ 145.554240362970859, -38.329670358499861 ], [ 145.554240362970859, -38.349447692009377 ], [ 145.511810576983464, -38.349447692009377 ] ] ] ] } },
 { "type": "Feature", "properties": { "id": "burdekin" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 147.521898649415476, -19.728861949887957 ], [ 147.521898649415476, -19.616742358013919 ], [ 147.625985792916424, -19.616742358013919 ], [ 147.625985792916424, -19.728861949887957 ], [ 147.521898649415476, -19.728861949887957 ] ] ] ] } },
 { "type": "Feature", "properties": { "id": "carnot" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 122.220209220116743, -17.197490231485581 ], [ 122.220209220116743, -17.101268193683236 ], [ 122.323961678268844, -17.101268193683236 ], [ 122.323961678268844, -17.197490231485581 ], [ 122.220209220116743, -17.197490231485581 ] ] ] ] } },
-{ "type": "Feature", "properties": { "id": "gulfcarpentaria4" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 136.315338389566364, -15.623326415467215 ], [ 136.315338389566364, -15.589869384809733 ], [ 136.363535790265161, -15.589869384809733 ], [ 136.363535790265161, -15.623326415467215 ], [ 136.315338389566364, -15.623326415467215 ] ] ] ] } }
+{ "type": "Feature", "properties": { "id": "gulfcarpentaria4" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 136.315338389566364, -15.623326415467215 ], [ 136.315338389566364, -15.589869384809733 ], [ 136.363535790265161, -15.589869384809733 ], [ 136.363535790265161, -15.623326415467215 ], [ 136.315338389566364, -15.623326415467215 ] ] ] ] } },
+{ "type": "Feature", "properties": { "id": "urangan" }, "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 152.913000413288557, -25.432524335826127 ], [ 152.985119628544425, -25.431438826573391 ], [ 152.988124595846813, -25.333702984459819 ], [ 152.917808360972259, -25.334789372635765 ], [ 152.897374583316434, -25.348911531348726 ], [ 152.899177563697833, -25.428182240140845 ], [ 152.913000413288557, -25.432524335826127 ] ] ] ] } }
 ]
 }

intertidal/elevation.py

@@ -101,6 +101,9 @@ def ds_to_flat(
     corr : xr.DataArray
         Correlation of NDWI pixel wetness with tide height.
     """
+    # Calculate clear count
+    clear = satellite_ds[index].notnull().sum(dim="time").rename("qa_count_clear")
+
     # Calculate frequency of wet per pixel, then threshold
     # to exclude always wet and always dry
     freq = (
@@ -126,6 +129,7 @@ def ds_to_flat(
         .dropna(dim="z", how="all")
     )
     freq = freq.stack(z=("y", "x"))
+    clear = clear.stack(z=("y", "x"))
 
     # Calculate correlations between NDWI water observations and tide
     # height. Because we are only interested in pixels with inundation
@@ -168,7 +172,7 @@ def ds_to_flat(
         f"{len(intertidal_candidates.z)} ({len(intertidal_candidates.z) * 100 / freq.count().item():.2f}%)"
     )
 
-    return flat_ds, freq, corr
+    return flat_ds, freq, corr, clear
 
 
 def rolling_tide_window(
@@ -881,10 +885,8 @@ def elevation(
         ds=satellite_ds,
         model=tide_model,
         directory=tide_model_dir,
-        ranking_points="data/raw/tide_correlations_2017-2019.geojson",
-        ensemble_top_n=3,
     )
-    
+
     # Set tide array pixels to nodata if the satellite data array pixels
     # contain nodata. This ensures that we ignore any tide observations
     # where we don't have matching satellite imagery
@@ -905,7 +907,7 @@ def elevation(
     )
     if valid_mask is not None:
         log.info(f"{run_id}: Applying valid data mask to constrain study area")
-    flat_ds, freq, corr = ds_to_flat(
+    flat_ds, freq, corr, clear = ds_to_flat(
         satellite_ds,
         min_freq=min_freq,
         max_freq=max_freq,
@@ -948,6 +950,7 @@ def elevation(
             flat_dem,  # DEM data
             freq,  # Frequency
             corr,  # Correlation
+            clear, # Clear count
         ],
     )
 
@@ -1015,9 +1018,9 @@ def elevation(
 @click.option(
     "--product_maturity",
     type=str,
-    default="provisional",
+    default="stable",
     help="Product maturity metadata to use for the output dataset. "
-    "Defaults to 'provisional', can also be 'stable'.",
+    "Defaults to 'stable', can also be 'provisional'.",
 )
 @click.option(
     "--dataset_maturity",

intertidal/exposure.py

@@ -445,8 +445,6 @@ def exposure(
         model=tide_model,
         times=time_range,
         directory=tide_model_dir,
-        ranking_points="data/raw/tide_correlations_2017-2019.geojson",
-        ensemble_top_n=3,
         resample=False,
     )
 

intertidal/io.py

@@ -790,11 +790,11 @@ def tidal_metadata(ds):
     metadata_dict["intertidal:tr_class"] = (
         "microtidal"
         if metadata_dict["intertidal:tr"] < 2
-        else "mesotidal"
-        if 2 <= metadata_dict["intertidal:tr"] <= 4
-        else "macrotidal"
-        if metadata_dict["intertidal:tr"] > 4
-        else np.nan
+        else (
+            "mesotidal"
+            if 2 <= metadata_dict["intertidal:tr"] <= 4
+            else "macrotidal" if metadata_dict["intertidal:tr"] > 4 else np.nan
+        )
     )
 
     return metadata_dict
@@ -822,9 +822,7 @@ def _ls_platform_instrument(year):
 
 def prepare_for_export(
     ds,
-    int_bands=None,
-    int_nodata=255,
-    int_dtype=np.uint8,
+    custom_dtypes=None,
     float_dtype=np.float32,
     output_location=None,
     overwrite=True,
@@ -839,16 +837,10 @@ def prepare_for_export(
     ----------
     ds : xarray.Dataset
         The dataset containing the bands to be exported.
-    int_bands : tuple or list, optional
-        A list of bands to export as integer datatype. If None, will use
-        the following list of bands: ("exposure", "extents",
-        "offset_hightide", "offset_lowtide", "spread")
-    int_nodata : int, optional
-        An integer that represents nodata values for integer bands
-        (default is 255).
-    int_dtype : string or numpy data type, optional
-        The data type to use for integer layers (default is
-        np.uint8).
+    custom_dtypes : dictionary, optional
+        An optional dictionary containing names of bands as keys,
+        and tuples in the form `(np.uint8, 255)` providing the 
+        dtype and nodata value to use for that band.
     float_dtype : string or numpy data type, optional
         The data type to use for floating point layers (default is
         np.float32).
@@ -866,12 +858,13 @@ def prepare_for_export(
     """
 
     def _prepare_band(
-        band, int_bands, int_nodata, int_dtype, float_dtype, output_location, overwrite
+        band, custom_dtypes, float_dtype, output_location, overwrite
     ):
         # Export specific bands as integer data types by first filling
         # NaN with nodata value before converting to int, then setting
         # nodata attribute on layer
-        if band.name in int_bands:
+        if band.name in custom_dtypes.keys():
+            int_dtype, int_nodata = custom_dtypes[band.name]
             band = band.fillna(int_nodata).astype(int_dtype)
             band.attrs["nodata"] = int_nodata
 
@@ -887,24 +880,25 @@ def _prepare_band(
 
         return band
 
-    # Use default list of bands to convert to integers if none provided
-    if int_bands is None:
-        int_bands = (
+    # Use default dictionary to convert to integers if none provided
+    if custom_dtypes is None:
+        custom_dtypes = {
             # Primary layers
-            "exposure",
-            "extents",
+            "exposure": (np.uint8, 255),
+            "extents": (np.uint8, 255),
             # Tide attribute layers
-            "ta_spread",
-            "ta_offset_high",
-            "ta_offset_low",
+            "ta_spread": (np.uint8, 255),
+            "ta_offset_high": (np.uint8, 255),
+            "ta_offset_low": (np.uint8, 255),
             # QA layers
-            "qa_ndwi_freq",
-        )
+            "qa_ndwi_freq": (np.uint8, 255),
+            "qa_count_clear": (np.int16, -999),
+        }
 
     # Apply to each array in the input `ds`
     return ds.apply(
         lambda x: _prepare_band(
-            x, int_bands, int_nodata, int_dtype, float_dtype, output_location, overwrite
+            x, custom_dtypes, float_dtype, output_location, overwrite
         )
     )
 

metadata/ga_s2ls_intertidal_cyear_3.odc-product.yaml

@@ -98,6 +98,13 @@ measurements:
     aliases:
       - ndwi_freq
 
+  - name: qa_count_clear
+    dtype: int16
+    units: "count"
+    nodata: -999
+    aliases:
+      - count_clear
+
 load:
   crs: "EPSG:3577"
   resolution: