diff --git a/src/pygeogrids/grids.py b/src/pygeogrids/grids.py
index 3a04fa6..81f8e37 100644
--- a/src/pygeogrids/grids.py
+++ b/src/pygeogrids/grids.py
@@ -391,11 +391,13 @@ def find_nearest_gpi(self, lon, lat, max_dist=np.Inf):
         Returns
         -------
         gpi : long
-            Grid point index.
+            Grid point index. If no point was found within the maximum
+            distance to consider, an empty array is returned.
         distance : float
             Distance of gpi to given lon, lat.
             At the moment not on a great circle but in spherical
-            cartesian coordinates.
+            cartesian coordinates. If no point was found within the maximum
+            distance to consider, an empty array is returned.
         """
         gpi, distance = self.find_k_nearest_gpi(lon, lat, max_dist=max_dist, k=1)
 
@@ -435,7 +437,7 @@ def find_k_nearest_gpi(self, lon, lat, max_dist=np.Inf, k=1):
         distance, ind = self.kdTree.find_nearest_index(
             lon, lat, max_dist=max_dist, k=k)
 
-        if self.gpidirect and self.allpoints:
+        if self.gpidirect and self.allpoints or len(ind) == 0:
             gpi = ind
         else:
             gpi = self.activegpis[ind]
diff --git a/tests/test_grid.py b/tests/test_grid.py
index a5c84d1..8becac3 100644
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -151,6 +151,13 @@ def test_nearest_neighbor_max_dist(self):
         assert len(gpi) == 0
         assert len(dist) == 0
 
+        # test with custom gpi, see issue #68
+        grid = grids.BasicGrid(lon=[16,17], lat=[45,46], gpis=[100,200])
+        gpi, dist = grid.find_nearest_gpi(0,0, max_dist=1000)
+        assert len(gpi) == 0
+        assert len(dist) == 0
+                                                
+
 class TestCellGridNotGpiDirect(unittest.TestCase):
 
     """