Repository: DHI/mikeio
Pull request #553 (2023-03-12 07:58:45)

Grid 2d find index performance


Patch:
diff --git a/Makefile b/Makefile
index 87f5c4b9..52c1f483 100644
--- a/Makefile
+++ b/Makefile
@@ -10,6 +10,9 @@ doctest:
 	pytest mikeio/dfs*.py mikeio/dfsu/*.py mikeio/eum.py mikeio/pfs/*.py mikeio/spatial/grid_geometry.py --doctest-modules
 	rm -f *.dfs* # remove temporary files, created from doctests
 
+perftest:
+	pytest tests/performance/ --durations=0
+
 typecheck:
 	mypy mikeio/dataarray.py mikeio/dataset.py
 
diff --git a/mikeio/spatial/grid_geometry.py b/mikeio/spatial/grid_geometry.py
index 7ebe9bee..67bcc4e0 100644
--- a/mikeio/spatial/grid_geometry.py
+++ b/mikeio/spatial/grid_geometry.py
@@ -712,10 +712,12 @@ class Grid2D(_Geometry):
     ):
         """Find nearest index (i,j) of point(s)
 
-        -1 is returned if point is outside grid
-
         Parameters
         ----------
+        x : float
+            x-coordinate of point
+        y : float
+            y-coordinate of point
         coords : array(float)
             xy-coordinate of points given as n-by-2 array
         area : array(float)
@@ -725,12 +727,21 @@ class Grid2D(_Geometry):
         -------
         array(int), array(int)
             i- and j-index of nearest cell
+
+        Raises
+        ------
+        ValueError if x or y are not scalar values, use coords instead
+        OutsideModelDomainError if point is outside grid
         """
         if x is None and y is not None and not np.isscalar(y):
-            raise ValueError(f"{y} is not a scalar value")
+            raise ValueError(
+                f"{y=} is not a scalar value, use the coords argument instead"
+            )
 
         if y is None and x is not None and not np.isscalar(x):
-            raise ValueError(f"{x} is not a a scalar value")
+            raise ValueError(
+                f"{x=} is not a scalar value, use the coords argument instead"
+            )
 
         if x is not None and y is not None:
             if coords is not None:
@@ -756,12 +767,12 @@ class Grid2D(_Geometry):
         jj = (-999999999) * np.ones_like(y, dtype=int)
 
         inside = self.contains(xy)
-        for j, xyp in enumerate(xy):
-            if inside[j]:
-                ii[j] = (np.abs(self.x - xyp[0])).argmin()
-                jj[j] = (np.abs(self.y - xyp[1])).argmin()
-            else:
-                raise OutsideModelDomainError(x=x[j], y=y[j])
+        if np.any(~inside):
+            raise OutsideModelDomainError(x=x[~inside], y=y[~inside])
+
+        # get index in x, and y for points inside based on the grid spacing and origin
+        ii = np.floor((x - (self.x[0] - self.dx / 2)) / self.dx).astype(int)
+        jj = np.floor((y - (self.y[0] - self.dy / 2)) / self.dy).astype(int)
 
         return ii, jj
 


Test patch:
diff --git a/.github/workflows/perf_test.yml b/.github/workflows/perf_test.yml
index 57f909cf..1b39e73c 100644
--- a/.github/workflows/perf_test.yml
+++ b/.github/workflows/perf_test.yml
@@ -32,4 +32,4 @@ jobs:
         pip install .[test]
     - name: Test with pytest
       run: |
-        pytest tests/performance/ --durations=0
+        make perftest
diff --git a/tests/performance/test_grid2d.py b/tests/performance/test_grid2d.py
new file mode 100644
index 00000000..10765856
--- /dev/null
+++ b/tests/performance/test_grid2d.py
@@ -0,0 +1,21 @@
+import numpy as np
+import mikeio
+
+
+def test_find_index_inside_domain():
+
+    g = mikeio.Grid2D(bbox=[-1, -1, 1, 5], dx=0.2)
+
+    np.random.seed(0)
+
+    # all points are inside the grid
+    x = np.random.uniform(-1, 1, size=100_000)
+    y = np.random.uniform(-1, 5, size=100_000)
+
+    xy = np.column_stack([x, y])
+    assert xy.shape == (100_000, 2)
+
+    ii, jj = g.find_index(coords=xy)
+
+    assert (ii[0], jj[0]) == g.find_index(x=x[0], y=y[0])
+    assert (ii[-1], jj[-1]) == g.find_index(x=x[-1], y=y[-1])
diff --git a/tests/test_geometry_grid.py b/tests/test_geometry_grid.py
index 358fe8d6..56420424 100644
--- a/tests/test_geometry_grid.py
+++ b/tests/test_geometry_grid.py
@@ -9,8 +9,6 @@ from mikeio.spatial.geometry import GeometryUndefined
 from mikeio.exceptions import OutsideModelDomainError
 
 
-
-
 def test_create_nx_ny():
 
     g = Grid2D(x0=180.0, y0=-90.0, dx=0.25, dy=0.25, nx=1440, ny=721)
@@ -280,6 +278,45 @@ def test_find_index():
         g.find_index(coords=[(-0.1, 0.1), (-0.1, 0.1)])
 
 
+def test_find_index_grid_with_negative_origin():
+    # create a grid with negative origin
+    bbox = [-1, -1, 1, 5]
+    g = Grid2D(bbox=bbox, dx=0.5)
+
+    xy1 = [0.52, 1.52]
+    ii, jj = g.find_index(x=xy1[0], y=xy1[1])
+    assert ii == 3
+    assert jj == 5
+
+    # to the left of the first element center, and above the last element center
+    xy2 = [-0.9, 4.95]
+    ii, jj = g.find_index(x=xy2[0], y=xy2[1])
+    assert ii == 0
+    assert jj == 11
+
+    # test three points, two inside and one outside
+    xy3 = [-1.5, 0.5]
+
+    xy = np.vstack([xy1, xy2, xy3])
+    with pytest.raises(OutsideModelDomainError):
+        g.find_index(coords=xy)
+
+
+def test_find_index_y_array_not_possible():
+    # create a simple grid
+    bbox = [0, 0, 1, 5]
+    g = Grid2D(bbox=bbox, dx=0.5)
+
+    xy1 = [0.52, 1.52]
+    xy2 = [0.4, 1.2]
+    xy = np.vstack([xy1, xy2])
+
+    # supply y as array
+    with pytest.raises(ValueError) as excinfo:
+        g.find_index(y=xy[:, 1])
+    assert "y=" in str(excinfo.value)
+
+
 def test_to_geometryFM():
     nx = 5
     ny = 3