Update graph func to use centroid observations (#133)

* Update graph func to use centroid observations

* Comment

* Automatically update integration test validation results

* Use pragmatic filter on clear low and high tide obs

* Automatically update integration test validation results

---------

Co-authored-by: robbibt <robbibt@users.noreply.github.com>

Author: robbibt

intertidal/composites.py

@@ -42,10 +42,14 @@ def tidal_thresholds(
     # Calculate per-pixel integer rankings for each tide height
     rank_n = tides_highres.rank(dim="time")
 
-    # Calculate low and high ranking thresholds from total rankings.
-    # Low threshold needs to be rounded up ("ceil"), and high tide
-    # rounded down ("floor") to ensure we capture all matching values.
-    rank_max = rank_n.max(dim="time")
+    # Calculate pixel-based low and high ranking thresholds from
+    # max ranking. Max ranking needs to be rounded up to the nearest
+    # integer using "ceil" as xarray will give multiple observation
+    # an average rank (e.g. 50.5) value if they are both identical.
+    # Additionally: to ensure we capture all matching values, Low
+    # threshold needs to be rounded up ("ceil"), and high tide
+    # rounded down ("floor").
+    rank_max = np.ceil(rank_n.max(dim="time"))
     rank_thresh_low = np.ceil(rank_max * threshold_lowtide)
     rank_thresh_high = np.floor(rank_max * threshold_hightide)
 
@@ -169,13 +173,6 @@ def tidal_composites(
     log.info(
         f"{run_id}: Calculating low and high tide thresholds with minimum {min_obs} observations"
     )
-    # threshold_ds = xr_quantile(
-    #     src=tides_highres.to_dataset(),
-    #     quantiles=[threshold_lowtide, threshold_hightide],
-    #     nodata=np.nan,
-    # )
-    # low_threshold = threshold_ds.isel(quantile=0).tide_height.drop("quantile")
-    # high_threshold = threshold_ds.isel(quantile=-1).tide_height.drop("quantile")
     low_threshold, high_threshold = tidal_thresholds(
         tides_highres=tides_highres,
         threshold_lowtide=threshold_lowtide,
@@ -188,9 +185,9 @@ def tidal_composites(
     low_mask = tides_highres <= low_threshold
     high_mask = tides_highres >= high_threshold
 
-    # Keep only scenes with at least some valid data to speed up geomedian
-    low_keep = low_mask.any(dim=["x", "y"])
-    high_keep = high_mask.any(dim=["x", "y"])
+    # Keep only scenes with at least 1% valid data to speed up geomedian
+    low_keep = low_mask.mean(dim=["x", "y"]) >= 0.01
+    high_keep = high_mask.mean(dim=["x", "y"]) >= 0.01
     ds_low = satellite_ds.sel(time=low_keep)
     ds_high = satellite_ds.sel(time=high_keep)
 
@@ -236,7 +233,7 @@ def tidal_composites(
     ds_lowtide["qa_low_threshold"] = low_threshold
     ds_hightide["qa_high_threshold"] = high_threshold
 
-    return ds_lowtide, ds_hightide, low_keep, high_keep
+    return ds_lowtide, ds_hightide
 
 
 @click.command()
@@ -346,8 +343,7 @@ def tidal_composites(
     "--include_coastal_aerosol/--no-include_coastal_aerosol",
     type=bool,
     default=True,
-    help="Whether to include the coastal aerosol band. "
-    "Defaults to True",
+    help="Whether to include the coastal aerosol band. Defaults to True",
 )
 @click.option(
     "--eps",
@@ -493,7 +489,7 @@ def tidal_composites_cli(
 
             # Calculate high and low tide geomedian composites
             log.info(f"{run_id}: Running DEA Tidal Composites workflow")
-            ds_lowtide, ds_hightide, low_keep, high_keep = tidal_composites(
+            ds_lowtide, ds_hightide = tidal_composites(
                 satellite_ds=satellite_ds,
                 threshold_lowtide=threshold_lowtide,
                 threshold_hightide=threshold_hightide,
@@ -554,20 +550,12 @@ def tidal_composites_cli(
                 log=log,
             )
 
-            # Add new array to data to label high or low tide
-            # images selected for geomedian analysis.
-            label = xr.where(
-                high_keep | low_keep,
-                "Clear low and high tide images",
-                "All satellite observations",
-            )
-            satellite_ds = satellite_ds.assign(label=label)
-
             # Calculate additional tile-level tidal metadata and graph.
             metadata_dict, tide_graph_fig = tidal_metadata(
                 product_family="tidal_composites",
+                threshold_lowtide=threshold_lowtide,
+                threshold_hightide=threshold_hightide,
                 data=satellite_ds,
-                plot_var="label",
                 modelled_freq="30min",
                 model=tide_model,
                 directory=tide_model_dir,

intertidal/io.py

@@ -785,15 +785,27 @@ def _write_stac(
     return stac
 
 
-def tidal_metadata(product_family, **tide_stats_kwargs):
+def tidal_metadata(
+    product_family,
+    threshold_lowtide=0.15,
+    threshold_hightide=0.85,
+    **tide_stats_kwargs,
+):
     """
     Generate tidal statistics and tide bias plot for a given input tile.
     Tidal statistics are calculated based on the centroid of the tile.
 
+    For `product_family=="tidal_composites"`, observations matching the
+    low and high tide thresholds will be plotted in white.
+
     Parameters
     ----------
     product_family : string
         Either "intertidal" or "tidal_composites".
+    threshold_lowtide : float, optional
+        Quantile used to identify low tide observations, by default 0.15.
+    threshold_hightide : float, optional
+        Quantile used to identify high tide observations, by default 0.85.
     **tide_stats_kwargs :
         Any required parameters to pass to `eo_tides.stats.tide_stats`,
         e.g. `data`, `model`, `directory` etc.
@@ -831,36 +843,50 @@ def tidal_metadata(product_family, **tide_stats_kwargs):
     )
 
     # Update figure line and point colours
-    if product_family == "intertidal":
-        fig.axes[0].get_lines()[0].set_color("#90b7d8")
-        fig.axes[0].get_lines()[0].set_alpha(1.0)
-        fig.axes[0].get_lines()[1].set_color("black")
-        fig.axes[0].get_lines()[1].set_markersize(4)
-        fig.axes[0].get_lines()[1].set_markeredgecolor("none")
-
-        # HAT/LOT lines
-        fig.axes[0].get_lines()[2].set_color("none")
-        fig.axes[0].get_lines()[3].set_color("none")
-        fig.axes[0].get_lines()[4].set_color("none")
-        fig.axes[0].get_lines()[5].set_color("none")
-
-    elif product_family == "tidal_composites":
-        fig.axes[0].get_lines()[0].set_color("#90b7d8")
-        fig.axes[0].get_lines()[0].set_alpha(1.0)
-        fig.axes[0].get_lines()[1].set_markeredgecolor("none")
-        fig.axes[0].get_lines()[2].set_markeredgecolor("#343c47")
-        fig.axes[0].get_lines()[1].set_marker("o")
-        fig.axes[0].get_lines()[2].set_marker("o")
-        fig.axes[0].get_lines()[1].set_markersize(4)
-        fig.axes[0].get_lines()[2].set_markersize(5)
-        fig.axes[0].get_lines()[1].set_color("black")
-        fig.axes[0].get_lines()[2].set_color("white")
-
-        # HAT/LOT lines
-        fig.axes[0].get_lines()[3].set_color("none")
-        fig.axes[0].get_lines()[4].set_color("none")
-        fig.axes[0].get_lines()[5].set_color("none")
-        fig.axes[0].get_lines()[6].set_color("none")
+    modelled = fig.axes[0].get_lines()[0]
+    observed = fig.axes[0].get_lines()[1]
+    hat = fig.axes[0].get_lines()[2]
+    hot = fig.axes[0].get_lines()[3]
+    lot = fig.axes[0].get_lines()[4]
+    lat = fig.axes[0].get_lines()[5]
+
+    # Set styling
+    modelled.set_color("#90b7d8")
+    modelled.set_alpha(1.0)
+    observed.set_color("black")
+    observed.set_markersize(4)
+    observed.set_markeredgecolor("none")
+
+    # Remove HAT/LOT lines
+    hat.set_color("none")
+    hot.set_color("none")
+    lot.set_color("none")
+    lat.set_color("none")
+
+    # For Tidal Composites, manually set low and high
+    # tide observations to white
+    if product_family == "tidal_composites":
+
+        # Extract observed data from plot
+        xdata = observed.get_xdata()
+        ydata = observed.get_ydata()
+
+        # Calculate thresholds and plot subset of points in white
+        min_thresh, max_thresh = np.quantile(
+            ydata, [threshold_lowtide, threshold_hightide]
+        )
+        mask = (ydata <= min_thresh) | (ydata >= max_thresh)
+        fig.axes[0].plot(
+            xdata[mask],
+            ydata[mask],
+            marker="o",
+            linestyle="None",
+            color="white",
+            markersize=5,
+            markeredgecolor="#343c47",
+            markeredgewidth=0.8,
+            label="Low and high tide images",
+        )
 
     # Set background to transparent
     fig.patch.set_facecolor("#5d646c00")

tests/README.md

@@ -10,7 +10,7 @@ Integration tests
 
 This directory contains tests that are run to verify that DEA Intertidal code runs correctly. The ``test_intertidal.py`` file runs a small-scale full workflow analysis over an intertidal flat in the Gulf of Carpentaria using the DEA Intertidal [Command Line Interface (CLI) tools](../notebooks/Intertidal_CLI.ipynb), and compares these results against a LiDAR validation DEM to produce some simple accuracy metrics.
 
-The latest integration test completed at **2025-05-02 11:43**. Compared to the previous run, it had an:
+The latest integration test completed at **2025-05-02 17:14**. Compared to the previous run, it had an:
 - RMSE accuracy of **0.14 m ( :heavy_minus_sign: no change)**
 - MAE accuracy of **0.12 m ( :heavy_minus_sign: no change)**
 - Bias of **0.12 m ( :heavy_minus_sign: no change)**

tests/validation.csv

@@ -108,3 +108,5 @@ time,Correlation,RMSE,MAE,R-squared,Bias,Regression slope
 2025-04-28 02:09:27.190062+00:00,0.975,0.145,0.123,0.95,0.117,1.119
 2025-04-28 06:23:57.332550+00:00,0.975,0.145,0.123,0.95,0.117,1.119
 2025-05-02 01:43:13.002005+00:00,0.975,0.145,0.123,0.95,0.117,1.119
+2025-05-02 06:32:51.594612+00:00,0.975,0.145,0.123,0.95,0.117,1.119
+2025-05-02 07:14:19.583638+00:00,0.975,0.145,0.123,0.95,0.117,1.119