diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index 3d4cda5e85..84acb477b6 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -297,6 +297,8 @@ def test_export_array():
 
     try:  # test GeoTIFF export
         import rasterio
+        from distutils.version import LooseVersion
+        rasterio_post101 = LooseVersion(rasterio.__version__) >= '1.0.1'
         m.sr.export_array(os.path.join(tpth, 'fb.tif'), m.dis.top.array,
                           nodata=nodata)
         with rasterio.open(os.path.join(tpth, 'fb.tif')) as src:
@@ -304,8 +306,10 @@ def test_export_array():
         # array should have identical values
         np.testing.assert_equal(arr, m.dis.top.array)
         assert src.shape == (m.nrow, m.ncol)
-        # Note: rasterio doesn't properly represent bounds for rotated rasters
-        np.testing.assert_almost_equal(src.bounds[0:2], m.sr.bounds[0:2])
+        if rasterio_post101:
+            np.testing.assert_almost_equal(src.bounds, m.sr.bounds)
+        else:  # Note: older rasterio had incomplete bounds on rotated rasters
+            np.testing.assert_almost_equal(src.bounds[0:2], m.sr.bounds[0:2])
         np.testing.assert_almost_equal(
             src.affine,
             [176.7766952966369, 176.77669529663686, 619653.0,
@@ -343,8 +347,10 @@ def test_export_array():
             arr = src.read(1)
         np.testing.assert_equal(arr, m.dis.top.array)
         assert src.shape == (10, 5)
-        # Note: rasterio doesn't properly represent bounds for rotated rasters
-        np.testing.assert_almost_equal(src.bounds[0:2], sr.bounds[0:2])
+        if rasterio_post101:
+            np.testing.assert_almost_equal(src.bounds, sr.bounds)
+        else:  # Note: older rasterio had incomplete bounds on rotated rasters
+            np.testing.assert_almost_equal(src.bounds[0:2], sr.bounds[0:2])
         np.testing.assert_equal(
             src.affine,
             [250.0, 0.0, 286.8,
@@ -611,8 +617,7 @@ def test_dynamic_xll_yll():
                                        xll=xll, yll=yll, rotation=30)
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
-    assert sr1.dx == 250.0
-    assert sr1.dy == 500.0
+    assert sr1.resolution == (250.0, 500.0)
     np.testing.assert_almost_equal(
         sr1.get_extent(), (-2213.2, 1369.3317547, 29.03, 4984.1570189))
     np.testing.assert_almost_equal(
@@ -626,8 +631,8 @@ def test_dynamic_xll_yll():
     assert sr1.yll == yll
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
-    assert sr1.dx == 250.0
-    assert sr1.dy == 500.0
+    np.testing.assert_almost_equal(
+        sr1.resolution, (76.1962816, 152.3925632))
     np.testing.assert_almost_equal(
         sr1.get_extent(), (-475.1628162, 616.7395778, 29.03, 1539.2790152))
     np.testing.assert_almost_equal(
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 1c8f2211b8..c4eaaa663d 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -92,14 +92,12 @@ class SpatialReference(object):
         1D array of cell vertices for whole grid in C-style (row-major) order
         (same as np.ravel())
 
-    dx : float
-        cell column width, or median for non-uniform widths along rows (delr)
-
-    dy : float
-        cell row height, or median for non-uniform widths along columns (delc)
+    resolution : tuple
+        model grid resolution (dx, dy) with length multiplier applied
 
     geotransform : tuple
-        GDAL-ordered geotransform tuple for exporting rasters
+        GDAL-ordered geotransform tuple for exporting rasters with uniform
+        grid spacings. Raises ValueError if grid spacing is not uniform.
 
     Notes
     -----
@@ -1036,8 +1034,7 @@ def export_array(self, filename, a, nodata=-9999,
                 xll, yll = xmin, ymin
             else:
                 xll, yll = self.xll, self.yll
-                dx = self.dx * self.length_multiplier
-                dy = self.dy * self.length_multiplier
+                dx, dy = self.resolution
                 if dx == dy:
                     cellsize = dx
             fmt = kwargs.get('fmt', '%.18e')
@@ -1282,31 +1279,23 @@ def _set_vertices(self):
         """
 
     @property
-    def dx(self):
-        dx = self.delr[0]
+    def resolution(self):
+        """Returns model grid resolution with length multiplier applied"""
+        dx = self.delr[0] * self.length_multiplier
+        dy = self.delc[0] * self.length_multiplier
         if self.delr.min() != self.delr.max():
-            dx = np.median(self.delr)
-            warn('delr values range from {0} to {1}; '
-                 'using median of {2} for dx'
-                 .format(self.delr.min(), self.delr.max(), dx))
-        return dx
-
-    @property
-    def dy(self):
-        dy = self.delc[0]
+            dx = None
         if self.delc.min() != self.delc.max():
-            dy = np.median(self.delc)
-            warn('delc values range from {0} to {1}; '
-                 'using median of {2} for dy'
-                 .format(self.delc.min(), self.delc.max(), dy))
-        return dy
+            dy = None
+        return dx, dy
 
     @property
     def geotransform(self):
         """Returns GDAL-ordered geotransform tuple"""
         c, f = self.xul, self.yul
-        dx = self.dx * self.length_multiplier
-        dy = self.dy * self.length_multiplier
+        dx, dy = self.resolution
+        if not all([dx, dy]):
+            raise ValueError('model grid is not uniform')
         rotation = self.rotation
         if rotation != 0:
             cr, sr = cos_sin(rotation)