diff --git a/csep/core/forecasts.py b/csep/core/forecasts.py
index 323e1299..9d6599d9 100644
--- a/csep/core/forecasts.py
+++ b/csep/core/forecasts.py
@@ -437,7 +437,6 @@ def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plo
             axes: matplotlib.Axes.axes
         """
         # no mutable function arguments
-        dh = round(self.region.dh, 5)
         if self.start_time is None or self.end_time is None:
             time = 'forecast period'
         else:
@@ -451,12 +450,12 @@ def plot(self, ax=None, show=False, log=True, extent=None, set_global=False, plo
 
         # this call requires internet connection and basemap
         if log:
-            plot_args.setdefault('clabel', f'log10 M{self.min_magnitude}+ rate per {str(dh)}° x {str(dh)}° per {time}')
+            plot_args.setdefault('clabel', f'log10 M{self.min_magnitude}+ rate per cell per {time}')
             with numpy.errstate(divide='ignore'):
                 ax = plot_spatial_dataset(numpy.log10(self.spatial_counts(cartesian=True)), self.region, ax=ax,
                                           show=show, extent=extent, set_global=set_global, plot_args=plot_args)
         else:
-            plot_args.setdefault('clabel', f'M{self.min_magnitude}+ rate per {str(dh)}° x {str(dh)}° per {time}')
+            plot_args.setdefault('clabel', f'M{self.min_magnitude}+ rate per cell per {time}')
             ax = plot_spatial_dataset(self.spatial_counts(cartesian=True), self.region, ax=ax,show=show, extent=extent,
                                       set_global=set_global, plot_args=plot_args)
         return ax
diff --git a/csep/core/regions.py b/csep/core/regions.py
index edf82e23..48295b40 100644
--- a/csep/core/regions.py
+++ b/csep/core/regions.py
@@ -602,7 +602,7 @@ def get_cartesian(self, data):
 
     def get_bbox(self):
         """ Returns rectangular bounding box around region. """
-        return (self.xs.min(), self.xs.max(), self.ys.min(), self.ys.max())
+        return (self.xs.min(), self.xs.max()+self.dh, self.ys.min(), self.ys.max()+self.dh)
 
     def midpoints(self):
         """ Returns midpoints of rectangular polygons in region """
@@ -970,7 +970,6 @@ def __init__(self, polygons, quadkeys, bounds, name='QuadtreeGrid2d', mask=None)
         # self.xs = xs
         # self.ys = ys
         # self.idx_map = a
-        self.dh = 0.5 #Temporary use, until 'dh' is removed from plot_spatial_datasets() of forecast.plot.
 
     @property
     def num_nodes(self):
diff --git a/csep/utils/plots.py b/csep/utils/plots.py
index cbae3016..c8c20153 100644
--- a/csep/utils/plots.py
+++ b/csep/utils/plots.py
@@ -862,7 +862,7 @@ def plot_spatial_dataset(gridded, region, ax=None, show=False, extent=None, set_
     # Get spatial information for plotting
     bbox = region.get_bbox()
     if extent is None and not set_global:
-        extent = [bbox[0], bbox[1], bbox[2] + region.dh, bbox[3] + region.dh]
+        extent = [bbox[0], bbox[1], bbox[2], bbox[3]] 
 
     # Retrieve plot arguments
     plot_args = plot_args or {}
@@ -930,9 +930,7 @@ def plot_spatial_dataset(gridded, region, ax=None, show=False, extent=None, set_
     cmap = matplotlib.colors.ListedColormap(cmap_tup)
 
     ## Plot spatial dataset
-    lons, lats = numpy.meshgrid(numpy.append(region.xs, region.xs[-1] + region.dh),
-                                numpy.append(region.ys, region.ys[-1] + region.dh))
-
+    lons, lats = numpy.meshgrid(numpy.append(region.xs, bbox[1]), numpy.append(region.ys, bbox[3]))
     im = ax.pcolor(lons, lats, gridded, cmap=cmap, alpha=alpha, snap=True, transform=ccrs.PlateCarree())
     im.set_clim(clim)