diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 68d4ea5035..02eb9fc684 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') @@ -10,66 +10,205 @@ 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""" + 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) + # binary test + m.upw.hk.export(os.path.join(binot, 'array_3d_test'), fmt='vtk', + name='hk_points', point_scalars=True, binary=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') + + # binary test + m.rch.rech.export(os.path.join(binot, 'transient_2d_test'), fmt='vtk', + binary=True) + + 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') + # 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 + + +# 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') + + # binary test + m.export(os.path.join(binot, m.name), fmt='vtk', binary=True) -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) 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 = ['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_binary_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, '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 + 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 + 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 + 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_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'), + kstpkper=[(0, 0), (0, 1), (0, 2)], point_scalars=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_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_binary_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..f8e125bfbc --- /dev/null +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -0,0 +1,387 @@ +{ + "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": 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": "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": {}, + "source": [ + "# 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": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# 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", + "# 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": "markdown", + "metadata": {}, + "source": [ + "#### Export model bottoms" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# 3D Array export\n", + "# export model bottoms\n", + "model_bottom_output = os.path.join(ot_arrays_folder, 'BOTM')\n", + "ml.dis.botm.export(model_bottom_output, fmt='vtk')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Export transient array recharge" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "# 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": "markdown", + "metadata": {}, + "source": [ + "#### Export HK with point scalars\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "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": "markdown", + "metadata": {}, + "source": [ + "# 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": "markdown", + "metadata": {}, + "source": [ + "#### Export dis and upw package" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# 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": "markdown", + "metadata": {}, + "source": [ + "# 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": "markdown", + "metadata": {}, + "source": [ + "#### Export model as binary .vtu files" + ] + }, + { + "cell_type": "code", + "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": {}, + "source": [ + "# 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", + " . 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", + " . binary: True turns on binary vtk export, default is False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 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": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# 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": "markdown", + "metadata": {}, + "source": [ + "### Export heads to binary .vtu files with point head scalars" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "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", + "# 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)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "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", + " . 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", + " . binary: True turns on binary vtk export, default is False" + ] + }, + { + "cell_type": "code", + "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, kstpkper=[(0, 0), (0, 9), (0, 10), (0, 11)],\n", + " text=['CONSTANT HEAD', 'STORAGE'], point_scalars=True, binary=True)" + ] + }, + { + "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..f3a5c15f8d 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,20 @@ 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. + + 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 + + """ assert isinstance(ml, ModelInterface) @@ -391,13 +408,23 @@ 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) + 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, binary=binary, + nanval=nanval) + 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 +434,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 +445,19 @@ 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. + + 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 @@ -455,6 +497,17 @@ 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) + binary = kwargs.get('binary', False) + nanval = kwargs.get('nanval', -1e20) + vtk.export_package(pak.parent, pak.name, f, smooth=smooth, + point_scalars=point_scalars, binary=binary, + nanval=nanval) + + else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -655,7 +708,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 +717,18 @@ 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 + binary: bool + For vtk, if set to True will output binary .vtu files. Default + is False which exports standard text xml .vtu files. """ @@ -768,11 +828,19 @@ 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) + nanval = kwargs.get('nanval', -1e20) + vtk.export_transient(t2d.model, t2d.array, f, name, smooth=smooth, + point_scalars=point_scalars, array2d=True, + nanval=nanval) 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 +849,18 @@ 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 + binary: bool + For vtk, if set to Ture will output binary .vtu files. Default + is False which exports standard text xml .vtu files. """ @@ -903,11 +978,25 @@ 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) + 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, binary=binary, + nanval=nanval) + 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 +1005,18 @@ 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 + 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 " \ @@ -1013,6 +1109,18 @@ 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) + 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, + 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 87852d2d91..a44cc7b3a7 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -2,9 +2,137 @@ 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 +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: + """ + 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 + 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 s = indent_level * indent_char + tag indent_level += 1 f.write(s + '\n') @@ -12,113 +140,223 @@ 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 + 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, output_filename, model, verbose=None): + def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, + point_scalars=False): + + 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 = {} + self.smooth = smooth + self.point_scalars = point_scalars + + # check if structured grid, vtk only supports structured grid + assert (isinstance(self.modelgrid, StructuredGrid)) + + # cbd + self.cbd_on = False + + # get ibound + if self.modelgrid.idomain is None: + # ibound = None + 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'): + + 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() + return - def add_array(self, name, a): - assert a.shape == self.shape + def add_array(self, name, a, array2d=False): + + """ + + Adds an array to the vtk object + + 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 + + """ + + if name == 'ibound': + return + + # 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 + + 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 - def write(self, shared_vertex=False, ibound_filter=False, htop=None): + def write(self, output_file, timeval=None): """ + writes the stored arrays to vtk file + 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. + 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 """ - indent_level = 0 - if self.verbose: - print('writing vtk file') - f = open(self.output_filename, 'w') + # make sure file ends with vtu + assert output_file.lower().endswith(".vtu") - ibound = None - if ibound_filter: - ibound = self.modelgrid.idomain + # get the active data cells based on the data arrays and ibound + actwcells3d = self._configure_data_arrays() + actwcells = actwcells3d.ravel() - 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) + # 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 +369,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 +404,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 +416,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 +449,1106 @@ 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 write_binary(self, output_file): + """ - Write a numpy array to the vtk file + + outputs binary .vtu file + + Parameters + ---------- + + output_file : str + vtk output file """ - # header tag - s = ''.format(name) + # 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 arrays and active cells to find where active data + exists, and what cells to output. + """ + + # 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 in self.arrays: + array = self.arrays[name] + # 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): + + """ + + 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: + actwcells = self.ibound + + 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 + + @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): + """ + + Writes the data array to the output vtk file + + 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 + + """ + + 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_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 * ' ' + 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 + Builds the iverts based on the vertices being output + + Parameters + ---------- + verts : array + vertices being output + + Returns + ------- + + iverts : array + array of ivert values + + """ + ncells = len(verts) npoints = ncells * 8 - verts = np.empty((npoints, 3), dtype=np.float) 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 + 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 + + +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', nanval=-1e+20, + kstpkper=None, text=None, smooth=False, + point_scalars=False, binary=False): + """ + + Exports cell by cell file to vtk + + 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 + 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 + 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. + + """ + + mg = model.modelgrid + shape = (mg.nlay, mg.nrow, mg.ncol) + + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + # set up the pvd file to make the output files time enabled + pvdfile = open( + os.path.join(otfolder, '{}_Heads.pvd'.format(model.name)), + 'w') + + pvdfile.write(""" + + \n""") + + # load cbc + + cbb = bf.CellBudgetFile(cbcfile, precision=precision) + + # totim_dict = dict(zip(cbb.get_kstpkper(), model.dis.get_totim())) + + # get records + records = _get_names(cbb.get_unique_record_names()) + + # build imeth lookup + imeth_dict = {record: imeth for (record, imeth) in zip(records, + cbb.imethlist)} + # get list of packages to export + 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 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') + + else: + 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 + + vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars) - verts[ipoint, 0:2] = np.array(pt1) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + # export data + addarray = False + count = 1 + for kper in kperlist: + for kstp in kstplist: - verts[ipoint, 0:2] = np.array(pt2) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + 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: - verts[ipoint, 0:2] = np.array(pt0) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + try: + rec = cbb.get_data(kstpkper=(kstp, kper), text=name, + full3D=True) - verts[ipoint, 0:2] = np.array(pt3) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + if len(rec) > 0: + array = rec[0] # need to fix for multiple pak + addarray = True - z = top_botm[k, i, j] + except ValueError: - verts[ipoint, 0:2] = np.array(pt1) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + rec = cbb.get_data(kstpkper=(kstp, kper), text=name)[0] - verts[ipoint, 0:2] = np.array(pt2) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + 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) - verts[ipoint, 0:2] = np.array(pt0) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + array[lyr, row, col] = q - verts[ipoint, 0:2] = np.array(pt3) - verts[ipoint, 2] = z - ivert.append(ipoint) - ipoint += 1 + addarray = True + else: + raise Exception('Data type not currently supported ' + 'for cbc output') + # print('Data type not currently supported ' + # 'for cbc output') - iverts.append(ivert) + if addarray: - return verts, iverts + # 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 + if binary: + vtk.write_binary(otfile) + else: + vtk.write(otfile) + count += 1 + # finish writing the pvd file + pvdfile.write(""" +""") + + pvdfile.close() + return + + +def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None, + smooth=False, point_scalars=False, + binary=False): + """ + + 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 + 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 + 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. + + """ + + # 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 the heads + hds = HeadFile(hdsfile) + + 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') + + else: + 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) + + # 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)]) + if binary: + vtk.write_binary(otfile) + else: + 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, + binary=False): + + """ + + 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. + + """ + + 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)) + 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, + binary=False): + """ + + 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): + 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)) + if binary: + vtk.write_binary(ot_file) + else: + 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, + binary=False): + + """ + Exports package to vtk + + 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 + 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)) + if binary: + vtk.write_binary(ot_file) + else: + 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) + 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, binary=False): + """ + + 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) + 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, + binary=binary)