Merge pull request #78 from GeoscienceAustralia/seasonal_correction

Add seasonal correction to correlation calculation

intertidal/elevation.py

@@ -44,6 +44,7 @@ def ds_to_flat(
     max_freq=0.99,
     min_correlation=0.15,
     corr_method="pearson",
+    correct_seasonality=False,
     valid_mask=None,
 ):
     """
@@ -77,6 +78,13 @@ def ds_to_flat(
     corr_method : str, optional
         Correlation method to use. Defaults to "pearson", also supports
         "spearman".
+    correct_seasonality : bool, optional
+        If True, remove any seasonal signal from the tide height data
+        by subtracting monthly mean tide height from each value. This
+        can reduce false tide correlations in regions where tide heights
+        correlate with seasonal changes in surface water. Note that 
+        seasonally corrected tides are only used to identify potentially
+        tide influenced pixels - not for elevation modelling itself.
     valid_mask : xr.DataArray, optional
         A boolean mask used to optionally constrain the analysis area,
         with the same spatial dimensions as `satellite_ds`. For example,
@@ -123,13 +131,22 @@ def ds_to_flat(
     # correlation. This prevents small changes in NDWI beneath the water
     # surface from producing correlations with tide height.
     wet_dry = flat_ds[index] > ndwi_thresh
+
+    # Use either tides directly or correct to remove seasonal signal
+    if correct_seasonality:
+        print("Removing seasonal signal before calculating tide correlations")
+        gb = flat_ds.tide_m.groupby('time.month')
+        tide_array = (gb - gb.mean())
+    else:
+        tide_array = flat_ds.tide_m  
+
     if corr_method == "pearson":
-        corr = xr.corr(wet_dry, flat_ds.tide_m, dim="time").rename("qa_ndwi_corr")
+        corr = xr.corr(wet_dry, tide_array, dim="time").rename("qa_ndwi_corr")
     elif corr_method == "spearman":
         import xskillscore
 
         corr = xskillscore.spearman_r(
-            flat_ds[index], flat_ds.tide_m, dim="time", skipna=True, keep_attrs=True
+            flat_ds[index], tide_array, dim="time", skipna=True, keep_attrs=True
         ).rename("qa_ndwi_corr")
 
     # TODO: investigate alternative function from DEA Tools
@@ -752,6 +769,7 @@ def elevation(
     min_correlation=0.15,
     windows_n=100,
     window_prop_tide=0.15,
+    correct_seasonality=False,
     max_workers=None,
     tide_model="FES2014",
     tide_model_dir="/var/share/tide_models",
@@ -788,6 +806,14 @@ def elevation(
     window_prop_tide : float, optional
         Proportion of the tide range to use for each window radius in
         the per-pixel rolling median calculation, by default 0.15
+    correct_seasonality : bool, optional
+        If True, remove any seasonal signal from the tide height data
+        by subtracting monthly mean tide height from each value prior to
+        correlation calculations. This can reduce false tide correlations
+        in regions where tide heights correlate with seasonal changes in
+        surface water. Note that seasonally corrected tides are only used
+        to identify potentially tide influenced pixels - not for elevation
+        modelling itself.
     max_workers : int, optional
         Maximum number of worker processes to use for parallel execution
         in the per-pixel rolling median calculation. Defaults to None,
@@ -873,6 +899,7 @@ def elevation(
         min_freq=min_freq,
         max_freq=max_freq,
         min_correlation=min_correlation,
+        correct_seasonality=correct_seasonality,
         valid_mask=valid_mask,
     )
 
@@ -1172,6 +1199,7 @@ def intertidal_cli(
             min_correlation=min_correlation,
             windows_n=windows_n,
             window_prop_tide=window_prop_tide,
+            correct_seasonality=True,
             tide_model=tide_model,
             tide_model_dir=tide_model_dir,
             run_id=run_id,

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 **2024-03-14 12:50**. Compared to the previous run, it had an:
+The latest integration test completed at **2024-03-22 14:25**. 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

@@ -47,3 +47,4 @@ time,Correlation,RMSE,MAE,R-squared,Bias,Regression slope
 2024-03-13 00:54:01.731099+00:00,0.975,0.141,0.121,0.95,0.116,1.11
 2024-03-13 23:17:46.582372+00:00,0.975,0.141,0.121,0.95,0.116,1.11
 2024-03-14 01:50:20.512235+00:00,0.975,0.141,0.121,0.95,0.116,1.11
+2024-03-22 03:25:50.523558+00:00,0.975,0.141,0.121,0.95,0.116,1.11