diff --git a/intertidal/composites.py b/intertidal/composites.py index 88905d9..0f798f6 100644 --- a/intertidal/composites.py +++ b/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, diff --git a/intertidal/io.py b/intertidal/io.py index 334c5ab..ded53e9 100644 --- a/intertidal/io.py +++ b/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") diff --git a/tests/README.md b/tests/README.md index 1577473..30524c1 100644 --- a/tests/README.md +++ b/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)** diff --git a/tests/validation.csv b/tests/validation.csv index 59799c6..3bf3fec 100644 --- a/tests/validation.csv +++ b/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 diff --git a/tests/validation.jpg b/tests/validation.jpg index a3f741b..4bd4dac 100644 Binary files a/tests/validation.jpg and b/tests/validation.jpg differ