From 2f99268268bf448b2b778f17548c5ce8306fa139 Mon Sep 17 00:00:00 2001 From: martindjm <43825685+martindjm@users.noreply.github.com> Date: Fri, 8 Nov 2019 09:24:26 -0800 Subject: [PATCH 1/5] feat(vtk) updated vtk export capabilities in flopy. Added support for array 2d/3d export, package export, model export, heads, cell by cell file. export interface change added fmt option to specify 'vtk' added examples/Notebooks/flopy3_vtk_export.ipynb a notebook for vtk export created vtk tests in autotest/t050_test.py --- autotest/t050_test.py | 187 +++- examples/Notebooks/flopy3_vtk_export.ipynb | 348 +++++++ flopy/export/utils.py | 84 +- flopy/export/vtk.py | 1071 ++++++++++++++++---- 4 files changed, 1421 insertions(+), 269 deletions(-) create mode 100644 examples/Notebooks/flopy3_vtk_export.ipynb diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 68d4ea5035..a1cbdfb72a 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -2,7 +2,7 @@ import os import numpy as np import flopy -from flopy.export.vtk import Vtk +from flopy.export import vtk # create output directory cpth = os.path.join('temp', 't050') @@ -11,65 +11,150 @@ os.makedirs(cpth) -def test_vtkoutput(): - """Make vtk with ibound_filter""" - nlay = 3 - nrow = 3 - ncol = 3 - ml = flopy.modflow.Modflow() - dis = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, top=0, - botm=[-1., -2., -3.]) - ibound = np.ones((nlay, nrow, ncol), dtype=np.int) - ibound[0, 1, 1] = 0 - bas = flopy.modflow.ModflowBas(ml, ibound=ibound) - - fvtkout = os.path.join(cpth, 'test.vtu') - vtkfile = Vtk(fvtkout, ml) - a = np.arange(nlay * nrow * ncol).reshape((nlay, nrow, ncol)) - vtkfile.add_array('testarray', a) - vtkfile.write(shared_vertex=False, ibound_filter=True) +def test_vtk_export_array2d(): + """Export 2d array""" + mpath = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + namfile = 'freyberg.nam' + m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False) + m.dis.top.export(os.path.join(cpth, 'array_2d_test'), fmt='vtk') + # with smoothing + m.dis.top.export(os.path.join(cpth, 'array_2d_test'), fmt='vtk', + name='top_smooth', smooth=True) + + +def test_vtk_export_array3d(): + """Vtk export 3d array""" + mpath = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + namfile = 'freyberg.nam' + m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False) + m.upw.hk.export(os.path.join(cpth, 'array_3d_test'), fmt='vtk') + # with point scalars + m.upw.hk.export(os.path.join(cpth, 'array_3d_test'), fmt='vtk', + name='hk_points', point_scalars=True) + + return + + +def test_vtk_transient_array_2d(): + """VTK export transient 2d array""" + mpath = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + namfile = 'freyberg.nam' + m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False) + m.rch.rech.export(os.path.join(cpth, 'transient_2d_test'), fmt='vtk') + + return + + +def test_vtk_export_packages(): + """testing vtk package export""" + mpath = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + namfile = 'freyberg.nam' + m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False) + # test dis export + m.dis.export(os.path.join(cpth, 'DIS'), fmt='vtk') + # upw with point scalar output + m.upw.export(os.path.join(cpth, 'UPW'), fmt='vtk', point_scalars=True) + # bas with smoothing on + m.bas6.export(os.path.join(cpth, 'BAS'), fmt='vtk', smooth=True) + # transient package drain + m.drn.export(os.path.join(cpth, 'DRN'), fmt='vtk') + + return + + +# add mf2005 model exports +def test_export_mf2005_vtk(): + """test vtk model export mf2005""" + pth = os.path.join('..', 'examples', 'data', 'mf2005_test') + namfiles = [namfile for namfile in os.listdir(pth) if + namfile.endswith('.nam')] + skip = ['bcf2ss.nam'] + for namfile in namfiles: + if namfile in skip: + continue + print('testing namefile', namfile) + m = flopy.modflow.Modflow.load(namfile, model_ws=pth, verbose=False) + m.export(os.path.join(cpth, m.name), fmt='vtk') + return -def test_vtkoutput_noibound(): - """Make vtk without ibound_filter""" - nlay = 3 - nrow = 3 - ncol = 3 - ml = flopy.modflow.Modflow() - dis = flopy.modflow.ModflowDis(ml, nlay=nlay, nrow=nrow, ncol=ncol, top=0, - botm=[-1., -2., -3.]) - fvtkout = os.path.join(cpth, 'test.vtu') - vtkfile = Vtk(fvtkout, ml) - a = np.arange(nlay * nrow * ncol).reshape((nlay, nrow, ncol)) - vtkfile.add_array('testarray', a) - vtkfile.write(shared_vertex=False, ibound_filter=False) +def test_vtk_mf6(): + mf6expth = os.path.join('..', 'examples', 'data', 'mf6') + # test vtk mf6 export + mf6sims = ['test005_advgw_tidal', 'test045_lake1ss_table', + 'test036_twrihfb', 'test045_lake2tr', 'test006_2models_mvr'] + # mf6sims = ['test005_advgw_tidal'] + # mf6sims = ['test036_twrihfb'] + + for simnm in mf6sims: + print(simnm) + simpth = os.path.join(mf6expth, simnm) + loaded_sim = flopy.mf6.MFSimulation.load(simnm, 'mf6', 'mf6', + simpth) + sim_models = loaded_sim.model_names + print(sim_models) + for mname in sim_models: + print(mname) + m = loaded_sim.get_model(mname) + m.export(os.path.join(cpth, m.name), fmt='vtk') + return -def test_vtkoutput_mf6(): - """Make vtk with ibound_filter""" - nlay = 3 - nrow = 3 - ncol = 3 - sim = flopy.mf6.MFSimulation() - gwf = flopy.mf6.ModflowGwf(sim) - idomain = np.ones((nlay, nrow, ncol), dtype=np.int) - idomain[0, 1, 1] = 0 - dis = flopy.mf6.ModflowGwfdis(gwf, nlay=nlay, nrow=nrow, ncol=ncol, top=0, - botm=[-1., -2., -3.], idomain=idomain) - - fvtkout = os.path.join(cpth, 'test.vtu') - vtkfile = Vtk(fvtkout, gwf) - a = np.arange(nlay * nrow * ncol).reshape((nlay, nrow, ncol)) - vtkfile.add_array('testarray', a) - vtkfile.write(shared_vertex=False, ibound_filter=True) +def test_vtk_bindary_head_export(): + + """test vet export of heads""" + + freyberg_pth = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + + hdsfile = os.path.join(freyberg_pth, 'freyberg.hds') + + m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_pth, + verbose=False) + otfolder = os.path.join(cpth, 'freyberg_hds_test') + + vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, + 1090]) + # test with points + vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, + 1090], point_scalars=True) + + # test vtk export heads with smoothing and no point scalars + vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, + 1090], + point_scalars=False, smooth=True) + return +def test_vtk_cbc(): + # test mf 2005 freyberg + freyberg_cbc = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient', 'freyberg.cbc') + + freyberg_mpth = os.path.join('..', 'examples', 'data', + 'freyberg_multilayer_transient') + m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_mpth, + verbose=False) + + vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'), + kperlist=[1, 2, 3], point_scalars=True) + + return if __name__ == '__main__': - test_vtkoutput() - test_vtkoutput_noibound() - test_vtkoutput_mf6() \ No newline at end of file + test_vtk_export_array2d() + test_vtk_export_array3d() + test_vtk_transient_array_2d() + test_vtk_export_packages() + test_export_mf2005_vtk() + test_vtk_mf6() + test_vtk_bindary_head_export() + test_vtk_cbc() diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb new file mode 100644 index 0000000000..b2d9c705b5 --- /dev/null +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -0,0 +1,348 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "flopy is installed in C:\\Users\\domartin\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\flopy\n", + "3.7.3 (default, Mar 27 2019, 17:13:21) [MSC v.1915 64 bit (AMD64)]\n", + "flopy version: 3.2.12\n" + ] + } + ], + "source": [ + "import os\n", + "import sys\n", + "import matplotlib as mpl\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "\n", + "try:\n", + " import flopy\n", + "except:\n", + " fpth = os.path.abspath(os.path.join('..', '..'))\n", + " sys.path.append(fpth)\n", + " import flopy\n", + " \n", + "from flopy.export import vtk\n", + " \n", + "print(sys.version)\n", + "print('flopy version: {}'.format(flopy.__version__))" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# load model for examples\n", + "nam_file = \"freyberg.nam\"\n", + "model_ws = os.path.join(\"..\", \"data\", \"freyberg_multilayer_transient\")\n", + "ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# output folder to store the outputs from the notebook\n", + "workspace = os.path.join(\".\", \"VTK_EXAMPLE_OUTPUTS\")\n", + "if not os.path.exists(workspace):\n", + " os.mkdir(workspace)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# for all export calls a folder is give and the the fmt flag is set to 'vtk'\n", + "# to initiate vtk export. Outputs are writtent to the provided folder.\n", + "# if it does not exist one will be created." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# export arrays\n", + "# 2d array export\n", + "# export model top\n", + "# **kwargs: smooth: True creates a smooth surface\n", + "# point_scalars: True creates scalar data values on vertices as well as cells\n", + " # this automatically turns on smoothing\n", + "# name: the name will be used for the output file and the scalar variable name\n", + "ot_arrays_folder = os.path.join(workspace, 'arrays_test')\n", + "model_top_output = os.path.join(ot_arrays_folder, 'TOP_point')\n", + "ml.dis.top.export(model_top_output, fmt='vtk', **{'smooth': False, 'point_scalars': True}) " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# transient 2d array\n", + "# export recharge for each stress period\n", + "# **kwargs: smooth: True creates a smooth surface\n", + "# point_scalars: True creates scalar data values on vertices as well as cells\n", + " # this automatically turns on smoothing\n", + "# name: the name will be used for the output file and the scalar variable name\n", + "model_recharge_output = os.path.join(ot_arrays_folder, 'RECH')\n", + "ml.rch.rech.export(model_recharge_output, fmt='vtk')" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# 3D Array export\n", + "# export model bottoms\n", + "# **kwargs: smooth: True creates a smooth surface\n", + "# point_scalars: True creates scalar data values on vertices as well as cells\n", + " # this automatically turns on smoothing\n", + "# name: the name will be used for the output file and the scalar variable name\n", + "model_bottom_output = os.path.join(ot_arrays_folder, 'BOTM')\n", + "ml.dis.botm.export(model_bottom_output, smooth=True, fmt='vtk')\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# 3D Array export\n", + "# export model bottoms\n", + "model_hk_export = os.path.join(ot_arrays_folder, 'HK')\n", + "ml.upw.hk.export(model_hk_export, smooth=True, fmt='vtk', name='HK')" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "# package export\n", + "# **kwargs: smooth: True creates a smooth surface\n", + "# point_scalars: True creates scalar data values on vertices as well as cells\n", + " # this automatically turns on smoothing\n", + "package_output = os.path.join(workspace, 'package_output_test')\n", + "# export dis\n", + "dis_output = os.path.join(package_output, 'DIS')\n", + "ml.dis.export(dis_output, fmt='vtk')\n", + "#export upw with point scalars\n", + "upw_output = os.path.join(package_output, 'UPW')\n", + "ml.upw.export(upw_output, fmt='vtk', point_scalars=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "'.\\\\VTK_EXAMPLE_OUTPUTS\\\\model_output_test'" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\n", + "# model export\n", + "# **kwargs: smooth: True creates a smooth surface\n", + "# point_scalars: True creates scalar data values on vertices as well as cells\n", + " # this automatically turns on smoothing\n", + "model_output = os.path.join(workspace, 'model_output_test')\n", + "ml.export(model_output, fmt='vtk')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# export binary heads\n", + "# args\n", + "# export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None,\n", + " # smooth=False, point_scalars=False, nanval=-1e+20):\n", + "# can limit time steps and stress periods with kstplist and kperlist\n", + "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", + "heads_output_folder = os.path.join(workspace, 'heads_output_test')\n", + "vtk.export_heads(ml, heads_file, heads_output_folder)\n", + "# outputs a .pvd file that can be loaded into paraview to view the heads as time series" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# export cbc file\n", + "#export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None,\n", + "# kperlist=None, keylist=None, smooth=False, point_scalars=False,\n", + "# nanval=-1e+20):\n", + "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", + "cbc_output_folder = os.path.join(workspace, 'cbc_output_test')\n", + "vtk.export_cbc(ml, cbc_file, cbc_output_folder)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# export heads with parameters\n", + "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", + "heads_output_folder = os.path.join(workspace, 'heads_output_test_parameters')\n", + "# kstp - time step\n", + "# kperlist - list of stress periods to export\n", + "# point_scalars - output heads as point values\n", + "vtk.export_heads(ml, heads_file, heads_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15], point_scalars=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", + "cbc_output_folder = os.path.join(workspace, 'cbc_output_test_parameters')\n", + "# kstp - time step\n", + "# kperlist - list of stress periods to export\n", + "# keylist - list of cbc terms to export\n", + "vtk.export_cbc(ml, cbc_file, cbc_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15],\n", + " keylist=['CONSTANT HEAD', 'STORAGE'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "nam_file = \"mnw1.nam\"\n", + "model_ws = os.path.join(\"..\", \"data\", \"mf2005_test\")\n", + "ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "ml.mnw1.export(os.path.join(workspace, 'mnw1_test'), fmt='vtk')" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "nam_file = \"lakeex3.nam\"\n", + "model_ws = os.path.join(\"..\", \"data\", \"mf2005_test\")\n", + "ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)\n", + "ml.lak.export(os.path.join(workspace, 'lak_Test'), fmt='vtk')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/flopy/export/utils.py b/flopy/export/utils.py index b1afc94b97..672c130fae 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -8,6 +8,7 @@ from ..datbase import DataType, DataInterface, DataListInterface from . import NetCdf, netcdf from . import shapefile_utils +from . import vtk NC_PRECISION_TYPE = {np.float64: "f8", np.float32: "f4", np.int: "i4", @@ -343,7 +344,7 @@ def output_helper(f, ml, oudic, **kwargs): return f -def model_export(f, ml, **kwargs): +def model_export(f, ml, fmt=None, **kwargs): """ Method to export a model to a shapefile or netcdf file @@ -354,6 +355,8 @@ def model_export(f, ml, **kwargs): or dictionary of .... ml : flopy.modflow.mbase.ModelInterface object flopy model object + fmt: str + output format flag. set to 'vtk' to export to vtk .vtu file **kwargs : keyword arguments modelgrid: flopy.discretization.Grid @@ -364,6 +367,14 @@ def model_export(f, ml, **kwargs): prj : str prj file name + smooth: bool + For vtk, if set to True will output smooth surface + + point_scalars: bool + For vtk, if set to True will create point scalar values along + with cell values. + + """ assert isinstance(ml, ModelInterface) @@ -391,13 +402,20 @@ def model_export(f, ml, **kwargs): for pak in ml.packagelist: f = package_export(f, pak, **kwargs) + elif fmt == 'vtk': + # call vtk model export + smooth = kwargs.get('smooth', False) + point_scalars = kwargs.get('point_scalars', False) + vtk.export_model(ml, f, package_names=package_names, smooth=smooth, + point_scalars=point_scalars) + else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) return f -def package_export(f, pak, **kwargs): +def package_export(f, pak, fmt=None, **kwargs): """ Method to export a package to shapefile or netcdf @@ -407,6 +425,8 @@ def package_export(f, pak, **kwargs): output file name (ends in .shp for shapefile or .nc for netcdf) pak : flopy.pakbase.Package object package to export + fmt: str + output format flag. set to 'vtk' to export to vtk .vtu file ** kwargs : keword arguments modelgrid: flopy.discretization.Grid user supplied modelgrid object which will supercede the built @@ -416,6 +436,13 @@ def package_export(f, pak, **kwargs): prj : str prj file name + smooth: bool + For vtk, if set to True will output smooth surface + + point_scalars: bool + For vtk, if set to True will create point scalar values along + with cell values. + Returns ------- f : NetCdf object or None @@ -455,6 +482,14 @@ def package_export(f, pak, **kwargs): f = array3d_export(f, v, **kwargs) return f + elif fmt == 'vtk': + # call vtk array export to folder + smooth = kwargs.get('smooth', False) + point_scalars = kwargs.get('point_scalars', False) + vtk.export_package(pak.parent, pak.name, f, smooth=smooth, + point_scalars=point_scalars) + + else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -655,7 +690,7 @@ def mflist_export(f, mfl, **kwargs): raise NotImplementedError("unrecognized export argument:{0}".format(f)) -def transient2d_export(f, t2d, **kwargs): +def transient2d_export(f, t2d, fmt=None, **kwargs): """ export helper for Transient2d instances @@ -664,11 +699,15 @@ def transient2d_export(f, t2d, **kwargs): f : str filename or existing export instance type (NetCdf only for now) u2d : Transient2d instance + fmt: set to 'vtk' to export a .vtu file **kwargs : keyword arguments min_valid : minimum valid value max_valid : maximum valid value modelgrid : flopy.discretization.Grid model grid instance which will supercede the flopy.model.modelgrid + smooth: for vtk to output a smooth represenation of the model + point_scalars: for vtk to output point value scalars as well as cell + name: for vtk to set a specific name for array and output file """ @@ -768,11 +807,17 @@ def transient2d_export(f, t2d, **kwargs): raise Exception(estr) return f + elif fmt == 'vtk': + smooth = kwargs.get('smooth', False) + point_scalars = kwargs.get('point_scalars', False) + name = kwargs.get('name', t2d.name) + vtk.export_transient(t2d.model, t2d.array, f, name, smooth=smooth, + point_scalars=point_scalars, array2d=True) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) -def array3d_export(f, u3d, **kwargs): +def array3d_export(f, u3d, fmt=None, **kwargs): """ export helper for Transient2d instances @@ -781,11 +826,15 @@ def array3d_export(f, u3d, **kwargs): f : str filename or existing export instance type (NetCdf only for now) u2d : Util3d instance + fmt: set to 'vtk' to export a .vtu file **kwargs : keyword arguments min_valid : minimum valid value max_valid : maximum valid value modelgrid : flopy.discretization.Grid model grid instance which will supercede the flopy.model.modelgrid + smooth: for vtk to output a smooth represenation of the model + point_scalars: for vtk to output point value scalars as well as cell + name: for vtk to set a specific name for array and output file """ @@ -903,11 +952,23 @@ def array3d_export(f, u3d, **kwargs): raise Exception(estr) return f + elif fmt == 'vtk': + # call vtk array export to folder + smooth = kwargs.get('smooth', False) + point_scalars = kwargs.get('point_scalars', False) + name = kwargs.get('name', u3d.name) + + if isinstance(name, list) or isinstance(name, tuple): + name = name[0] + + vtk.export_array(u3d.model, u3d.array, f, name, smooth=smooth, + point_scalars=point_scalars) + else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) -def array2d_export(f, u2d, **kwargs): +def array2d_export(f, u2d, fmt=None, **kwargs): """ export helper for Util2d instances @@ -916,11 +977,15 @@ def array2d_export(f, u2d, **kwargs): f : str filename or existing export instance type (NetCdf only for now) u2d : Util2d instance + fmt: set to 'vtk' to export a .vtu file **kwargs : keyword arguments min_valid : minimum valid value max_valid : maximum valid value modelgrid : flopy.discretization.Grid model grid instance which will supercede the flopy.model.modelgrid + smooth: for vtk to output a smooth represenation of the model + point_scalars: for vtk to output point value scalars as well as cell + name: for vtk to set a specific name for array and output file """ assert isinstance(u2d, DataInterface), "util2d_helper only helps " \ @@ -1013,6 +1078,15 @@ def array2d_export(f, u2d, **kwargs): raise Exception(estr) return f + elif fmt == 'vtk': + + # call vtk array export to folder + smooth = kwargs.get('smooth', False) + point_scalars = kwargs.get('point_scalars', False) + name = kwargs.get('name', u2d.name) + vtk.export_array(u2d.model, u2d.array, f, name, smooth=smooth, + point_scalars=point_scalars, array2d=True) + else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 87852d2d91..2cce473a66 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -2,9 +2,18 @@ import os import numpy as np from ..discretization import StructuredGrid +from ..datbase import DataType, DataInterface +import flopy.utils.binaryfile as bf +from flopy.utils import HeadFile +import numpy.ma as ma + +""" +Module for exporting vtk from flopy +""" def start_tag(f, tag, indent_level, indent_char=' '): + # starts xml tag s = indent_level * indent_char + tag indent_level += 1 f.write(s + '\n') @@ -12,113 +21,201 @@ def start_tag(f, tag, indent_level, indent_char=' '): def end_tag(f, tag, indent_level, indent_char=' '): + # ends xml tag indent_level -= 1 s = indent_level * indent_char + tag f.write(s + '\n') return indent_level +class _Array(object): + # class to store array and tell if array is 2d + def __init__(self, array, array2d): + self.array = array + self.array2d = array2d + + +def _get_basic_modeltime(perlen_list): + modeltim = 0 + totim = [] + for tim in perlen_list: + totim.append(modeltim) + modeltim += tim + return totim + + class Vtk(object): - """ - Support for writing a model to a vtk file - """ + def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, + point_scalars=False): + + """ + Make Vtk object for exporting flopy vtk + + :param model: flopy model instance + :param verbose: if True, stdout is verbose + :param nanval: no data value + :param smooth: if True will create a smooth output surface + :param point_scalars: if True will output poin scalar values, + this will set smooth to True. + """ - def __init__(self, output_filename, model, verbose=None): + if point_scalars: + smooth = True - assert output_filename.lower().endswith(".vtu") if verbose is None: verbose = model.verbose self.verbose = verbose - if os.path.exists(output_filename): - if self.verbose: - print('removing existing vtk file: ' + output_filename) - os.remove(output_filename) - self.output_filename = output_filename - + # set up variables self.model = model self.modelgrid = model.modelgrid + self.arrays = {} self.shape = (self.modelgrid.nlay, self.modelgrid.nrow, self.modelgrid.ncol) - + self.shape2d = (self.shape[1], self.shape[2]) + self.nlay = self.modelgrid.nlay + self.nrow = self.modelgrid.nrow + self.ncol = self.modelgrid.ncol + self.ncol = self.modelgrid.ncol + self.nanval = nanval + + self.cell_type = 11 self.arrays = {} - return + self.smooth = smooth + self.point_scalars = point_scalars + + # check if structured grid, vtk only supports structured grid + assert (isinstance(self.modelgrid, StructuredGrid)) + + self.cbd_on = False + + # get ibound + if self.modelgrid.idomain is None: + # ibound = None + ibound = np.ones(self.shape) + else: + ibound = self.modelgrid.idomain + + if ibound is not None and hasattr(self.model, 'dis') and \ + hasattr(self.model.dis, 'laycbd'): + + self.cbd = np.where(self.model.dis.laycbd.array > 0) + ibound = np.insert(ibound, self.cbd[0] + 1, ibound[self.cbd[ + 0], :, :], + axis=0) + self.cbd_on = True + + self.ibound = ibound + + # build the model vertices + self.verts, self.iverts, self.zverts = self.get_3d_vertex_connectivity() - def add_array(self, name, a): - assert a.shape == self.shape - self.arrays[name] = a return - def write(self, shared_vertex=False, ibound_filter=False, htop=None): + def add_array(self, name, a, array2d=False): + """ - Parameters - ---------- - shared_vertex : bool - Make a smoothed representation of model layers by calculating - an interpolated z value for the cell corner vertices. - - ibound_filter : bool - Use the ibound array in the basic package of the model that - was passed in to exclude cells in the vtk file that have an - ibound value of zero. - - htop : ndarray - This array must of shape (nlay, nrow, ncol). If htop is passed - then these htop values will be used to set the z elevation for the - cell tops. This makes it possible to show cells based on the - saturated thickness. htop should be calculated by the user as the - minimum of the cell top and the head and the maximum of the cell - bottom and the head. + Adds an array to the vtk object + :param name: name of the array + :param a: the array to be added to the vtk object + :param array2d: True if the array is 2d """ - indent_level = 0 - if self.verbose: - print('writing vtk file') - f = open(self.output_filename, 'w') + if name == 'ibound': + return - ibound = None - if ibound_filter: - ibound = self.modelgrid.idomain + # limited recarray support + if type(a) == np.recarray: + matrix = np.zeros(self.shape) + for row in a: + matrix[a.k, a.i, a.j ] = 1 + a = matrix - dis = self.model.dis - z = np.vstack([dis.top.array.reshape(1, self.modelgrid.nrow, - self.modelgrid.ncol), - dis.botm.array]) - if shared_vertex: - verts, iverts = self.get_3d_shared_vertex_connectivity( - self.modelgrid) - #verts, iverts = dis.sr.get_3d_shared_vertex_connectivity(nlay, - # z, - # ibound=ibound) else: - top = z[:-1] - bot = z[1:] - if htop is not None: - top = htop - verts, iverts = self.get_3d_vertex_connectivity(self.modelgrid) - #verts, iverts = dis.sr.get_3d_vertex_connectivity(nlay, top, - # bot, - # ibound=ibound) + + # if array is 2d reformat to 3d array + if array2d: + assert a.shape == self.shape2d + array = np.full(self.shape, self.nanval) + array[0, :, :] = a + a = array + + try: + + assert a.shape == self.shape + except AssertionError: + return + # raise AssertionError + if a.dtype == int: + a = a.astype(float) + # add array to self.arrays + self.arrays[name] = a + return + + def write(self, output_file, timeval=None): + """ + writes the arrays from self.arrays to vtk file + + :param output_file: output file name to write the vtk data + :param timeval: model time value to be stored in the time section of + the vtk file + """ + + # make sure file ends with vtu + assert output_file.lower().endswith(".vtu") + + # get the active data cells based on the data arrays and ibound + actwcells3d = self._configure_data_arrays() + actwcells = actwcells3d.ravel() + + # get the indexes of the active cells + idxs = np.argwhere(actwcells != 0).ravel() + + # get the verts and iverts to be output + verts = [self.verts[idx] for idx in idxs] + iverts = self._build_iverts(verts) + + # get the total number of cells and vertices ncells = len(iverts) - npoints = verts.shape[0] + npoints = ncells * 8 + if self.verbose: - s = 'Number of point is {}\n ' \ - 'Number of cells is {}\n'.format(npoints, ncells) - print(s) + print('Writing vtk file: ' + output_file) + print('Number of point is {}, Number of cells is {}\n'.format( + npoints, ncells)) - # xml + # open output file for writing + f = open(output_file, 'w') + + # write xml + indent_level = 0 s = '' f.write(s + '\n') indent_level = start_tag(f, '', indent_level) - # unstructured grid indent_level = start_tag(f, '', indent_level) + # if time value write time section + if timeval: + + indent_level = start_tag(f, '', indent_level) + + s = '' + indent_level = start_tag(f, s, indent_level) + + f.write(indent_level * ' ' + '{}\n'.format(timeval)) + + indent_level = end_tag(f, '', indent_level) + + indent_level = end_tag(f, '', indent_level) + # piece s = ''.format(npoints, ncells) @@ -131,18 +228,10 @@ def write(self, shared_vertex=False, ibound_filter=False, htop=None): s = '' indent_level = start_tag(f, s, indent_level) assert (isinstance(self.modelgrid, StructuredGrid)) - z = np.vstack([self.modelgrid.top.reshape(1, self.modelgrid.nrow, - self.modelgrid.ncol), - self.modelgrid.botm]) - if shared_vertex: - verts, iverts = self.get_3d_shared_vertex_connectivity( - self.modelgrid) - else: - verts, iverts = self.get_3d_vertex_connectivity(self.modelgrid) - - for row in verts: - s = indent_level * ' ' + '{} {} {} \n'.format(*row) - f.write(s) + for cell in verts: + for row in cell: + s = indent_level * ' ' + '{} {} {} \n'.format(*row) + f.write(s) s = '' indent_level = end_tag(f, s, indent_level) @@ -174,7 +263,7 @@ def write(self, shared_vertex=False, ibound_filter=False, htop=None): s = '' indent_level = start_tag(f, s, indent_level) for row in iverts: - s = indent_level * ' ' + '{} \n'.format(11) + s = indent_level * ' ' + '{} \n'.format(self.cell_type) f.write(s) s = '' indent_level = end_tag(f, s, indent_level) @@ -186,14 +275,28 @@ def write(self, shared_vertex=False, ibound_filter=False, htop=None): s = '' indent_level = start_tag(f, s, indent_level) - self._write_data_array(f, indent_level, 'top', z[0:-1], ibound) - - for name, a in self.arrays.items(): - self._write_data_array(f, indent_level, name, a, ibound) + # write data arrays to file + for arrayName, arrayValues in self.arrays.items(): + self.write_data_array(f, indent_level, arrayName, + arrayValues, actwcells3d) s = '' indent_level = end_tag(f, s, indent_level) + # write arrays as point data if point scalars is set to True + if self.point_scalars: + s = '' + indent_level = start_tag(f, s, indent_level) + for array_name, array_values in self.arrays.items(): + self.write_point_value(f, indent_level, array_values, + array_name, actwcells3d) + + s = '' + indent_level = end_tag(f, s, indent_level) + + else: + pass + # end piece indent_level = end_tag(f, '', indent_level) @@ -205,172 +308,714 @@ def write(self, shared_vertex=False, ibound_filter=False, htop=None): # end file f.close() + self.arrays.clear() return - def _write_data_array(self, f, indent_level, name, a, ibound): + def _configure_data_arrays(self): + """ + Compares all the stored arrays int the vtk class along with the + ibound to figure out what cells and points need to be written to the + output vtk file. """ - Write a numpy array to the vtk file + # get 1d shape + shape1d = self.shape[0] * self.shape[1] * self.shape[2] + + # build index array + ot_idx_array = np.zeros(shape1d, dtype=np.int) + # loop through arrays + for name, array in self.arrays.items(): + # make array 1d + a = array.ravel() + # find where no data + where_nan = np.isnan(a) + # where no data set to the class nan val + a[where_nan] = self.nanval + # get the indexes where there is data + idxs = np.argwhere(a != self.nanval) + # set the active array to 1 + ot_idx_array[idxs] = 1 + + # reset the shape of the active data array + ot_idx_array = ot_idx_array.reshape(self.shape) + # where the ibound is 0 set the active array to 0 + ot_idx_array[self.ibound == 0] = 0 + + return ot_idx_array + + def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): """ - # header tag - s = ''.format(name) + :param actwcells: array of where data exists + :param zvalues: array of values to be used instead of the zvalues of + the vertices. This allows point scalars to be interpolated. + :return: dictionaries of verts, iverts, and zvalues + """ + # set up active cells + if actwcells is not None: + ncells = int(actwcells.sum()) + else: + ncells = self.nlay * self.nrow * self.ncol + actwcells = self.ibound + npoints = ncells * 8 + ipoint = 0 + + vertsdict = {} + ivertsdict = {} + zvertsdict = {} + + # if smoothing interpolate the z values + if self.smooth: + if zvalues is not None: + # use the given data array values + zVertices = self.extendedDataArray(zvalues) + else: + zVertices = self.extendedDataArray(self.modelgrid.top_botm) + else: + zVertices = None + + # model cellid based on 1darray + # build the vertices + cellid = -1 + for k in range(self.nlay): + for i in range(self.nrow): + for j in range(self.ncol): + cellid += 1 + if actwcells[k, i, j] == 0: + continue + verts = [] + ivert = [] + zverts = [] + pts = self.modelgrid._cell_vert_list(i, j) + pt0, pt1, pt2, pt3, pt0 = pts + + if not self.smooth: + cellBot = self.modelgrid.top_botm[k + 1, i, j] + cellTop = self.modelgrid.top_botm[k, i, j] + celElev = [cellBot, cellTop] + for elev in celElev: + # verts[ipoint, :] = np.append(pt1,elev) + verts.append([pt1[0], pt1[1], elev]) + # verts[ipoint+1, :] = np.append(pt2,elev) + verts.append([pt2[0], pt2[1], elev]) + # verts[ipoint+2, :] = np.append(pt0,elev) + verts.append([pt0[0], pt0[1], elev]) + # verts[ipoint+3, :] = np.append(pt3,elev) + verts.append([pt3[0], pt3[1], elev]) + ivert.extend([ipoint, ipoint+1, ipoint+2, ipoint+3]) + zverts.extend([elev, elev, elev, elev]) + ipoint += 4 + vertsdict[cellid] = verts + ivertsdict[cellid] = ivert + zvertsdict[cellid] = zverts + else: + layers = [k+1, k] + for lay in layers: + verts.append([pt1[0], pt1[1], zVertices[lay, i+1, + j]]) + verts.append([pt2[0], pt2[1], zVertices[lay, i+1, + j+1]]) + + verts.append([pt0[0], pt0[1], zVertices[lay, i, j]]) + + verts.append([pt3[0], pt3[1], zVertices[lay, i, + j+1]]) + ivert.extend([ipoint, ipoint+1, ipoint+2, ipoint+3]) + zverts.extend([zVertices[lay, i+1, j], zVertices[ + lay, i+1, j+1], + zVertices[lay, i, j], zVertices[lay, i, + j+1]]) + ipoint += 4 + vertsdict[cellid] = verts + ivertsdict[cellid] = ivert + zvertsdict[cellid] = zverts + return vertsdict, ivertsdict, zvertsdict + + def extendedDataArray(self, dataArray): + if dataArray.shape[0] == self.nlay+1: + dataArray = dataArray + else: + listArray = [dataArray[0]] + for lay in range(dataArray.shape[0]): + listArray.append(dataArray[lay]) + dataArray = np.stack(listArray) + + matrix = np.zeros([self.nlay+1, self.nrow+1, self.ncol+1]) + for lay in range(self.nlay+1): + for row in range(self.nrow+1): + for col in range(self.ncol+1): + + indexList = [[row-1, col-1], [row-1, col], [row, col-1], + [row, col]] + neighList = [] + for index in indexList: + if index[0] in range(self.nrow) and index[1] in \ + range(self.ncol): + neighList.append(dataArray[lay, index[0], index[1]]) + neighList = np.array(neighList) + if neighList[neighList > self.nanval].shape[0] > 0: + headMean = neighList[neighList > self.nanval].mean() + else: + headMean = self.nanval + matrix[lay, row, col] = headMean + return matrix + + def write_data_array(self, f, indent_level, arrayName, arrayValues, + actWCells): + """ + + Writes the data array to the output vtk file + + :param f: output vtk file + :param indent_level: current indent of the xml + :param arrayName: name of the output array + :param arrayValues: the data array being output + :param actWCells: array of the active cells + + """ + + s = ''.format( + arrayName) indent_level = start_tag(f, s, indent_level) # data - nlay = a.shape[0] - - # combine ibound with laycbd when model supports laycbd - if ibound is not None and hasattr(self.model, 'dis') and \ - hasattr(self.model.dis, 'laycbd'): - cbd = np.where(self.model.dis.laycbd.array > 0) - ibound = np.insert(ibound, cbd[0] + 1, ibound[cbd[0], :, :], - axis=0) + nlay = arrayValues.shape[0] - for k in range(nlay): + for lay in range(nlay): s = indent_level * ' ' f.write(s) - if ibound is None: - ak = a[k].flatten() - else: - idx = (ibound[k] != 0) - ak = a[k][idx].flatten() - for v in ak: - s = ' {}'.format(v) + idx = (actWCells[lay] != 0) + arrayValuesLay = arrayValues[lay][idx].flatten() + for layValues in arrayValuesLay: + s = ' {}'.format(layValues) f.write(s) f.write('\n') + s = '' + indent_level = end_tag(f, s, indent_level) + return + + def write_point_value(self, f, indent_level, data_array, array_name, + actwcells): + """ + Writes the data array to the output vtk file as point scalars + """ + # header tag + s = ''.format( + array_name) + indent_level = start_tag(f, s, indent_level) + + # data + verts, secus, zverts = self.get_3d_vertex_connectivity( + actwcells=actwcells, zvalues=data_array) + + for cellid in sorted(zverts): + for z in zverts[cellid]: + s = indent_level * ' ' + f.write(s) + s = ' {}'.format(z) + f.write(s) + f.write('\n') + # ending tag s = '' indent_level = end_tag(f, s, indent_level) return @staticmethod - def get_3d_shared_vertex_connectivity(mg): - - # get the x and y points for the grid - x, y = mg.xyzvertices[0:2] - x = x.flatten() - y = y.flatten() - - # set the size of the vertex grid - nrowvert = mg.nrow + 1 - ncolvert = mg.ncol + 1 - nlayvert = mg.nlay + 1 - nrvncv = nrowvert * ncolvert - npoints = nrvncv * nlayvert - - # create and fill a 3d points array for the grid - verts = np.empty((npoints, 3), dtype=np.float) - verts[:, 0] = np.tile(x, nlayvert) - verts[:, 1] = np.tile(y, nlayvert) - istart = 0 - istop = nrvncv - top_botm = mg.top_botm - for k in range(mg.nlay + 1): - verts[istart:istop, 2] = mg.interpolate(top_botm[k], - verts[istart:istop, :2], - method='linear') - istart = istop - istop = istart + nrvncv - - # create the list of points comprising each cell. points must be - # listed a specific way according to vtk requirements. - iverts = [] - for k in range(mg.nlay): - koffset = k * nrvncv - for i in range(mg.nrow): - for j in range(mg.ncol): - if mg._idomain is not None: - if mg._idomain[k, i, j] == 0: - continue - iv1 = i * ncolvert + j + koffset - iv2 = iv1 + 1 - iv4 = (i + 1) * ncolvert + j + koffset - iv3 = iv4 + 1 - iverts.append([iv4 + nrvncv, iv3 + nrvncv, - iv1 + nrvncv, iv2 + nrvncv, - iv4, iv3, iv1, iv2]) - - return verts, iverts + def _build_iverts(verts): + """ - @staticmethod - def get_3d_vertex_connectivity(mg): - if mg.idomain is None: - ncells = mg.nlay * mg.nrow * mg.ncol - ibound = np.ones((mg.nlay, mg.nrow, mg.ncol), dtype=np.int) - else: - ncells = (mg.idomain != 0).sum() - ibound = mg.idomain - npoints = ncells * 8 - verts = np.empty((npoints, 3), dtype=np.float) + Builds the iverts based on the vertices being output + + :param verts: vertices being output + :return: iverts for output + """ + count = 0 iverts = [] - ipoint = 0 - top_botm = mg.top_botm - for k in range(mg.nlay): - for i in range(mg.nrow): - for j in range(mg.ncol): - if ibound[k, i, j] == 0: - continue + for a in verts: + ivert = [] + for i in range(len(a)): + ivert.append(count) + count += 1 + iverts.append(ivert) + iverts = np.array(iverts) + return iverts + + +def _get_names(in_list): + ot_list = [] + for x in in_list: + if isinstance(x, bytes): + ot_list.append(str(x.decode('UTF-8'))) + else: + ot_list.append(x) + return ot_list - ivert = [] - pts = mg._cell_vert_list(i, j) - pt0, pt1, pt2, pt3, pt0 = pts - z = top_botm[k + 1, i, j] +def export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None, + kperlist=None, keylist=None, smooth=False, point_scalars=False, + nanval=-1e+20): + """ + + Exports cell by cell file to vtk - verts[ipoint, 0:2] = np.array(pt1) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + :param model: the flopy model instance + :param cbcfile: the cell by cell file + :param otfolder: output folder to write the data to + :param precision: bindary file precision + :param kstplist: list of timesteps to be written + :param kperlist: list of stress periods to be writeen + :param keylist: list of flow term names to be output + :param smooth: If true a smooth surface will be output + :param point_scalars: If True point scalar values will be written + :param nanval: the no data value + + """ - verts[ipoint, 0:2] = np.array(pt2) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + mg = model.modelgrid + shape = (mg.nlay, mg.nrow, mg.ncol) - verts[ipoint, 0:2] = np.array(pt0) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + if not os.path.exists(otfolder): + os.mkdir(otfolder) - verts[ipoint, 0:2] = np.array(pt3) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + # set up the pvd file to make the output files time enabled + pvdfile = open( + os.path.join(otfolder, '{}_Heads.pvd'.format(model.name)), + 'w') - z = top_botm[k, i, j] + pvdfile.write(""" + + \n""") - verts[ipoint, 0:2] = np.array(pt1) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + # load cbc - verts[ipoint, 0:2] = np.array(pt2) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + cbb = bf.CellBudgetFile(cbcfile, precision=precision) - verts[ipoint, 0:2] = np.array(pt0) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + # totim_dict = dict(zip(cbb.get_kstpkper(), model.dis.get_totim())) - verts[ipoint, 0:2] = np.array(pt3) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + # get records + records = _get_names(cbb.get_unique_record_names()) - iverts.append(ivert) + # build imeth lookup + imeth_dict = {record: imeth for (record, imeth) in zip(records, + cbb.imethlist)} + # get list of packages to export + if not keylist: + keylist = records + + if not kperlist: + kperlist = list(set([x[1] for x in cbb.get_kstpkper() if x[1] > -1])) + else: + kperlist = [kper - 1 for kper in kperlist] + + if not kstplist: + kstplist = list(set([x[0] for x in cbb.get_kstpkper() if x[0] > -1])) + else: + kstplist = [kstp - 1 for kstp in kstplist] + + # get model name + model_name = model.name + + vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) + + # export data + addarray = False + count = 1 + for kper in kperlist: + for kstp in kstplist: + + ot_base = '{}_CBC_KPER{}_KSTP{}.vtu'.format( + model_name, kper + 1, kstp + 1) + otfile = os.path.join(otfolder, ot_base) + pvdfile.write("""\n""".format(count, ot_base)) + for name in keylist: + + try: + rec = cbb.get_data(kstpkper=(kstp, kper), text=name, + full3D=True) + + if len(rec) > 0: + array = rec[0] # need to fix for multiple pak + addarray = True + + except ValueError: + + rec = cbb.get_data(kstpkper=(kstp, kper), text=name)[0] + + if imeth_dict[name] == 6: + array = np.full(shape, nanval) + # rec array + for [node, q] in zip(rec['node'], rec['q']): + lyr, row, col = np.unravel_index(node - 1, shape) + + array[lyr, row, col] = q + + addarray = True + else: + raise Exception('Data type not currenlty supported ' + 'for cbc output') + + if addarray: + + # set the data to no data value + if ma.is_masked(array): + array = np.where(array.mask, nanval, array) + + # add array to vtk + vtk.add_array(name.strip(), array) # need to adjust for + + # write the vtk data to the output file + vtk.write(otfile) + count += 1 + # finish writing the pvd file + pvdfile.write(""" +""") + + pvdfile.close() + return + + +def export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None, + smooth=False, point_scalars=False, nanval=-1e+20): + """ + Exports heads to vtk files by timestep and stressperiod + + :param model: the model instance + :param hdsfile: the binary heads file + :param otfolder: the output folder to write the .vtu files + :param kstplist: list of time steps to output + :param kperlist: list of stress periods to output + :param smooth: If set to True a smooth surface will be output + :param point_scalars: If set to True the heads will be written to the + vertices as point scalars as well as cell values + :param nanval: The no data value + :return: Heads will be written to files named by stress period and + timestep to the otfolder + """ - return verts, iverts + # setup output folder + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + # start writing the pvd file to make the data time aware + pvdfile = open(os.path.join(otfolder, '{}_Heads.pvd'.format(model.name)), + 'w') + + pvdfile.write(""" + + \n""") + # get teh ehads + hds = HeadFile(hdsfile) + + # get the model time + # totim_dict = dict(zip(hds.get_kstpkper(), model.dis.get_totim())) + + # set up time step and stress periods for output + if not kperlist: + kperlist = list(set([x[1] for x in hds.get_kstpkper() if x[1] > -1])) + else: + kperlist = [kper - 1 for kper in kperlist] + + if not kstplist: + kstplist = list(set([x[0] for x in hds.get_kstpkper() if x[0] > -1])) + else: + kstplist = [kstp - 1 for kstp in kstplist] + + # set upt the vtk + vtk = Vtk(model, smooth=smooth, point_scalars=point_scalars, nanval=nanval) + + # output data + count = 0 + for kper in kperlist: + for kstp in kstplist: + hdarr = hds.get_data((kstp, kper)) + vtk.add_array('head', hdarr) + ot_base = '{}_Heads_KPER{}_KSTP{}.vtu'.format( + model.name, kper + 1, kstp + 1) + otfile = os.path.join(otfolder, ot_base) + # vtk.write(otfile, timeval=totim_dict[(kstp, kper)]) + vtk.write(otfile) + pvdfile.write("""\n""".format(count, ot_base)) + count += 1 + + pvdfile.write(""" +""") + + pvdfile.close() + + +def export_array(model, array, output_folder, name, nanval=-1e+20, + array2d=False, smooth=False, point_scalars=False): + + """ + + :param model: array the model belongs to + :param array: array to be exported + :param arrayname: name of the array to be used in vtk file + :param output_folder: output folder to store the output .vtu file + :param array2d: True if array is 2d + :nanval nan value + :return: outputs a .vtu file of array + + """ + + if not os.path.exists(output_folder): + os.mkdir(output_folder) + + vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) + vtk.add_array(name, array, array2d=array2d) + otfile = os.path.join(output_folder, '{}.vtu'.format(name)) + vtk.write(otfile) + + return + + +def export_transient(model, array, output_folder, name, nanval=-1e+20, + array2d=False, smooth=False, point_scalars=False): + + """ + + :param model: model of transient array + :param array: transient array to export + :param name: name of the data + :param output_folder: output folder location + :param array2d: True if array is 2d + :param nanval: nan value + :return: ouputs .vtu files of transient data to output folder + """ + + if not os.path.exists(output_folder): + os.mkdir(output_folder) + + to_tim = model.dis.get_totim() + + vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) + + if array2d: + + for kper in range(array.shape[0]): + + t2d_array_kper = array[kper] + t2d_array_kper_shape = t2d_array_kper.shape + t2d_array_input = t2d_array_kper.reshape(t2d_array_kper_shape[1], + t2d_array_kper_shape[2]) + + vtk.add_array(name, t2d_array_input, array2d=True) + + ot_name = '{}_0{}'.format(name, kper + 1) + ot_file = os.path.join(output_folder, '{}.vtu'.format(ot_name)) + vtk.write(ot_file, timeval=to_tim[kper]) + + else: + + for kper in range(array.shape[0]): + vtk.add_array(name, array[kper]) + + ot_name = '{}_0{}'.format(name, kper + 1) + ot_file = os.path.join(output_folder, '{}.vtu'.format(ot_name)) + vtk.write(ot_file, timeval=to_tim[kper]) + + return + + +def trans_dict(in_dict, name, trans_array, array2d=False): + """ + Builds or adds to dictionary trans_array + """ + if not in_dict: + in_dict = {} + for kper in range(trans_array.shape[0]): + if kper not in in_dict: + in_dict[kper] = {} + in_dict[kper][name] = _Array(trans_array[kper], array2d=array2d) + + return in_dict + + +def export_package(pak_model, pak_name, ot_folder, vtkobj=None, + nanval=-1e+20, smooth=False, point_scalars=False): + + """ + Exports package to vtk + + :param pak_model: the model of the package + :param pak_name: the name of the package + :param ot_folder: the folder to output the package data + :param vtkobj: a vtk object (allows export_package to be called from + export_model) + :param nanval: no data value + :param smooth: If True the output will be a smooth represenation + :param point_scalars: If True the package data will be written as point + values as well as cell values. + """ + + # see if there is vtk object being supplied by export_model + if not vtkobj: + # if not build one + vtk = Vtk(pak_model, nanval=nanval, smooth=smooth, + point_scalars=point_scalars) + else: + # otherwise use the vtk object that was supplied + vtk = vtkobj + + if not os.path.exists(ot_folder): + os.mkdir(ot_folder) + + # is there output data + has_output = False + # is there output transient data + vtk_trans_dict = None + + # get package + if isinstance(pak_name, list): + pak_name = pak_name[0] + + pak = pak_model.get_package(pak_name) + + shape_check_3d = (pak_model.modelgrid.nlay, pak_model.modelgrid.nrow, + pak_model.modelgrid.ncol) + shape_check_2d = (shape_check_3d[1], shape_check_3d[2]) + + # loop through the items in the package + for item, value in pak.__dict__.items(): + + if value is None or not hasattr(value, 'data_type'): + continue + + if isinstance(value, list): + raise NotImplementedError("LIST") + + elif isinstance(value, DataInterface): + # if transiet list data add to the vtk_trans_dict for later output + if value.data_type == DataType.transientlist: + + try: + list(value.masked_4D_arrays_itr()) + except AttributeError: + + continue + except ValueError: + continue + has_output = True + for name, array in value.masked_4D_arrays_itr(): + + vtk_trans_dict = trans_dict(vtk_trans_dict, name, array) + + elif value.data_type == DataType.array3d: + # if 3d array add array to the vtk and set has_output to True + if value.array is not None: + has_output = True + + vtk.add_array(item, value.array) + + elif value.data_type == DataType.array2d and value.array.shape ==\ + shape_check_2d: + # if 2d array add array to vtk object and turn on has output + if value.array is not None: + has_output = True + vtk.add_array(item, value.array, array2d=True) + + elif value.data_type == DataType.transient2d: + # if transient data add data to vtk_trans_dict for later output + if value.array is not None: + has_output = True + vtk_trans_dict = trans_dict(vtk_trans_dict, item, + value.array, array2d=True) + + elif value.data_type == DataType.list: + # this data type is not being output + if value.array is not None: + has_output = True + if isinstance(value.array, np.recarray): + pass + + else: + raise Exception('Data type not understond in data list') + + elif value.data_type == DataType.transient3d: + # add to transient dictionary for output + if value.array is not None: + has_output = True + # vtk_trans_dict = _export_transient_3d(vtk, value.array, + # vtkdict=vtk_trans_dict) + vtk_trans_dict = trans_dict(vtk_trans_dict, item, + value.array) + + else: + pass + else: + pass + + if not has_output: + # there is no data to output + pass + + else: + + # write out data + # write array data + if len(vtk.arrays) > 0: + ot_file = os.path.join(ot_folder, '{}.vtu'.format(pak_name)) + vtk.write(ot_file) + + # write transient data + if vtk_trans_dict: + + # get model time + # to_tim = _get_basic_modeltime(pak_model.modeltime.perlen) + # loop through the transient array data that was stored in the + # trans_array_dict and output + for kper, array_dict in vtk_trans_dict.items(): + # if to_tim: + # time = to_tim[kper] + # else: + # time = None + # set up output file + ot_file = os.path.join(ot_folder, '{} _0{}.vtu'.format( + pak_name, kper + 1)) + for name, array in sorted(array_dict.items()): + if array.array2d: + array_shape = array.array.shape + a = array.array.reshape(array_shape[1], array_shape[2]) + else: + a = array.array + vtk.add_array(name, a, array.array2d) + # vtk.write(ot_file, timeval=time) + vtk.write(ot_file) + return + + +def export_model(model, ot_folder, package_names=None, nanval=-1e+20, + smooth=False, point_scalars=False): + """ + + :param model: flopy model instance + :param ot_folder: output folder + :param package_names: list ofpackage names to be exported + :param nanval: no data value + :param smooth: smoothing + :param point_scalars: If True array data will be written as point scalars + as well as cell scalars + + """ + vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) + if package_names is not None: + if not isinstance(package_names, list): + package_names = [package_names] + else: + package_names = [pak.name[0] for pak in ml.packagelist] -if __name__ == '__main__': - import flopy + if not os.path.exists(ot_folder): + os.mkdir(ot_folder) - ml = flopy.modflow.Modflow() - dis = flopy.modflow.ModflowDis(ml, nlay=3, nrow=3, ncol=3, top=0, - botm=[-1., -2., -3.]) - vtkfile = Vtk('test.vtu', ml) - vtkfile.write() + for pak_name in package_names: + export_package(model, pak_name, ot_folder, vtkobj=vtk, nanval=nanval, + smooth=smooth, point_scalars=point_scalars) From 79b22b595efa8f7408de3f6dc19408b25fac45ec Mon Sep 17 00:00:00 2001 From: martindjm <43825685+martindjm@users.noreply.github.com> Date: Fri, 8 Nov 2019 11:42:27 -0800 Subject: [PATCH 2/5] fix(vtk) fix pep8 errors --- examples/Notebooks/flopy3_vtk_export.ipynb | 63 ++++++---------- flopy/export/vtk.py | 83 +++++++++++----------- 2 files changed, 65 insertions(+), 81 deletions(-) diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb index b2d9c705b5..5a87cfff66 100644 --- a/examples/Notebooks/flopy3_vtk_export.ipynb +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -2,19 +2,9 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "flopy is installed in C:\\Users\\domartin\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\flopy\n", - "3.7.3 (default, Mar 27 2019, 17:13:21) [MSC v.1915 64 bit (AMD64)]\n", - "flopy version: 3.2.12\n" - ] - } - ], + "outputs": [], "source": [ "import os\n", "import sys\n", @@ -38,7 +28,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -50,7 +40,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -62,7 +52,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -73,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -85,13 +75,15 @@ " # this automatically turns on smoothing\n", "# name: the name will be used for the output file and the scalar variable name\n", "ot_arrays_folder = os.path.join(workspace, 'arrays_test')\n", + "if not os.path.exists(ot_arrays_folder):\n", + " os.mkdir(ot_arrays_folder)\n", "model_top_output = os.path.join(ot_arrays_folder, 'TOP_point')\n", "ml.dis.top.export(model_top_output, fmt='vtk', **{'smooth': False, 'point_scalars': True}) " ] }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -107,7 +99,7 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -123,7 +115,7 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -135,7 +127,7 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -144,6 +136,8 @@ "# point_scalars: True creates scalar data values on vertices as well as cells\n", " # this automatically turns on smoothing\n", "package_output = os.path.join(workspace, 'package_output_test')\n", + "if not os.path.exists(package_output):\n", + " os.mkdir(package_output)\n", "# export dis\n", "dis_output = os.path.join(package_output, 'DIS')\n", "ml.dis.export(dis_output, fmt='vtk')\n", @@ -154,22 +148,11 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": { "scrolled": true }, - "outputs": [ - { - "data": { - "text/plain": [ - "'.\\\\VTK_EXAMPLE_OUTPUTS\\\\model_output_test'" - ] - }, - "execution_count": 10, - "metadata": {}, - "output_type": "execute_result" - } - ], + "outputs": [], "source": [ "\n", "# model export\n", @@ -189,7 +172,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -206,7 +189,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -221,7 +204,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -236,7 +219,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -258,7 +241,7 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -269,7 +252,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -278,7 +261,7 @@ }, { "cell_type": "code", - "execution_count": 17, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 2cce473a66..3beb941ec9 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -7,9 +7,7 @@ from flopy.utils import HeadFile import numpy.ma as ma -""" -Module for exporting vtk from flopy -""" +# Module for exporting vtk from flopy def start_tag(f, tag, indent_level, indent_char=' '): @@ -110,7 +108,8 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, self.ibound = ibound # build the model vertices - self.verts, self.iverts, self.zverts = self.get_3d_vertex_connectivity() + self.verts, self.iverts, self.zverts = \ + self.get_3d_vertex_connectivity() return @@ -128,21 +127,12 @@ def add_array(self, name, a, array2d=False): if name == 'ibound': return - # limited recarray support - if type(a) == np.recarray: - matrix = np.zeros(self.shape) - for row in a: - matrix[a.k, a.i, a.j ] = 1 - a = matrix - - else: - - # if array is 2d reformat to 3d array - if array2d: - assert a.shape == self.shape2d - array = np.full(self.shape, self.nanval) - array[0, :, :] = a - a = array + # if array is 2d reformat to 3d array + if array2d: + assert a.shape == self.shape2d + array = np.full(self.shape, self.nanval) + array[0, :, :] = a + a = array try: @@ -150,8 +140,8 @@ def add_array(self, name, a, array2d=False): except AssertionError: return # raise AssertionError - if a.dtype == int: - a = a.astype(float) + # if a.dtype == int: + a = a.astype(float) # add array to self.arrays self.arrays[name] = a return @@ -324,7 +314,8 @@ def _configure_data_arrays(self): # build index array ot_idx_array = np.zeros(shape1d, dtype=np.int) # loop through arrays - for name, array in self.arrays.items(): + for name in self.arrays: + array = self.arrays[name] # make array 1d a = array.ravel() # find where no data @@ -352,12 +343,9 @@ def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): :return: dictionaries of verts, iverts, and zvalues """ # set up active cells - if actwcells is not None: - ncells = int(actwcells.sum()) - else: - ncells = self.nlay * self.nrow * self.ncol + if actwcells is None: actwcells = self.ibound - npoints = ncells * 8 + ipoint = 0 vertsdict = {} @@ -402,7 +390,8 @@ def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): verts.append([pt0[0], pt0[1], elev]) # verts[ipoint+3, :] = np.append(pt3,elev) verts.append([pt3[0], pt3[1], elev]) - ivert.extend([ipoint, ipoint+1, ipoint+2, ipoint+3]) + ivert.extend([ipoint, ipoint+1, ipoint+2, + ipoint+3]) zverts.extend([elev, elev, elev, elev]) ipoint += 4 vertsdict[cellid] = verts @@ -416,11 +405,13 @@ def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): verts.append([pt2[0], pt2[1], zVertices[lay, i+1, j+1]]) - verts.append([pt0[0], pt0[1], zVertices[lay, i, j]]) + verts.append([pt0[0], pt0[1], zVertices[lay, i, + j]]) verts.append([pt3[0], pt3[1], zVertices[lay, i, j+1]]) - ivert.extend([ipoint, ipoint+1, ipoint+2, ipoint+3]) + ivert.extend([ipoint, ipoint+1, ipoint+2, + ipoint+3]) zverts.extend([zVertices[lay, i+1, j], zVertices[ lay, i+1, j+1], zVertices[lay, i, j], zVertices[lay, i, @@ -451,7 +442,8 @@ def extendedDataArray(self, dataArray): for index in indexList: if index[0] in range(self.nrow) and index[1] in \ range(self.ncol): - neighList.append(dataArray[lay, index[0], index[1]]) + neighList.append(dataArray[lay, index[0], + index[1]]) neighList = np.array(neighList) if neighList[neighList > self.nanval].shape[0] > 0: headMean = neighList[neighList > self.nanval].mean() @@ -460,7 +452,8 @@ def extendedDataArray(self, dataArray): matrix[lay, row, col] = headMean return matrix - def write_data_array(self, f, indent_level, arrayName, arrayValues, + @staticmethod + def write_data_array(f, indent_level, arrayName, arrayValues, actWCells): """ @@ -506,9 +499,11 @@ def write_point_value(self, f, indent_level, data_array, array_name, indent_level = start_tag(f, s, indent_level) # data - verts, secus, zverts = self.get_3d_vertex_connectivity( + verts_info = self.get_3d_vertex_connectivity( actwcells=actwcells, zvalues=data_array) + zverts = verts_info[2] + for cellid in sorted(zverts): for z in zverts[cellid]: s = indent_level * ' ' @@ -531,15 +526,20 @@ def _build_iverts(verts): :param verts: vertices being output :return: iverts for output """ - count = 0 + ncells = len(verts) + npoints = ncells * 8 iverts = [] - for a in verts: - ivert = [] - for i in range(len(a)): - ivert.append(count) - count += 1 - iverts.append(ivert) + ivert = [] + count = 1 + for i in range(npoints): + ivert.append(i) + if count == 8: + iverts.append(ivert) + ivert = [] + count = 0 + count += 1 iverts = np.array(iverts) + return iverts @@ -937,7 +937,8 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None, pass else: - raise Exception('Data type not understond in data list') + raise Exception('Data type not understond in data ' + 'list') elif value.data_type == DataType.transient3d: # add to transient dictionary for output From 5c24c8698ba7f4c1df4b6906576664455e050a2a Mon Sep 17 00:00:00 2001 From: martindjm <43825685+martindjm@users.noreply.github.com> Date: Thu, 14 Nov 2019 14:00:10 -0800 Subject: [PATCH 3/5] feat(vtk) added support for binary output --- autotest/t050_test.py | 43 ++- examples/Notebooks/flopy3_vtk_export.ipynb | 341 ++++++++++------- flopy/export/utils.py | 46 ++- flopy/export/vtk.py | 425 +++++++++++++++++++-- 4 files changed, 673 insertions(+), 182 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index a1cbdfb72a..15611f33ef 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -10,6 +10,12 @@ shutil.rmtree(cpth) os.makedirs(cpth) +# binary output directory +binot = os.path.join(cpth, 'bin') +if os.path.isdir(binot): + shutil.rmtree(binot) +os.makedirs(binot) + def test_vtk_export_array2d(): """Export 2d array""" @@ -33,6 +39,9 @@ def test_vtk_export_array3d(): # with point scalars m.upw.hk.export(os.path.join(cpth, 'array_3d_test'), fmt='vtk', name='hk_points', point_scalars=True) + # binary test + m.upw.hk.export(os.path.join(binot, 'array_3d_test'), fmt='vtk', + name='hk_points', point_scalars=True, binary=True) return @@ -45,6 +54,10 @@ def test_vtk_transient_array_2d(): m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False) m.rch.rech.export(os.path.join(cpth, 'transient_2d_test'), fmt='vtk') + # binary test + m.rch.rech.export(os.path.join(binot, 'transient_2d_test'), fmt='vtk', + binary=True) + return @@ -62,6 +75,11 @@ def test_vtk_export_packages(): m.bas6.export(os.path.join(cpth, 'BAS'), fmt='vtk', smooth=True) # transient package drain m.drn.export(os.path.join(cpth, 'DRN'), fmt='vtk') + # binary test + m.dis.export(os.path.join(binot, 'DIS'), fmt='vtk', binary=True) + # upw with point scalar output + m.upw.export(os.path.join(binot, 'UPW'), fmt='vtk', point_scalars=True, + binary=True) return @@ -80,13 +98,16 @@ def test_export_mf2005_vtk(): m = flopy.modflow.Modflow.load(namfile, model_ws=pth, verbose=False) m.export(os.path.join(cpth, m.name), fmt='vtk') + # binary test + m.export(os.path.join(binot, m.name), fmt='vtk', binary=True) + return def test_vtk_mf6(): mf6expth = os.path.join('..', 'examples', 'data', 'mf6') # test vtk mf6 export - mf6sims = ['test005_advgw_tidal', 'test045_lake1ss_table', + mf6sims = ['test045_lake1ss_table', 'test036_twrihfb', 'test045_lake2tr', 'test006_2models_mvr'] # mf6sims = ['test005_advgw_tidal'] # mf6sims = ['test036_twrihfb'] @@ -119,16 +140,25 @@ def test_vtk_bindary_head_export(): verbose=False) otfolder = os.path.join(cpth, 'freyberg_hds_test') - vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, - 1090]) + vtk.export_heads(m, hdsfile, otfolder, nanval=-999.99, kperlist=[1, 200, + 355, + 455, + 1090]) # test with points vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, - 1090], point_scalars=True) + 1090], + point_scalars=True, nanval=-999.99) # test vtk export heads with smoothing and no point scalars vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, 1090], - point_scalars=False, smooth=True) + point_scalars=False, smooth=True, nanval=-999.99) + + # test binary output + vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, + 1090], + point_scalars=False, smooth=True, binary=True, + nanval=-999.99) return @@ -147,6 +177,9 @@ def test_vtk_cbc(): vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'), kperlist=[1, 2, 3], point_scalars=True) + vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'), + kperlist=[1, 2, 3], point_scalars=True, binary=True) + return if __name__ == '__main__': diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb index 5a87cfff66..ed00998a07 100644 --- a/examples/Notebooks/flopy3_vtk_export.ipynb +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -1,10 +1,31 @@ { "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FloPy VTK Export Demo\n", + " This notebook demonstrates how to use FloPy to export to vtk (.vtu) files. This example will cover:\n", + " \n", + " . basic exporting of information for a model, individual package, or array\n", + " . exporting heads and model output data" + ] + }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "flopy is installed in C:\\Users\\domartin\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\flopy\n", + "3.7.3 (default, Mar 27 2019, 17:13:21) [MSC v.1915 64 bit (AMD64)]\n", + "flopy version: 3.2.12\n" + ] + } + ], "source": [ "import os\n", "import sys\n", @@ -28,7 +49,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -40,7 +61,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, "outputs": [], "source": [ @@ -51,253 +72,287 @@ ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exporting to VTK\n", + "\n", + " . for all export calls a folder is provided and the the fmt flag is set to 'vtk'\n" + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# for all export calls a folder is give and the the fmt flag is set to 'vtk'\n", - "# to initiate vtk export. Outputs are writtent to the provided folder.\n", - "# if it does not exist one will be created." + "# Exporting FloPy arrays to .vtu files\n", + "\n", + " All array exports have the following kwargs:\n", + " \n", + " . smooth: True creates a smooth surface, default is False\n", + " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", + " Outputing point scalars will set smooth to True\n", + " . name: A name can be specified to use for the output filename and array scalar name,\n", + " by default the FloPy array name is used\n", + " . binary: True turns on binary vtk export, deault is False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Export model top" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "# export arrays\n", - "# 2d array export\n", - "# export model top\n", - "# **kwargs: smooth: True creates a smooth surface\n", - "# point_scalars: True creates scalar data values on vertices as well as cells\n", - " # this automatically turns on smoothing\n", - "# name: the name will be used for the output file and the scalar variable name\n", + "# create output folder\n", "ot_arrays_folder = os.path.join(workspace, 'arrays_test')\n", "if not os.path.exists(ot_arrays_folder):\n", " os.mkdir(ot_arrays_folder)\n", - "model_top_output = os.path.join(ot_arrays_folder, 'TOP_point')\n", - "ml.dis.top.export(model_top_output, fmt='vtk', **{'smooth': False, 'point_scalars': True}) " + "# export model top\n", + "model_top_output = os.path.join(ot_arrays_folder, 'TOP')\n", + "ml.dis.top.export(model_top_output, fmt='vtk') " ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# transient 2d array\n", - "# export recharge for each stress period\n", - "# **kwargs: smooth: True creates a smooth surface\n", - "# point_scalars: True creates scalar data values on vertices as well as cells\n", - " # this automatically turns on smoothing\n", - "# name: the name will be used for the output file and the scalar variable name\n", - "model_recharge_output = os.path.join(ot_arrays_folder, 'RECH')\n", - "ml.rch.rech.export(model_recharge_output, fmt='vtk')" + "#### Export model bottoms" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# 3D Array export\n", "# export model bottoms\n", - "# **kwargs: smooth: True creates a smooth surface\n", - "# point_scalars: True creates scalar data values on vertices as well as cells\n", - " # this automatically turns on smoothing\n", - "# name: the name will be used for the output file and the scalar variable name\n", "model_bottom_output = os.path.join(ot_arrays_folder, 'BOTM')\n", - "ml.dis.botm.export(model_bottom_output, smooth=True, fmt='vtk')\n" + "ml.dis.botm.export(model_bottom_output, fmt='vtk')" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# 3D Array export\n", - "# export model bottoms\n", - "model_hk_export = os.path.join(ot_arrays_folder, 'HK')\n", - "ml.upw.hk.export(model_hk_export, smooth=True, fmt='vtk', name='HK')" + "#### Export transient array recharge" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, "outputs": [], "source": [ - "# package export\n", - "# **kwargs: smooth: True creates a smooth surface\n", - "# point_scalars: True creates scalar data values on vertices as well as cells\n", - " # this automatically turns on smoothing\n", - "package_output = os.path.join(workspace, 'package_output_test')\n", - "if not os.path.exists(package_output):\n", - " os.mkdir(package_output)\n", - "# export dis\n", - "dis_output = os.path.join(package_output, 'DIS')\n", - "ml.dis.export(dis_output, fmt='vtk')\n", - "#export upw with point scalars\n", - "upw_output = os.path.join(package_output, 'UPW')\n", - "ml.upw.export(upw_output, fmt='vtk', point_scalars=True)" + "# transient 2d array\n", + "# export recharge\n", + "model_recharge_output = os.path.join(ot_arrays_folder, 'RECH')\n", + "ml.rch.rech.export(model_recharge_output, fmt='vtk')" ] }, { - "cell_type": "code", - "execution_count": null, - "metadata": { - "scrolled": true - }, - "outputs": [], + "cell_type": "markdown", + "metadata": {}, "source": [ - "\n", - "# model export\n", - "# **kwargs: smooth: True creates a smooth surface\n", - "# point_scalars: True creates scalar data values on vertices as well as cells\n", - " # this automatically turns on smoothing\n", - "model_output = os.path.join(workspace, 'model_output_test')\n", - "ml.export(model_output, fmt='vtk')\n" + "#### Export HK with point scalars\n" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# 3D Array export\n", + "# hk export, with points\n", + "model_hk_export = os.path.join(ot_arrays_folder, 'HK')\n", + "ml.upw.hk.export(model_hk_export, smooth=True, fmt='vtk', name='HK', point_scalars=True)" + ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# export binary heads\n", - "# args\n", - "# export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None,\n", - " # smooth=False, point_scalars=False, nanval=-1e+20):\n", - "# can limit time steps and stress periods with kstplist and kperlist\n", - "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", - "heads_output_folder = os.path.join(workspace, 'heads_output_test')\n", - "vtk.export_heads(ml, heads_file, heads_output_folder)\n", - "# outputs a .pvd file that can be loaded into paraview to view the heads as time series" + "# Package export to .vtu files\n", + "\n", + " Package export has the following kwargs:\n", + " \n", + " . smooth: True creates a smooth surface, default is False\n", + " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", + " Outputing point scalars will set smooth to True\n", + " . binary: True turns on binary vtk export, default is False" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "# export cbc file\n", - "#export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None,\n", - "# kperlist=None, keylist=None, smooth=False, point_scalars=False,\n", - "# nanval=-1e+20):\n", - "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", - "cbc_output_folder = os.path.join(workspace, 'cbc_output_test')\n", - "vtk.export_cbc(ml, cbc_file, cbc_output_folder)" + "#### Export dis and upw package" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 8, "metadata": {}, "outputs": [], "source": [ - "# export heads with parameters\n", - "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", - "heads_output_folder = os.path.join(workspace, 'heads_output_test_parameters')\n", - "# kstp - time step\n", - "# kperlist - list of stress periods to export\n", - "# point_scalars - output heads as point values\n", - "vtk.export_heads(ml, heads_file, heads_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15], point_scalars=True)" + "# package export\n", + "# set up package export folder\n", + "package_output = os.path.join(workspace, 'package_output_test')\n", + "if not os.path.exists(package_output):\n", + " os.mkdir(package_output)\n", + "# export dis\n", + "dis_output = os.path.join(package_output, 'DIS')\n", + "ml.dis.export(dis_output, fmt='vtk')\n", + "#export upw with point scalars\n", + "upw_output = os.path.join(package_output, 'UPW')\n", + "ml.upw.export(upw_output, fmt='vtk', point_scalars=True)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", - "cbc_output_folder = os.path.join(workspace, 'cbc_output_test_parameters')\n", - "# kstp - time step\n", - "# kperlist - list of stress periods to export\n", - "# keylist - list of cbc terms to export\n", - "vtk.export_cbc(ml, cbc_file, cbc_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15],\n", - " keylist=['CONSTANT HEAD', 'STORAGE'])" + "# Model export to .vtu files\n", + "\n", + " Package export has the following kwargs:\n", + " \n", + " . package_names: list of package names to export\n", + " . smooth: True creates a smooth surface, default is False\n", + " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", + " Outputing point scalars will set smooth to True\n", + " . binary: True turns on binary vtk export, default is False" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "#### Export model as binary .vtu files" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "data": { + "text/plain": [ + "'.\\\\VTK_EXAMPLE_OUTPUTS\\\\model_output_test'" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "\n", + "model_output = os.path.join(workspace, 'model_output_test')\n", + "ml.export(model_output, fmt='vtk', binary=True)\n" + ] + }, + { + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "nam_file = \"mnw1.nam\"\n", - "model_ws = os.path.join(\"..\", \"data\", \"mf2005_test\")\n", - "ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)" + "# Export heads\n", + "\n", + " vtk.export_heads(model, headsfile, output_folder)\n", + "\n", + " Heads export has the following arguments:\n", + "\n", + " . nanval: No data value default is -1e20\n", + " . kstplist: Time steps list\n", + " . kperlist: Stress period list\n", + " . smooth: True creates a smooth surface, default is False\n", + " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", + " Outputing point scalars will set smooth to True\n", + " . binary: True turns on binary vtk export, default is False" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "ml.mnw1.export(os.path.join(workspace, 'mnw1_test'), fmt='vtk')" + "### Export heads to binary .vtu files\n", + "\n", + "#### Exports a .pvd file that can be loaded into Paraview to view heads as time series" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": {}, "outputs": [], "source": [ - "nam_file = \"lakeex3.nam\"\n", - "model_ws = os.path.join(\"..\", \"data\", \"mf2005_test\")\n", - "ml = flopy.modflow.Modflow.load(nam_file, model_ws=model_ws, check=False)\n", - "ml.lak.export(os.path.join(workspace, 'lak_Test'), fmt='vtk')" + "# export binary heads\n", + "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", + "heads_output_folder = os.path.join(workspace, 'heads_output_test')\n", + "vtk.export_heads(ml, heads_file, heads_output_folder, binary=True, nanval=-999.99)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "### Export heads to binary .vtu files with point head scalars" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 11, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# export heads with parameters\n", + "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", + "heads_output_folder = os.path.join(workspace, 'heads_output_test_parameters')\n", + "vtk.export_heads(ml, heads_file, heads_output_folder, kstplist=[1], kperlist=[1, 50, 100, 1000], \n", + " point_scalars=True, nanval=-999.99)" + ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "# Export output cell by cell file to .vtu\n", + "\n", + " vtk.export_heads(model, cellbycellfile, outputfolder)\n", + "\n", + " Heads export has the following arguments:\n", + "\n", + " . precision: Default is single\n", + " . nanval: No data value default is -1e20\n", + " . kstplist: Time steps list\n", + " . kperlist: Stress period list\n", + " . keylist: Flow terms to output\n", + " . smooth: True creates a smooth surface, default is False\n", + " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", + " Outputing point scalars will set smooth to True\n", + " . binary: True turns on binary vtk export, default is False" + ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 13, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", + "cbc_output_folder = os.path.join(workspace, 'cbc_output_test_parameters')\n", + "vtk.export_cbc(ml, cbc_file, cbc_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15],\n", + " keylist=['CONSTANT HEAD', 'STORAGE'], point_scalars=True, binary=True)" + ] }, { "cell_type": "code", diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 672c130fae..6f082911fa 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -374,6 +374,12 @@ def model_export(f, ml, fmt=None, **kwargs): For vtk, if set to True will create point scalar values along with cell values. + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. + + nanval: For vtk, no data value + """ @@ -406,8 +412,11 @@ def model_export(f, ml, fmt=None, **kwargs): # call vtk model export smooth = kwargs.get('smooth', False) point_scalars = kwargs.get('point_scalars', False) + binary = kwargs.get('binary', False) + nanval = kwargs.get('nanval', -1e20) vtk.export_model(ml, f, package_names=package_names, smooth=smooth, - point_scalars=point_scalars) + point_scalars=point_scalars, binary=binary, + nanval=nanval) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -443,6 +452,12 @@ def package_export(f, pak, fmt=None, **kwargs): For vtk, if set to True will create point scalar values along with cell values. + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. + + nanval: For vtk, no data value + Returns ------- f : NetCdf object or None @@ -486,8 +501,11 @@ def package_export(f, pak, fmt=None, **kwargs): # call vtk array export to folder smooth = kwargs.get('smooth', False) point_scalars = kwargs.get('point_scalars', False) + binary = kwargs.get('binary', False) + nanval = kwargs.get('nanval', -1e20) vtk.export_package(pak.parent, pak.name, f, smooth=smooth, - point_scalars=point_scalars) + point_scalars=point_scalars, binary=binary, + nanval=nanval) else: @@ -708,6 +726,9 @@ def transient2d_export(f, t2d, fmt=None, **kwargs): smooth: for vtk to output a smooth represenation of the model point_scalars: for vtk to output point value scalars as well as cell name: for vtk to set a specific name for array and output file + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. """ @@ -811,8 +832,10 @@ def transient2d_export(f, t2d, fmt=None, **kwargs): smooth = kwargs.get('smooth', False) point_scalars = kwargs.get('point_scalars', False) name = kwargs.get('name', t2d.name) + nanval = kwargs.get('nanval', -1e20) vtk.export_transient(t2d.model, t2d.array, f, name, smooth=smooth, - point_scalars=point_scalars, array2d=True) + point_scalars=point_scalars, array2d=True, + nanval=nanval) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -835,6 +858,9 @@ def array3d_export(f, u3d, fmt=None, **kwargs): smooth: for vtk to output a smooth represenation of the model point_scalars: for vtk to output point value scalars as well as cell name: for vtk to set a specific name for array and output file + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. """ @@ -957,12 +983,14 @@ def array3d_export(f, u3d, fmt=None, **kwargs): smooth = kwargs.get('smooth', False) point_scalars = kwargs.get('point_scalars', False) name = kwargs.get('name', u3d.name) - + binary = kwargs.get('binary', False) + nanval = kwargs.get('nanval', -1e20) if isinstance(name, list) or isinstance(name, tuple): name = name[0] vtk.export_array(u3d.model, u3d.array, f, name, smooth=smooth, - point_scalars=point_scalars) + point_scalars=point_scalars, binary=binary, + nanval=nanval) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -986,6 +1014,9 @@ def array2d_export(f, u2d, fmt=None, **kwargs): smooth: for vtk to output a smooth represenation of the model point_scalars: for vtk to output point value scalars as well as cell name: for vtk to set a specific name for array and output file + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. """ assert isinstance(u2d, DataInterface), "util2d_helper only helps " \ @@ -1084,8 +1115,11 @@ def array2d_export(f, u2d, fmt=None, **kwargs): smooth = kwargs.get('smooth', False) point_scalars = kwargs.get('point_scalars', False) name = kwargs.get('name', u2d.name) + binary = kwargs.get('binary', False) + nanval = kwargs.get('nanval', -1e20) vtk.export_array(u2d.model, u2d.array, f, name, smooth=smooth, - point_scalars=point_scalars, array2d=True) + point_scalars=point_scalars, array2d=True, + binary=binary, nanval=nanval) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 3beb941ec9..9609db54bb 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -6,9 +6,120 @@ import flopy.utils.binaryfile as bf from flopy.utils import HeadFile import numpy.ma as ma +import struct +import sys # Module for exporting vtk from flopy +# BINARY ********************************************* +np_to_struct = {'int8': 'b', + 'uint8': 'B', + 'int16': 'h', + 'uint16': 'H', + 'int32': 'i', + 'uint32': 'I', + 'int64': 'q', + 'uint64': 'Q', + 'float32': 'f', + 'float64': 'd'} + + +class BinaryXml: + def __init__(self, file_path): + self.stream = open(file_path, "wb") + self.open_tag = False + self.current = [] + self.stream.write(b'') + if sys.byteorder == "little": + self.byte_order = '<' + else: + self.byte_order = '>' + + def write_size(self, block_size): + # size is a 64 bit unsigned integer + byte_order = self.byte_order + 'Q' + block_size = struct.pack(byte_order, block_size) + self.stream.write(block_size) + + def write_array(self, data): + assert (data.flags['C_CONTIGUOUS'] or data.flags['F_CONTIGUOUS']) + + # ravel in fortran order + dd = np.ravel(data, order='F') + + data_format = self.byte_order + str(data.size) + \ + np_to_struct[data.dtype.name] + binary_data = struct.pack(data_format, *dd) + self.stream.write(binary_data) + + def write_coord_arrays(self, x, y, z): + # check that arrays are the same shape and data type + assert (x.size == y.size == z.size) + assert (x.dtype.itemsize == y.dtype.itemsize == z.dtype.itemsize) + + # check if arrays are contiguous + assert (x.flags['C_CONTIGUOUS'] or x.flags['F_CONTIGUOUS']) + assert (y.flags['C_CONTIGUOUS'] or y.flags['F_CONTIGUOUS']) + assert (z.flags['C_CONTIGUOUS'] or z.flags['F_CONTIGUOUS']) + + data_format = self.byte_order + str(1) + \ + np_to_struct[x.dtype.name] + + xrav = np.ravel(x, order='F') + yrav = np.ravel(y, order='F') + zrav = np.ravel(z, order='F') + + for idx in range(x.size): + bx = struct.pack(data_format, xrav[idx]) + by = struct.pack(data_format, yrav[idx]) + bz = struct.pack(data_format, zrav[idx]) + self.stream.write(bx) + self.stream.write(by) + self.stream.write(bz) + + def close(self): + assert (not self.open_tag) + self.stream.close() + + def open_element(self, tag): + if self.open_tag: + self.stream.write(b">") + tag_string = "\n<%s" % tag + self.stream.write(str.encode(tag_string)) + self.open_tag = True + self.current.append(tag) + return self + + def close_element(self, tag=None): + if tag: + assert (self.current.pop() == tag) + if self.open_tag: + self.stream.write(b">") + self.open_tag = False + string = "\n" % tag + self.stream.write(str.encode(string)) + else: + self.stream.write(b"/>") + self.open_tag = False + self.current.pop() + return self + + def add_text(self, text): + if self.open_tag: + self.stream.write(b">\n") + self.open_tag = False + self.stream.write(str.encode(text)) + return self + + def add_attributes(self, **kwargs): + assert self.open_tag + for key in kwargs: + st = ' %s="%s"' % (key, kwargs[key]) + self.stream.write(str.encode(st)) + return self + +# END BINARY ********************************************* + def start_tag(f, tag, indent_level, indent_char=' '): # starts xml tag @@ -87,6 +198,7 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, # check if structured grid, vtk only supports structured grid assert (isinstance(self.modelgrid, StructuredGrid)) + # cbd self.cbd_on = False # get ibound @@ -95,7 +207,7 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, ibound = np.ones(self.shape) else: ibound = self.modelgrid.idomain - + # build cbd ibound if ibound is not None and hasattr(self.model, 'dis') and \ hasattr(self.model.dis, 'laycbd'): @@ -139,9 +251,11 @@ def add_array(self, name, a, array2d=False): assert a.shape == self.shape except AssertionError: return - # raise AssertionError - # if a.dtype == int: + + a = np.where(a == self.nanval, np.nan, a) a = a.astype(float) + # idxs = np.argwhere(a == self.nanval) + # add array to self.arrays self.arrays[name] = a return @@ -301,6 +415,217 @@ def write(self, output_file, timeval=None): self.arrays.clear() return + def write_binary(self, output_file): + + """ + + :param output_file: vtk output file + :return: outputs binary .vtu file + + # timeval not yet supported in binary + + """ + + # make sure file ends with vtu + assert output_file.lower().endswith(".vtu") + + if self.verbose: + print('writing binary vtk file') + + xml = BinaryXml(output_file) + offset = 0 + grid_type = 'UnstructuredGrid' + + # get the active data cells based on the data arrays and ibound + actwcells3d = self._configure_data_arrays() + actwcells = actwcells3d.ravel() + + # get the indexes of the active cells + idxs = np.argwhere(actwcells != 0).ravel() + + # get the verts and iverts to be output + verts = [self.verts[idx] for idx in idxs] + iverts = self._build_iverts(verts) + + # check if there is data to be written out + if len(verts) == 0: + # if not cannot write binary .vtu file + return + + # get the total number of cells and vertices + ncells = len(iverts) + npoints = ncells * 8 + + if self.verbose: + print('Writing vtk file: ' + output_file) + print('Number of point is {}, Number of cells is {}\n'.format( + npoints, ncells)) + + # format verts and iverts + verts = np.array(verts) + verts.reshape(npoints, 3) + iverts = np.ascontiguousarray(iverts, np.float64) + + # write xml file info + xml.open_element("VTKFile"). \ + add_attributes(type=grid_type, version="1.0", + byte_order=self._get_byte_order(), + header_type="UInt64") + # unstructured grid + xml.open_element(grid_type) + + # piece + xml.open_element('Piece') + xml.add_attributes(NumberOfPoints=npoints, NumberOfCells=ncells) + + # points + xml.open_element('Points') + + xml.open_element('DataArray') + xml.add_attributes(Name='points', NumberOfComponents='3', + type='Float64', + format='appended', offset=offset) + + # calculate the offset of the start of the next piece of data + # offset is calculated from beginning of data section + points_size = verts.size * verts[0].dtype.itemsize + offset += points_size + 8 + + xml.close_element('DataArray') + + xml.close_element('Points') + + # cells + xml.open_element('Cells') + + # connectivity + xml.open_element('DataArray') + xml.add_attributes(Name='connectivity', NumberOfComponents='1', + type='Float64', + format='appended', offset=offset) + conn_size = iverts.size * iverts[0].dtype.itemsize + offset += conn_size + 8 + + xml.close_element('DataArray') + + xml.open_element('DataArray') + xml.add_attributes(Name='offsets', NumberOfComponents='1', + type='Float64', + format='appended', offset=offset) + offsets_size = iverts.shape[0] * iverts[0].dtype.itemsize + offset += offsets_size + 8 + + xml.close_element('DataArray') + + xml.open_element('DataArray') + xml.add_attributes(Name='types', NumberOfComponents='1', + type='Float64', + format='appended', offset=offset) + types_size = iverts.shape[0] * iverts[0].dtype.itemsize + offset += types_size + 8 + + xml.close_element('DataArray') + + xml.close_element('Cells') + + xml.open_element('CellData') + xml.add_attributes(Scalars='scalars') + + # format data arrays and store for later output + processed_arrays = [] + for name, a in self.arrays.items(): + a = a.ravel()[idxs] + xml.open_element('DataArray') + xml.add_attributes(Name=name, NumberOfComponents='1', + type='Float64', + format='appended', offset=offset) + a = np.ascontiguousarray(a, np.float64) + processed_arrays.append([a, a.size * a[0].dtype.itemsize]) + offset += processed_arrays[-1][-1] + 8 + xml.close_element('DataArray') + + xml.close_element('CellData') + + # for data array point scalars + if self.point_scalars: + + xml.open_element('PointData') + xml.add_attributes(Scalars='scalars') + + # get output point arrays + # loop through stored arrays + for name, a in self.arrays.items(): + # get the array values onto vertices + verts_info = self.get_3d_vertex_connectivity( + actwcells=actwcells3d, zvalues=a) + # get values + point_values_dict = verts_info[2] + a = np.array([point_values_dict[cellid] for cellid in + sorted(point_values_dict.keys())]).ravel() + + xml.open_element('DataArray') + xml.add_attributes(Name=name, NumberOfComponents='1', + type='Float64', + format='appended', offset=offset) + a = np.ascontiguousarray(a, np.float64) + processed_arrays.append([a, a.size * a[0].dtype.itemsize]) + offset += processed_arrays[-1][-1] + 8 + + xml.close_element('DataArray') + xml.close_element('PointData') + + # end piece + xml.close_element('Piece') + + # end unstructured grid + xml.close_element('UnstructuredGrid') + + # build data section + xml.open_element("AppendedData").add_attributes( + encoding="raw").add_text("_") + + xml.write_size(points_size) + # format verts for output + verts_x = np.ascontiguousarray(np.ravel(verts[:, :, 0]), + np.float64) + verts_y = np.ascontiguousarray(np.ravel(verts[:, :, 1]), + np.float64) + verts_z = np.ascontiguousarray(np.ravel(verts[:, :, 2]), + np.float64) + # write coordinates + xml.write_coord_arrays(verts_x, verts_y, verts_z) + + # write iverts + xml.write_size(conn_size) + rav_iverts = np.ascontiguousarray(np.ravel(iverts), np.float64) + xml.write_array(rav_iverts) + + xml.write_size(offsets_size) + data = np.empty((iverts.shape[0]), np.float64) + icount = 0 + for index, row in enumerate(iverts): + icount += len(row) + data[index] = icount + xml.write_array(data) + + # write cell types (11) + xml.write_size(types_size) + data = np.empty((iverts.shape[0]), np.float64) + data.fill(self.cell_type) + xml.write_array(data) + + # write out the array scalars and array point scalars + for a, block_size in processed_arrays: + xml.write_size(block_size) + xml.write_array(a) + + # end xml + xml.close_element("AppendedData") + xml.close_element('VTKFile') + xml.close() + # clear arrays + self.arrays.clear() + def _configure_data_arrays(self): """ Compares all the stored arrays int the vtk class along with the @@ -423,6 +748,7 @@ def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): return vertsdict, ivertsdict, zvertsdict def extendedDataArray(self, dataArray): + if dataArray.shape[0] == self.nlay+1: dataArray = dataArray else: @@ -445,13 +771,20 @@ def extendedDataArray(self, dataArray): neighList.append(dataArray[lay, index[0], index[1]]) neighList = np.array(neighList) - if neighList[neighList > self.nanval].shape[0] > 0: - headMean = neighList[neighList > self.nanval].mean() + if neighList[neighList != self.nanval].shape[0] > 0: + headMean = neighList[neighList != self.nanval].mean() else: headMean = self.nanval matrix[lay, row, col] = headMean return matrix + @staticmethod + def _get_byte_order(): + if sys.byteorder == "little": + return "LittleEndian" + else: + return "BigEndian" + @staticmethod def write_data_array(f, indent_level, arrayName, arrayValues, actWCells): @@ -553,9 +886,9 @@ def _get_names(in_list): return ot_list -def export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None, - kperlist=None, keylist=None, smooth=False, point_scalars=False, - nanval=-1e+20): +def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, + kstplist=None, kperlist=None, keylist=None, smooth=False, + point_scalars=False, binary=False): """ Exports cell by cell file to vtk @@ -564,12 +897,14 @@ def export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None, :param cbcfile: the cell by cell file :param otfolder: output folder to write the data to :param precision: bindary file precision + :param nanval: the no data value :param kstplist: list of timesteps to be written :param kperlist: list of stress periods to be writeen :param keylist: list of flow term names to be output :param smooth: If true a smooth surface will be output :param point_scalars: If True point scalar values will be written - :param nanval: the no data value + :param binary: if True the output .vtu file will be binary, defualt is + False. """ @@ -669,7 +1004,10 @@ def export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None, vtk.add_array(name.strip(), array) # need to adjust for # write the vtk data to the output file - vtk.write(otfile) + if binary: + vtk.write_binary(otfile) + else: + vtk.write(otfile) count += 1 # finish writing the pvd file pvdfile.write(""" @@ -679,22 +1017,26 @@ def export_cbc(model, cbcfile, otfolder, precision='single', kstplist=None, return -def export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None, - smooth=False, point_scalars=False, nanval=-1e+20): +def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstplist=None, + kperlist=None, smooth=False, point_scalars=False, + binary=False): """ Exports heads to vtk files by timestep and stressperiod :param model: the model instance :param hdsfile: the binary heads file :param otfolder: the output folder to write the .vtu files + :param nanval: The no data value :param kstplist: list of time steps to output :param kperlist: list of stress periods to output :param smooth: If set to True a smooth surface will be output :param point_scalars: If set to True the heads will be written to the vertices as point scalars as well as cell values - :param nanval: The no data value + :param binary: if True the output .vtu file will be binary, defualt is + False. :return: Heads will be written to files named by stress period and timestep to the otfolder + """ # setup output folder @@ -740,7 +1082,10 @@ def export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None, model.name, kper + 1, kstp + 1) otfile = os.path.join(otfolder, ot_base) # vtk.write(otfile, timeval=totim_dict[(kstp, kper)]) - vtk.write(otfile) + if binary: + vtk.write_binary(otfile) + else: + vtk.write(otfile) pvdfile.write("""\n""".format(count, ot_base)) count += 1 @@ -752,16 +1097,21 @@ def export_heads(model, hdsfile, otfolder, kstplist=None, kperlist=None, def export_array(model, array, output_folder, name, nanval=-1e+20, - array2d=False, smooth=False, point_scalars=False): + array2d=False, smooth=False, point_scalars=False, + binary=False): """ :param model: array the model belongs to :param array: array to be exported - :param arrayname: name of the array to be used in vtk file + :param name: name of the array to be used in vtk file :param output_folder: output folder to store the output .vtu file :param array2d: True if array is 2d - :nanval nan value + :param nanval nan value + :param smooth: If set to True a smooth surface will be output + :param point_scalars: If set to True the heads will be written to the + vertices as point scalars as well as cell values + :param binary: if Ture vtk will export as binary :return: outputs a .vtu file of array """ @@ -772,22 +1122,29 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) vtk.add_array(name, array, array2d=array2d) otfile = os.path.join(output_folder, '{}.vtu'.format(name)) - vtk.write(otfile) + if binary: + vtk.write_binary(otfile) + else: + vtk.write(otfile) return def export_transient(model, array, output_folder, name, nanval=-1e+20, - array2d=False, smooth=False, point_scalars=False): + array2d=False, smooth=False, point_scalars=False, + binary=False): """ - :param model: model of transient array :param array: transient array to export :param name: name of the data + :param nanval: nan value :param output_folder: output folder location :param array2d: True if array is 2d - :param nanval: nan value + :param smooth: If set to True a smooth surface will be output + :param point_scalars: If set to True the heads will be written to the + vertices as point scalars as well as cell values + :param binary: if Ture vtk will export as binary :return: ouputs .vtu files of transient data to output folder """ @@ -820,8 +1177,10 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20, ot_name = '{}_0{}'.format(name, kper + 1) ot_file = os.path.join(output_folder, '{}.vtu'.format(ot_name)) - vtk.write(ot_file, timeval=to_tim[kper]) - + if binary: + vtk.write_binary(ot_file) + else: + vtk.write(ot_file, timeval=to_tim[kper]) return @@ -840,7 +1199,8 @@ def trans_dict(in_dict, name, trans_array, array2d=False): def export_package(pak_model, pak_name, ot_folder, vtkobj=None, - nanval=-1e+20, smooth=False, point_scalars=False): + nanval=-1e+20, smooth=False, point_scalars=False, + binary=False): """ Exports package to vtk @@ -854,6 +1214,7 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None, :param smooth: If True the output will be a smooth represenation :param point_scalars: If True the package data will be written as point values as well as cell values. + :param binary: if Ture vtk will export as binary """ # see if there is vtk object being supplied by export_model @@ -964,7 +1325,10 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None, # write array data if len(vtk.arrays) > 0: ot_file = os.path.join(ot_folder, '{}.vtu'.format(pak_name)) - vtk.write(ot_file) + if binary: + vtk.write_binary(ot_file) + else: + vtk.write(ot_file) # write transient data if vtk_trans_dict: @@ -989,12 +1353,15 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None, a = array.array vtk.add_array(name, a, array.array2d) # vtk.write(ot_file, timeval=time) - vtk.write(ot_file) + if binary: + vtk.write_binary(ot_file) + else: + vtk.write(ot_file) return def export_model(model, ot_folder, package_names=None, nanval=-1e+20, - smooth=False, point_scalars=False): + smooth=False, point_scalars=False, binary=False): """ :param model: flopy model instance @@ -1003,6 +1370,7 @@ def export_model(model, ot_folder, package_names=None, nanval=-1e+20, :param nanval: no data value :param smooth: smoothing :param point_scalars: If True array data will be written as point scalars + :param binary: if Ture vtk will export as binary as well as cell scalars """ @@ -1019,4 +1387,5 @@ def export_model(model, ot_folder, package_names=None, nanval=-1e+20, for pak_name in package_names: export_package(model, pak_name, ot_folder, vtkobj=vtk, nanval=nanval, - smooth=smooth, point_scalars=point_scalars) + smooth=smooth, point_scalars=point_scalars, + binary=binary) From 257e4efe125f82fc011149f7e311c9cb789a878b Mon Sep 17 00:00:00 2001 From: martindjm <43825685+martindjm@users.noreply.github.com> Date: Thu, 14 Nov 2019 18:39:26 -0800 Subject: [PATCH 4/5] fix(vtk): updated docstrings to flopy format --- examples/Notebooks/flopy3_vtk_export.ipynb | 2 +- flopy/export/utils.py | 2 +- flopy/export/vtk.py | 357 +++++++++++++++------ 3 files changed, 253 insertions(+), 108 deletions(-) diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb index ed00998a07..04e5f95094 100644 --- a/examples/Notebooks/flopy3_vtk_export.ipynb +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -77,7 +77,7 @@ "source": [ "# Exporting to VTK\n", "\n", - " . for all export calls a folder is provided and the the fmt flag is set to 'vtk'\n" + " for all export calls a folder is provided and the the fmt flag is set to 'vtk'\n" ] }, { diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 6f082911fa..f3a5c15f8d 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -727,7 +727,7 @@ def transient2d_export(f, t2d, fmt=None, **kwargs): point_scalars: for vtk to output point value scalars as well as cell name: for vtk to set a specific name for array and output file binary: bool - For vtk, if set to Ture will output binary .vtu files. Default + For vtk, if set to True will output binary .vtu files. Default is False which exports standard text xml .vtu files. """ diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 9609db54bb..6e288588ae 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -25,6 +25,16 @@ class BinaryXml: + """ + Helps write binary vtk files + + Parameters + ---------- + + file_path : str + output file path + + """ def __init__(self, file_path): self.stream = open(file_path, "wb") self.open_tag = False @@ -47,8 +57,8 @@ def write_array(self, data): # ravel in fortran order dd = np.ravel(data, order='F') - data_format = self.byte_order + str(data.size) + \ - np_to_struct[data.dtype.name] + data_format = self.byte_order + str(data.size) + np_to_struct[ + data.dtype.name] binary_data = struct.pack(data_format, *dd) self.stream.write(binary_data) @@ -154,21 +164,34 @@ def _get_basic_modeltime(perlen_list): class Vtk(object): + """ + Class to build VTK object for exporting flopy vtk + + Parameters + ---------- + + model : MFModel + flopy model instance + verbose : bool + if True, stdout is verbose + nanval : float + no data value, default is -1e20 + smooth : bool + If True will create smooth output surface + point_scalars : bool + if True will output point scalar values, this will set smooth to True. + + Attributes + ---------- + + arrays : dict + Stores data arrays added to VTK object + + """ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, point_scalars=False): - """ - Make Vtk object for exporting flopy vtk - - :param model: flopy model instance - :param verbose: if True, stdout is verbose - :param nanval: no data value - :param smooth: if True will create a smooth output surface - :param point_scalars: if True will output poin scalar values, - this will set smooth to True. - """ - if point_scalars: smooth = True @@ -230,9 +253,16 @@ def add_array(self, name, a, array2d=False): """ Adds an array to the vtk object - :param name: name of the array - :param a: the array to be added to the vtk object - :param array2d: True if the array is 2d + + Parameters + ---------- + + name : str + name of the array + a : flopy array + the array to be added to the vtk object + array2d : bool + True if the array is 2d """ @@ -262,11 +292,18 @@ def add_array(self, name, a, array2d=False): def write(self, output_file, timeval=None): """ - writes the arrays from self.arrays to vtk file - :param output_file: output file name to write the vtk data - :param timeval: model time value to be stored in the time section of - the vtk file + writes the stored arrays to vtk file + + Parameters + ---------- + + output_file : str + output file name to write the vtk data + + timeval : scalar + model time value to be stored in the time section of the vtk + file, default is None """ # make sure file ends with vtu @@ -419,10 +456,13 @@ def write_binary(self, output_file): """ - :param output_file: vtk output file - :return: outputs binary .vtu file + outputs binary .vtu file - # timeval not yet supported in binary + Parameters + ---------- + + output_file : str + vtk output file """ @@ -628,9 +668,8 @@ def write_binary(self, output_file): def _configure_data_arrays(self): """ - Compares all the stored arrays int the vtk class along with the - ibound to figure out what cells and points need to be written to the - output vtk file. + Compares arrays and active cells to find where active data + exists, and what cells to output. """ # get 1d shape @@ -660,12 +699,28 @@ def _configure_data_arrays(self): return ot_idx_array def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None): + """ - :param actwcells: array of where data exists - :param zvalues: array of values to be used instead of the zvalues of - the vertices. This allows point scalars to be interpolated. - :return: dictionaries of verts, iverts, and zvalues + Builds x,y,z vertices + + Parameters + ---------- + actwcells : array + array of where data exists + zvalues: array + array of values to be used instead of the zvalues of + the vertices. This allows point scalars to be interpolated. + + Returns + ------- + vertsdict : dict + dictionary of verts + ivertsdict : dict + dictionary of iverts + zvertsdict : dict + dictionary of zverts + """ # set up active cells if actwcells is None: @@ -792,11 +847,18 @@ def write_data_array(f, indent_level, arrayName, arrayValues, Writes the data array to the output vtk file - :param f: output vtk file - :param indent_level: current indent of the xml - :param arrayName: name of the output array - :param arrayValues: the data array being output - :param actWCells: array of the active cells + Parameters + ---------- + f : file object + output vtk file + indent_level : int + current indent of the xml + arrayName : str + name of the output array + arrayValues : array + the data array being output + actWCells : array + array of the active cells """ @@ -856,8 +918,17 @@ def _build_iverts(verts): Builds the iverts based on the vertices being output - :param verts: vertices being output - :return: iverts for output + Parameters + ---------- + verts : array + vertices being output + + Returns + ------- + + iverts : array + array of ivert values + """ ncells = len(verts) npoints = ncells * 8 @@ -893,18 +964,32 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, Exports cell by cell file to vtk - :param model: the flopy model instance - :param cbcfile: the cell by cell file - :param otfolder: output folder to write the data to - :param precision: bindary file precision - :param nanval: the no data value - :param kstplist: list of timesteps to be written - :param kperlist: list of stress periods to be writeen - :param keylist: list of flow term names to be output - :param smooth: If true a smooth surface will be output - :param point_scalars: If True point scalar values will be written - :param binary: if True the output .vtu file will be binary, defualt is - False. + Parameters + ---------- + + model : flopy model instance + the flopy model instance + cbcfile : str + the cell by cell file + otfolder : str + output folder to write the data + precision : str: + binary file precision, default is 'single' + nanval : scalar + no data value + kstplist : list + list of timesteps + kperlist : list + list of stress periods + keylist : list + list of flow term names + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. """ @@ -991,8 +1076,10 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, addarray = True else: - raise Exception('Data type not currenlty supported ' + raise Exception('Data type not currently supported ' 'for cbc output') + # print('Data type not currently supported ' + # 'for cbc output') if addarray: @@ -1021,21 +1108,31 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstplist=None, kperlist=None, smooth=False, point_scalars=False, binary=False): """ - Exports heads to vtk files by timestep and stressperiod - - :param model: the model instance - :param hdsfile: the binary heads file - :param otfolder: the output folder to write the .vtu files - :param nanval: The no data value - :param kstplist: list of time steps to output - :param kperlist: list of stress periods to output - :param smooth: If set to True a smooth surface will be output - :param point_scalars: If set to True the heads will be written to the - vertices as point scalars as well as cell values - :param binary: if True the output .vtu file will be binary, defualt is - False. - :return: Heads will be written to files named by stress period and - timestep to the otfolder + + Exports binary head file to vtk + + Parameters + ---------- + + model : MFModel + the flopy model instance + hdsfile : str + binary heads file + otfolder : str + output folder to write the data + nanval : scalar + no data value, default value is -1e20 + kstplist : list + list of timesteps + kperlist : list + list of stress periods + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. """ @@ -1102,17 +1199,30 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, """ - :param model: array the model belongs to - :param array: array to be exported - :param name: name of the array to be used in vtk file - :param output_folder: output folder to store the output .vtu file - :param array2d: True if array is 2d - :param nanval nan value - :param smooth: If set to True a smooth surface will be output - :param point_scalars: If set to True the heads will be written to the - vertices as point scalars as well as cell values - :param binary: if Ture vtk will export as binary - :return: outputs a .vtu file of array + Export array to vtk + + Parameters + ---------- + + model : flopy model instance + the flopy model instance + array : flopy array + flopy 2d or 3d array + output_folder : str + output folder to write the data + name : str + name of array + nanval : scalar + no data value, default value is -1e20 + array2d : bool + True if array is 2d, default is False + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. """ @@ -1133,19 +1243,33 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, def export_transient(model, array, output_folder, name, nanval=-1e+20, array2d=False, smooth=False, point_scalars=False, binary=False): - """ - :param model: model of transient array - :param array: transient array to export - :param name: name of the data - :param nanval: nan value - :param output_folder: output folder location - :param array2d: True if array is 2d - :param smooth: If set to True a smooth surface will be output - :param point_scalars: If set to True the heads will be written to the - vertices as point scalars as well as cell values - :param binary: if Ture vtk will export as binary - :return: ouputs .vtu files of transient data to output folder + + Export transient 2d array to vtk + + Parameters + ---------- + + model : MFModel + the flopy model instance + array : Transient instance + flopy transient array + output_folder : str + output folder to write the data + name : str + name of array + nanval : scalar + no data value, default value is -1e20 + array2d : bool + True if array is 2d, default is False + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. + """ if not os.path.exists(output_folder): @@ -1205,16 +1329,28 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None, """ Exports package to vtk - :param pak_model: the model of the package - :param pak_name: the name of the package - :param ot_folder: the folder to output the package data - :param vtkobj: a vtk object (allows export_package to be called from - export_model) - :param nanval: no data value - :param smooth: If True the output will be a smooth represenation - :param point_scalars: If True the package data will be written as point - values as well as cell values. - :param binary: if Ture vtk will export as binary + Parameters + ---------- + + pak_model : flopy model instance + the model of the package + pak_name : str + the name of the package + ot_folder : str + output folder to write the data + vtkobj : VTK instance + a vtk object (allows export_package to be called from + export_model) + nanval : scalar + no data value, default value is -1e20 + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. + """ # see if there is vtk object being supplied by export_model @@ -1364,14 +1500,23 @@ def export_model(model, ot_folder, package_names=None, nanval=-1e+20, smooth=False, point_scalars=False, binary=False): """ - :param model: flopy model instance - :param ot_folder: output folder - :param package_names: list ofpackage names to be exported - :param nanval: no data value - :param smooth: smoothing - :param point_scalars: If True array data will be written as point scalars - :param binary: if Ture vtk will export as binary - as well as cell scalars + model : flopy model instance + flopy model + ot_folder : str + output folder + package_names : list + list of package names to be exported + nanval : scalar + no data value, default value is -1e20 + array2d : bool + True if array is 2d, default is False + smooth : bool + If true a smooth surface will be output, default is False + point_scalars : bool + If True point scalar values will be written, default is False + binary : bool + if True the output .vtu file will be binary, default is + False. """ vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) From 5d235d3ac03669a07179cdacb2de1fb10bfcddaf Mon Sep 17 00:00:00 2001 From: martindjm <43825685+martindjm@users.noreply.github.com> Date: Fri, 15 Nov 2019 12:29:11 -0800 Subject: [PATCH 5/5] fix(vtk): updated kstpkper and text parameter inputs for export functions to match flopy format --- autotest/t050_test.py | 59 +++++++++++----- examples/Notebooks/flopy3_vtk_export.ipynb | 21 +++--- flopy/export/vtk.py | 82 +++++++++++++--------- 3 files changed, 101 insertions(+), 61 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 15611f33ef..02eb9fc684 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -127,7 +127,7 @@ def test_vtk_mf6(): return -def test_vtk_bindary_head_export(): +def test_vtk_binary_head_export(): """test vet export of heads""" @@ -138,50 +138,71 @@ def test_vtk_bindary_head_export(): m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_pth, verbose=False) - otfolder = os.path.join(cpth, 'freyberg_hds_test') - - vtk.export_heads(m, hdsfile, otfolder, nanval=-999.99, kperlist=[1, 200, - 355, - 455, - 1090]) + otfolder = os.path.join(cpth, 'heads_test') + + vtk.export_heads(m, hdsfile, otfolder, nanval=-999.99, kstpkper=[(0, 0), + (0, 199), + (0, 354), + (0, 454), + (0, + 1089)]) # test with points - vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, - 1090], + otfolder = os.path.join(cpth, 'heads_test_1') + vtk.export_heads(m, hdsfile, otfolder, + kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, + 1089)], point_scalars=True, nanval=-999.99) # test vtk export heads with smoothing and no point scalars - vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, - 1090], + otfolder = os.path.join(cpth, 'heads_test_2') + vtk.export_heads(m, hdsfile, otfolder, + kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, + 1089)], point_scalars=False, smooth=True, nanval=-999.99) # test binary output - vtk.export_heads(m, hdsfile, otfolder, kperlist=[1, 200, 355, 455, - 1090], + otfolder = os.path.join(cpth, 'heads_test_3') + vtk.export_heads(m, hdsfile, otfolder, + kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, + 1089)], point_scalars=False, smooth=True, binary=True, nanval=-999.99) + otfolder = os.path.join(cpth, 'heads_test_4') + vtk.export_heads(m, hdsfile, otfolder, kstpkper=(0, 0), + point_scalars=False, + smooth=True, binary=True, nanval=-999.99) + return def test_vtk_cbc(): # test mf 2005 freyberg freyberg_cbc = os.path.join('..', 'examples', 'data', - 'freyberg_multilayer_transient', 'freyberg.cbc') + 'freyberg_multilayer_transient', + 'freyberg.cbc') freyberg_mpth = os.path.join('..', 'examples', 'data', - 'freyberg_multilayer_transient') + 'freyberg_multilayer_transient') m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_mpth, verbose=False) vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'), - kperlist=[1, 2, 3], point_scalars=True) + kstpkper=[(0, 0), (0, 1), (0, 2)], point_scalars=True) - vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'), - kperlist=[1, 2, 3], point_scalars=True, binary=True) + vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST_bin'), + kstpkper=[(0, 0), (0, 1), (0, 2)], point_scalars=True, + binary=True) + + vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, + 'freyberg_CBCTEST_bin2'), + kstpkper=(0, 0), text='CONSTANT HEAD', + point_scalars=True, binary=True) return + if __name__ == '__main__': test_vtk_export_array2d() test_vtk_export_array3d() @@ -189,5 +210,5 @@ def test_vtk_cbc(): test_vtk_export_packages() test_export_mf2005_vtk() test_vtk_mf6() - test_vtk_bindary_head_export() + test_vtk_binary_head_export() test_vtk_cbc() diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb index 04e5f95094..f8e125bfbc 100644 --- a/examples/Notebooks/flopy3_vtk_export.ipynb +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -272,8 +272,8 @@ " Heads export has the following arguments:\n", "\n", " . nanval: No data value default is -1e20\n", - " . kstplist: Time steps list\n", - " . kperlist: Stress period list\n", + " . kstpkper: A tuple containing the time step and stress period (kstp, kper) or\n", + " list of (kstp, kper) tuples. The kstp and kper values are zero based.\n", " . smooth: True creates a smooth surface, default is False\n", " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", " Outputing point scalars will set smooth to True\n", @@ -291,7 +291,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 10, "metadata": {}, "outputs": [], "source": [ @@ -317,7 +317,8 @@ "# export heads with parameters\n", "heads_file = os.path.join(model_ws, 'freyberg.hds')\n", "heads_output_folder = os.path.join(workspace, 'heads_output_test_parameters')\n", - "vtk.export_heads(ml, heads_file, heads_output_folder, kstplist=[1], kperlist=[1, 50, 100, 1000], \n", + "# export heads for time step 1, stress periods 1, 50, 100, 1000\n", + "vtk.export_heads(ml, heads_file, heads_output_folder, kstpkper=[(0,0), (0, 49), (0, 99), (0, 999)],\n", " point_scalars=True, nanval=-999.99)" ] }, @@ -333,9 +334,9 @@ "\n", " . precision: Default is single\n", " . nanval: No data value default is -1e20\n", - " . kstplist: Time steps list\n", - " . kperlist: Stress period list\n", - " . keylist: Flow terms to output\n", + " . kstpkper: A tuple containing the time step and stress period (kstp, kper) or\n", + " list of (kstp, kper) tuples. The kstp and kper values are zero based.\n", + " . text: The text identifier for the record. Can be a str or list of strings.\n", " . smooth: True creates a smooth surface, default is False\n", " . point_scalars: True outputs point scalar values as well as cell values, default is False.\n", " Outputing point scalars will set smooth to True\n", @@ -344,14 +345,14 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "cbc_file = os.path.join(model_ws, 'freyberg.cbc')\n", "cbc_output_folder = os.path.join(workspace, 'cbc_output_test_parameters')\n", - "vtk.export_cbc(ml, cbc_file, cbc_output_folder, kstplist=[1], kperlist=[10, 11, 12, 13, 14, 15],\n", - " keylist=['CONSTANT HEAD', 'STORAGE'], point_scalars=True, binary=True)" + "vtk.export_cbc(ml, cbc_file, cbc_output_folder, kstpkper=[(0, 0), (0, 9), (0, 10), (0, 11)],\n", + " text=['CONSTANT HEAD', 'STORAGE'], point_scalars=True, binary=True)" ] }, { diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 6e288588ae..a44cc7b3a7 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -958,7 +958,7 @@ def _get_names(in_list): def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, - kstplist=None, kperlist=None, keylist=None, smooth=False, + kstpkper=None, text=None, smooth=False, point_scalars=False, binary=False): """ @@ -977,12 +977,12 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, binary file precision, default is 'single' nanval : scalar no data value - kstplist : list - list of timesteps - kperlist : list - list of stress periods - keylist : list - list of flow term names + kstpkper : tuple of ints or list of tuple of ints + A tuple containing the time step and stress period (kstp, kper). + The kstp and kper values are zero based. + text : str or list of str + The text identifier for the record. Examples include + 'RIVER LEAKAGE', 'STORAGE', 'FLOW RIGHT FACE', etc. smooth : bool If true a smooth surface will be output, default is False point_scalars : bool @@ -1023,18 +1023,33 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, imeth_dict = {record: imeth for (record, imeth) in zip(records, cbb.imethlist)} # get list of packages to export - if not keylist: + if text is not None: + # build keylist + if isinstance(text, str): + keylist = [text] + elif isinstance(text, list): + keylist = text + else: + raise Exception('text must be type str or list of str') + else: keylist = records - if not kperlist: - kperlist = list(set([x[1] for x in cbb.get_kstpkper() if x[1] > -1])) - else: - kperlist = [kper - 1 for kper in kperlist] + if kstpkper is not None: + if isinstance(kstpkper, tuple): + kstplist = [kstpkper[0]] + kperlist = [kstpkper[1]] + elif isinstance(kstpkper, list): + kstpkper_list = list(map(list, zip(*kstpkper))) + kstplist = kstpkper_list[0] + kperlist = kstpkper_list[1] + + else: + raise Exception('kstpkper must be tuple of (kstp, kper) or list ' + 'of tuples') - if not kstplist: - kstplist = list(set([x[0] for x in cbb.get_kstpkper() if x[0] > -1])) else: - kstplist = [kstp - 1 for kstp in kstplist] + kperlist = list(set([x[1] for x in cbb.get_kstpkper() if x[1] > -1])) + kstplist = list(set([x[0] for x in cbb.get_kstpkper() if x[0] > -1])) # get model name model_name = model.name @@ -1104,8 +1119,8 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, return -def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstplist=None, - kperlist=None, smooth=False, point_scalars=False, +def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None, + smooth=False, point_scalars=False, binary=False): """ @@ -1122,10 +1137,9 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstplist=None, output folder to write the data nanval : scalar no data value, default value is -1e20 - kstplist : list - list of timesteps - kperlist : list - list of stress periods + kstpkper : tuple of ints or list of tuple of ints + A tuple containing the time step and stress period (kstp, kper). + The kstp and kper values are zero based. smooth : bool If true a smooth surface will be output, default is False point_scalars : bool @@ -1149,22 +1163,26 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstplist=None, byte_order="LittleEndian" compressor="vtkZLibDataCompressor"> \n""") - # get teh ehads + # get the heads hds = HeadFile(hdsfile) - # get the model time - # totim_dict = dict(zip(hds.get_kstpkper(), model.dis.get_totim())) + if kstpkper is not None: + if isinstance(kstpkper, tuple): + kstplist = [kstpkper[0]] + kperlist = [kstpkper[1]] - # set up time step and stress periods for output - if not kperlist: - kperlist = list(set([x[1] for x in hds.get_kstpkper() if x[1] > -1])) - else: - kperlist = [kper - 1 for kper in kperlist] + elif isinstance(kstpkper, list): + kstpkper_list = list(map(list, zip(*kstpkper))) + kstplist = kstpkper_list[0] + kperlist = kstpkper_list[1] + + else: + raise Exception('kstpkper must be tuple of (kstp, kper) or list ' + 'of tuples') - if not kstplist: - kstplist = list(set([x[0] for x in hds.get_kstpkper() if x[0] > -1])) else: - kstplist = [kstp - 1 for kstp in kstplist] + kperlist = list(set([x[1] for x in hds.get_kstpkper() if x[1] > -1])) + kstplist = list(set([x[0] for x in hds.get_kstpkper() if x[0] > -1])) # set upt the vtk vtk = Vtk(model, smooth=smooth, point_scalars=point_scalars, nanval=nanval)