Add nearest interp method

Author: AntSimi

CHANGELOG.md

@@ -11,6 +11,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Acces at dataset variable like attribute, and lifetime/age are available for all observations
 - Add **EddyInfos** application to get general information about eddies dataset
 - Add method to inspect contour rejection (which are not in eddies)
+- Grid interp could be "nearest" or "bilinear" 
 ### Changed
 
 ### Removed

src/py_eddy_tracker/dataset/grid.py

@@ -1842,13 +1842,14 @@ def regrid(self, other, grid_name, new_name=None):
         # self.variables_description[new_name]['infos'] = False
         # self.variables_description[new_name]['kwargs']['dimensions'] = ...
 
-    def interp(self, grid_name, lons, lats):
+    def interp(self, grid_name, lons, lats, method="bilinear"):
         """
         Compute z over lons, lats
 
         :param str grid_name: Grid to be interpolated
         :param lons: new x
         :param lats: new y
+        :param str method: Could be 'bilinear' or 'nearest'
 
         :return: new z
         """
@@ -1857,7 +1858,9 @@ def interp(self, grid_name, lons, lats):
             m = g.mask
         else:
             m = ones(g.shape) if g.mask else zeros(g.shape)
-        return interp2d_geo(self.x_c, self.y_c, g, m, lons, lats)
+        return interp2d_geo(
+            self.x_c, self.y_c, g, m, lons, lats, nearest=method == "nearest"
+        )
 
 
 @njit(cache=True, fastmath=True)

src/py_eddy_tracker/generic.py

@@ -177,7 +177,7 @@ def cumsum_by_track(field, track):
 
 
 @njit(cache=True, fastmath=True)
-def interp2d_geo(x_g, y_g, z_g, m_g, x, y):
+def interp2d_geo(x_g, y_g, z_g, m_g, x, y, nearest=False):
     """
     For geographic grid, test of cicularity.
 
@@ -187,6 +187,7 @@ def interp2d_geo(x_g, y_g, z_g, m_g, x, y):
     :param array m_g: Boolean grid, True if value is masked
     :param array x: coordinate where interpolate z
     :param array y: coordinate where interpolate z
+    :param bool nearest: if true we will take nearest pixel
     :return: z interpolated
     :rtype: array
     """
@@ -223,9 +224,21 @@ def interp2d_geo(x_g, y_g, z_g, m_g, x, y):
         if m_g[i0, j0] or m_g[i0, j1] or m_g[i1, j0] or m_g[i1, j1]:
             z[i] = nan
         else:
-            z[i] = (z00 * (1 - xd) + (z10 * xd)) * (1 - yd) + (
-                z01 * (1 - xd) + z11 * xd
-            ) * yd
+            if nearest:
+                if xd <= 0.5:
+                    if yd <= 0.5:
+                        z[i] = z00
+                    else:
+                        z[i] = z01
+                else:
+                    if yd <= 0.5:
+                        z[i] = z10
+                    else:
+                        z[i] = z11
+            else:
+                z[i] = (z00 * (1 - xd) + (z10 * xd)) * (1 - yd) + (
+                    z01 * (1 - xd) + z11 * xd
+                ) * yd
     return z
 
 

src/py_eddy_tracker/observations/observation.py

@@ -1934,12 +1934,12 @@ def interp_grid(
         :param grid_object: Handler of grid to interp
         :type grid_object: py_eddy_tracker.dataset.grid.RegularGridDataset
         :param str varname: Name of variable to use
-        :param str method: 'center', 'mean', 'max', 'min'
+        :param str method: 'center', 'mean', 'max', 'min', 'nearest'
         :param str dtype: if None we use var dtype
         :param bool intern: Use extern or intern contour
         """
-        if method == "center":
-            return grid_object.interp(varname, self.longitude, self.latitude)
+        if method in ("center", "nearest"):
+            return grid_object.interp(varname, self.longitude, self.latitude, method)
         elif method in ("min", "max", "mean", "count"):
             x0 = grid_object.x_bounds[0]
             x_name, y_name = self.intern(False if intern is None else intern)

tests/test_grid.py

@@ -1,4 +1,5 @@
 from matplotlib.path import Path
+from numpy import array, ma
 from pytest import approx
 
 from py_eddy_tracker.data import get_path
@@ -53,3 +54,24 @@ def test_bounds():
     x0, x1, y0, y1 = G.bounds
     assert x0 == -1 / 120.0 and x1 == 360 - 1 / 120
     assert y0 == approx(-90 - 1 / 120.0) and y1 == approx(90 - 1 / 120)
+
+
+def test_interp():
+    # Fake grid
+    g = RegularGridDataset.with_array(
+        coordinates=("x", "y"),
+        datas=dict(
+            z=ma.array(((0, 1), (2, 3)), dtype="f4"),
+            x=array((0, 20)),
+            y=array((0, 10)),
+        ),
+        centered=True,
+    )
+    x0, y0 = array((10,)), array((5,))
+    x1, y1 = array((15,)), array((5,))
+    # Interp nearest
+    assert g.interp("z", x0, y0, method="nearest") == 0
+    assert g.interp("z", x1, y1, method="nearest") == 2
+    # Interp bilinear
+    assert g.interp("z", x0, y0) == 1.5
+    assert g.interp("z", x1, y1) == 2