diff --git a/CHANGES.rst b/CHANGES.rst
index 1b99131..1bcad83 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -5,10 +5,13 @@ Changelog
 v0.2.4
 ======
 
-- Add option to subset a grid with a shape file (OGRGeometry) in get_shp_grid_points.
+- Add option to subset a grid with a shape file (OGRGeometry) in
+  get_shp_grid_points.
 - Add shapefile module for reading shapefiles from
   http://biogeo.ucdavis.edu/data/gadm2.8/gadm28_levels.shp.zip by Global
   Administrative Level
+- Ensure that `get_bbox_grid_points` returns points while taking cell order into
+  account.
 
 v0.2.3
 ======
diff --git a/pygeogrids/grids.py b/pygeogrids/grids.py
index ed9d442..e5a4e66 100755
--- a/pygeogrids/grids.py
+++ b/pygeogrids/grids.py
@@ -199,7 +199,7 @@ def __init__(self, lon, lat, gpis=None, geodatum='WGS84', subset=None,
         self.geodatum = GeodeticDatum(geodatum)
 
         if gpis is None:
-            self.gpis = np.arange(self.n_gpi)
+            self.gpis = np.arange(self.n_gpi, dtype=int)
             self.gpidirect = True
         else:
             if lat.shape != gpis.shape:
@@ -607,7 +607,7 @@ def get_bbox_grid_points(self, latmin=-90, latmax=90, lonmin=-180,
         """
 
         gp_info = np.array(list(self.grid_points()))
-        gpis = gp_info[:, 0]
+        gpis = gp_info[:, 0].astype(int)
         lons = gp_info[:, 1]
         lats = gp_info[:, 2]
         index = np.where((lats <= latmax) &
diff --git a/tests/test_grid.py b/tests/test_grid.py
index e62e0d4..c08c507 100755
--- a/tests/test_grid.py
+++ b/tests/test_grid.py
@@ -156,6 +156,14 @@ def test_gpi2cell_iterable(self):
         cell = self.cellgrid.gpi2cell(gpi)
         assert np.all(cell == [1043, 2015])
 
+    def test_gpi2cell_numpy(self):
+        """
+        test if gpi to cell lookup works correctly
+        """
+        gpi = np.array([200, 255])
+        cell = self.cellgrid.gpi2cell(gpi)
+        assert np.all(cell == [1043, 2015])
+
     def test_gpi2cell_numpy_single(self):
         """
         test if gpi to row column lookup works correctly
@@ -412,6 +420,52 @@ def test_subgrid_from_gpis(self):
             lons_should, lats_should, cells_should, gpis=gpis)
         assert subgrid == subgrid_should
 
+    def test_get_bbox_grid_points(self):
+        gpis = self.cellgrid.get_bbox_grid_points(latmin=-10,
+                                                  latmax=-5,
+                                                  lonmin=-10,
+                                                  lonmax=-5)
+        nptest.assert_allclose(gpis,
+                               np.array([5684, 5685, 5828, 5829,
+                                         5540, 5541, 5686, 5830, 5542]))
+        # gpis should come back sorted by cells
+        nptest.assert_allclose(self.cellgrid.gpi2cell(gpis),
+                               np.array([1240, 1240, 1240, 1240,
+                                         1241, 1241, 1276, 1276, 1277]))
+        lats, lons = self.cellgrid.get_bbox_grid_points(latmin=-10,
+                                                        latmax=-5,
+                                                        lonmin=-10,
+                                                        lonmax=-5,
+                                                        coords=True)
+        lats_should = np.array([-7.5, -7.5, -10., -10.,
+                                -5., -5., -7.5, -10., -5.])
+        lons_should = np.array([-10.,  -7.5, -10.,  -7.5,
+                                -10.,  -7.5,  -5.,  -5.,  -5.])
+        nptest.assert_allclose(lats,
+                               lats_should)
+        nptest.assert_allclose(lons,
+                               lons_should)
+        gpis, lats, lons = self.cellgrid.get_bbox_grid_points(latmin=-10,
+                                                              latmax=-5,
+                                                              lonmin=-10,
+                                                              lonmax=-5,
+                                                              both=True)
+        lats_should = np.array([-7.5, -7.5, -10., -10.,
+                                -5., -5., -7.5, -10., -5.])
+        lons_should = np.array([-10.,  -7.5, -10.,  -7.5,
+                                -10.,  -7.5,  -5.,  -5.,  -5.])
+        nptest.assert_allclose(lats,
+                               lats_should)
+        nptest.assert_allclose(lons,
+                               lons_should)
+        nptest.assert_allclose(gpis,
+                               np.array([5684, 5685, 5828, 5829,
+                                         5540, 5541, 5686, 5830, 5542]))
+        # gpis should come back sorted by cells
+        nptest.assert_allclose(self.cellgrid.gpi2cell(gpis),
+                               np.array([1240, 1240, 1240, 1240,
+                                         1241, 1241, 1276, 1276, 1277]))
+
 
 def test_setup_grid_with_lists():
 
@@ -586,5 +640,6 @@ def test_shpgrid(self):
         assert subgrid.activearrlon == 15
         assert subgrid.activearrlat == 46
 
+
 if __name__ == "__main__":
     unittest.main()