diff --git a/csep/core/forecasts.py b/csep/core/forecasts.py
index 9d6599d9..68cf27e7 100644
--- a/csep/core/forecasts.py
+++ b/csep/core/forecasts.py
@@ -397,22 +397,23 @@ def load_ascii(cls, ascii_fname, start_date=None, end_date=None, name=None, swap
         data = numpy.loadtxt(ascii_fname)
         # this is very ugly, but since unique returns a sorted list, we want to get the index, sort that and then return
         # from the original array. same for magnitudes below.
-        all_polys = data[:,:4]
-        all_poly_mask = data[:,-1]
+        all_polys = data[:, :4]
+        all_poly_mask = data[:, -1]
         sorted_idx = numpy.sort(numpy.unique(all_polys, return_index=True, axis=0)[1], kind='stable')
         unique_poly = all_polys[sorted_idx]
         # gives the flag for a spatial cell in the order it was presented in the file
         poly_mask = all_poly_mask[sorted_idx]
         # create magnitudes bins using Mag_0, ignoring Mag_1 bc they are regular until last bin. we dont want binary search for this
-        all_mws = data[:,-4]
+        all_mws = data[:, -4]
         sorted_idx = numpy.sort(numpy.unique(all_mws, return_index=True)[1], kind='stable')
         mws = all_mws[sorted_idx]
         # csep1 stores the lat lons as min values and not (x,y) tuples
-        bboxes = [tuple(itertools.product(bbox[:2], bbox[2:])) for bbox in unique_poly]
         if swap_latlon:
-            bboxes = [tuple(itertools.product(bbox[2:], bbox[:2])) for bbox in unique_poly]
+            bboxes = [((i[2], i[0]), (i[3], i[0]), (i[3], i[1]), (i[2], i[1])) for i in unique_poly]
+        else:
+            bboxes = [((i[0], i[2]), (i[0], i[3]), (i[1], i[3]), (i[1], i[2])) for i in unique_poly]
         # the spatial cells are arranged fast in latitude, so this only works for the specific csep1 file format
-        dh = float(unique_poly[0,3] - unique_poly[0,2])
+        dh = float(unique_poly[0, 3] - unique_poly[0, 2])
         # create CarteisanGrid of points
         region = CartesianGrid2D([Polygon(bbox) for bbox in bboxes], dh, mask=poly_mask)
         # get dims of 2d np.array
diff --git a/csep/core/regions.py b/csep/core/regions.py
index 73aa9123..e9551f03 100644
--- a/csep/core/regions.py
+++ b/csep/core/regions.py
@@ -8,6 +8,8 @@
 import numpy
 import numpy as np
 import mercantile
+from shapely import geometry
+from shapely.ops import unary_union
 
 # PyCSEP imports
 from csep.utils.calc import bin1d_vec, cleaner_range, first_nonnan, last_nonnan
@@ -806,35 +808,19 @@ def _build_bitmask_vec(self):
 
         return a, xs, ys
 
-    def tight_bbox(self):
-        # creates tight bounding box around the region, probably a faster way to do this.
-        ny, nx = self.idx_map.shape
-        asc = []
-        desc = []
-        for j in range(ny):
-            row = self.idx_map[j, :]
-            argmin = first_nonnan(row)
-            argmax = last_nonnan(row)
-            # points are stored clockwise
-            poly_min = self.polygons[int(row[argmin])].points
-            asc.insert(0, poly_min[0])
-            asc.insert(0, poly_min[1])
-            poly_max = self.polygons[int(row[argmax])].points
-            lat_0 = poly_max[2][1]
-            lat_1 = poly_max[3][1]
-            # last two points are 'right hand side of polygon'
-            if lat_0 < lat_1:
-                desc.append(poly_max[2])
-                desc.append(poly_max[3])
-            else:
-                desc.append(poly_max[3])
-                desc.append(poly_max[2])
-        # close the loop
-        poly = np.array(asc + desc)
+    def tight_bbox(self, precision=4):
+        # creates tight bounding box around the region
+        poly = np.array([i.points for i in self.polygons])
+
         sorted_idx = np.sort(np.unique(poly, return_index=True, axis=0)[1], kind='stable')
         unique_poly = poly[sorted_idx]
-        unique_poly = np.append(unique_poly, [unique_poly[0, :]], axis=0)
-        return unique_poly
+
+        # merges all the cell polygons into one
+        polygons = [geometry.Polygon(np.round(i, precision)) for i in unique_poly]
+        joined_poly = unary_union(polygons)
+        bounds = np.array([i for i in joined_poly.boundary.xy]).T
+
+        return bounds
 
     def get_cell_area(self):
         """ Compute the area of each polygon in sq. kilometers.
diff --git a/requirements.txt b/requirements.txt
index e5f26286..78dfc73b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -9,6 +9,7 @@ python-dateutil
 pytest
 vcrpy
 pytest-cov
+shapely
 sphinx
 sphinx-gallery
 sphinx-rtd-theme
diff --git a/requirements.yml b/requirements.yml
index 75a1c52a..85779ee0 100644
--- a/requirements.yml
+++ b/requirements.yml
@@ -13,6 +13,7 @@ dependencies:
   - python-dateutil
   - pytest
   - cartopy
+  - shapely
   - sphinx
   - sphinx-gallery
   - sphinx_rtd_theme
diff --git a/setup.py b/setup.py
index ae8d3fa3..f11037a0 100644
--- a/setup.py
+++ b/setup.py
@@ -34,7 +34,8 @@ def get_version():
         'obspy',
         'pyproj',
         'python-dateutil',
-        'mercantile'
+        'mercantile',
+        'shapely'
     ],
     extras_require = {
         'test': [