From 799012564db50f305cba59bdeaf468f9d469996c Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Thu, 12 Jul 2018 21:28:08 +1200
Subject: [PATCH 1/6] add xlength and ylength properties to SpatialReference

add missing autotests for get_extent
---
 autotest/t007_test.py    | 18 ++++++++++---
 flopy/utils/reference.py | 56 ++++++++++++++++++++++++++--------------
 2 files changed, 50 insertions(+), 24 deletions(-)

diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index 0511177cab..e3549ad7f0 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -553,10 +553,20 @@ def test_dynamic_xll_yll():
     sr1 = flopy.utils.SpatialReference(delr=ms2.dis.delr.array,
                                        delc=ms2.dis.delc.array, lenuni=2,
                                        xll=xll, yll=yll, rotation=30)
+    assert sr1.xlength == 1250.0
+    assert sr1.ylength == 5000.0
+    np.testing.assert_almost_equal(
+        sr1.get_extent(), [-2213.2, 1369.3317547, 29.03, 4984.1570189])
+
     xul, yul = sr1.xul, sr1.yul
     sr1.length_multiplier = 1.0 / 3.281
     assert sr1.xll == xll
     assert sr1.yll == yll
+    assert sr1.xlength == 1250.0
+    assert sr1.ylength == 5000.0
+    np.testing.assert_almost_equal(
+        sr1.get_extent(), [-475.1628162, 616.7395778, 29.03, 1539.2790152])
+
     sr2 = flopy.utils.SpatialReference(delr=ms2.dis.delr.array,
                                        delc=ms2.dis.delc.array, lenuni=2,
                                        xul=xul, yul=yul, rotation=30)
@@ -585,12 +595,12 @@ def test_dynamic_xll_yll():
     assert np.abs(sr4.xul - (xul + 10.)) < 1e-6 # these shouldn't move because ul has priority
     assert np.abs(sr4.yul - (yul + 21.)) < 1e-6
     assert np.abs(sr4.xll - sr4.xul) < 1e-6
-    assert np.abs(sr4.yll - (sr4.yul - sr4.yedge[0])) < 1e-6
+    assert np.abs(sr4.yll - (sr4.yul - sr4.ylength)) < 1e-6
     sr4.xll = 0.
     sr4.yll = 10.
     assert sr4.origin_loc == 'll'
     assert sr4.xul == 0.
-    assert sr4.yul == sr4.yedge[0] + 10.
+    assert sr4.yul == sr4.ylength + 10.
     sr4.xul = xul
     sr4.yul = yul
     assert sr4.origin_loc == 'ul'
@@ -604,10 +614,10 @@ def test_dynamic_xll_yll():
                                        rotation=0, epsg=26915)
     sr5.lenuni = 1
     assert sr5.length_multiplier == .3048
-    assert sr5.yul == sr5.yll + sr5.yedge[0] * sr5.length_multiplier
+    assert sr5.yul == sr5.yll + sr5.ylength * sr5.length_multiplier
     sr5.lenuni = 2
     assert sr5.length_multiplier == 1.
-    assert sr5.yul == sr5.yll + sr5.yedge[0]
+    assert sr5.yul == sr5.yll + sr5.ylength
     sr5.proj4_str = '+proj=utm +zone=16 +datum=NAD83 +units=us-ft +no_defs'
     assert sr5.units == 'feet'
     assert sr5.length_multiplier == 1/.3048
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 80a4291f8f..415168e95e 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -57,6 +57,12 @@ class SpatialReference(object):
 
     Attributes
     ----------
+    xlength : float
+        horizontal length of model across columns (x-direction)
+
+    ylength : float
+        horizontal length of model across rows (y-direction)
+
     xedge : ndarray
         array of column edges
 
@@ -161,7 +167,7 @@ def xll(self):
             xll = self._xll if self._xll is not None else 0.
         elif self.origin_loc == 'ul':
             # calculate coords for lower left corner
