From 4679d8fbe6844ce3ec8e6e94af2af4753922de23 Mon Sep 17 00:00:00 2001 From: Samuel Scherrer Date: Fri, 27 Nov 2020 10:05:22 +0100 Subject: [PATCH] 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 | 2 +- src/pygeogrids/nearest_neighbor.py | 24 +++++++++++++++++++----- tests/test_grid.py | 13 +++++++++++++ 3 files changed, 33 insertions(+), 6 deletions(-) diff --git a/src/pygeogrids/grids.py b/src/pygeogrids/grids.py index 24eabb7..3a04fa6 100644 --- a/src/pygeogrids/grids.py +++ b/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] diff --git a/src/pygeogrids/nearest_neighbor.py b/src/pygeogrids/nearest_neighbor.py index e088e95..45e4c3d 100644 --- a/src/pygeogrids/nearest_neighbor.py +++ b/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) diff --git a/tests/test_grid.py b/tests/test_grid.py index 08ad010..a5c84d1 100644 --- a/tests/test_grid.py +++ b/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):