Skip to content

Commit

Permalink
fix(_plot_util3d_helper, export_array): (#895)
Browse files Browse the repository at this point in the history
* fix naming convention for mf6 plotting in _plot_util3d_helper
* update Affine translation to use xul and yul in export_array
* code and docstring cleanups in utils.py
  • Loading branch information
jlarsen-usgs committed Jun 1, 2020
1 parent 8e70c92 commit 617b98d
Show file tree
Hide file tree
Showing 3 changed files with 46 additions and 31 deletions.
Binary file modified examples/data/mf6/test003_gwfs_disv/test003_gwfs_disv.dbf
Binary file not shown.
73 changes: 44 additions & 29 deletions flopy/export/utils.py
Expand Up @@ -228,8 +228,6 @@ def _add_output_nc_zonebudget_variable(f, array, var_name, flux,
Parameters
----------
f : NetCdf object
times : list
list of times
array : np.ndarray
zonebudget output budget group array
var_name : str
Expand Down Expand Up @@ -453,7 +451,6 @@ def output_helper(f, ml, oudic, **kwargs):
plotarray = np.atleast_3d(out_obj.get_alldata()
.transpose()).transpose()


for per in range(plotarray.shape[0]):
for k in range(plotarray.shape[1]):
if kper is not None and per != kper:
Expand Down Expand Up @@ -496,7 +493,7 @@ def output_helper(f, ml, oudic, **kwargs):
if logger:
logger.lraise("unrecognized export argument:{0}".format(f))
else:
raise NotImplementedError("unrecognized export argument" + \
raise NotImplementedError("unrecognized export argument" +
":{0}".format(f))
return f

Expand Down Expand Up @@ -620,9 +617,9 @@ def package_export(f, pak, fmt=None, **kwargs):
"error adding {0} as variable".format(a.name))
elif a.data_type == DataType.array3d:
f = array3d_export(f, a, **kwargs)
elif a.data_type == DataType.transient2d:
elif a.data_type == DataType.transient2d:
f = transient2d_export(f, a, **kwargs)
elif a.data_type == DataType.transientlist:
elif a.data_type == DataType.transientlist:
f = mflist_export(f, a, **kwargs)
elif isinstance(a, list):
for v in a:
Expand Down Expand Up @@ -729,9 +726,10 @@ def mflist_export(f, mfl, **kwargs):
model grid instance which will supercede the flopy.model.modelgrid
"""
assert isinstance(mfl, DataListInterface) and \
isinstance(mfl, DataInterface) \
, "mflist_helper only helps instances that support DataListInterface"
if not isinstance(mfl, (DataListInterface, DataInterface)):
err = "mflist_helper only helps instances that support " \
"DataListInterface"
raise AssertionError(err)

modelgrid = mfl.model.modelgrid
if "modelgrid" in kwargs:
Expand Down Expand Up @@ -770,8 +768,7 @@ def mflist_export(f, mfl, **kwargs):
else:
from ..export.shapefile_utils import recarray2shp
from ..utils.geometry import Polygon
#if np.isscalar(kper):
# kper = [kper]

df = mfl.get_dataframe(squeeze=squeeze)
if 'kper' in kwargs or df is None:
ra = mfl[kper]
Expand Down Expand Up @@ -853,7 +850,7 @@ def transient2d_export(f, t2d, fmt=None, **kwargs):
-----------
f : str
filename or existing export instance type (NetCdf only for now)
u2d : Transient2d instance
t2d : Transient2d instance
fmt : str
output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
Expand All @@ -865,9 +862,10 @@ def transient2d_export(f, t2d, fmt=None, **kwargs):
"""

assert isinstance(t2d, DataInterface) \
, "transient2d_helper only helps instances that support " \
"DataInterface"
if not isinstance(t2d, DataInterface):
err = "transient2d_helper only helps instances that support " \
"DataInterface"
raise AssertionError(err)

min_valid = kwargs.get("min_valid", -1.0e+9)
max_valid = kwargs.get("max_valid", 1.0e+9)
Expand Down Expand Up @@ -986,7 +984,7 @@ def array3d_export(f, u3d, fmt=None, **kwargs):
-----------
f : str
filename or existing export instance type (NetCdf only for now)
u2d : Util3d instance
u3d : Util3d instance
fmt : str
output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
Expand Down Expand Up @@ -1204,9 +1202,9 @@ def array2d_export(f, u2d, fmt=None, **kwargs):
array[array <= min_valid] = netcdf.FILLVALUE
array[array >= max_valid] = netcdf.FILLVALUE
if modelgrid.idomain is not None and \
"ibound" not in u2d.name.lower() and \
"idomain" not in u2d.name.lower() and \
"icbund" not in u2d.name.lower():
"ibound" not in u2d.name.lower() and \
"idomain" not in u2d.name.lower() and \
"icbund" not in u2d.name.lower():
array[modelgrid.idomain[0, :, :] == 0] = \
netcdf.FILLVALUE
var_name = u2d.name
Expand Down Expand Up @@ -1275,6 +1273,8 @@ def export_array(modelgrid, filename, a, nodata=-9999,
Parameters
----------
modelgrid : flopy.discretization.StructuredGrid object
model grid
filename : str
Path of output file. Export format is determined by
file extention.
Expand Down Expand Up @@ -1308,27 +1308,31 @@ def export_array(modelgrid, filename, a, nodata=-9999,
"""

if filename.lower().endswith(".asc"):
if len(np.unique(modelgrid.delr)) != len(np.unique(modelgrid.delc)) != 1 \
if len(np.unique(modelgrid.delr)) != \
len(np.unique(modelgrid.delc)) != 1 \
or modelgrid.delr[0] != modelgrid.delc[0]:
raise ValueError('Arc ascii arrays require a uniform grid.')

xoffset, yoffset = modelgrid.xoffset, modelgrid.yoffset
cellsize = modelgrid.delr[0] # * self.length_multiplier
cellsize = modelgrid.delr[0]
fmt = kwargs.get('fmt', '%.18e')
a = a.copy()
a[np.isnan(a)] = nodata
if modelgrid.angrot != 0:
try:
from scipy.ndimage import rotate
except ImportError:
rotate = None
print('scipy package required to export rotated grid.')

if rotate is not None:
a = rotate(a, modelgrid.angrot, cval=nodata)
height_rot, width_rot = a.shape
xmin, ymin, xmax, ymax = modelgrid.extent
dx = (xmax - xmin) / width_rot
dy = (ymax - ymin) / height_rot
cellsize = np.max((dx, dy))
xoffset, yoffset = xmin, ymin
except ImportError:
print('scipy package required to export rotated grid.')

filename = '.'.join(
filename.split('.')[:-1]) + '.asc' # enforce .asc ending
Expand All @@ -1348,7 +1352,8 @@ def export_array(modelgrid, filename, a, nodata=-9999,
print('wrote {}'.format(filename))

elif filename.lower().endswith(".tif"):
if len(np.unique(modelgrid.delr)) != len(np.unique(modelgrid.delc)) != 1 \
if len(np.unique(modelgrid.delr)) != \
len(np.unique(modelgrid.delc)) != 1 \
or modelgrid.delr[0] != modelgrid.delc[0]:
raise ValueError('GeoTIFF export require a uniform grid.')
try:
Expand All @@ -1357,10 +1362,14 @@ def export_array(modelgrid, filename, a, nodata=-9999,
except ImportError:
print('GeoTIFF export requires the rasterio package.')
return
dxdy = modelgrid.delc[0] # * self.length_multiplier
trans = Affine.translation(modelgrid.xoffset, modelgrid.yoffset) * \
Affine.rotation(modelgrid.angrot) * \
Affine.scale(dxdy, -dxdy)
dxdy = modelgrid.delc[0]
# because this is only implemented for a structured grid,
# we can get the xul and yul coordinate from modelgrid.xvertices(0, 0)
verts = modelgrid.get_cell_vertices(0, 0)
xul, yul = verts[0]
trans = Affine.translation(xul, yul) * \
Affine.rotation(modelgrid.angrot) * \
Affine.scale(dxdy, -dxdy)

# third dimension is the number of bands
a = a.copy()
Expand Down Expand Up @@ -1499,7 +1508,7 @@ def export_contourf(filename, contours, fieldname='level', epsg=None,

try:
from shapely import geometry
except:
except (ImportError, ModuleNotFoundError):
raise ImportError('export_contourf requires python shapely package')

from ..utils.geometry import Polygon
Expand Down Expand Up @@ -1563,12 +1572,16 @@ def export_array_contours(modelgrid, filename, a,
Parameters
----------
modelgrid : flopy.discretization.Grid object
model grid object
filename : str
Path of output file with '.shp' extention.
a : 2D numpy array
Array to contour
fieldname : str
gis field name
interval : float
interval to calculate levels from
levels : list
list of contour levels
maxlevels : int
Expand Down Expand Up @@ -1609,6 +1622,8 @@ def contour_array(modelgrid, ax, a, **kwargs):
Parameters
----------
modelgrid : flopy.discretization.Grid object
modelgrid object
ax : matplotlib.axes.Axes
ax to add the contours
Expand Down
4 changes: 2 additions & 2 deletions flopy/plot/plotutil.py
Expand Up @@ -977,9 +977,9 @@ def _plot_util3d_helper(util3d, filename_base=None,
else:
i0 = 0
i1 = array.shape[0]
# build filenames
# build filenames, use local "name" variable (flopy6 adaptation)
filenames = ['{}_{}_Layer{}.{}'.format(
filename_base, util3d.name[k],
filename_base, name[k],
k + 1, fext)
for k in range(i0, i1)]

Expand Down

0 comments on commit 617b98d

Please sign in to comment.