tests for max_dist in find_nearest_gpi

The tests revealed that d=inf and gpi=len(gpis)+1 is returned in case no gpi is
found within max_dist.
I added a check for infinite distances in nearest_neighbor.find_nearest_index.

src/pygeogrids/grids.py

@@ -399,7 +399,7 @@ def find_nearest_gpi(self, lon, lat, max_dist=np.Inf):
         """
         gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=max_dist, k=1)
 
-        if not _element_iterable(lon):
+        if not _element_iterable(lon) and len(gpi) > 0:
             gpi = gpi[0]
             distance = distance[0]
 

src/pygeogrids/nearest_neighbor.py

@@ -178,12 +178,22 @@ def find_nearest_index(self, lon, lat, max_dist=np.Inf, k=1):
             great circle distance at the moment. This should be OK for
             most applications that look for the nearest neighbor which
             should not be hundreds of kilometers away.
+            If no point was found within the maximum distance to consider, an
+            empty array is returned.
         ind : int, numpy.array
-            indices of nearest neighbor
+            If ``self.grid`` is ``False`` indices of nearest neighbor.
+            If no point was found within the maximum distance to consider, an
+            empty array is returned.
         index_lon : numpy.array, optional
-            if self.grid is True then return index into lon array of grid definition
+            If ``self.grid`` is ``True`` then return index into lon array of
+            grid definition.
+            If no point was found within the maximum distance to consider, an
+            empty array is returned.
         index_lat : numpy.array, optional
-            if self.grid is True then return index into lat array of grid definition
+            If ``self.grid`` is ``True`` then return index into lat array of
+            grid definition.
+            If no point was found within the maximum distance to consider, an
+            empty array is returned.
         """
         if self.kdtree is None:
             self._build_kdtree()
@@ -193,11 +203,15 @@ def find_nearest_index(self, lon, lat, max_dist=np.Inf, k=1):
         d, ind = self.kdtree.query(
             query_coords, distance_upper_bound=max_dist, k=k)
 
+        # if no point was found, d == inf
+        if not np.all(np.isfinite(d)):
+            d, ind = np.array([]), np.array([])
+
         if not self.grid:
             return d, ind
         else:
-            # calculate index position in grid definition arrays assuming row-major
-            # flattening of arrays after numpy.meshgrid
+            # calculate index position in grid definition arrays assuming
+            # row-major flattening of arrays after numpy.meshgrid
             index_lat = ind / self.lon_size
             index_lon = ind % self.lon_size
             return d, index_lon.astype(np.int32), index_lat.astype(np.int32)

tests/test_grid.py

@@ -137,6 +137,19 @@ def test_k_nearest_neighbor_list(self):
         assert gpi[1, 0] == 38430
         assert gpi[1, 1] == 38429
 
+    def test_nearest_neighbor_max_dist(self):
+        # test with maxdist higher than nearest point
+        gpi, dist = self.grid.find_nearest_gpi(14.3, 18.5, max_dist=100e3)
+        assert gpi == 25754
+        assert len([dist]) == 1
+        lon, lat = self.grid.gpi2lonlat(gpi)
+        assert lon == 14.5
+        assert lat == 18.5
+
+        # test with maxdist lower than nearest point
+        gpi, dist = self.grid.find_nearest_gpi(14.3, 18.5, max_dist=10e3)
+        assert len(gpi) == 0
+        assert len(dist) == 0
 
 class TestCellGridNotGpiDirect(unittest.TestCase):