diff --git a/scenes/__init__.py b/scenes/__init__.py
index f6d872852e951d5c48c1697849d0d835a2b5cdac..7daf55f709d093503d01829558a58d71d23bb32e 100644
--- a/scenes/__init__.py
+++ b/scenes/__init__.py
@@ -7,3 +7,4 @@ __version__ = "1.0.1"
 from .core import load_scenes, save_scenes  # noqa: 401
 from .indexation import Index  # noqa: 401
 from .download import TheiaDownloader  # noqa: 401
+from .bbox import BoundingBox  # noqa: 401
diff --git a/scenes/bbox.py b/scenes/bbox.py
new file mode 100644
index 0000000000000000000000000000000000000000..16acde3114ddd0f08ee407a40eee7166b2f1dd38
--- /dev/null
+++ b/scenes/bbox.py
@@ -0,0 +1,37 @@
+"""
+This module provides bounding box class
+"""
+
+
+class BoundingBox:
+    """
+    The bounding box class
+    """
+    def __init__(self, xmin, xmax, ymin, ymax):
+        """
+        Args:
+            xmin: Lower value on the x-axis
+            xmax: Higher value on the x-axis
+            ymin: Lower value on the y-axis
+            ymax: Higher value on the y-axis
+        """
+        self.xmin = xmin
+        self.xmax = xmax
+        self.ymin = ymin
+        self.ymax = ymax
+
+    def union(self, other):
+        """
+        Return a new bounding box resulting in the union of self and other
+
+        Args:
+            other: another bounding box
+
+        Returns:
+            a new bounding box
+
+        """
+        return BoundingBox(xmin=min(self.xmin, other.xmin),
+                           xmax=max(self.xmax, other.xmax),
+                           ymin=min(self.ymin, other.ymin),
+                           ymax=max(self.ymax, other.ymax))
diff --git a/scenes/download.py b/scenes/download.py
index eb55009a5f1ed502bea55cf77715cac157fbb81f..cf950ffaf8b91fc1dd1faef4bd6e275e70abd1a3 100644
--- a/scenes/download.py
+++ b/scenes/download.py
@@ -12,6 +12,15 @@ import pycurl
 import tqdm
 
 
+def bbox2str(bbox):
+    """
+    Return a str containing the bounding box
+    :param bbox: the bounding box (BoundingBox instance)
+    :return: a str
+    """
+    return '{},{},{},{}'.format(bbox.ymin, bbox.xmin, bbox.ymax, bbox.xmax)
+
+
 def compute_md5(fname):
     """Compute md5sum of a file
 
@@ -244,10 +253,8 @@ class TheiaDownloader:
             # TODO
         """
         start_date, end_date = dates_range
