diff --git a/autotest/t007_test.py b/autotest/t007_test.py index 0511177cab..84acb477b6 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,64 +237,125 @@ def test_write_shapefile(): def test_export_array(): - try: - from scipy.ndimage import rotate - except: - 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: - pass - if rasterio is not None: + 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: 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 + 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, + 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) + 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, + 0.0, -500.0, 5029.03, 0.0, 0.0, 1.0]) + def test_mbase_sr(): import numpy as np @@ -553,10 +615,31 @@ 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 + 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( + sr1.geotransform, + (-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 assert sr1.xll == xll assert sr1.yll == yll + assert sr1.xlength == 1250.0 + assert sr1.ylength == 5000.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( + sr1.geotransform, + (-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, xul=xul, yul=yul, rotation=30) @@ -585,12 +668,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 +687,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..c4eaaa663d 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): @@ -57,6 +58,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 @@ -85,6 +92,12 @@ class SpatialReference(object): 1D array of cell vertices for whole grid in C-style (row-major) order (same as np.ravel()) + resolution : tuple + model grid resolution (dx, dy) with length multiplier applied + + geotransform : tuple + GDAL-ordered geotransform tuple for exporting rasters with uniform + grid spacings. Raises ValueError if grid spacing is not uniform. Notes ----- @@ -158,44 +171,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.yedge[0] * - 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.yedge[0] * - 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.yedge[0] * - 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.yedge[0] * - 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 @@ -478,18 +488,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 @@ -635,11 +639,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): @@ -651,30 +654,13 @@ 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.yedge[0] * - self.length_multiplier) - self.yul = self.yll + (np.cos(self.theta) * self.yedge[0] * - 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.length_multiplier) - self.yll = self.yul - (np.cos(self.theta) * self.yedge[0] * - self.length_multiplier) - self._reset() - return + @property + def xlength(self): + return self.delr.sum() @property - def theta(self): - return -self.rotation * np.pi / 180. + def ylength(self): + return self.delc.sum() @property def xedge(self): @@ -732,17 +718,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): @@ -774,15 +755,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 +789,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 +810,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 +864,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 @@ -988,16 +971,16 @@ 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 ---------- filename : str Path of output file. Export format is determined by - file extention. + file extension. '.asc' Arc Ascii grid - '.tif' GeoTIFF (requries rasterio package) + '.tif' GeoTIFF (requires rasterio package) '.shp' Shapefile a : 2D numpy.ndarray Array to export @@ -1005,7 +988,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) @@ -1013,78 +996,81 @@ 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, dy = self.resolution + 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: + with open(filename, 'ab') as output: # for Python 3 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 - 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 @@ -1105,14 +1091,14 @@ 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: 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) @@ -1121,6 +1107,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, @@ -1289,6 +1278,38 @@ def _set_vertices(self): self._vertices = np.array(map(list, zip(ij, i1j, i1j1, ij1, ij))) """ + @property + 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 = None + if self.delc.min() != self.delc.max(): + dy = None + return dx, dy + + @property + def geotransform(self): + """Returns GDAL-ordered geotransform tuple""" + c, f = self.xul, self.yul + 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) + 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 @@ -2128,3 +2149,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)