Merge pull request #73 from GeoscienceAustralia/update_interpolation

Use IDW interpolation for ensemble tides

intertidal/tide_modelling.py

@@ -2,16 +2,17 @@
 import geopandas as gpd
 
 from dea_tools.coastal import pixel_tides, _pixel_tides_resample
-from dea_tools.spatial import interpolate_2d
+from dea_tools.spatial import xr_interpolate
 
 
 def pixel_tides_ensemble(
     ds,
     ancillary_points,
     model="ensemble",
     top_n=3,
-    interp_method="nearest",
+    interp_method="idw",
     reduce_method="mean",
+    ancillary_valid_perc=0.02,
     **pixel_tides_kwargs,
 ):
     """
@@ -52,6 +53,9 @@ def pixel_tides_ensemble(
     reduce_method : str, optional
         Method used to reduce values from the `top_n` tide models into
         a single enemble output. Defaults to "mean", supports "median".
+    ancillary_valid_perc : float, optional
+        The minimum valid percent used to filter input point correlations.
+        Defaults to 0.02.
     **pixel_tides_kwargs
         All other optional keyword arguments to provide to the underlying
         `pixel_tides` function.
@@ -78,10 +82,10 @@ def pixel_tides_ensemble(
     # ensemble tide modelling
     else:
         print("Running ensemble tide modelling")
-        # Extract the `resample` param if it exists so we can run 
+        # Extract the `resample` param if it exists so we can run
         # `pixel_tides` with `resample=False`, and then resample later
         resample_param = pixel_tides_kwargs.pop("resample", True)
-        
+
         # Run `pixel_tides` on all tide models and return low-res output
         ensemble_models = [
             "FES2014",
@@ -99,36 +103,25 @@ def pixel_tides_ensemble(
             **pixel_tides_kwargs,
         )
 
-    # Load ancillary points from file, reproject to match satellite
-    # data, and drop empty points
+    # Load ancillary points from file, filter by minimum valid data perc
+    # and drop any empty points/unnecessary columns
     print("Generating ensemble tide model from point inputs")
     corr_gdf = (
-        gpd.read_file(ancillary_points)[ensemble_models + ["geometry"]]
-        .to_crs(ds.odc.crs)
-        .dropna()
+        gpd.read_file(ancillary_points)
+        .query(f"valid_perc > {ancillary_valid_perc}")
+        .dropna()[ensemble_models + ["geometry"]]
     )
 
-    # Loop through each model, interpolating correlations into
-    # low-res tide grid
-    out_list = []
-
-    for model in ensemble_models:
-        out = interpolate_2d(
-            tides_lowres,
-            x_coords=corr_gdf.geometry.x,
-            y_coords=corr_gdf.geometry.y,
-            z_coords=corr_gdf[model],
-            method=interp_method,
-        ).expand_dims({"tide_model": [model]})
-
-        out_list.append(out)
+    # Spatially interpolate each tide model
+    print(f"Interpolating model weights using '{interp_method}' interpolation")
+    weights_ds = xr_interpolate(
+        tides_lowres, gdf=corr_gdf, columns=ensemble_models, method=interp_method
+    ).to_array(dim="tide_model")
 
-    # Combine along tide model dimension into a single xarray.Dataset
-    weights_ds = xr.concat(out_list, dim="tide_model")
-
     # Print models in order of correlation
     print(
-        weights_ds.mean(dim=["x", "y"])
+        weights_ds.drop("spatial_ref")
+        .mean(dim=["x", "y"])
         .to_dataframe("weights")
         .sort_values("weights", ascending=False)
     )
@@ -141,21 +134,21 @@ def pixel_tides_ensemble(
     # Reduce remaining models to produce a single ensemble output
     # for each time/x/y
     if reduce_method == "median":
-        print("Reducing multiple model outputs using 'median'")
+        print("Reducing multiple models into single ensemble model using 'median'")
         tides_lowres_ensemble = tides_top_n.median("tide_model")
     elif reduce_method == "mean":
-        print("Reducing multiple model outputs using 'mean'")
+        print("Reducing multiple models into single ensemble model using 'mean'")
         tides_lowres_ensemble = tides_top_n.mean("tide_model")
 
-    # Optionally resample/reproject ensemble tides to match high-res 
+    # Optionally resample/reproject ensemble tides to match high-res
     # satellite data
     if resample_param:
-        print("Reprojecting tides into original array")
+        print("Reprojecting ensemble tides into original array")
         tides_highres_ensemble, _ = _pixel_tides_resample(
             tides_lowres=tides_lowres_ensemble, ds=ds
         )
         return tides_highres_ensemble, tides_lowres_ensemble
 
     else:
-        print("Returning low resolution tide array")
+        print("Returning low resolution ensemble tide array")
         return tides_lowres_ensemble

requirements.in

@@ -3,7 +3,7 @@ aiohttp
 botocore
 click==8.1.3
 datacube[s3,performance]==1.8.13
-dea-tools==0.3.1
+dea-tools==0.3.2.dev49
 eodatasets3==0.30.2
 geopandas==0.13.2
 matplotlib==3.7.1

tests/README.md

@@ -10,11 +10,11 @@ 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-08 15:13**. Compared to the previous run, it had an:
-- RMSE accuracy of **0.14 m ( :heavy_check_mark: improved by 0.006)**
-- MAE accuracy of **0.12 m ( :heavy_check_mark: improved by 0.005)**
-- Bias of **0.12 m ( :heavy_check_mark: improved by 0.005)**
-- Pearson correlation of **0.975 ( :heavy_exclamation_mark: worsened by 0.002)**
+The latest integration test completed at **2024-03-12 16:16**. 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)**
+- Pearson correlation of **0.975 ( :heavy_minus_sign: no change)**
 
 
 <img src="validation.jpg" width="950"/>

tests/validation.csv

@@ -43,3 +43,4 @@ time,Correlation,RMSE,MAE,R-squared,Bias,Regression slope
 2024-03-05 05:33:39.493036+00:00,0.976,0.141,0.121,0.952,0.117,1.109
 2024-03-07 07:52:01.276183+00:00,0.977,0.147,0.126,0.955,0.121,1.125
 2024-03-08 04:13:32.992746+00:00,0.975,0.141,0.121,0.95,0.116,1.11
+2024-03-12 05:16:34.997844+00:00,0.975,0.141,0.121,0.95,0.116,1.11