-            xll = self._xul - (np.sin(self.theta) * self.yedge[0] *
+            xll = self._xul - (np.sin(self.theta) * self.ylength *
                                self.length_multiplier)
         return xll
 
@@ -171,7 +177,7 @@ def yll(self):
             yll = self._yll if self._yll is not None else 0.
         elif self.origin_loc == 'ul':
             # calculate coords for lower left corner
-            yll = self._yul - (np.cos(self.theta) * self.yedge[0] *
+            yll = self._yul - (np.cos(self.theta) * self.ylength *
                                self.length_multiplier)
         return yll
 
@@ -179,7 +185,7 @@ def yll(self):
     def xul(self):
         if self.origin_loc == 'll':
             # calculate coords for upper left corner
-            xul = self._xll + (np.sin(self.theta) * self.yedge[0] *
+            xul = self._xll + (np.sin(self.theta) * self.ylength *
                                self.length_multiplier)
         if self.origin_loc == 'ul':
             # calculate coords for lower left corner
@@ -190,7 +196,7 @@ def xul(self):
     def yul(self):
         if self.origin_loc == 'll':
             # calculate coords for upper left corner
-            yul = self._yll + (np.cos(self.theta) * self.yedge[0] *
+            yul = self._yll + (np.cos(self.theta) * self.ylength *
                                self.length_multiplier)
 
         if self.origin_loc == 'ul':
@@ -656,18 +662,18 @@ def set_origin(self, xul=None, yul=None, xll=None, yll=None):
             # calculate coords for upper left corner
             self._xll = xll if xll is not None else 0.
             self.yll = yll if yll is not None else 0.
-            self.xul = self._xll + (np.sin(self.theta) * self.yedge[0] *
+            self.xul = self._xll + (np.sin(self.theta) * self.ylength *
                                     self.length_multiplier)
-            self.yul = self.yll + (np.cos(self.theta) * self.yedge[0] *
+            self.yul = self.yll + (np.cos(self.theta) * self.ylength *
                                    self.length_multiplier)
 
         if self.origin_loc == 'ul':
             # calculate coords for lower left corner
             self.xul = xul if xul is not None else 0.
             self.yul = yul if yul is not None else 0.
-            self._xll = self.xul - (np.sin(self.theta) * self.yedge[0] *
+            self._xll = self.xul - (np.sin(self.theta) * self.ylength *
                                     self.length_multiplier)
-            self.yll = self.yul - (np.cos(self.theta) * self.yedge[0] *
+            self.yll = self.yul - (np.cos(self.theta) * self.ylength *
                                    self.length_multiplier)
         self._reset()
         return
@@ -676,6 +682,14 @@ def set_origin(self, xul=None, yul=None, xll=None, yll=None):
     def theta(self):
         return -self.rotation * np.pi / 180.
 
+    @property
+    def xlength(self):
+        return self.delr.sum()
+
+    @property
+    def ylength(self):
+        return self.delc.sum()
+
     @property
     def xedge(self):
         return self.get_xedge_array()
@@ -774,15 +788,15 @@ def transform(self, x, y, inverse=False):
 
     def get_extent(self):
         """
-        Get the extent of the rotated and offset grid
+        Get the extent of the scaled, rotated and offset grid
 
         Return (xmin, xmax, ymin, ymax)
 
         """
-        x0 = self.xedge[0]
-        x1 = self.xedge[-1]
-        y0 = self.yedge[0]
-        y1 = self.yedge[-1]
+        x0 = 0.0
+        x1 = self.xlength
+        y0 = self.ylength
+        y1 = 0.0
 
         # upper left point
         x0r, y0r = self.transform(x0, y0)
@@ -808,14 +822,16 @@ def get_grid_lines(self):
             Get the grid lines as a list
 
         """
-        xmin = self.xedge[0]
-        xmax = self.xedge[-1]
-        ymin = self.yedge[-1]
-        ymax = self.yedge[0]
+        xedge = self.xedge
+        yedge = self.yedge
+        xmin = xedge[0]
+        xmax = xedge[-1]
+        ymin = yedge[-1]
+        ymax = yedge[0]
         lines = []
         # Vertical lines
         for j in range(self.ncol + 1):
-            x0 = self.xedge[j]
+            x0 = xedge[j]
             x1 = x0
             y0 = ymin
             y1 = ymax
@@ -827,7 +843,7 @@ def get_grid_lines(self):
         for i in range(self.nrow + 1):
             x0 = xmin
             x1 = xmax
-            y0 = self.yedge[i]
+            y0 = yedge[i]
             y1 = y0
             x0r, y0r = self.transform(x0, y0)
             x1r, y1r = self.transform(x1, y1)
@@ -881,7 +897,7 @@ def get_yedge_array(self):
         rotated. Array is of size (nrow + 1)
 
         """
-        length_y = np.add.reduce(self.delc)
+        length_y = self.ylength
         yedge = np.concatenate(([length_y], length_y -
                                 np.add.accumulate(self.delc)))
         return yedge

From 354cd933238904a73c02a913365628c9ea4c3a7d Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Thu, 12 Jul 2018 23:04:16 +1200
Subject: [PATCH 2/6] tidy trigonometry for SpatialReference  - remove unused
 set_origin (assumed to not be needed?)  - remove theta property, which was
 used in equations that assumed CW    rotations (we use CCW elsewhere)  -
 reverse sign in xll and xul to use CCW rotation convention  - add generic
 cos_sin(deg), which snaps to right angles  - use shorthand logic "x or y"
 instead of "x if x is not None else y"

---
 flopy/utils/reference.py | 95 +++++++++++++++-------------------------
 1 file changed, 35 insertions(+), 60 deletions(-)

diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 415168e95e..505c23b94c 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -164,44 +164,41 @@ def __init__(self, delr=np.array([]), delc=np.array([]), lenuni=2,
     @property
     def xll(self):
         if self.origin_loc == 'll':
-            xll = self._xll if self._xll is not None else 0.
+            xll = self._xll or 0.
         elif self.origin_loc == 'ul':
             # calculate coords for lower left corner
-            xll = self._xul - (np.sin(self.theta) * self.ylength *
-                               self.length_multiplier)
+            _, sin_t = cos_sin(self.rotation)
+            xll = self._xul + (sin_t * self.ylength * self.length_multiplier)
         return xll
 
     @property
     def yll(self):
         if self.origin_loc == 'll':
-            yll = self._yll if self._yll is not None else 0.
+            yll = self._yll or 0.
         elif self.origin_loc == 'ul':
             # calculate coords for lower left corner
-            yll = self._yul - (np.cos(self.theta) * self.ylength *
-                               self.length_multiplier)
+            cos_t, _ = cos_sin(self.rotation)
+            yll = self._yul - (cos_t * self.ylength * self.length_multiplier)
         return yll
 
     @property
     def xul(self):
         if self.origin_loc == 'll':
             # calculate coords for upper left corner
-            xul = self._xll + (np.sin(self.theta) * self.ylength *
-                               self.length_multiplier)
+            _, sin_t = cos_sin(self.rotation)
+            xul = self._xll - (sin_t * self.ylength * self.length_multiplier)
         if self.origin_loc == 'ul':
-            # calculate coords for lower left corner
-            xul = self._xul if self._xul is not None else 0.
+            xul = self._xul or 0.
         return xul
 
     @property
     def yul(self):
         if self.origin_loc == 'll':
             # calculate coords for upper left corner
-            yul = self._yll + (np.cos(self.theta) * self.ylength *
-                               self.length_multiplier)
-
+            cos_t, _ = cos_sin(self.rotation)
+            yul = self._yll + (cos_t * self.ylength * self.length_multiplier)
         if self.origin_loc == 'ul':
-            # calculate coords for lower left corner
-            yul = self._yul if self._yul is not None else 0.
+            yul = self._yul or 0.
         return yul
 
     @property
@@ -484,18 +481,12 @@ def __setattr__(self, key, value):
         elif key == "length_multiplier":
             super(SpatialReference, self). \
                 __setattr__("_length_multiplier", float(value))
-            # self.set_origin(xul=self.xul, yul=self.yul, xll=self.xll,
-            #                yll=self.yll)
         elif key == "rotation":
             super(SpatialReference, self). \
                 __setattr__("rotation", float(value))
-            # self.set_origin(xul=self.xul, yul=self.yul, xll=self.xll,
-            #                yll=self.yll)
         elif key == "lenuni":
             super(SpatialReference, self). \
                 __setattr__("_lenuni", int(value))
-            # self.set_origin(xul=self.xul, yul=self.yul, xll=self.xll,
-            #                yll=self.yll)
         elif key == "units":
             value = value.lower()
             assert value in self.supported_units
@@ -641,11 +632,10 @@ def set_spatialreference(self, xul=None, yul=None, xll=None, yll=None,
             self.origin_loc = 'ul'
 
         self.rotation = rotation
-        self._xll = xll if xll is not None else 0.
-        self._yll = yll if yll is not None else 0.
-        self._xul = xul if xul is not None else 0.
-        self._yul = yul if yul is not None else 0.
-        # self.set_origin(xul, yul, xll, yll)
+        self._xll = xll or 0.
+        self._yll = yll or 0.
+        self._xul = xul or 0.
+        self._yul = yul or 0.
         return
 
     def __repr__(self):
@@ -657,31 +647,6 @@ def __repr__(self):
         s += "length_multiplier:{}".format(self.length_multiplier)
         return s
 
-    def set_origin(self, xul=None, yul=None, xll=None, yll=None):
-        if self.origin_loc == 'll':
-            # calculate coords for upper left corner
-            self._xll = xll if xll is not None else 0.
-            self.yll = yll if yll is not None else 0.
-            self.xul = self._xll + (np.sin(self.theta) * self.ylength *
-                                    self.length_multiplier)
-            self.yul = self.yll + (np.cos(self.theta) * self.ylength *
-                                   self.length_multiplier)
-
-        if self.origin_loc == 'ul':
-            # calculate coords for lower left corner
-            self.xul = xul if xul is not None else 0.
-            self.yul = yul if yul is not None else 0.
-            self._xll = self.xul - (np.sin(self.theta) * self.ylength *
-                                    self.length_multiplier)
-            self.yll = self.yul - (np.cos(self.theta) * self.ylength *
-                                   self.length_multiplier)
-        self._reset()
-        return
-
-    @property
-    def theta(self):
-        return -self.rotation * np.pi / 180.
-
     @property
     def xlength(self):
         return self.delr.sum()
@@ -746,17 +711,12 @@ def rotate(x, y, theta, xorigin=0., yorigin=0.):
         """
         Given x and y array-like values calculate the rotation about an
         arbitrary origin and then return the rotated coordinates.  theta is in
-        degrees.
+        degrees, positive CCW.
 
         """
-        # jwhite changed on Oct 11 2016 - rotation is now positive CCW
-        # theta = -theta * np.pi / 180.
-        theta = theta * np.pi / 180.
-
-        xrot = xorigin + np.cos(theta) * (x - xorigin) - np.sin(theta) * \
-               (y - yorigin)
-        yrot = yorigin + np.sin(theta) * (x - xorigin) + np.cos(theta) * \
-               (y - yorigin)
+        cos_t, sin_t = cos_sin(theta)
+        xrot = xorigin + cos_t * (x - xorigin) - sin_t * (y - yorigin)
+        yrot = yorigin + sin_t * (x - xorigin) + cos_t * (y - yorigin)
         return xrot, yrot
 
     def transform(self, x, y, inverse=False):
@@ -2144,3 +2104,18 @@ def getproj4(epsg):
         text for a projection (*.prj) file.
     """
     return get_spatialreference(epsg, text='proj4')
+
+
+def cos_sin(deg):
+    """Return cos and sin for a given angle in degrees"""
+    deg = deg % 360.0
+    if deg == 0.0:
+        return 1.0, 0.0
+    elif deg == 90.0:
+        return 0.0, 1.0
+    elif deg == 180.0:
+        return -1.0, 0.0
+    elif deg == 270.0:
+        return 0.0, -1.0
+    rad = np.radians(deg)
+    return np.cos(rad), np.sin(rad)

From b851fe5e61b8ecb0bb14a50979bd5e252fbe077a Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Sun, 15 Jul 2018 19:47:09 +1200
Subject: [PATCH 3/6] Add geotransform property to return GDAL-style affine
 transformation

Use this for export_array with .tif, since it is more robust with
older versions of affine (before v.2, rotation was clock-wise, and after
it is counter clock-wise. E.g. the added autotest previously failed with
affine 1.2.0 packaged with Ubuntu).
---
 autotest/t007_test.py    |  9 ++++++--
 flopy/utils/reference.py | 49 ++++++++++++++++++++++++++++++++--------
 2 files changed, 47 insertions(+), 11 deletions(-)

diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index e3549ad7f0..1345586541 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -238,7 +238,7 @@ def test_export_array():
 
     try:
         from scipy.ndimage import rotate
-    except:
+    except ImportError:
         rotate = False
         pass
 
@@ -284,7 +284,7 @@ def test_export_array():
     try:
         # test GeoTIFF export
         import rasterio
-    except:
+    except ImportError:
         pass
     if rasterio is not None:
         m.sr.export_array(os.path.join(tpth, 'fb.tif'), m.dis.top.array,
@@ -294,6 +294,11 @@ def test_export_array():
         assert src.shape == (m.nrow, m.ncol)
         assert np.abs(src.bounds[0] - m.sr.bounds[0]) < 1e-6
         assert np.abs(src.bounds[1] - m.sr.bounds[1]) < 1e-6
+        np.testing.assert_almost_equal(
+            src.affine,
+            [176.7766952966369, 176.77669529663686, 619653.0,
+             176.77669529663686, -176.7766952966369, 3353277.0, 0.0, 0.0, 1.0])
+
 
 def test_mbase_sr():
     import numpy as np
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 505c23b94c..6fa2ed2360 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -5,6 +5,7 @@
 import sys
 import os
 import numpy as np
+from warnings import warn
 
 
 class SpatialReference(object):
@@ -91,6 +92,10 @@ class SpatialReference(object):
         1D array of cell vertices for whole grid in C-style (row-major) order
         (same as np.ravel())
 
+    geotransform : tuple
+        GDAL-ordered geotransform tuple for exporting rasters. If delc or delr
+        are not constant, model cell size is determined from the median, and a
+        warning is shown.
 
     Notes
     -----
@@ -964,7 +969,7 @@ def plot_array(self, a, ax=None, **kwargs):
     def export_array(self, filename, a, nodata=-9999,
                      fieldname='value',
                      **kwargs):
-        """Write a numpy array to Arc Ascii grid
+        """Write a numpy array to Arc Ascii grid, GeoTIFF,
         or shapefile with the model reference.
 
         Parameters
@@ -973,7 +978,7 @@ def export_array(self, filename, a, nodata=-9999,
             Path of output file. Export format is determined by
             file extention.
             '.asc'  Arc Ascii grid
-            '.tif'  GeoTIFF (requries rasterio package)
+            '.tif'  GeoTIFF (requires rasterio package)
             '.shp'  Shapefile
         a : 2D numpy.ndarray
             Array to export
@@ -1049,18 +1054,14 @@ def export_array(self, filename, a, nodata=-9999,
             try:
                 import rasterio
                 from rasterio import Affine
-            except:
+            except ImportError:
                 print('GeoTIFF export requires the rasterio package.')
                 return
-            dxdy = self.delc[0] * self.length_multiplier
-            trans = Affine.translation(self.xul, self.yul) * \
-                    Affine.rotation(self.rotation) * \
-                    Affine.scale(dxdy, -dxdy)
 
             # third dimension is the number of bands
             a = a.copy()
             if len(a.shape) == 2:
-                a = np.reshape(a, (1, a.shape[0], a.shape[1]))
+                a.shape = (1,) + a.shape
             if a.dtype.name == 'int64':
                 a = a.astype('int32')
                 dtype = rasterio.int32
@@ -1081,7 +1082,7 @@ def export_array(self, filename, a, nodata=-9999,
                     'dtype': dtype,
                     'driver': 'GTiff',
                     'crs': self.proj4_str,
-                    'transform': trans
+                    'transform': Affine.from_gdal(*self.geotransform)
                     }
             meta.update(kwargs)
             with rasterio.open(filename, 'w', **meta) as dst:
@@ -1265,6 +1266,36 @@ def _set_vertices(self):
         self._vertices = np.array(map(list, zip(ij, i1j, i1j1, ij1, ij)))
         """
 
+    @property
+    def geotransform(self):
+        """Returns GDAL-ordered geotransform tuple"""
+        c, f = self.xul, self.yul
+        dx = self.delc[0] * self.length_multiplier
+        if self.delc.min() != self.delc.max():
+            dx = np.median(self.delc) * self.length_multiplier
+            warn('delc values range from {0} to {1}; '
+                 'using median of {2} for raster dx'
+                 .format(self.delc.min(), self.delc.max(), dx))
+        dy = self.delr[0] * self.length_multiplier
+        if self.delr.min() != self.delr.max():
+            dy = np.median(self.delr) * self.length_multiplier
+            warn('delr values range from {0} to {1}; '
+                 'using median of {2} for raster dy'
+                 .format(self.delr.min(), self.delr.max(), dy))
+        rotation = self.rotation
+        if rotation != 0:
+            cr, sr = cos_sin(rotation)
+            a = cr * dx
+            b = -sr * -dy
+            d = sr * dx
+            e = cr * -dy
+        else:
+            a = dx
+            e = -dy
+            b = d = 0.0
+        # GDAL order of affine transformation coefficients
+        return c, a, b, f, d, e
+
     def interpolate(self, a, xi, method='nearest'):
         """
         Use the griddata method to interpolate values from an array onto the

From 9c316c2574e33b7fc2b987829b96cfdc9ff103b9 Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Wed, 25 Jul 2018 22:30:01 +1200
Subject: [PATCH 4/6] Modify Arc Ascii grid export:  - Remove restriction that
 a uniform grid is needed; use dx and dy    fields, as is supported by GDAL
 and Surfer, but perhaps not others.  - Raise a warning when rotation is used,
 as it currently modifies    grid values via spline interpolation. Add TODO to
 consider removing    this feature, perhaps making an example in the notes or
 elsewhere.

Rewrite test_export_array:
 - First do some simpler tests on unrotated exports.
 - More compact reading/testing of .asc
 - Use np.testing.assert_almost_equal list/tuple/array checks
   that are not 100% identical.

Other changes:
 - Move dx/dy to properties, which raise warnings if not constant.
 - Rewrite part of 'Notes' section for export_array
 - Raise ValueError if filename extension not supported
 - Add some other tests for dx/dy in test_dynamic_xll_yll
---
 autotest/t007_test.py    | 155 +++++++++++++++++++++++++++------------
 flopy/utils/reference.py | 148 +++++++++++++++++++++----------------
 2 files changed, 195 insertions(+), 108 deletions(-)

diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index 1345586541..bb37545630 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -7,6 +7,7 @@
 import shutil
 import numpy as np
 import flopy
+import warnings
 
 pth = os.path.join('..', 'examples', 'data', 'mf2005_test')
 namfiles = [namfile for namfile in os.listdir(pth) if namfile.endswith('.nam')]
@@ -236,68 +237,118 @@ def test_write_shapefile():
 
 def test_export_array():
 
-    try:
-        from scipy.ndimage import rotate
-    except ImportError:
-        rotate = False
-        pass
-
     namfile = 'freyberg.nam'
     model_ws = '../examples/data/freyberg_multilayer_transient/'
     m = flopy.modflow.Modflow.load(namfile, model_ws=model_ws, verbose=False,
                                    load_only=['DIS', 'BAS6'])
-    m.sr.rotation = 45.
     nodata = -9999
-    m.sr.export_array(os.path.join(tpth, 'fb.asc'), m.dis.top.array, nodata=nodata)
-    arr = np.loadtxt(os.path.join(tpth, 'fb.asc'), skiprows=6)
-
-    m.sr.write_shapefile(os.path.join(tpth, 'grid.shp'))
-    # check bounds
-    with open(os.path.join(tpth, 'fb.asc')) as src:
-        for line in src:
-            if 'xllcorner' in line.lower():
-                val = float(line.strip().split()[-1])
-                # ascii grid origin will differ if it was unrotated
-                if rotate:
-                    assert np.abs(val - m.sr.bounds[0]) < 1e-6
-                else:
-                    assert np.abs(val - m.sr.xll) < 1e-6
-            if 'yllcorner' in line.lower():
-                val = float(line.strip().split()[-1])
-                if rotate:
-                    assert np.abs(val - m.sr.bounds[1]) < 1e-6
-                else:
-                    assert np.abs(val - m.sr.yll) < 1e-6
-            if 'cellsize' in line.lower():
-                val = float(line.strip().split()[-1])
-                rot_cellsize = np.cos(np.radians(m.sr.rotation)) * m.sr.delr[0] * m.sr.length_multiplier
-                #assert np.abs(val - rot_cellsize) < 1e-6
-                break
-    rotate = False
-    rasterio = None
-    if rotate:
-        rotated = rotate(m.dis.top.array, m.sr.rotation, cval=nodata)
-
-    if rotate:
-        assert rotated.shape == arr.shape
+    # write a simple unrotated Arc Ascii grid file first
+    m.sr.rotation = 0.
+    assert m.sr.geotransform == (619653.0, 250.0, 0.0, 3353277.0, 0.0, -250.0)
+    assert m.sr.bounds == (619653.0, 3343277.0, 624653.0, 3353277.0)
+    m.sr.export_array(os.path.join(tpth, 'fb-unrotated.asc'), m.dis.top.array,
+                      nodata=nodata)
+    with open(os.path.join(tpth, 'fb-unrotated.asc')) as src:
+        split_lines = list(src.readline().split() for _ in range(6))
+        header = {x[0]: int(x[1]) for x in split_lines[0:2]}
+        header.update({x[0]: float(x[1]) for x in split_lines[2:6]})
+        arr = np.loadtxt(src).reshape((m.nrow, m.ncol))
+    assert header == {'ncols': 20, 'nrows': 40,
+                      'xllcorner': 619653.0, 'yllcorner': 3343277.0,
+                      'cellsize': 250.0, 'NODATA_value': -9999.0}
+    np.testing.assert_almost_equal(arr, m.dis.top.array)
+
+    # modify rotation
+    m.sr.rotation = 45.
+    np.testing.assert_almost_equal(
+        m.sr.geotransform,
+        (619653.0, 176.7766952966369, 176.77669529663686,
+         3353277.0, 176.77669529663686, -176.7766952966369))
+    np.testing.assert_almost_equal(
+        m.sr.bounds,
+        (619653.0, 3346205.9321881346, 630259.6017177983, 3356812.5339059327))
+    m.sr.write_shapefile(os.path.join(tpth, 'grid.shp'))  # TODO: now what?
+    try:  # test Arc Ascii grid export on rotated grid
+        # TODO: consider removing this feature (e.g. see summary stats)
+        from scipy.ndimage import rotate
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            m.sr.export_array(os.path.join(tpth, 'fb.asc'), m.dis.top.array,
+                              nodata=nodata)
+        with open(os.path.join(tpth, 'fb.asc')) as src:
+            split_lines = list(src.readline().split() for _ in range(6))
+            header = {x[0]: int(x[1]) for x in split_lines[0:2]}
+            header.update({x[0]: float(x[1]) for x in split_lines[2:6]})
+            vals = np.ma.masked_equal(np.loadtxt(src), nodata)
+        expected = {
+            'ncols': 42, 'nrows': 42, 'xllcorner': 619653.0,
+            'yllcorner': 3346205.932188, 'cellsize': 252.53813613805409,
+            'NODATA_value': -9999.0}
+        assert set(header.keys()) == set(expected.keys())
+        for key in header.keys():
+            np.testing.assert_approx_equal(
+                header[key], expected[key], err_msg=key)
+        # Note: these are different than original values due to spline interp
+        np.testing.assert_almost_equal(vals.min(), -2.6153745651245117)
+        np.testing.assert_almost_equal(vals.mean(), 4.645456585243616)
+        np.testing.assert_almost_equal(vals.max(), 19.73752784729004)
+    except ImportError:  # scipy not available
+        pass
 
-    try:
-        # test GeoTIFF export
+    try:  # test GeoTIFF export
         import rasterio
-    except ImportError:
-        pass
-    if rasterio is not None:
         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:
             arr = src.read(1)
+        # array should have identical values
+        np.testing.assert_equal(arr, m.dis.top.array)
         assert src.shape == (m.nrow, m.ncol)
-        assert np.abs(src.bounds[0] - m.sr.bounds[0]) < 1e-6
-        assert np.abs(src.bounds[1] - m.sr.bounds[1]) < 1e-6
+        # Note: rasterio doesn't properly represent bounds for 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,
              176.77669529663686, -176.7766952966369, 3353277.0, 0.0, 0.0, 1.0])
+    except ImportError:  # rasterio not available
+        rasterio = False
+
+    # export arrays with different grid spacings in x and y directions
+    nlay, nrow, ncol = 1, 10, 5
+    delr, delc = 250, 500
+    xll, yll = 286.80, 29.03
+    m = flopy.modflow.Modflow()
+    dis = flopy.modflow.ModflowDis(m, nlay=nlay, nrow=nrow, ncol=ncol,
+                                   delr=delr, delc=delc)
+    sr = flopy.utils.SpatialReference(delr=m.dis.delr.array,
+                                      delc=m.dis.delc.array, lenuni=2,
+                                      xll=xll, yll=yll)
+    np.testing.assert_almost_equal(
+        sr.geotransform, (286.8, 250.0, 0.0, 5029.03, 0.0, -500.0))
+    sr.export_array(os.path.join(tpth, 'dxdy.asc'), m.dis.top.array,
+                    nodata=nodata)
+    with open(os.path.join(tpth, 'dxdy.asc')) as src:
+        split_lines = list(src.readline().split() for _ in range(7))
+        header = {x[0]: int(x[1]) for x in split_lines[0:2]}
+        header.update({x[0]: float(x[1]) for x in split_lines[2:7]})
+        arr = np.loadtxt(src).reshape((m.nrow, m.ncol))
+    assert header == {'ncols': 5, 'nrows': 10,
+                      'xllcorner': 286.8, 'yllcorner': 29.03,
+                      'dx': 250.0, 'dy': 500.0, 'NODATA_value': -9999.0}
+    np.testing.assert_almost_equal(arr, m.dis.top.array)
+    if rasterio:
+        sr.export_array(os.path.join(tpth, 'dxdy.tif'), m.dis.top.array,
+                        nodata=nodata)
+        with rasterio.open(os.path.join(tpth, 'dxdy.tif')) as src:
+            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])
+        np.testing.assert_equal(
+            src.affine,
+            [250.0, 0.0, 286.8,
+             0.0, -500.0, 5029.03, 0.0, 0.0, 1.0])
 
 
 def test_mbase_sr():
@@ -560,8 +611,14 @@ def test_dynamic_xll_yll():
                                        xll=xll, yll=yll, rotation=30)
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
+    assert sr1.dx == 500.0
+    assert sr1.dy == 250.0
     np.testing.assert_almost_equal(
         sr1.get_extent(), [-2213.2, 1369.3317547, 29.03, 4984.1570189])
+    np.testing.assert_almost_equal(
+        sr1.geotransform,
+        (-2213.2, 433.01270189221935, 125.0,
+         4359.157018922193, 250.0, -216.50635094610968))
 
     xul, yul = sr1.xul, sr1.yul
     sr1.length_multiplier = 1.0 / 3.281
@@ -569,8 +626,14 @@ def test_dynamic_xll_yll():
     assert sr1.yll == yll
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
+    assert sr1.dx == 500.0
+    assert sr1.dy == 250.0
     np.testing.assert_almost_equal(
         sr1.get_extent(), [-475.1628162, 616.7395778, 29.03, 1539.2790152])
+    np.testing.assert_almost_equal(
+        sr1.geotransform,
+        (-475.1628162145685, 131.97583111618997, 38.098140810728424,
+         1348.7883111618996, 76.19628162145685, -65.98791555809498))
 
     sr2 = flopy.utils.SpatialReference(delr=ms2.dis.delr.array,
                                        delc=ms2.dis.delc.array, lenuni=2,
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 6fa2ed2360..3a720284b3 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -92,10 +92,14 @@ 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)
+
     geotransform : tuple
-        GDAL-ordered geotransform tuple for exporting rasters. If delc or delr
-        are not constant, model cell size is determined from the median, and a
-        warning is shown.
+        GDAL-ordered geotransform tuple for exporting rasters
 
     Notes
     -----
@@ -976,7 +980,7 @@ def export_array(self, filename, a, nodata=-9999,
         ----------
         filename : str
             Path of output file. Export format is determined by
-            file extention.
+            file extension.
             '.asc'  Arc Ascii grid
             '.tif'  GeoTIFF (requires rasterio package)
             '.shp'  Shapefile
@@ -986,7 +990,7 @@ def export_array(self, filename, a, nodata=-9999,
             Value to assign to np.nan entries (default -9999)
         fieldname : str
             Attribute field name for array values (shapefile export only).
-            (default 'values')
+            (default 'value')
         kwargs: 
             keyword arguments to np.savetxt (ascii) 
             rasterio.open (GeoTIFF)
@@ -994,63 +998,70 @@ def export_array(self, filename, a, nodata=-9999,
             
         Notes
         -----
-        Rotated grids will be either be unrotated prior to export,
-        using scipy.ndimage.rotate (Arc Ascii format) or rotation will be
-        included in their transform property (GeoTiff format). In either case
-        the pixels will be displayed in the (unrotated) projected geographic coordinate system,
-        so the pixels will no longer align exactly with the model grid
-        (as displayed from a shapefile, for example). A key difference between
-        Arc Ascii and GeoTiff (besides disk usage) is that the
-        unrotated Arc Ascii will have a different grid size, whereas the GeoTiff
-        will have the same number of rows and pixels as the original.
+        Raster formats normally expect a uniform grid. If the spacing of rows
+        and columns are different, the Arc Ascii format will use dx and dy
+        instead of cellsize, as supported by GDAL and Golden Surfer (but
+        possibly not other software). Furthermore, if variable spacings are
+        used for either rows or columns, dx and/or dy are determined from
+        the median spacing distance, as described by a warning message.
+
+        Rotated grids are supported by GeoTIFF, and the output will preserve
+        the same values and dimension as the array. However, Arc Ascii formats
+        do no support rotation, and will be resampled using a spline
+        interpolation using scipy.ndimage.rotate, if available. As a result,
+        rotated Arc Ascii exports will have different values and dimensions,
+        and a warning message will be shown.
         """
-
-        if filename.lower().endswith(".asc"):
-            if len(np.unique(self.delr)) != len(np.unique(self.delc)) != 1 \
-                    or self.delr[0] != self.delc[0]:
-                raise ValueError('Arc ascii arrays require a uniform grid.')
-
-            xll, yll = self.xll, self.yll
-            cellsize = self.delr[0] * self.length_multiplier
-            fmt = kwargs.get('fmt', '%.18e')
-            a = a.copy()
-            a[np.isnan(a)] = nodata
+        ext = os.path.splitext(filename)[1].lower()
+        if ext == ".asc":
+            cellsize = None
+            rotate = False
             if self.rotation != 0:
                 try:
                     from scipy.ndimage import rotate
-                    a = rotate(a, self.rotation, cval=nodata)
-                    height_rot, width_rot = a.shape
-                    xmin, ymin, xmax, ymax = self.bounds
-                    dx = (xmax - xmin) / width_rot
-                    dy = (ymax - ymin) / height_rot
-                    cellsize = np.max((dx, dy))
-                    # cellsize = np.cos(np.radians(self.rotation)) * cellsize
-                    xll, yll = xmin, ymin
                 except ImportError:
-                    print('scipy package required to export rotated grid.')
-                    pass
-
-            filename = '.'.join(
-                filename.split('.')[:-1]) + '.asc'  # enforce .asc ending
-            nrow, ncol = a.shape
+                    print('scipy package required to export rotated '
+                          'Arc Ascii grid.')
+            if rotate:
+                # TODO: consider removing this feature
+                warn('using spline interpolation to potentially modify '
+                     'raster values in exported Arc Ascii file')
+                a = rotate(a, self.rotation, cval=nodata)
+                height_rot, width_rot = a.shape
+                xmin, ymin, xmax, ymax = self.bounds
+                dx = (xmax - xmin) / width_rot
+                dy = (ymax - ymin) / height_rot
+                cellsize = np.max((dx, dy))
+                # cellsize = np.cos(np.radians(self.rotation)) * cellsize
+                xll, yll = xmin, ymin
+            else:
+                xll, yll = self.xll, self.yll
+                dx = self.dx * self.length_multiplier
+                dy = self.dy * self.length_multiplier
+                if dx == dy:
+                    cellsize = dx
+            fmt = kwargs.get('fmt', '%.18e')
+            a = a.copy()
             a[np.isnan(a)] = nodata
-            txt = 'ncols  {:d}\n'.format(ncol)
-            txt += 'nrows  {:d}\n'.format(nrow)
-            txt += 'xllcorner  {:f}\n'.format(xll)
-            txt += 'yllcorner  {:f}\n'.format(yll)
-            txt += 'cellsize  {}\n'.format(cellsize)
+            nrow, ncol = a.shape
+            txt = (
+                'ncols  {:d}\n'
+                'nrows  {:d}\n'
+                'xllcorner  {:f}\n'
+                'yllcorner  {:f}\n'
+                .format(ncol, nrow, xll, yll))
+            if cellsize:
+                txt += 'cellsize  {}\n'.format(cellsize)
+            else:  # supported by GDAL, Surfer, but possibly not others
+                txt += 'dx  {}\ndy  {}\n'.format(dx, dy)
             # ensure that nodata fmt consistent w values
             txt += 'NODATA_value  {}\n'.format(fmt) % (nodata)
             with open(filename, 'w') as output:
                 output.write(txt)
-            with open(filename, 'ab') as output:
                 np.savetxt(output, a, **kwargs)
             print('wrote {}'.format(filename))
 
-        elif filename.lower().endswith(".tif"):
-            if len(np.unique(self.delr)) != len(np.unique(self.delc)) != 1 \
-                    or self.delr[0] != self.delc[0]:
-                raise ValueError('GeoTIFF export require a uniform grid.')
+        elif ext == ".tif":
             try:
                 import rasterio
                 from rasterio import Affine
@@ -1089,7 +1100,7 @@ def export_array(self, filename, a, nodata=-9999,
                 dst.write(a)
             print('wrote {}'.format(filename))
 
-        elif filename.lower().endswith(".shp"):
+        elif ext == ".shp":
             from ..export.shapefile_utils import write_grid_shapefile2
             epsg = kwargs.get('epsg', None)
             prj = kwargs.get('prj', None)
@@ -1098,6 +1109,9 @@ def export_array(self, filename, a, nodata=-9999,
             write_grid_shapefile2(filename, self, array_dict={fieldname: a},
                                   nan_val=nodata,
                                   epsg=epsg, prj=prj)
+        else:
+            raise ValueError('export_array: filename extension {!r} not '
+                             'supported'.format(ext))
 
     def export_contours(self, filename, contours,
                         fieldname='level', epsg=None, prj=None,
@@ -1266,22 +1280,32 @@ def _set_vertices(self):
         self._vertices = np.array(map(list, zip(ij, i1j, i1j1, ij1, ij)))
         """
 
+    @property
+    def dx(self):
+        dx = self.delr[0]
+        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]
+        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
+
     @property
     def geotransform(self):
         """Returns GDAL-ordered geotransform tuple"""
         c, f = self.xul, self.yul
-        dx = self.delc[0] * self.length_multiplier
-        if self.delc.min() != self.delc.max():
-            dx = np.median(self.delc) * self.length_multiplier
-            warn('delc values range from {0} to {1}; '
-                 'using median of {2} for raster dx'
-                 .format(self.delc.min(), self.delc.max(), dx))
-        dy = self.delr[0] * self.length_multiplier
-        if self.delr.min() != self.delr.max():
-            dy = np.median(self.delr) * self.length_multiplier
-            warn('delr values range from {0} to {1}; '
-                 'using median of {2} for raster dy'
-                 .format(self.delr.min(), self.delr.max(), dy))
+        dx = self.dx * self.length_multiplier
+        dy = self.dy * self.length_multiplier
         rotation = self.rotation
         if rotation != 0:
             cr, sr = cos_sin(rotation)

From 51f05b808ef42c57a08394f6597c83437ba30ffb Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Wed, 25 Jul 2018 23:32:17 +1200
Subject: [PATCH 5/6] corrections for a few unit tests

---
 autotest/t007_test.py    | 20 ++++++++++----------
 flopy/utils/reference.py |  1 +
 2 files changed, 11 insertions(+), 10 deletions(-)

diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index bb37545630..3d4cda5e85 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -611,14 +611,14 @@ def test_dynamic_xll_yll():
                                        xll=xll, yll=yll, rotation=30)
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
-    assert sr1.dx == 500.0
-    assert sr1.dy == 250.0
+    assert sr1.dx == 250.0
+    assert sr1.dy == 500.0
     np.testing.assert_almost_equal(
-        sr1.get_extent(), [-2213.2, 1369.3317547, 29.03, 4984.1570189])
+        sr1.get_extent(), (-2213.2, 1369.3317547, 29.03, 4984.1570189))
     np.testing.assert_almost_equal(
         sr1.geotransform,
-        (-2213.2, 433.01270189221935, 125.0,
-         4359.157018922193, 250.0, -216.50635094610968))
+        (-2213.2, 216.50635094610968, 250.0,
+         4359.157018922193, 125.0, -433.01270189221935))
 
     xul, yul = sr1.xul, sr1.yul
     sr1.length_multiplier = 1.0 / 3.281
@@ -626,14 +626,14 @@ def test_dynamic_xll_yll():
     assert sr1.yll == yll
     assert sr1.xlength == 1250.0
     assert sr1.ylength == 5000.0
-    assert sr1.dx == 500.0
-    assert sr1.dy == 250.0
+    assert sr1.dx == 250.0
+    assert sr1.dy == 500.0
     np.testing.assert_almost_equal(
-        sr1.get_extent(), [-475.1628162, 616.7395778, 29.03, 1539.2790152])
+        sr1.get_extent(), (-475.1628162, 616.7395778, 29.03, 1539.2790152))
     np.testing.assert_almost_equal(
         sr1.geotransform,
-        (-475.1628162145685, 131.97583111618997, 38.098140810728424,
-         1348.7883111618996, 76.19628162145685, -65.98791555809498))
+        (-475.1628162145685, 65.987915558094983, 76.196281621456848,
+         1348.7883111618996, 38.098140810728424, -131.97583111618997))
 
     sr2 = flopy.utils.SpatialReference(delr=ms2.dis.delr.array,
                                        delc=ms2.dis.delc.array, lenuni=2,
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 3a720284b3..1c8f2211b8 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -1058,6 +1058,7 @@ def export_array(self, filename, a, nodata=-9999,
             txt += 'NODATA_value  {}\n'.format(fmt) % (nodata)
             with open(filename, 'w') as output:
                 output.write(txt)
+            with open(filename, 'ab') as output:  # for Python 3
                 np.savetxt(output, a, **kwargs)
             print('wrote {}'.format(filename))
 

From 3f434b3cd53db53c66440b0e4b5668d0604481c6 Mon Sep 17 00:00:00 2001
From: Mike Toews <mwtoews@gmail.com>
Date: Wed, 8 Aug 2018 21:59:11 +1200
Subject: [PATCH 6/6] Refine properties after discussion

 - replace dx/dy with resolution with length multiplier applied
 - return None rather than median of delr/delc
 - geotransform raises ValueError for non-uniform grids
 - check full bounds for recent rasterio only
---
 autotest/t007_test.py    | 21 ++++++++++++--------
 flopy/utils/reference.py | 41 +++++++++++++++-------------------------
 2 files changed, 28 insertions(+), 34 deletions(-)

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)