again fixing scipy verision and updating testcase

requirements.txt

@@ -1,4 +1,4 @@
-scipy
+scipy<1.10.0
 numpy
 matplotlib
 pandas

setup.cfg

@@ -27,7 +27,8 @@ classifiers =
 [options]
 python_requires = >=3.8
 install_requires = 
-	scipy
+	#install scipy<1.10.0 to avoid da.interp() issue: https://github.com/Deltares/dfm_tools/issues/287 and https://github.com/pydata/xarray/issues/7701
+	scipy<1.10.0 
 	numpy
 	matplotlib
 	#install pandas<2.0.0 in case of py38 to avoid conflict with xarray<2023.3.0: https://github.com/Deltares/xugrid/issues/78#issuecomment-1597723955 (xarray 2023.1.0 is the latest py38 release, and py38 is still supported by dfm_tools)

tests/test_external_packages.py

@@ -60,50 +60,6 @@ def test_xarray_pandas_resample():
     ds.resample(time='D')
 
 
-@pytest.mark.unittest
-def test_xarray_interp_to_newdim():
-    """
-    this one fails with scipy>=1.10.0: https://github.com/pydata/xarray/issues/7701
-    """
-    ds = xr.Dataset()
-    so_np = np.array([[[35.819576, 35.82568 , 35.82873 ],
-                       [35.819576, 35.824154, 35.831783],
-                       [35.822628, 35.824154, 35.82873 ]],
-
-                      [[35.802788, 35.80584 , 35.815   ],
-                       [35.815   , 35.810417, 35.821102],
-                       [35.824154, 35.813473, 35.81805 ]],
-
-                      [[35.786003, 35.789055, np.nan],
-                       [35.807365, 35.796684, np.nan],
-                       [35.824154, 35.80584 , np.nan]],
-
-                      [[35.776848, np.nan,    np.nan],
-                       [35.792107, np.nan,    np.nan],
-                       [35.822628, np.nan,    np.nan]],
-
-                      [[35.781425, np.nan,    np.nan],
-                       [35.792107, np.nan,    np.nan],
-                       [35.789055, np.nan,    np.nan]]])
-    ds['so'] = xr.DataArray(so_np,dims=('depth','latitude','longitude'))
-    ds['longitude'] = xr.DataArray([-9.6, -9.5, -9.4], dims=('longitude'))
-    ds['latitude'] = xr.DataArray([42.9, 43.0, 43.1], dims=('latitude'))
-
-    x_xr = xr.DataArray([-9.5],dims=('plipoints'))
-    y_xr = xr.DataArray([43],dims=('plipoints'))
-
-    interp_with_floats = ds.interp(longitude=x_xr[0], latitude=y_xr[0], method='linear').so #selecting one value from the da drops the new plipoints dimension
-    interp_with_da_existing = ds.interp(longitude=x_xr.values, latitude=y_xr.values, method='linear').so.isel(longitude=0,latitude=0) #using the DataArray values keeps lat/lon dimenions, gives the same interp result
-    interp_with_da_newdim = ds.interp(longitude=x_xr, latitude=y_xr, method='linear').so.isel(plipoints=0) #using the DataArray introduces a plipoints dimension, which gives different interp result
-    print(interp_with_floats.to_numpy())
-    print(interp_with_da_existing.to_numpy())
-    print(interp_with_da_newdim.to_numpy())
-    print(xr.__version__)
-
-    assert (interp_with_floats.isnull()==interp_with_da_existing.isnull()).all() #success
-    assert (interp_with_floats.isnull()==interp_with_da_newdim.isnull()).all() #fails with scipy>=1.10.0
-
-
 @pytest.mark.unittest
 def test_xarray_decode_default_fillvals():
     """

tests/test_interpolate_grid2bnd.py

@@ -9,6 +9,7 @@
 import dfm_tools as dfmt
 import numpy as np
 import datetime as dt
+import xarray as xr
 
 
 @pytest.mark.unittest
@@ -55,6 +56,89 @@ def test_components_translate_upper():
     assert comp_list == ['M2','EPSILON2','EPSILON2']
 
 
+@pytest.mark.unittest
+def test_xarray_interp_to_newdim():
+    """
+    Linear interpolation to a new dimension in dfmt.interp_regularnc_to_plipoints() 
+    resulted in unexpected nan values since scipy 1.10.0.
+    More info in https://github.com/pydata/xarray/issues/7701 and related issues.
+    Since this interpn is way faster than interp1d on a grid or all separate xy points, 
+    we sticked to interpn but filled nan values with a nearest interpolation. 
+    This only happens if there are non-nan values in the surrounding cells.
+    In this test we check whether this combination of interpolation results in the same values as interp1d.
+    This is mainly relevant when interpolating to existing lat/lon values, so that is the case we test here.
+    This method gives as much valid values as possible given the input dataset, but does not fill nans where it should not do that.
+    """
+
+    ds = xr.Dataset()
+    so_np = np.array([[[35.819576, 35.82568 , 35.82873 ],
+                       [35.819576, 35.824154, 35.831783],
+                       [35.822628, 35.824154, 35.82873 ]],
+
+                      [[35.802788, 35.80584 , 35.815   ],
+                       [35.815   , 35.810417, 35.821102],
+                       [35.824154, 35.813473, 35.81805 ]],
+
+                      [[35.786003, 35.789055, np.nan],
+                       [35.807365, 35.796684, np.nan],
+                       [35.824154, 35.80584 , np.nan]],
+
+                      [[35.776848, np.nan,    np.nan],
+                       [35.792107, np.nan,    np.nan],
+                       [35.822628, np.nan,    np.nan]],
+
+                      [[35.781425, np.nan,    np.nan],
+                       [35.792107, np.nan,    np.nan],
+                       [35.789055, np.nan,    np.nan]]])
+    ds['so'] = xr.DataArray(so_np,dims=('depth','latitude','longitude'))
+    lons = [-9.6, -9.5, -9.4]
+    lats = [42.9, 43.0, 43.1]
+    ds['longitude'] = xr.DataArray(lons, dims=('longitude'))
+    ds['latitude'] = xr.DataArray(lats, dims=('latitude'))
+
+    for ipoint in range(3):
+        x_xr = xr.DataArray([lons[ipoint]],dims=('plipoints'))
+        y_xr = xr.DataArray([lats[ipoint]],dims=('plipoints'))
+
+        # interp1d # these are actually irrelevant now, are not used for testing
+        interp_with_floats = ds.interp(longitude=x_xr[0], latitude=y_xr[0], method='linear').so #selecting one value from the da drops the new plipoints dimension
+        interp_with_da_existing = ds.interp(longitude=x_xr.values, latitude=y_xr.values, method='linear').so.isel(longitude=0,latitude=0) #using the DataArray values keeps lat/lon dimenions, gives the same interp result
+
+        # combination of linear and nearest interpn
+        interp_with_da_newdim_lin = ds.interp(longitude=x_xr, latitude=y_xr, method='linear').so.isel(plipoints=0) #using the DataArray introduces a plipoints dimension, which gives different interp result
+        interp_with_da_newdim_near = ds.interp(longitude=x_xr, latitude=y_xr, method='nearest').so.isel(plipoints=0) #using the DataArray introduces a plipoints dimension, which gives different interp result
+        interp_with_da_newdim = interp_with_da_newdim_lin.combine_first(interp_with_da_newdim_near)
+        interp_da_expected = xr.DataArray(so_np[:,ipoint,ipoint],dims=('depth'))
+
+        print(ipoint)
+        print(interp_with_floats.to_numpy())
+        print(interp_with_da_existing.to_numpy())
+        print(interp_with_da_newdim.to_numpy())
+        print(interp_da_expected.to_numpy())
+
+        assert (interp_with_floats.isnull()==interp_with_da_existing.isnull()).all()
+        assert (interp_da_expected.isnull()==interp_with_da_newdim.isnull()).all()
+
+    """
+    prints with scipy 1.11.3:
+    0
+    [35.819576 35.802788 35.786003       nan       nan]
+    [35.819576 35.802788 35.786003       nan       nan]
+    [35.819576 35.802788 35.786003 35.776848 35.781425]
+    [35.819576 35.802788 35.786003 35.776848 35.781425]
+    1
+    [35.824154 35.810417 35.796684       nan       nan]
+    [35.824154 35.810417 35.796684       nan       nan]
+    [35.824154 35.810417 35.796684       nan       nan]
+    [35.824154 35.810417 35.796684       nan       nan]
+    2
+    [35.82873 35.81805      nan      nan      nan]
+    [35.82873 35.81805      nan      nan      nan]
+    [35.82873 35.81805      nan      nan      nan]
+    [35.82873 35.81805      nan      nan      nan]
+    """
+
+
 @pytest.mark.systemtest
 @pytest.mark.requireslocaldata
 def test_interpolate_tide_to_plipoints():