-        # lonmin, latmin, lonmax, latmax
-        box = '{},{},{},{}'.format(bbox_wgs84[2], bbox_wgs84[0], bbox_wgs84[3], bbox_wgs84[1])
         dict_query = {
-            "box": box,
+            "box": bbox2str(bbox_wgs84),  # lonmin, latmin, lonmax, latmax
             "startDate": start_date.strftime("%Y-%m-%d"),
             "completionDate": end_date.strftime("%Y-%m-%d"),
             "maxRecords": self.max_records,
@@ -281,9 +288,7 @@ class TheiaDownloader:
         ndays_seek = datetime.timedelta(days=17)  # temporal range to check for monthly synthesis
 
         # Query products
-        # box is [lonmin, latmin, lonmax, latmax]
-        box = '{},{},{},{}'.format(bbox_wgs84[2], bbox_wgs84[0], bbox_wgs84[3], bbox_wgs84[1])
-        dict_query = {'box': box}
+        dict_query = {'box': bbox2str(bbox_wgs84)}  # lonmin, latmin, lonmax, latmax
         start_date = acq_date - ndays_seek
         end_date = acq_date + ndays_seek
 
diff --git a/scenes/indexation.py b/scenes/indexation.py
index 6843a177cfa63a727bd45dd199fe61c483ce8f87..c10cd0e42b87ac675c1d27cc18b5b9c98769a7ff 100644
--- a/scenes/indexation.py
+++ b/scenes/indexation.py
@@ -39,7 +39,7 @@ def bbox(bbox_wgs84, dt_min, dt_max):
     """Return a bounding box in the domain (lat, lon, time)
 
     Args:
-        bbox_wgs84: Bounding box (in WGS84)
+        bbox_wgs84: The bounding box in WGS84 (BoundingBox instance)
         dt_min: date min (datetime.datetime)
         dt_max: date max (datetime.datetime)
 
@@ -47,8 +47,8 @@ def bbox(bbox_wgs84, dt_min, dt_max):
         item for rtree
 
     """
-    (xmin, xmax, ymin, ymax) = bbox_wgs84
-    return xmin, ymin, get_timestamp(dt_min), xmax, ymax, get_timestamp(dt_max)
+    return (bbox_wgs84.xmin, bbox_wgs84.ymin, get_timestamp(dt_min),
+            bbox_wgs84.xmax, bbox_wgs84.ymax, get_timestamp(dt_max))
 
 
 class Index:
@@ -72,7 +72,7 @@ class Index:
         """Search the intersecting elements, and return their indices
 
         Args:
-            bbox_wgs84: bounding box (WGS84)
+            bbox_wgs84: The bounding box in WGS84 (BoundingBox instance)
             dt_min: date min (datetime.datetime) (Default value = None)
             dt_max: date max (datetime.datetime) (Default value = None)
 
@@ -91,7 +91,7 @@ class Index:
         """Search the intersecting elements, and return them
 
         Args:
-            bbox_wgs84: bounding box (WGS84)
+            bbox_wgs84: The bounding box in WGS84 (BoundingBox instance)
             dt_min: date min (datetime.datetime) (Default value = None)
             dt_max: date max (datetime.datetime) (Default value = None)
 
diff --git a/scenes/utils.py b/scenes/utils.py
index d15d88b52c75a16428c141726890ca08c64f0b06..0f7a544c643639799d93f76f5d702feab0b24c02 100644
--- a/scenes/utils.py
+++ b/scenes/utils.py
@@ -9,6 +9,7 @@ import zipfile
 import re
 import fnmatch
 from osgeo import osr, ogr, gdal
+from scenes.bbox import BoundingBox
 
 
 def epsg2srs(epsg):
@@ -106,7 +107,7 @@ def extent_to_bbox(extent):
         extent: extent
 
     Returns:
-        bounding box (xmin, xmax, ymin, ymax)
+        BoundingBox instance
 
     """
     xmin, xmax = math.inf, -math.inf
@@ -117,7 +118,7 @@ def extent_to_bbox(extent):
         ymin = min(ymin, y)
         ymax = max(ymax, y)
 
-    return xmin, xmax, ymin, ymax
+    return BoundingBox(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
 
 
 def get_bbox_wgs84_from_gdal_ds(gdal_ds):
@@ -127,7 +128,7 @@ def get_bbox_wgs84_from_gdal_ds(gdal_ds):
         gdal_ds: GDAL dataset
 
     Returns:
-        bounding box (xmin, xmax, ymin, ymax)
+        bounding box in WGS84 CRS (BoundingBox instance)
 
     """
     extend_wgs84 = get_extent_wgs84(gdal_ds)
@@ -249,7 +250,7 @@ def get_bbox_wgs84_from_vector(vector_file):
         vector_file: vector data filename
 
     Returns:
-        bounding box in WGS84 CRS
+        bounding box in WGS84 CRS (BoundingBox instance)
 
     """
     poly_ds = open_vector_layer(vector_file=vector_file)
@@ -261,7 +262,18 @@ def get_bbox_wgs84_from_vector(vector_file):
     [(xmin, ymin), (xmax, ymax)] = reproject_coords(coords=coords,  # pylint: disable=unbalanced-tuple-unpacking
                                                     src_srs=src_srs,
                                                     tgt_srs=tgt_srs)
-    return xmin, xmax, ymin, ymax
+    return BoundingBox(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
+
+
+def get_bbox_wgs84_from_raster(raster_file):
+    """
+    Returns the bounding box in WGS84 CRS from a raster
+    :param raster_file: raster file
+    :return: bounding box in WGS84 CRS (BoundingBox instance)
+    """
+    gdal_ds = gdal.Open(raster_file)
+
+    return get_bbox_wgs84_from_gdal_ds(gdal_ds=gdal_ds)
 
 
 def find_files_in_all_subdirs(pth, pattern, case_sensitive=True):
diff --git a/test/download_test.py b/test/download_test.py
index ab5f9b11402127472aa1c8989e65b79622b183ec..5c84a86d2e20572887c0eb4ab975632a7bece915 100644
--- a/test/download_test.py
+++ b/test/download_test.py
@@ -13,7 +13,7 @@ parser.add_argument("--day", type=int, default=1)
 params = parser.parse_args()
 
 # Get all scenes in the root_dir
-_, _, bbox = utils.get_epsg_extent_bbox(params.refimage)
+bbox = utils.get_bbox_wgs84_from_raster(params.refimage)
 acq_date = datetime.datetime(year=params.year, month=params.month, day=params.day)
 downloader = TheiaDownloader(config_file=params.theia_cfg)
 downloader.download_closest(bbox_wgs84=bbox, acq_date=acq_date, download_dir=params.download_dir)
diff --git a/test/drs_search.py b/test/drs_search.py
index 94febe5b0c89b53efb3284d08b93e5e4e4c554fa..eb11015ef2bdab4c4e0b6d807a5afd44b1fecaff 100644
--- a/test/drs_search.py
+++ b/test/drs_search.py
@@ -1,5 +1,4 @@
 import argparse
-import gdal
 from scenes import utils, Index
 
 # Arguments
@@ -17,12 +16,9 @@ idx = Index(scenes)
 
 # search
 print("search roi")
-if params.roi.lower().endswith(".tif"):
-    gdal_ds = gdal.Open(params.roi)
-    bbox = utils.get_bbox_wgs84_from_gdal_ds(gdal_ds=gdal_ds)
-else:
-    bbox = utils.get_bbox_wgs84_from_vector(params.roi)
-matches = idx.find(bbox_wgs84=bbox)
+bbox_fn = utils.get_bbox_wgs84_from_raster if params.roi.lower().endswith(".tif") \
+    else utils.get_bbox_wgs84_from_vector
+matches = idx.find(bbox_wgs84=bbox_fn(params.roi))
 print("{} scenes found.".format(len(matches)))
 #for scene_match in matches:
 #    print(scene_match)
diff --git a/test/drs_stack.py b/test/drs_stack.py
index 7901aa1d521447968a9389ec5ebd20e8e6e9f5f7..7deaaba3aa2a0beda3497aaaaaf1fa2f90029b7d 100644
--- a/test/drs_stack.py
+++ b/test/drs_stack.py
@@ -1,5 +1,4 @@
 import argparse
-import gdal
 import pyotb
 from scenes import utils, Index
 from scenes.spot import get_spot67_scenes
@@ -20,8 +19,7 @@ idx = Index(scenes)
 
 # search
 print("search roi")
-gdal_ds = gdal.Open(params.roi)
-bbox = utils.get_extent_wgs84(gdal_ds)
+bbox = utils.get_bbox_wgs84_from_raster(params.roi)
 matches = idx.find(bbox_wgs84=bbox)
 print("{} scenes found.".format(len(matches)))
 
diff --git a/test/indexation_test.py b/test/indexation_test.py
index 5630f4c20ac6114c5ebf146006b01961ec796283..f8733ee9a9e71d08dd4e000d38eb89d6c27a2ec1 100644
--- a/test/indexation_test.py
+++ b/test/indexation_test.py
@@ -1,14 +1,14 @@
 # -*- coding: utf-8 -*-
 from scenes_test_base import ScenesTestBase
-from scenes import Index, utils
+from scenes import Index, utils, BoundingBox
 import tests_data
 
 class ImageryTest(ScenesTestBase):
 
     def test_scene1_indexation(self):
         index = Index(scenes_list=[tests_data.SCENE1])
-        self.assertTrue(index.find(bbox_wgs84=(43.706, 43.708, 4.317, 4.420)))
-        self.assertFalse(index.find(bbox_wgs84=(43.000, 43.001, 3.000, 3.001)))
+        self.assertTrue(index.find(bbox_wgs84=BoundingBox(43.706, 43.708, 4.317, 4.420)))
+        self.assertFalse(index.find(bbox_wgs84=BoundingBox(43.000, 43.001, 3.000, 3.001)))
         self.assertTrue(index.find(bbox_wgs84=utils.get_bbox_wgs84_from_vector(tests_data.ROI_MTP_4326)))
         self.assertTrue(index.find(bbox_wgs84=utils.get_bbox_wgs84_from_vector(tests_data.ROI_MTP_2154)))