From 8df4d211e3146efbb6b0227c4a5eaf1c6a2bdb5f Mon Sep 17 00:00:00 2001
From: Yuri <yuri.shendryk@gmail.com>
Date: Fri, 14 Aug 2026 10:28:06 +1000
Subject: [PATCH] Fix ensemble_tides() re-interpolating duplicate locations

Ran into out-of-memory (OOM) errors during ensemble modelling. x/y
passed to idw() come straight from tide_df's full index, so each
location gets interpolated once per timestep and model instead of once.
Dedupe before calling idw(), not after. Memory now scales with unique
locations, not locations x times x models.
---
 eo_tides/model.py   | 18 ++++++++++-----
 tests/test_model.py | 56 +++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 68 insertions(+), 6 deletions(-)

diff --git a/eo_tides/model.py b/eo_tides/model.py
index 8495493..7fff480 100644
--- a/eo_tides/model.py
+++ b/eo_tides/model.py
@@ -397,25 +397,31 @@ def ensemble_tides(
         """
         raise Exception(textwrap.dedent(error_msg).strip()) from None
 
+    # x/y repeat once per timestep and tide model since they come from tide_df's full
+    # index, but rankings only depend on location. Interpolating the full, repeated
+    # coordinates re-runs the same IDW computation for no benefit, and can raise a
+    # MemoryError for large study areas or long time series. Interpolating only the
+    # unique locations avoids this. The join back onto tide_df by (tide_model, x, y)
+    # below already broadcasts each location's ranking to every matching row.
+    unique_xy = pd.DataFrame({"x": x, "y": y}).drop_duplicates()
+
     # Use points to interpolate model rankings into requested x and y
     id_kwargs_str = "" if idw_kwargs == {} else idw_kwargs
     print(f"Interpolating model rankings using IDW interpolation {id_kwargs_str}")
     ensemble_ranks_df = (
-        # Run IDW interpolation on subset of ranking columns
+        # Run IDW interpolation on subset of ranking columns, once per unique location
         pd.DataFrame(
             idw(
                 input_z=model_ranks_gdf[model_ranking_cols],
                 input_x=model_ranks_gdf.geometry.x,
                 input_y=model_ranks_gdf.geometry.y,
-                output_x=x,
-                output_y=y,
+                output_x=unique_xy.x,
+                output_y=unique_xy.y,
                 **idw_kwargs,
             ),
             columns=model_ranking_cols,
         )
-        .assign(x=x, y=y)
-        # Drop any duplicates then melt columns into long format
-        .drop_duplicates()
+        .assign(x=unique_xy.x.to_numpy(), y=unique_xy.y.to_numpy())
         .melt(id_vars=["x", "y"], var_name="tide_model", value_name="rank")
         # Remove "rank_" prefix to get plain model names
         .replace({"^rank_": ""}, regex=True)
diff --git a/tests/test_model.py b/tests/test_model.py
index 10ec405..d9a8f03 100644
--- a/tests/test_model.py
+++ b/tests/test_model.py
@@ -1,3 +1,4 @@
+import geopandas as gpd
 import numpy as np
 import pandas as pd
 import pytest
@@ -7,6 +8,7 @@
     _parallel_splits,
     _set_directory,
     ensemble_tides,
+    idw,
     model_phases,
     model_tides,
 )
@@ -519,6 +521,60 @@ def test_model_tides_ensemble():
     }
 
 
+def test_ensemble_tides_deduplicates_ranking_queries(monkeypatch, tmp_path):
+    """`ensemble_tides` should interpolate model rankings once per unique (x, y)
+    location, not once per (time, x, y, tide_model) row. `tide_df` repeats each
+    location once per timestep and tide model, but rankings are static per
+    location, so passing the full repeated coordinates into `idw` performs the
+    same interpolation over and over for no benefit. For a large number of
+    timesteps and/or models this scales query size unnecessarily and can raise
+    a MemoryError for large study areas or long time series."""
+    locations = [(0.0, 0.0), (10.0, 0.0), (5.0, 8.0)]
+    times = pd.date_range("2021-01-01", periods=5, freq="6h")
+    models = ["EOT20", "HAMTIDE11"]
+
+    rng = np.random.default_rng(0)
+    rows = [
+        {"time": t, "x": x, "y": y, "tide_model": m, "tide_height": rng.uniform(-1, 1)}
+        for t in times
+        for x, y in locations
+        for m in models
+    ]
+    tide_df = pd.DataFrame(rows).set_index(["time", "x", "y"])
+
+    ranking_gdf = gpd.GeoDataFrame(
+        {"rank_EOT20": [1, 2], "rank_HAMTIDE11": [2, 1], "valid_perc": [0.5, 0.5]},
+        geometry=gpd.points_from_xy([-2, 12], [-2, -1]),
+        crs="EPSG:3857",
+    )
+    ranking_path = tmp_path / "ranking.gpkg"
+    ranking_gdf.to_file(ranking_path, driver="GPKG")
+
+    # Spy on `idw` to record how many query points it was actually called with
+    query_sizes = []
+    original_idw = idw
+
+    def spy_idw(*args, **kwargs):
+        query_sizes.append(len(np.atleast_1d(kwargs["output_x"])))
+        return original_idw(*args, **kwargs)
+
+    monkeypatch.setattr("eo_tides.model.idw", spy_idw)
+
+    ensemble_tides(
+        tide_df,
+        crs="EPSG:3857",
+        ensemble_models=models,
+        ranking_points=str(ranking_path),
+        ranking_valid_perc=0.0,
+        k=2,
+    )
+
+    n_unique_locations = len(locations)
+    n_total_rows = len(locations) * len(times) * len(models)
+    assert n_total_rows > n_unique_locations  # sanity check the scenario is meaningful
+    assert query_sizes == [n_unique_locations]
+
+
 # Test ensemble dtype is set correctly
 @pytest.mark.parametrize(
     "dtype",