diff --git a/idefix/utils.py b/idefix/utils.py
index 6fde6e7..701d0f8 100644
--- a/idefix/utils.py
+++ b/idefix/utils.py
@@ -55,5 +55,42 @@ def bbox(data):
 
         [[xmin, ymin, zmin],
          [xmax, ymax, zmax]]
+
+    See Also
+    --------
+    fit_bbox : Return a bounding box fit on fixed coordinates.
     """
     return np.array((np.min(data, axis=0), np.max(data, axis=0)))
+
+def fit_bbox(data, decimals=0):
+    """Return a bounding box fit on fixed coordinates.
+
+    - Round $x$ and $y$ coordinates to match most orthoimagery tiles.
+    - Ceil and floor $z$ coordinates to include all the point on the vertical axis.
+
+    Parameters
+    ----------
+    data : ndarray (N, 3)
+        Bbox or point cloud data of shape (N, 3), i.e. (x,y,z).
+    decimals : int
+        The precision for the rounding, ceiling and flooring operations.
+
+    Returns
+    -------
+    bbox : ndarray
+        Lower and upper points describing the bounding box such as::
+
+        [[xmin, ymin, zmin],
+         [xmax, ymax, zmax]]
+
+    See Also
+    --------
+    bbox : Returns a raw bounding box on the data.
+    """
+    bbox = bbox(data)
+
+    nbbox = np.round(bbox, decimals)
+    nbbox[0,2] = np.floor(bbox * 10 ** decimals)[0,2] / 10 ** decimals
+    nbbox[1,2] = np.ceil(bbox * 10 ** decimals)[1,2] / 10 ** decimals
+
+    return nbbox
diff --git a/idefix/vxl.py b/idefix/vxl.py
index 44f18cc..bfdce3b 100644
--- a/idefix/vxl.py
+++ b/idefix/vxl.py
@@ -11,7 +11,7 @@ General functions to transform point clouds to voxels compatible with numpy.
 import logging
 import numpy as np
 import humanize
-from .utils import bbox
+from .utils import bbox, fit_bbox
 
 log = logging.getLogger(__name__)
 
@@ -38,8 +38,8 @@ def _ui_step(step, spatial):
 def get_grid(spatial, step):
     '''Return grid bins.
 
-    Compute the grid bins of a point cloud or the corresponding bounding box
-    according to given step (or steps for anisotropic grid).
+    Compute the grid bins of a point cloud or the corresponding bounding
+    box according to given step (or steps for anisotropic grid).
 
     Parameters
     ----------
@@ -54,8 +54,26 @@ def get_grid(spatial, step):
     Returns
     -------
     grid : array of array (n,)
-        Grid of spatial given step. Return three arrays (not necessarily of the
-        same size) defining the bins of axis `x`, `y` and `z`.
+        Grid of spatial given step. Return three arrays (not necessarily
+        of the same size) defining the bins of axis `x`, `y` and `z`.
+
+    Notes
+    -----
+    The grid is built considering the half-open interval
+    $[min_of_the_axis, max_of_the_axis)$. If the positions of the points
+    are directly used as spatial parameter, the points at the upper
+    limits will be excluded from further processing.
+
+    You can define a more precise bounding box to take into account all
+    the points of your dataset (e.g. by adding a small distance on the
+    upper limits).
+
+    See Also
+    --------
+    bbox : Returns the raw bounding box of the point cloud (excluding
+        points on upper limit).
+    fit_bbox : Returns a bounding box on rounded coordinates (can
+        include all the points).
     '''
     spatial = np.array(spatial)
     bb = bbox(spatial)
@@ -65,7 +83,7 @@ def get_grid(spatial, step):
     for a_min, a_max, a_s in zip(bb[0], bb[1], step):
         # Beware of float underflow
         if a_s:
-            bins = np.trunc((a_max - a_min) / a_s).astype(int) + 1
+            bins = np.trunc((a_max - a_min) / a_s).astype(int)
             grid += [np.linspace(a_min, a_min + bins * a_s, bins + 1)]
         else:
             grid += [np.array((a_min, a_max + 1))]
@@ -347,7 +365,7 @@ def squash(voxel_grid, method='top', axis=-1):
         return voxel_grid.max(axis=axis)
     elif method == 'std':
         return voxel_grid.std(axis=axis)
-    
+
     raise NotImplementedError('Method \'{}\' does not exist.'.format(method))
 
 def plot(voxel_grid, vmin=None, vmax=None):
diff --git a/test/test_vxl.py b/test/test_vxl.py
index 99c01cf..43f0e30 100644
--- a/test/test_vxl.py
+++ b/test/test_vxl.py
@@ -86,7 +86,14 @@ def test_get_grid(datadir, set_id, step, grid_id):
     assert res is not None, 'Test data empty, test function is broken!'
     assert len(res) == 3, 'Test data malformed, test function is broken!'
 
-    test = vxl.get_grid(spatial, step)
+    bbox = vxl.bbox(spatial)
+    try:
+        bbox[1,:] += step * 1.01
+    except:
+        # Workaround for the None on last axis
+        bbox[1,:] += step[0]
+
+    test = vxl.get_grid(bbox, step)
 
     assert test is not None, 'Function did not return anything :('
     assert len(test) == 3, 'Function doesn\'t give right number of axis'
@@ -178,7 +185,7 @@ def test__geo_to_np_coordinate():
     raster_truth[0, -1] = 25
     raster_truth[-1, 2] = 7
 
-    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 2D raters' 
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 2D raters'
 
     raster = np.zeros((5, 5, 3), dtype=np.uint8)
     raster[0, 0, 0] = 42
@@ -190,7 +197,7 @@ def test__geo_to_np_coordinate():
     raster_truth[0, -1, 1] = 25
     raster_truth[-1, 2, 2] = 7
 
-    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters' 
+    assert (raster_truth == vxl._geo_to_np_coordinate(raster)).all(), 'Missmatch between 3D raters'
 
 @pytest.mark.parametrize('set_id, grid_id, axis, method', [
     (1, 1, 2, 'top'),
diff --git a/test/test_vxl/pc0_grid_s1-1-n.txt b/test/test_vxl/pc0_grid_s1-1-n.txt
index 1382f54..4a2c07a 100644
--- a/test/test_vxl/pc0_grid_s1-1-n.txt
+++ b/test/test_vxl/pc0_grid_s1-1-n.txt
@@ -1,3 +1,3 @@
 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0
-1.0 11.0
+1.0 12.0