From fa27c4a6753c81901d4ee330dc012aa221d3df4e Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 28 Sep 2021 12:30:28 -0700 Subject: [PATCH 1/5] update(flopy3_vtk_export.ipynb): updates for new VTK class and routines --- examples/Notebooks/flopy3_vtk_export.ipynb | 674 +++++++++++++++++---- 1 file changed, 561 insertions(+), 113 deletions(-) diff --git a/examples/Notebooks/flopy3_vtk_export.ipynb b/examples/Notebooks/flopy3_vtk_export.ipynb index a8ca0815d..880049bd8 100644 --- a/examples/Notebooks/flopy3_vtk_export.ipynb +++ b/examples/Notebooks/flopy3_vtk_export.ipynb @@ -5,10 +5,13 @@ "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" + "\n", + "The `Vtk()` class in FloPy allows users to export Structured, Vertex, and Unstructured Grid based models to Visualization ToolKit files for display. 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 to `Vtk()`\n", + " - example usage of the `Vtk()` class object to output data \n", + " - exporting heads and model output data\n", + " - exporting modpath pathlines to `Vtk()`" ] }, { @@ -20,9 +23,8 @@ "name": "stdout", "output_type": "stream", "text": [ - "3.8.10 (default, May 19 2021, 11:01:55) \n", - "[Clang 10.0.0 ]\n", - "flopy version: 3.3.4\n" + "3.7.4 (default, Aug 9 2019, 18:34:13) [MSC v.1915 64 bit (AMD64)]\n", + "flopy version: 3.3.5\n" ] } ], @@ -51,16 +53,7 @@ "cell_type": "code", "execution_count": 2, "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/Users/jdhughes/Documents/Development/flopy_git/flopy_fork/flopy/mbase.py:354: DeprecationWarning: xul/yul have been deprecated. Use xll/yll instead.\n", - " warnings.warn(\n" - ] - } - ], + "outputs": [], "source": [ "# load model for examples\n", "nam_file = \"freyberg.nam\"\n", @@ -84,32 +77,36 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Exporting to VTK\n", + "# Basic exporting to VTK using the `.export()` method\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 path must be provided** provided and the the `fmt` flag should be set to 'vtk'\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Exporting FloPy arrays to .vtu files\n", + "## 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" + "All array exports have the following optional keyword arguments:\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", + " - `name`: A name can be specified to use for the output filename and array scalar name, by default the FloPy array name is used\n", + " - `binary`: argument that can be specified to switch between binary and ASCII, default is True\n", + " - `xml`: True will write an xml base vtk file, default is False\n", + " - `masked_values`: list or tuple of values to mask (set to nan) when writing a array\n", + " - `vertical_exageration`: floating point value that can be used to scale the vertical exageration of the vtk points. Default is 1.\n", + " \n", + "Tranient type array exports (\"stress_period_data\"; ex. recharge data, well flux, etc ...) have additional optional keyword arguments:\n", + " - `pvd`: True will write a paraview data file with simulation time for animations. Default is False\n", + " - `kper`: a list, tuple, or integer value of specific stess periods to output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export model top" + "### Export model top" ] }, { @@ -119,19 +116,20 @@ "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", + "output_dir = os.path.join(workspace, 'arrays_test')\n", + "if not os.path.exists(output_dir):\n", + " os.mkdir(output_dir)\n", + " \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') " + "model_top_dir = os.path.join(output_dir, 'TOP')\n", + "ml.dis.top.export(model_top_dir, fmt='vtk') " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export model bottoms" + "### Export model bottoms" ] }, { @@ -142,34 +140,42 @@ "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')" + "model_bottom_dir = os.path.join(output_dir, 'BOTM')\n", + "ml.dis.botm.export(model_bottom_dir, fmt='vtk')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export transient array recharge" + "### Export transient array recharge" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Switching to xml, ASCII and standard binary are not supported by Paraview's PVD reader\n" + ] + } + ], "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')" + "model_recharge_dir = os.path.join(output_dir, 'RECH')\n", + "ml.rch.rech.export(model_recharge_dir, fmt='vtk', pvd=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export HK with point scalars\n" + "### Export HK with point scalars\n" ] }, { @@ -180,29 +186,33 @@ "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)" + "model_hk_dir = os.path.join(output_dir, 'HK')\n", + "ml.upw.hk.export(model_hk_dir, smooth=True, fmt='vtk', name='HK', point_scalars=True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Package export to .vtu files\n", + "## 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" + "Package export has the following keyword arguments:\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", + " - `name`: A name can be specified to use for the output filename and array scalar name, by default the FloPy array name is used\n", + " - `binary`: argument that can be specified to switch between binary and ASCII, default is True\n", + " - `xml`: True will write an xml base vtk file, default is False\n", + " - `masked_values`: list or tuple of values to mask (set to nan) when writing a array\n", + " - `vertical_exageration`: floating point value that can be used to scale the vertical exageration of the vtk points. Default is 1.\n", + " - `pvd`: True will write a paraview data file with simulation time for animations. Default is False\n", + " - `kper`: a list, tuple, or integer value of specific stess periods to output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export dis and upw package" + "### Export dis and upw package" ] }, { @@ -213,15 +223,17 @@ "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", + "output_dir = os.path.join(workspace, 'package_output_test')\n", + "if not os.path.exists(output_dir):\n", + " os.mkdir(output_dir)\n", + " \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)" + "dis_output_dir = os.path.join(output_dir, 'DIS')\n", + "ml.dis.export(dis_output_dir, fmt='vtk')\n", + "\n", + "#export upw with point scalars as a binary xml based vtk file\n", + "upw_output_dir = os.path.join(output_dir, 'UPW')\n", + "ml.upw.export(upw_output_dir, fmt='vtk', point_scalars=True, xml=True)" ] }, { @@ -230,20 +242,25 @@ "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" + "Model export has the following optional keyword arguments:\n", + " \n", + " - `package_names`: a list of package names to export, default is None and will export all packages in the model.\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", + " - `name`: A name can be specified to use for the output filename and array scalar name, by default the FloPy array name is used\n", + " - `binary`: argument that can be specified to switch between binary and ASCII, default is True\n", + " - `xml`: True will write an xml base vtk file, default is False\n", + " - `masked_values`: list or tuple of values to mask (set to nan) when writing a array\n", + " - `vertical_exageration`: floating point value that can be used to scale the vertical exageration of the vtk points. Default is 1.\n", + " - `pvd`: True will write a paraview data file with simulation time for animations. Default is False\n", + " - `kper`: a list, tuple, or integer value of specific stess periods to output" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "#### Export model as binary .vtu files" + "### Export model as a binary unstructured vtk file" ] }, { @@ -256,7 +273,7 @@ { "data": { "text/plain": [ - "'./temp/VTK_EXAMPLE_OUTPUTS/model_output_test'" + "'.\\\\temp\\\\VTK_EXAMPLE_OUTPUTS\\\\model_output_test'" ] }, "execution_count": 9, @@ -265,103 +282,534 @@ } ], "source": [ + "model_output_dir = os.path.join(workspace, 'model_output_test')\n", + "ml.export(model_output_dir, fmt='vtk')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Exporting using the `Vtk()` class object \n", + "\n", + "To export custom arrays, or choose a custom combination of model inputs to view, the user first needs to instantiate a new `Vtk()` object. The `Vtk()` object has a single required parameter and a number of optional parameters that the user can take advantage of. These parameters are as follows:\n", "\n", - "model_output = os.path.join(workspace, 'model_output_test')\n", - "ml.export(model_output, fmt='vtk', binary=True)\n" + " - `model`: any flopy model object can be supplied to create the vtk geometry. Either the model (recommended!) or modelgrid parameter must be supplied to the Vtk() object.\n", + " - `modelgrid`: any flopy modelgrid object (StructuredGrid, VertexGrid, UnstructuredGrid) can be supplied, in leiu of a model object, to create the vtk geometery. \n", + " - `vertical_exageration`: floating point value that can be used to scale the vertical exageration of the vtk points. Default is 1.\n", + " - `binary`: boolean flag to switch between binary and ASCII vtk files. Default is True.\n", + " - `xml`: boolean flag to write xml based vtk files. Default is False\n", + " - `pvd`: boolean flag to write a paraview data file for transient series of vtu files. This file relates model time to vtu file for animations. Default is False. If set to True Vtk() will automatically write xml based vtu files. \n", + " - `shared_points`: boolean flag to write shared vertices within the vtk file. Default is False.\n", + " - `smooth`: boolean flag to interpolate vertex elevations using IDW based on shared cell elevations. Default is False.\n", + " - `point_scalars`: boolean flag to write interpolated data at each point." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# create a binary XML VTK object and enable PVD file writing\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, vertical_exageration=10)\n", + "vtkobj" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Export heads\n", + "## Adding array data to the `Vtk()` object\n", "\n", - " vtk.export_heads(model, headsfile, output_folder)\n", + "The `Vtk()` object has an `add_array()` method that lets the user add array data to the Field data section of the VTK file.\n", "\n", - " Heads export has the following arguments:\n", + "`add_array()` has a few parameters for the user:\n", + " - `array` : numpy array that has a size equal to the number of cells in the model (modelgrid.nnodes).\n", + " - `name` : array name (string)\n", + " - `masked_values` : list of array values to mask/set to NaN " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a vtk object\n", + "vtkobj = vtk.Vtk(ml, vertical_exageration=10)\n", + "\n", + "## create some random array data\n", + "r_array = np.random.random(ml.modelgrid.nnodes) * 100\n", + "\n", + "## add random data to the VTK object\n", + "vtkobj.add_array(r_array, \"random_data\")\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" + "## add the model botom data to the VTK object\n", + "vtkobj.add_array(ml.dis.botm.array, \"botm\")\n", + "\n", + "## write the vtk object to file\n", + "vtkobj.write(os.path.join(output_dir, \"Array_example\", \"model.vtu\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Export heads to binary .vtu files\n", + "## Adding transient array data to the `Vtk()` object\n", + "\n", + "The `Vtk()` class has an `add_transient_array()` method that allows the user to create a series of time varying `Vtk()` files that can be used for animation in VTK viewers.\n", "\n", - "#### Exports a .pvd file that can be loaded into Paraview to view heads as time series" + "The `add_transient_array()` method accepts a dictionary of array2d, array3d, or numpy array objects. Parameters include:\n", + " - `d`: dictionary of array2d, array3d, or numpy array objects\n", + " - `name`: parameter name, required when user provides a dictionary of numpy arrays\n", + " - `masked_values`: optional list of values to set equal to NaN." ] }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 12, "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)" + "# create a vtk object\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, vertical_exageration=10)\n", + "\n", + "## add recharge to the VTK object\n", + "recharge = ml.rch.rech.transient_2ds\n", + "vtkobj.add_transient_array(recharge, \"recharge\", masked_values=[0,])\n", + "\n", + "## write vtk files\n", + "vtkobj.write(os.path.join(output_dir, \"tr_array_example\", \"recharge.vtu\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Export heads to binary .vtu files with point head scalars" + "## Adding transient list data to the `Vtk()` object\n", + "\n", + "The `Vtk()` class has an `add_transient_list()` method that allows the user to create a series of time varying `Vtk()` files that can be used for animation in VTK viewers.\n", + "\n", + "The `add_transient_list()` method accepts a FloPy mflist (transient list) type object. Parameters include:\n", + " - `mflist`: flopy transient list object\n", + " - `masked_values`: list of values to set equal to NaN" ] }, { "cell_type": "code", - "execution_count": 11, + "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", - "# 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)" + "# create the vtk object\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, vertical_exageration=10)\n", + "\n", + "## add well fluxes to the VTK object\n", + "spd = ml.wel.stress_period_data\n", + "vtkobj.add_transient_list(spd, masked_values=[0,])\n", + "\n", + "## write vtk files\n", + "vtkobj.write(os.path.join(output_dir, \"tr_list_example\", \"wel_flux.vtu\"))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Export output cell by cell file to .vtu\n", + "## Adding packages to the `Vtk` object\n", + "\n", + "The `Vtk()` class has a method for adding package data to a VTK file as Field Data. The `add_package()` method allows the user to add packages for subsequent export. `add_package()` takes the following parameters:\n", + "\n", + " - `pkg`: flopy package object\n", + " - `masked_values`: optional list of values to set to NaN.\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the following example, a HFB package is added to the existing freyberg model and then exported with the WEL package." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [], + "source": [ + "# create a HFB package for the example\n", + "hfb_data = []\n", + "for k in range(3):\n", + " for i in range(20):\n", + " rec = [k, i, 6, i, 7, 1e-06]\n", + " hfb_data.append(rec)\n", + " \n", + "hfb = flopy.modflow.ModflowHfb(ml, hfb_data=hfb_data)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "# export HFB and WEL packages using Vtk()\n", + "vtkobj = vtk.Vtk(ml, vertical_exageration=10)\n", "\n", - " vtk.export_heads(model, cellbycellfile, outputfolder)\n", + "vtkobj.add_package(hfb)\n", + "vtkobj.add_package(ml.wel)\n", "\n", - " Heads export has the following arguments:\n", + "vtkobj.write(os.path.join(output_dir, \"package_example\", \"package_export.vtu\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Export heads to binary .vtu files\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" + "Once a `Vtk()` object is instantiated (see above), the `Vtk()` class has an `add_heads()` method. This method has a few parameters:\n", + " - `hds`: a flopy FormattedHeadFile or HeadFile object. This method also accepts DrawdownFile, and ConcentrationFile objects.\n", + " - `kstpkper`: optional list of zero based (timestep, stress period) tuples to output. Default is None and will output all data to a series of vtu files\n", + " - `masked_values`: optional list of values to set to NaN, default is None.\n" ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 16, + "metadata": {}, + "outputs": [], + "source": [ + "# import the HeadFile reader and read in the head file\n", + "from flopy.utils import HeadFile\n", + "\n", + "head_file = os.path.join(model_ws, 'freyberg.hds')\n", + "hds = HeadFile(head_file)\n", + "\n", + "\n", + "# create the vtk object and export heads\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, vertical_exageration=10)\n", + "vtkobj.add_heads(hds)\n", + "vtkobj.write(os.path.join(workspace, 'heads_output_test', 'freyberg_head.vtu'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Export heads as point scalar arrays" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [], + "source": [ + "# export heads as point scalars\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, point_scalars=True, vertical_exageration=10)\n", + "\n", + "# export heads for time step 1, stress periods 1, 50, 100, 1000\n", + "vtkobj.add_heads(hds, kstpkper=[(0,0), (0, 49), (0, 99), (0, 999)])\n", + "vtkobj.write(os.path.join(workspace, 'heads_output_test_parameters', 'freyberg_head.vtu'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Export cell budget information\n", + "\n", + "Once a `Vtk()` object is instantiated (see above), the `Vtk()` class has an `add_cell_budget()` method. This method has a few parameters:\n", + " - `cbc`: flopy CellBudgetFile object\n", + " - `text`: Optional text identifier for a record type. Examples include 'RIVER LEAKAGE', 'STORAGE', etc... Default is None and will export all cell budget information to vtk files\n", + " - `kstpkper`: optional list of zero based (timestep, stress period) tuples to output. Default is None and will output all data to a series of vtu files\n", + " - `masked_values`: optional list of values to set to NaN, default is None." + ] + }, + { + "cell_type": "code", + "execution_count": 18, "metadata": {}, "outputs": [], "source": [ + "# import the CellBudgetFile reader and read the CBC file\n", + "from flopy.utils import CellBudgetFile\n", + "\n", "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)" + "cbc = CellBudgetFile(cbc_file)\n", + "\n", + "# export the cbc file to a series of Vtu files with a PVD file for animation\n", + "vtkobj = vtk.Vtk(ml, xml=True, pvd=True, vertical_exageration=10)\n", + "vtkobj.add_cell_budget(cbc, kstpkper=[(0, 0), (0, 9), (0, 10), (0, 11)])\n", + "vtkobj.write(os.path.join(workspace, 'cbc_output_test_parameters', 'freyberg_cbc.vtu'))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Exporting vectors from the Cell Budget File with the `Vtk()` class\n", + "\n", + "The `Vtk()` class has an `add_vector()` method that allows the user to write vector information to VTK files. This method can be used to export information such as cell centered specific discharge.\n", + "\n", + "The `add_vector()` method accepts a numpy array of vector data. The array size must be 3 * the number of model cells (3 * modelgrid.nnodes). Parameters include:\n", + " - `vector`: numpy array of size 3 * nnodes\n", + " - `name`: name of the vector \n", + " - `masked_values`: list of values to set equal to NaN\n" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "# get frf, fff, flf from the Cell Budget file (or SPDIS when using MF6)\n", + "from flopy.utils import postprocessing\n", + "\n", + "frf = cbc.get_data(text=\"FLOW RIGHT FACE\", kstpkper=(0, 9), full3D=True)[0]\n", + "fff = cbc.get_data(text=\"FLOW FRONT FACE\", kstpkper=(0, 9), full3D=True)[0]\n", + "flf = cbc.get_data(text=\"FLOW LOWER FACE\", kstpkper=(0, 9), full3D=True)[0]\n", + "\n", + "spdis = postprocessing.get_specific_discharge((frf, fff, flf), ml)\n", + "\n", + "# create the Vtk() object\n", + "vtkobj = vtk.Vtk(ml, vertical_exageration=10)\n", + "\n", + "# add the vector\n", + "vtkobj.add_vector(spdis, name='spdis')\n", + "\n", + "# write to file\n", + "vtkobj.write(os.path.join(output_dir, \"vector_example\", \"spdis_vector.vtu\"))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Adding Modpath timeseries or pathline data to `Vtk()`\n", + "\n", + "The `Vtk()` class support writing MODPATH pathline/timeseries data to a Vtk file. To start the example, let's first load and run a MODPATH simulation (see flopy3_modpath7_unstructured_example for details) and then add the output to a `Vtk()` object." + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "WARNING: Unable to resolve dimension of ('gwf6', 'disv', 'cell2d', 'cell2d', 'icvert') based on shape \"ncvert\".\n", + "writing simulation...\n", + " writing simulation name file...\n", + " writing simulation tdis package...\n", + " writing ims package ims...\n", + " writing model mp7p2...\n", + " writing model name file...\n", + " writing package disv...\n", + " writing package ic...\n", + " writing package npf...\n", + " writing package wel_0...\n", + "INFORMATION: maxbound in ('gwf6', 'wel', 'dimensions') changed to 1 based on size of stress_period_data\n", + " writing package rcha_0...\n", + " writing package riv_0...\n", + "INFORMATION: maxbound in ('gwf6', 'riv', 'dimensions') changed to 21 based on size of stress_period_data\n", + " writing package oc...\n", + "FloPy is using the following executable to run the model: .\\mf6.EXE\n", + " MODFLOW 6\n", + " U.S. GEOLOGICAL SURVEY MODULAR HYDROLOGIC MODEL\n", + " VERSION 6.2.1 02/18/2021\n", + "\n", + " MODFLOW 6 compiled Feb 18 2021 21:14:51 with IFORT compiler (ver. 19.10.3)\n", + "\n", + "This software has been approved for release by the U.S. Geological \n", + "Survey (USGS). Although the software has been subjected to rigorous \n", + "review, the USGS reserves the right to update the software as needed \n", + "pursuant to further analysis and review. No warranty, expressed or \n", + "implied, is made by the USGS or the U.S. Government as to the \n", + "functionality of the software and related material nor shall the \n", + "fact of release constitute any such warranty. Furthermore, the \n", + "software is released on condition that neither the USGS nor the U.S. \n", + "Government shall be held liable for any damages resulting from its \n", + "authorized or unauthorized use. Also refer to the USGS Water \n", + "Resources Software User Rights Notice for complete use, copyright, \n", + "and distribution information.\n", + "\n", + " \n", + " Run start date and time (yyyy/mm/dd hh:mm:ss): 2021/09/27 18:11:17\n", + " \n", + " Writing simulation list file: mfsim.lst\n", + " Using Simulation name file: mfsim.nam\n", + " \n", + " Solving: Stress period: 1 Time step: 1\n", + " \n", + " Run end date and time (yyyy/mm/dd hh:mm:ss): 2021/09/27 18:11:17\n", + " Elapsed run time: 0.180 Seconds\n", + " \n", + "\n", + "WARNING REPORT:\n", + "\n", + " 1. NONLINEAR BLOCK VARIABLE 'OUTER_HCLOSE' IN FILE 'mp7p2.ims' WAS\n", + " DEPRECATED IN VERSION 6.1.1. SETTING OUTER_DVCLOSE TO OUTER_HCLOSE VALUE.\n", + " 2. LINEAR BLOCK VARIABLE 'INNER_HCLOSE' IN FILE 'mp7p2.ims' WAS DEPRECATED\n", + " IN VERSION 6.1.1. SETTING INNER_DVCLOSE TO INNER_HCLOSE VALUE.\n", + " Normal termination of simulation.\n", + "FloPy is using the following executable to run the model: .\\mp7.EXE\n", + "\n", + "MODPATH Version 7.2.001 \n", + "Program compiled Feb 18 2021 21:19:24 with IFORT compiler (ver. 19.10.3) \n", + " \n", + " \n", + "Run particle tracking simulation ...\n", + "Processing Time Step 1 Period 1. Time = 1.00000E+03 Steady-state flow \n", + "\n", + "Particle Summary:\n", + " 0 particles are pending release.\n", + " 0 particles remain active.\n", + " 16 particles terminated at boundary faces.\n", + " 0 particles terminated at weak sink cells.\n", + " 0 particles terminated at weak source cells.\n", + " 0 particles terminated at strong source/sink cells.\n", + " 0 particles terminated in cells with a specified zone number.\n", + " 0 particles were stranded in inactive or dry cells.\n", + " 0 particles were unreleased.\n", + " 0 particles have an unknown status.\n", + " \n", + "Normal termination. \n", + "FloPy is using the following executable to run the model: .\\mp7.EXE\n", + "\n", + "MODPATH Version 7.2.001 \n", + "Program compiled Feb 18 2021 21:19:24 with IFORT compiler (ver. 19.10.3) \n", + " \n", + " \n", + "Run particle tracking simulation ...\n", + "Processing Time Step 1 Period 1. Time = 1.00000E+03 Steady-state flow \n", + "\n", + "Particle Summary:\n", + " 0 particles are pending release.\n", + " 0 particles remain active.\n", + " 416 particles terminated at boundary faces.\n", + " 0 particles terminated at weak sink cells.\n", + " 0 particles terminated at weak source cells.\n", + " 0 particles terminated at strong source/sink cells.\n", + " 0 particles terminated in cells with a specified zone number.\n", + " 0 particles were stranded in inactive or dry cells.\n", + " 0 particles were unreleased.\n", + " 0 particles have an unknown status.\n", + " \n", + "Normal termination. \n", + "Output file located: mp7p2.hds\n", + "Output file located: mp7p2.cbb\n" + ] + } + ], + "source": [ + "# load and run the vertex grid model and modpath7\n", + "sys.path.append(os.path.join(\"..\", \"common\"))\n", + "import setup_pmv_demo\n", + "setup_pmv_demo.run()\n", + "\n", + "# check if model ran properly\n", + "modelpth = os.path.join(\"data\", \"mp7_ex2\", \"mf6\")\n", + "files = ['mp7p2.hds', 'mp7p2.cbb']\n", + "for f in files:\n", + " if os.path.isfile(os.path.join(modelpth, f)):\n", + " msg = 'Output file located: {}'.format(f)\n", + " print (msg)\n", + " else:\n", + " errmsg = 'Error. Output file cannot be found: {}'.format(f)\n", + " print (errmsg)" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "loading simulation...\n", + " loading simulation name file...\n", + " loading tdis package...\n", + " loading model gwf6...\n", + " loading package disv...\n", + "WARNING: Unable to resolve dimension of ('gwf6', 'disv', 'cell2d', 'cell2d', 'icvert') based on shape \"ncvert\".\n", + " loading package ic...\n", + " loading package npf...\n", + " loading package wel...\n", + " loading package rch...\n", + " loading package riv...\n", + " loading package oc...\n", + " loading ims package mp7p2...\n", + "[[0, 0, 1, 2, 3, 0, None], [1, 1, 4, 5, 2, 1, None], [2, 4, 6, 7, 5, 4, None], [3, 6, 8, 9, 7, 6, None], [4, 8, 10, 11, 9, 8, None], [5, 10, 12, 13, 11, 10, None], [6, 12, 14, 15, 13, 12, None], [7, 14, 16, 17, 15, 14, None], [8, 16, 18, 19, 17, 16, None], [9, 18, 20, 21, 19, 18, None], [10, 20, 22, 23, 21, 20, None], [11, 22, 24, 25, 23, 22, None], [12, 24, 26, 27, 25, 24, None], [13, 26, 28, 29, 27, 26, None], [14, 28, 30, 31, 29, 28, None], [15, 30, 32, 33, 31, 30, None], [16, 32, 34, 35, 33, 32, None], [17, 34, 36, 37, 35, 34, None], [18, 36, 38, 39, 37, 36, None], [19, 38, 40, 41, 39, 38, None], [20, 3, 2, 42, 43, 3, None], [21, 2, 5, 44, 42, 2, None], [22, 5, 7, 45, 44, 5, None], [23, 7, 9, 46, 45, 7, None], [24, 9, 11, 47, 46, 9, None], [25, 11, 13, 48, 47, 11, None], [26, 13, 15, 49, 48, 13, None], [27, 15, 17, 50, 49, 15, None], [28, 17, 19, 51, 50, 17, None], [29, 19, 21, 52, 51, 19, None], [30, 21, 23, 53, 52, 21, None], [31, 23, 25, 54, 53, 23, None], [32, 25, 27, 55, 54, 25, None], [33, 27, 29, 56, 55, 27, None], [34, 29, 31, 57, 56, 29, None], [35, 31, 33, 58, 57, 31, None], [36, 33, 35, 59, 58, 33, None], [37, 35, 37, 60, 59, 35, None], [38, 37, 39, 61, 60, 37, None], [39, 39, 41, 62, 61, 39, None], [40, 43, 42, 63, 64, 43, None], [41, 42, 44, 65, 63, 42, None], [42, 44, 45, 66, 65, 44, None], [43, 45, 46, 67, 66, 45, None], [44, 46, 47, 68, 67, 46, None], [45, 47, 48, 69, 68, 47, None], [46, 48, 49, 70, 69, 48, None], [47, 49, 50, 71, 70, 49, None], [48, 50, 51, 72, 71, 50, None], [49, 51, 52, 73, 72, 51, None], [50, 52, 53, 74, 73, 52, None], [51, 53, 54, 75, 74, 53, None], [52, 54, 55, 76, 75, 54, None], [53, 55, 56, 77, 76, 55, None], [54, 56, 57, 78, 77, 56, None], [55, 57, 58, 79, 78, 57, None], [56, 58, 59, 80, 79, 58, None], [57, 59, 60, 81, 80, 59, None], [58, 60, 61, 82, 81, 60, None], [59, 61, 62, 83, 82, 61, None], [60, 64, 63, 84, 85, 64, None], [61, 63, 65, 86, 84, 63, None], [62, 65, 66, 87, 86, 65, None], [63, 66, 67, 88, 87, 66, None], [64, 67, 68, 89, 88, 67, None], [65, 68, 69, 90, 89, 68, None], [66, 69, 70, 91, 90, 69, None], [67, 70, 71, 92, 91, 70, None], [68, 71, 72, 93, 92, 71, None], [69, 72, 73, 94, 93, 72, None], [70, 73, 74, 95, 94, 73, None], [71, 74, 75, 96, 95, 74, None], [72, 75, 76, 97, 96, 75, None], [73, 76, 77, 98, 97, 76, None], [74, 77, 78, 99, 98, 77, None], [75, 78, 79, 100, 99, 78, None], [76, 79, 80, 101, 100, 79, None], [77, 80, 81, 102, 101, 80, None], [78, 81, 82, 103, 102, 81, None], [79, 82, 83, 104, 103, 82, None], [80, 85, 84, 105, 106, 85, None], [81, 84, 86, 107, 105, 84, None], [82, 86, 87, 108, 107, 86, None], [83, 87, 88, 109, 108, 87, None], [84, 88, 89, 110, 109, 88, None], [85, 89, 90, 111, 110, 89, None], [86, 90, 91, 112, 111, 90, None], [87, 91, 92, 113, 112, 91, None], [88, 92, 93, 114, 113, 92, None], [89, 93, 94, 115, 114, 93, None], [90, 94, 95, 116, 115, 94, None], [91, 95, 96, 117, 116, 95, None], [92, 96, 97, 118, 117, 96, None], [93, 97, 98, 119, 118, 97, None], [94, 98, 99, 120, 119, 98, None], [95, 99, 100, 121, 120, 99, None], [96, 100, 101, 122, 121, 100, None], [97, 101, 102, 123, 122, 101, None], [98, 102, 103, 124, 123, 102, None], [99, 103, 104, 125, 124, 103, None], [100, 106, 105, 126, 127, 106, None], [101, 105, 107, 128, 126, 105, None], [102, 107, 108, 129, 128, 107, None], [103, 108, 109, 130, 129, 108, None], [104, 109, 110, 131, 130, 109, None], [105, 110, 111, 132, 131, 110, None], [106, 111, 112, 133, 132, 111, None], [107, 112, 113, 134, 133, 112, None], [108, 113, 114, 135, 134, 113, None], [109, 114, 115, 136, 135, 114, None], [110, 115, 116, 137, 136, 115, None], [111, 116, 117, 138, 137, 116, None], [112, 117, 118, 139, 138, 117, None], [113, 118, 119, 140, 139, 118, None], [114, 119, 120, 141, 140, 119, None], [115, 120, 121, 142, 141, 120, None], [116, 121, 122, 143, 142, 121, None], [117, 122, 123, 144, 143, 122, None], [118, 123, 124, 145, 144, 123, None], [119, 124, 125, 146, 145, 124, None], [120, 127, 126, 147, 148, 127, None], [121, 126, 128, 149, 147, 126, None], [122, 128, 129, 150, 149, 128, None], [123, 129, 130, 151, 150, 129, None], [124, 130, 131, 152, 151, 130, None], [125, 131, 132, 153, 152, 131, None], [126, 132, 133, 154, 153, 132, None], [127, 133, 134, 155, 154, 133, None], [128, 134, 135, 156, 155, 134, None], [129, 135, 136, 157, 156, 135, None], [130, 136, 137, 158, 157, 136, None], [131, 137, 138, 159, 158, 137, None], [132, 138, 139, 160, 159, 138, None], [133, 139, 140, 161, 160, 139, None], [134, 140, 141, 162, 161, 140, None], [135, 141, 142, 163, 162, 141, None], [136, 142, 143, 164, 163, 142, None], [137, 143, 144, 165, 164, 143, None], [138, 144, 145, 166, 165, 144, None], [139, 145, 146, 167, 166, 145, None], [140, 148, 147, 168, 169, 148, None], [141, 147, 149, 170, 168, 147, None], [142, 149, 150, 171, 170, 149, None], [143, 150, 151, 172, 171, 150, None], [144, 151, 152, 173, 172, 151, None], [145, 152, 153, 174, 173, 152, None], [146, 153, 154, 175, 174, 153, None], [147, 154, 155, 176, 197, 175, 154], [148, 155, 156, 177, 203, 176, 155], [149, 156, 157, 178, 208, 177, 156], [150, 157, 158, 179, 213, 178, 157], [151, 158, 159, 180, 218, 179, 158], [152, 159, 160, 181, 180, 159, None], [153, 160, 161, 182, 181, 160, None], [154, 161, 162, 183, 182, 161, None], [155, 162, 163, 184, 183, 162, None], [156, 163, 164, 185, 184, 163, None], [157, 164, 165, 186, 185, 164, None], [158, 165, 166, 187, 186, 165, None], [159, 166, 167, 188, 187, 166, None], [160, 169, 168, 189, 190, 169, None], [161, 168, 170, 191, 189, 168, None], [162, 170, 171, 192, 191, 170, None], [163, 171, 172, 193, 192, 171, None], [164, 172, 173, 194, 193, 172, None], [165, 173, 174, 195, 194, 173, None], [166, 174, 175, 199, 196, 195, 174], [167, 175, 197, 198, 199, 175, None], [168, 197, 176, 200, 198, 197, None], [169, 199, 198, 201, 196, 199, None], [170, 198, 200, 202, 201, 198, None], [171, 176, 203, 204, 200, 176, None], [172, 203, 177, 205, 204, 203, None], [173, 200, 204, 206, 244, 202, 200], [174, 204, 205, 207, 250, 206, 204], [175, 177, 208, 209, 205, 177, None], [176, 208, 178, 210, 209, 208, None], [177, 205, 209, 211, 264, 207, 205], [178, 209, 210, 212, 269, 211, 209], [179, 178, 213, 214, 210, 178, None], [180, 213, 179, 215, 214, 213, None], [181, 210, 214, 216, 282, 212, 210], [182, 214, 215, 217, 287, 216, 214], [183, 179, 218, 219, 215, 179, None], [184, 218, 180, 220, 219, 218, None], [185, 215, 219, 221, 217, 215, None], [186, 219, 220, 222, 221, 219, None], [187, 180, 181, 223, 222, 220, 180], [188, 181, 182, 224, 223, 181, None], [189, 182, 183, 225, 224, 182, None], [190, 183, 184, 226, 225, 183, None], [191, 184, 185, 227, 226, 184, None], [192, 185, 186, 228, 227, 185, None], [193, 186, 187, 229, 228, 186, None], [194, 187, 188, 230, 229, 187, None], [195, 190, 189, 231, 232, 190, None], [196, 189, 191, 233, 231, 189, None], [197, 191, 192, 234, 233, 191, None], [198, 192, 193, 235, 234, 192, None], [199, 193, 194, 236, 235, 193, None], [200, 194, 195, 237, 236, 194, None], [201, 195, 196, 240, 238, 237, 195], [202, 196, 201, 239, 240, 196, None], [203, 201, 202, 246, 241, 239, 201], [204, 240, 239, 242, 238, 240, None], [205, 239, 241, 256, 243, 242, 239], [206, 202, 244, 245, 246, 202, None], [207, 244, 206, 247, 245, 244, None], [208, 246, 245, 248, 241, 246, None], [209, 245, 247, 249, 248, 245, None], [210, 206, 250, 251, 247, 206, None], [211, 250, 207, 252, 251, 250, None], [212, 247, 251, 253, 249, 247, None], [213, 251, 252, 254, 253, 251, None], [214, 241, 248, 255, 256, 241, None], [215, 248, 249, 257, 255, 248, None], [216, 256, 255, 258, 243, 256, None], [217, 255, 257, 259, 258, 255, None], [218, 249, 253, 260, 257, 249, None], [219, 253, 254, 261, 260, 253, None], [220, 257, 260, 262, 259, 257, None], [221, 260, 261, 263, 262, 260, None], [222, 207, 264, 265, 252, 207, None], [223, 264, 211, 266, 265, 264, None], [224, 252, 265, 267, 254, 252, None], [225, 265, 266, 268, 267, 265, None], [226, 211, 269, 270, 266, 211, None], [227, 269, 212, 271, 270, 269, None], [228, 266, 270, 272, 268, 266, None], [229, 270, 271, 273, 272, 270, None], [230, 254, 267, 274, 261, 254, None], [231, 267, 268, 275, 274, 267, None], [232, 261, 274, 276, 343, 263, 261], [233, 274, 275, 277, 349, 276, 274], [234, 268, 272, 278, 275, 268, None], [235, 272, 273, 279, 278, 272, None], [236, 275, 278, 280, 363, 277, 275], [237, 278, 279, 281, 368, 280, 278], [238, 212, 282, 283, 271, 212, None], [239, 282, 216, 284, 283, 282, None], [240, 271, 283, 285, 273, 271, None], [241, 283, 284, 286, 285, 283, None], [242, 216, 287, 288, 284, 216, None], [243, 287, 217, 289, 288, 287, None], [244, 284, 288, 290, 286, 284, None], [245, 288, 289, 291, 290, 288, None], [246, 273, 285, 292, 279, 273, None], [247, 285, 286, 293, 292, 285, None], [248, 279, 292, 294, 281, 279, None], [249, 292, 293, 295, 294, 292, None], [250, 286, 290, 296, 293, 286, None], [251, 290, 291, 297, 296, 290, None], [252, 293, 296, 298, 295, 293, None], [253, 296, 297, 299, 298, 296, None], [254, 217, 221, 300, 291, 289, 217], [255, 221, 222, 301, 300, 221, None], [256, 291, 300, 302, 299, 297, 291], [257, 300, 301, 303, 302, 300, None], [258, 222, 223, 304, 303, 301, 222], [259, 223, 224, 305, 304, 223, None], [260, 224, 225, 306, 305, 224, None], [261, 225, 226, 307, 306, 225, None], [262, 226, 227, 308, 307, 226, None], [263, 227, 228, 309, 308, 227, None], [264, 228, 229, 310, 309, 228, None], [265, 229, 230, 311, 310, 229, None], [266, 232, 231, 312, 313, 232, None], [267, 231, 233, 314, 312, 231, None], [268, 233, 234, 315, 314, 233, None], [269, 234, 235, 316, 315, 234, None], [270, 235, 236, 317, 316, 235, None], [271, 236, 237, 318, 317, 236, None], [272, 237, 238, 321, 319, 318, 237], [273, 238, 242, 320, 321, 238, None], [274, 242, 243, 326, 322, 320, 242], [275, 321, 320, 323, 319, 321, None], [276, 320, 322, 335, 324, 323, 320], [277, 243, 258, 325, 326, 243, None], [278, 258, 259, 327, 325, 258, None], [279, 326, 325, 328, 322, 326, None], [280, 325, 327, 329, 328, 325, None], [281, 259, 262, 330, 327, 259, None], [282, 262, 263, 345, 331, 330, 262], [283, 327, 330, 332, 329, 327, None], [284, 330, 331, 355, 333, 332, 330], [285, 322, 328, 334, 335, 322, None], [286, 328, 329, 336, 334, 328, None], [287, 335, 334, 337, 324, 335, None], [288, 334, 336, 338, 337, 334, None], [289, 329, 332, 339, 336, 329, None], [290, 332, 333, 382, 340, 339, 332], [291, 336, 339, 341, 338, 336, None], [292, 339, 340, 391, 342, 341, 339], [293, 263, 343, 344, 345, 263, None], [294, 343, 276, 346, 344, 343, None], [295, 345, 344, 347, 331, 345, None], [296, 344, 346, 348, 347, 344, None], [297, 276, 349, 350, 346, 276, None], [298, 349, 277, 351, 350, 349, None], [299, 346, 350, 352, 348, 346, None], [300, 350, 351, 353, 352, 350, None], [301, 331, 347, 354, 355, 331, None], [302, 347, 348, 356, 354, 347, None], [303, 355, 354, 357, 333, 355, None], [304, 354, 356, 358, 357, 354, None], [305, 348, 352, 359, 356, 348, None], [306, 352, 353, 360, 359, 352, None], [307, 356, 359, 361, 358, 356, None], [308, 359, 360, 362, 361, 359, None], [309, 277, 363, 364, 351, 277, None], [310, 363, 280, 365, 364, 363, None], [311, 351, 364, 366, 353, 351, None], [312, 364, 365, 367, 366, 364, None], [313, 280, 368, 369, 365, 280, None], [314, 368, 281, 370, 369, 368, None], [315, 365, 369, 371, 367, 365, None], [316, 369, 370, 372, 371, 369, None], [317, 353, 366, 373, 360, 353, None], [318, 366, 367, 374, 373, 366, None], [319, 360, 373, 375, 362, 360, None], [320, 373, 374, 376, 375, 373, None], [321, 367, 371, 377, 374, 367, None], [322, 371, 372, 378, 377, 371, None], [323, 374, 377, 379, 376, 374, None], [324, 377, 378, 380, 379, 377, None], [325, 333, 357, 381, 382, 333, None], [326, 357, 358, 383, 381, 357, None], [327, 382, 381, 384, 340, 382, None], [328, 381, 383, 385, 384, 381, None], [329, 358, 361, 386, 383, 358, None], [330, 361, 362, 387, 386, 361, None], [331, 383, 386, 388, 385, 383, None], [332, 386, 387, 389, 388, 386, None], [333, 340, 384, 390, 391, 340, None], [334, 384, 385, 392, 390, 384, None], [335, 391, 390, 393, 342, 391, None], [336, 390, 392, 394, 393, 390, None], [337, 385, 388, 395, 392, 385, None], [338, 388, 389, 396, 395, 388, None], [339, 392, 395, 397, 394, 392, None], [340, 395, 396, 398, 397, 395, None], [341, 362, 375, 399, 387, 362, None], [342, 375, 376, 400, 399, 375, None], [343, 387, 399, 401, 389, 387, None], [344, 399, 400, 402, 401, 399, None], [345, 376, 379, 403, 400, 376, None], [346, 379, 380, 404, 403, 379, None], [347, 400, 403, 405, 402, 400, None], [348, 403, 404, 406, 405, 403, None], [349, 389, 401, 407, 396, 389, None], [350, 401, 402, 408, 407, 401, None], [351, 396, 407, 409, 398, 396, None], [352, 407, 408, 410, 409, 407, None], [353, 402, 405, 411, 408, 402, None], [354, 405, 406, 412, 411, 405, None], [355, 408, 411, 413, 410, 408, None], [356, 411, 412, 414, 413, 411, None], [357, 281, 294, 415, 372, 370, 281], [358, 294, 295, 416, 415, 294, None], [359, 372, 415, 417, 380, 378, 372], [360, 415, 416, 418, 417, 415, None], [361, 295, 298, 419, 416, 295, None], [362, 298, 299, 420, 419, 298, None], [363, 416, 419, 421, 418, 416, None], [364, 419, 420, 422, 421, 419, None], [365, 380, 417, 423, 406, 404, 380], [366, 417, 418, 424, 423, 417, None], [367, 406, 423, 425, 414, 412, 406], [368, 423, 424, 426, 425, 423, None], [369, 418, 421, 427, 424, 418, None], [370, 421, 422, 428, 427, 421, None], [371, 424, 427, 429, 426, 424, None], [372, 427, 428, 430, 429, 427, None], [373, 299, 302, 431, 422, 420, 299], [374, 302, 303, 432, 431, 302, None], [375, 422, 431, 433, 430, 428, 422], [376, 431, 432, 434, 433, 431, None], [377, 303, 304, 435, 434, 432, 303], [378, 304, 305, 436, 435, 304, None], [379, 305, 306, 437, 436, 305, None], [380, 306, 307, 438, 437, 306, None], [381, 307, 308, 439, 438, 307, None], [382, 308, 309, 440, 439, 308, None], [383, 309, 310, 441, 440, 309, None], [384, 310, 311, 442, 441, 310, None], [385, 313, 312, 443, 444, 313, None], [386, 312, 314, 445, 443, 312, None], [387, 314, 315, 446, 445, 314, None], [388, 315, 316, 447, 446, 315, None], [389, 316, 317, 448, 447, 316, None], [390, 317, 318, 449, 448, 317, None], [391, 318, 319, 452, 450, 449, 318], [392, 319, 323, 451, 452, 319, None], [393, 323, 324, 457, 453, 451, 323], [394, 452, 451, 454, 450, 452, None], [395, 451, 453, 466, 455, 454, 451], [396, 324, 337, 456, 457, 324, None], [397, 337, 338, 458, 456, 337, None], [398, 457, 456, 459, 453, 457, None], [399, 456, 458, 460, 459, 456, None], [400, 338, 341, 461, 458, 338, None], [401, 341, 342, 462, 461, 341, None], [402, 458, 461, 463, 460, 458, None], [403, 461, 462, 464, 463, 461, None], [404, 453, 459, 465, 466, 453, None], [405, 459, 460, 467, 465, 459, None], [406, 466, 465, 468, 455, 466, None], [407, 465, 467, 469, 468, 465, None], [408, 460, 463, 470, 467, 460, None], [409, 463, 464, 471, 470, 463, None], [410, 467, 470, 472, 469, 467, None], [411, 470, 471, 473, 472, 470, None], [412, 342, 393, 394, 474, 462, 342], [413, 394, 397, 398, 475, 474, 394], [414, 462, 474, 476, 464, 462, None], [415, 474, 475, 477, 476, 474, None], [416, 398, 409, 410, 478, 475, 398], [417, 410, 413, 414, 479, 478, 410], [418, 475, 478, 480, 477, 475, None], [419, 478, 479, 481, 480, 478, None], [420, 464, 476, 482, 471, 464, None], [421, 476, 477, 483, 482, 476, None], [422, 471, 482, 484, 473, 471, None], [423, 482, 483, 485, 484, 482, None], [424, 477, 480, 486, 483, 477, None], [425, 480, 481, 487, 486, 480, None], [426, 483, 486, 488, 485, 483, None], [427, 486, 487, 489, 488, 486, None], [428, 414, 425, 490, 479, 414, None], [429, 425, 426, 491, 490, 425, None], [430, 479, 490, 492, 481, 479, None], [431, 490, 491, 493, 492, 490, None], [432, 426, 429, 494, 491, 426, None], [433, 429, 430, 495, 494, 429, None], [434, 491, 494, 496, 493, 491, None], [435, 494, 495, 497, 496, 494, None], [436, 481, 492, 498, 487, 481, None], [437, 492, 493, 499, 498, 492, None], [438, 487, 498, 500, 489, 487, None], [439, 498, 499, 501, 500, 498, None], [440, 493, 496, 502, 499, 493, None], [441, 496, 497, 503, 502, 496, None], [442, 499, 502, 504, 501, 499, None], [443, 502, 503, 505, 504, 502, None], [444, 430, 433, 506, 497, 495, 430], [445, 433, 434, 507, 506, 433, None], [446, 497, 506, 508, 505, 503, 497], [447, 506, 507, 509, 508, 506, None], [448, 434, 435, 510, 509, 507, 434], [449, 435, 436, 511, 510, 435, None], [450, 436, 437, 512, 511, 436, None], [451, 437, 438, 513, 512, 437, None], [452, 438, 439, 514, 513, 438, None], [453, 439, 440, 515, 514, 439, None], [454, 440, 441, 516, 515, 440, None], [455, 441, 442, 517, 516, 441, None], [456, 444, 443, 518, 519, 444, None], [457, 443, 445, 520, 518, 443, None], [458, 445, 446, 521, 520, 445, None], [459, 446, 447, 522, 521, 446, None], [460, 447, 448, 523, 522, 447, None], [461, 448, 449, 524, 523, 448, None], [462, 449, 450, 527, 525, 524, 449], [463, 450, 454, 526, 527, 450, None], [464, 454, 455, 528, 526, 454, None], [465, 527, 526, 529, 525, 527, None], [466, 526, 528, 530, 529, 526, None], [467, 455, 468, 469, 531, 528, 455], [468, 469, 472, 473, 532, 531, 469], [469, 528, 531, 533, 530, 528, None], [470, 531, 532, 534, 533, 531, None], [471, 473, 484, 485, 535, 532, 473], [472, 485, 488, 489, 536, 535, 485], [473, 532, 535, 537, 534, 532, None], [474, 535, 536, 538, 537, 535, None], [475, 489, 500, 501, 539, 536, 489], [476, 501, 504, 505, 540, 539, 501], [477, 536, 539, 541, 538, 536, None], [478, 539, 540, 542, 541, 539, None], [479, 505, 508, 543, 540, 505, None], [480, 508, 509, 544, 543, 508, None], [481, 540, 543, 545, 542, 540, None], [482, 543, 544, 546, 545, 543, None], [483, 509, 510, 547, 546, 544, 509], [484, 510, 511, 548, 547, 510, None], [485, 511, 512, 549, 548, 511, None], [486, 512, 513, 550, 549, 512, None], [487, 513, 514, 551, 550, 513, None], [488, 514, 515, 552, 551, 514, None], [489, 515, 516, 553, 552, 515, None], [490, 516, 517, 554, 553, 516, None], [491, 519, 518, 555, 556, 519, None], [492, 518, 520, 557, 555, 518, None], [493, 520, 521, 558, 557, 520, None], [494, 521, 522, 559, 558, 521, None], [495, 522, 523, 560, 559, 522, None], [496, 523, 524, 561, 560, 523, None], [497, 524, 525, 562, 561, 524, None], [498, 525, 529, 530, 563, 562, 525], [499, 530, 533, 534, 564, 563, 530], [500, 534, 537, 538, 565, 564, 534], [501, 538, 541, 542, 566, 565, 538], [502, 542, 545, 546, 567, 566, 542], [503, 546, 547, 568, 567, 546, None], [504, 547, 548, 569, 568, 547, None], [505, 548, 549, 570, 569, 548, None], [506, 549, 550, 571, 570, 549, None], [507, 550, 551, 572, 571, 550, None], [508, 551, 552, 573, 572, 551, None], [509, 552, 553, 574, 573, 552, None], [510, 553, 554, 575, 574, 553, None], [511, 556, 555, 576, 577, 556, None], [512, 555, 557, 578, 576, 555, None], [513, 557, 558, 579, 578, 557, None], [514, 558, 559, 580, 579, 558, None], [515, 559, 560, 581, 580, 559, None], [516, 560, 561, 582, 581, 560, None], [517, 561, 562, 583, 582, 561, None], [518, 562, 563, 584, 583, 562, None], [519, 563, 564, 585, 584, 563, None], [520, 564, 565, 586, 585, 564, None], [521, 565, 566, 587, 586, 565, None], [522, 566, 567, 588, 587, 566, None], [523, 567, 568, 589, 588, 567, None], [524, 568, 569, 590, 589, 568, None], [525, 569, 570, 591, 590, 569, None], [526, 570, 571, 592, 591, 570, None], [527, 571, 572, 593, 592, 571, None], [528, 572, 573, 594, 593, 572, None], [529, 573, 574, 595, 594, 573, None], [530, 574, 575, 596, 595, 574, None], [531, 577, 576, 597, 598, 577, None], [532, 576, 578, 599, 597, 576, None], [533, 578, 579, 600, 599, 578, None], [534, 579, 580, 601, 600, 579, None], [535, 580, 581, 602, 601, 580, None], [536, 581, 582, 603, 602, 581, None], [537, 582, 583, 604, 603, 582, None], [538, 583, 584, 605, 604, 583, None], [539, 584, 585, 606, 605, 584, None], [540, 585, 586, 607, 606, 585, None], [541, 586, 587, 608, 607, 586, None], [542, 587, 588, 609, 608, 587, None], [543, 588, 589, 610, 609, 588, None], [544, 589, 590, 611, 610, 589, None], [545, 590, 591, 612, 611, 590, None], [546, 591, 592, 613, 612, 591, None], [547, 592, 593, 614, 613, 592, None], [548, 593, 594, 615, 614, 593, None], [549, 594, 595, 616, 615, 594, None], [550, 595, 596, 617, 616, 595, None], [551, 598, 597, 618, 619, 598, None], [552, 597, 599, 620, 618, 597, None], [553, 599, 600, 621, 620, 599, None], [554, 600, 601, 622, 621, 600, None], [555, 601, 602, 623, 622, 601, None], [556, 602, 603, 624, 623, 602, None], [557, 603, 604, 625, 624, 603, None], [558, 604, 605, 626, 625, 604, None], [559, 605, 606, 627, 626, 605, None], [560, 606, 607, 628, 627, 606, None], [561, 607, 608, 629, 628, 607, None], [562, 608, 609, 630, 629, 608, None], [563, 609, 610, 631, 630, 609, None], [564, 610, 611, 632, 631, 610, None], [565, 611, 612, 633, 632, 611, None], [566, 612, 613, 634, 633, 612, None], [567, 613, 614, 635, 634, 613, None], [568, 614, 615, 636, 635, 614, None], [569, 615, 616, 637, 636, 615, None], [570, 616, 617, 638, 637, 616, None], [571, 619, 618, 639, 640, 619, None], [572, 618, 620, 641, 639, 618, None], [573, 620, 621, 642, 641, 620, None], [574, 621, 622, 643, 642, 621, None], [575, 622, 623, 644, 643, 622, None], [576, 623, 624, 645, 644, 623, None], [577, 624, 625, 646, 645, 624, None], [578, 625, 626, 647, 646, 625, None], [579, 626, 627, 648, 647, 626, None], [580, 627, 628, 649, 648, 627, None], [581, 628, 629, 650, 649, 628, None], [582, 629, 630, 651, 650, 629, None], [583, 630, 631, 652, 651, 630, None], [584, 631, 632, 653, 652, 631, None], [585, 632, 633, 654, 653, 632, None], [586, 633, 634, 655, 654, 633, None], [587, 634, 635, 656, 655, 634, None], [588, 635, 636, 657, 656, 635, None], [589, 636, 637, 658, 657, 636, None], [590, 637, 638, 659, 658, 637, None], [591, 640, 639, 660, 661, 640, None], [592, 639, 641, 662, 660, 639, None], [593, 641, 642, 663, 662, 641, None], [594, 642, 643, 664, 663, 642, None], [595, 643, 644, 665, 664, 643, None], [596, 644, 645, 666, 665, 644, None], [597, 645, 646, 667, 666, 645, None], [598, 646, 647, 668, 667, 646, None], [599, 647, 648, 669, 668, 647, None], [600, 648, 649, 670, 669, 648, None], [601, 649, 650, 671, 670, 649, None], [602, 650, 651, 672, 671, 650, None], [603, 651, 652, 673, 672, 651, None], [604, 652, 653, 674, 673, 652, None], [605, 653, 654, 675, 674, 653, None], [606, 654, 655, 676, 675, 654, None], [607, 655, 656, 677, 676, 655, None], [608, 656, 657, 678, 677, 656, None], [609, 657, 658, 679, 678, 657, None], [610, 658, 659, 680, 679, 658, None], [611, 661, 660, 681, 682, 661, None], [612, 660, 662, 683, 681, 660, None], [613, 662, 663, 684, 683, 662, None], [614, 663, 664, 685, 684, 663, None], [615, 664, 665, 686, 685, 664, None], [616, 665, 666, 687, 686, 665, None], [617, 666, 667, 688, 687, 666, None], [618, 667, 668, 689, 688, 667, None], [619, 668, 669, 690, 689, 668, None], [620, 669, 670, 691, 690, 669, None], [621, 670, 671, 692, 691, 670, None], [622, 671, 672, 693, 692, 671, None], [623, 672, 673, 694, 693, 672, None], [624, 673, 674, 695, 694, 673, None], [625, 674, 675, 696, 695, 674, None], [626, 675, 676, 697, 696, 675, None], [627, 676, 677, 698, 697, 676, None], [628, 677, 678, 699, 698, 677, None], [629, 678, 679, 700, 699, 678, None], [630, 679, 680, 701, 700, 679, None], [631, 682, 681, 702, 703, 682, None], [632, 681, 683, 704, 702, 681, None], [633, 683, 684, 705, 704, 683, None], [634, 684, 685, 706, 705, 684, None], [635, 685, 686, 707, 706, 685, None], [636, 686, 687, 708, 707, 686, None], [637, 687, 688, 709, 708, 687, None], [638, 688, 689, 710, 709, 688, None], [639, 689, 690, 711, 710, 689, None], [640, 690, 691, 712, 711, 690, None], [641, 691, 692, 713, 712, 691, None], [642, 692, 693, 714, 713, 692, None], [643, 693, 694, 715, 714, 693, None], [644, 694, 695, 716, 715, 694, None], [645, 695, 696, 717, 716, 695, None], [646, 696, 697, 718, 717, 696, None], [647, 697, 698, 719, 718, 697, None], [648, 698, 699, 720, 719, 698, None], [649, 699, 700, 721, 720, 699, None], [650, 700, 701, 722, 721, 700, None]]\n" + ] + } + ], + "source": [ + "# load the simulation and get the model\n", + "vertex_sim_name = \"mfsim.nam\"\n", + "vertex_sim = flopy.mf6.MFSimulation.load(sim_name=vertex_sim_name, exe_name=\"mf6\", \n", + " sim_ws=modelpth)\n", + "vertex_ml6 = vertex_sim.get_model(\"mp7p2\")\n", + "\n", + "# load the MODPATH-7 results\n", + "mp_namea = 'mp7p2a_mp'\n", + "fpth = os.path.join(modelpth, mp_namea + \".mppth\")\n", + "p = flopy.utils.PathlineFile(fpth)\n", + "p0 = p.get_alldata()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Create the `Vtk()` object and add all of the model data to it" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [], + "source": [ + "vtkobj = vtk.Vtk(vertex_ml6, xml=True, pvd=True, vertical_exageration=10)\n", + "vtkobj.add_model(vertex_ml6)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Add modpath data to the `Vtk()` object.\n", + "\n", + "*Note: this will create a second vtk file that has the file signiture (myfilename)_pathline.vtu*" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [], + "source": [ + "vtkobj.add_pathline_points(p0, timeseries=False)\n", + "vtkobj.write(os.path.join(output_dir, \"mp7_vertex_example\", \"vertex_ex.vtu\"))" ] }, { @@ -388,7 +836,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.10" + "version": "3.7.4" } }, "nbformat": 4, From 6b9769e203691f1b4f8ad2f497a316a0a40d37c1 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 28 Sep 2021 12:39:25 -0700 Subject: [PATCH 2/5] update(vtk.py): rewrite of Vtk() class and functions for Vertex and Unstructured Grids, MODPATH * New dependency added (python vtk package: pip install vtk) * Added support for HFB package * Added support for MODPATH pathline and timeseries output * Updated the support for head exports and cell budget files --- autotest/t050_test.py | 1037 ++++------- etc/environment.yml | 1 + flopy/export/utils.py | 141 +- flopy/export/vtk.py | 4079 +++++++++++++++++++---------------------- 4 files changed, 2366 insertions(+), 2892 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 2fa3909d9..9d9501063 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -2,13 +2,9 @@ import numpy as np import os import flopy -from flopy.export import vtk +from flopy.export.vtk import Vtk # Test vtk export -# Note: initially thought about asserting that exported file size in bytes is -# unchanged, but this seems to be sensitive to the running environment. -# Thus, only asserting that the number of lines is unchanged. -# Still keeping the file size check commented for development purposes. # create output directory cpth = os.path.join("temp", "t050") @@ -29,6 +25,10 @@ def count_lines_in_file(filepath, binary=False): def test_vtk_export_array2d(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg mpath = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" @@ -40,25 +40,24 @@ def test_vtk_export_array2d(): output_dir = os.path.join(cpth, "array_2d_test") # export and check - m.dis.top.export(output_dir, name="top", fmt="vtk") - filetocheck = os.path.join(output_dir, "top.vtu") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==351846) + m.dis.top.export(output_dir, name="top", fmt="vtk", binary=False) + filetocheck = os.path.join(output_dir, "top.vtk") nlines = count_lines_in_file(filetocheck) - assert nlines == 2846 + assert nlines == 17615 # with smoothing - m.dis.top.export(output_dir, fmt="vtk", name="top_smooth", smooth=True) - filetocheck = os.path.join(output_dir, "top_smooth.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==351715) + m.dis.top.export(output_dir, fmt="vtk", name="top_smooth", + binary=False, smooth=True) + filetocheck = os.path.join(output_dir, "top_smooth.vtk") nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 2846 - - return + assert nlines1 == 17615 def test_vtk_export_array3d(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg mpath = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" @@ -73,22 +72,18 @@ def test_vtk_export_array3d(): output_dir = os.path.join(cpth, "array_3d_test") # export and check - m.upw.hk.export(output_dir, fmt="vtk", name="hk") - filetocheck = os.path.join(output_dir, "hk.vtu") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==992036) + m.upw.hk.export(output_dir, fmt="vtk", name="hk", binary=False) + filetocheck = os.path.join(output_dir, "hk.vtk") nlines = count_lines_in_file(filetocheck) - assert nlines == 8486 + assert nlines == 17615 # with point scalars m.upw.hk.export( - output_dir, fmt="vtk", name="hk_points", point_scalars=True + output_dir, fmt="vtk", name="hk_points", point_scalars=True, binary=False ) - filetocheck = os.path.join(output_dir, "hk_points.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==1320666) + filetocheck = os.path.join(output_dir, "hk_points.vtk") nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 10605 + assert nlines1 == 19482 # with point scalars and binary m.upw.hk.export( @@ -96,19 +91,16 @@ def test_vtk_export_array3d(): fmt="vtk", name="hk_points_bin", point_scalars=True, - binary=True, ) - filetocheck = os.path.join(output_dir, "hk_points_bin.vtu") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==629401) - # nlines2 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines2==2105) + filetocheck = os.path.join(output_dir, "hk_points_bin.vtk") assert os.path.exists(filetocheck) - return - def test_vtk_transient_array_2d(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg mpath = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" @@ -125,37 +117,27 @@ def test_vtk_transient_array_2d(): kpers = [0, 1, 1096] # export and check - m.rch.rech.export(output_dir, fmt="vtk", kpers=kpers) - filetocheck = os.path.join(output_dir, "rech_01.vtu") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==355144) + m.rch.rech.export(output_dir, fmt="vtk", kpers=kpers, binary=False, xml=True) + filetocheck = os.path.join(output_dir, "rech_000001.vtk") nlines = count_lines_in_file(filetocheck) - assert nlines == 2851 - filetocheck = os.path.join(output_dir, "rech_01097.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==354442) + assert nlines == 26837 + filetocheck = os.path.join(output_dir, "rech_001096.vtk") nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 2851 + assert nlines1 == 26837 # with binary m.rch.rech.export(output_dir_bin, fmt="vtk", binary=True, kpers=kpers) - filetocheck = os.path.join(output_dir_bin, "rech_01.vtu") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==168339) - # nlines2 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines2==846) + filetocheck = os.path.join(output_dir_bin, "rech_000001.vtk") assert os.path.exists(filetocheck) - filetocheck = os.path.join(output_dir_bin, "rech_01097.vtu") - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==168339) - # nlines3 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines3==846) + filetocheck = os.path.join(output_dir_bin, "rech_001096.vtk") assert os.path.exists(filetocheck) - return - def test_vtk_export_packages(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg mpath = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" @@ -170,70 +152,59 @@ def test_vtk_export_packages(): # dis export and check output_dir = os.path.join(cpth, "DIS") - m.dis.export(output_dir, fmt="vtk") - filetocheck = os.path.join(output_dir, "DIS.vtu") + # todo: pakbase.export() for vtk!!!! + m.dis.export(output_dir, fmt="vtk", xml=True, binary=False) + filetocheck = os.path.join(output_dir, "DIS.vtk") # totalbytes = os.path.getsize(filetocheck) # assert(totalbytes==1019857) nlines = count_lines_in_file(filetocheck) - assert nlines == 8491, f"nlines ({nlines}) not equal to 8491" + assert nlines == 27239, f"nlines ({nlines}) not equal to 27239" # upw with point scalar output output_dir = os.path.join(cpth, "UPW") - m.upw.export(output_dir, fmt="vtk", point_scalars=True) - filetocheck = os.path.join(output_dir, "UPW.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==2559173) + m.upw.export( + output_dir, fmt="vtk", xml=True, binary=False, point_scalars=True + ) + filetocheck = os.path.join(output_dir, "UPW.vtk") nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 21215, f"nlines ({nlines}) not equal to 21215" + assert nlines1 == 42445, f"nlines ({nlines}) not equal to 42445" # bas with smoothing on output_dir = os.path.join(cpth, "BAS") - m.bas6.export(output_dir, fmt="vtk", smooth=True) - filetocheck = os.path.join(output_dir, "BAS6.vtu") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==1001580) + m.bas6.export(output_dir, fmt="vtk", binary=False, smooth=True) + filetocheck = os.path.join(output_dir, "BAS6.vtk") nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 8491 + assert nlines2 == 17883 # transient package drain kpers = [0, 1, 1096] output_dir = os.path.join(cpth, "DRN") - m.drn.export(output_dir, fmt="vtk", kpers=kpers) - filetocheck = os.path.join(output_dir, "DRN_01.vtu") - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==20670) + m.drn.export(output_dir, fmt="vtk", binary=False, xml=True, kpers=kpers, pvd=True) + filetocheck = os.path.join(output_dir, "DRN_000001.vtu") nlines3 = count_lines_in_file(filetocheck) - assert nlines3 == 191 - filetocheck = os.path.join(output_dir, "DRN_01097.vtu") - # totalbytes4 = os.path.getsize(filetocheck) - # assert(totalbytes4==20670) + assert nlines3 == 27239 + filetocheck = os.path.join(output_dir, "DRN_001096.vtu") nlines4 = count_lines_in_file(filetocheck) - assert nlines4 == 191 + assert nlines4 == 27239 # dis with binary output_dir = os.path.join(cpth, "DIS_bin") m.dis.export(output_dir, fmt="vtk", binary=True) - filetocheck = os.path.join(output_dir, "DIS.vtu") - # totalbytes5 = os.path.getsize(filetocheck) - # assert(totalbytes5==519516) - # nlines5 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines5==1797) + filetocheck = os.path.join(output_dir, "DIS.vtk") assert os.path.exists(filetocheck) # upw with point scalars and binary output_dir = os.path.join(cpth, "UPW_bin") m.upw.export(output_dir, fmt="vtk", point_scalars=True, binary=True) - filetocheck = os.path.join(output_dir, "UPW.vtu") - # totalbytes6 = os.path.getsize(filetocheck) - # assert(totalbytes6==1349801) - # nlines6 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines6==4240) + filetocheck = os.path.join(output_dir, "UPW.vtk") assert os.path.exists(filetocheck) - return - def test_vtk_mf6(): + try: + import vtk + except ImportError: + return # test mf6 mf6expth = os.path.join("..", "examples", "data", "mf6") mf6sims = [ @@ -252,183 +223,117 @@ def test_vtk_mf6(): for mname in sim_models: print(mname) m = loaded_sim.get_model(mname) - m.export(os.path.join(cpth, m.name), fmt="vtk") + m.export(os.path.join(cpth, m.name), fmt="vtk", binary=False) # check one - filetocheck = os.path.join(cpth, "twrihfb2015", "npf.vtr") + filetocheck = os.path.join(cpth, "twrihfb2015", "twrihfb2015_000000.vtk") # totalbytes = os.path.getsize(filetocheck) # assert(totalbytes==21609) nlines = count_lines_in_file(filetocheck) - assert nlines == 76 - - return + assert nlines == 9537 def test_vtk_binary_head_export(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg + from flopy.utils import HeadFile + mpth = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" ) namfile = "freyberg.nam" hdsfile = os.path.join(mpth, "freyberg.hds") + heads = HeadFile(hdsfile) m = flopy.modflow.Modflow.load( namfile, model_ws=mpth, verbose=False, load_only=["dis", "bas6"] ) - filenametocheck = "freyberg_head_KPER455_KSTP1.vtu" + filenametocheck = "freyberg_head_000003.vtu" # export and check 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)], + vtkobj = Vtk(m, pvd=True, xml=True) + vtkobj.add_heads( + heads, kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, 1089)] ) + vtkobj.write(os.path.join(otfolder, "freyberg_head")) + filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==993215) nlines = count_lines_in_file(filetocheck) - assert nlines == 8486 + assert nlines == 34 # with point scalars 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, + vtkobj = Vtk(m, pvd=True, xml=True, point_scalars=True) + vtkobj.add_heads( + heads, kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, 1089)] ) + vtkobj.write(os.path.join(otfolder, "freyberg_head")) + filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==1331858) nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 10605 + assert nlines1 == 34 # with smoothing 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)], - smooth=True, - nanval=-999.99, + vtkobj = Vtk(m, pvd=True, xml=True, smooth=True) + vtkobj.add_heads( + heads, kstpkper=[(0, 0), (0, 199), (0, 354), (0, 454), (0, 1089)] ) + vtkobj.write(os.path.join(otfolder, "freyberg_head")) + filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==993077) nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 8486 - - # with smoothing and binary - 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)], - smooth=True, - binary=True, - nanval=-999.99, - ) - filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==493853) - # nlines3 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines3==1781) - assert os.path.exists(filetocheck) - - # with smoothing and binary, single time - 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, - ) - filetocheck = os.path.join(otfolder, "freyberg_head_KPER1_KSTP1.vtu") - # totalbytes4 = os.path.getsize(filetocheck) - # assert(totalbytes4==493853) - # nlines4 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines4==1787) - assert os.path.exists(filetocheck) - - return + assert nlines2 == 34 def test_vtk_cbc(): + try: + import vtk + except ImportError: + return # test mf 2005 freyberg + from flopy.utils import CellBudgetFile + mpth = os.path.join( "..", "examples", "data", "freyberg_multilayer_transient" ) namfile = "freyberg.nam" cbcfile = os.path.join(mpth, "freyberg.cbc") + cbc = CellBudgetFile(cbcfile) m = flopy.modflow.Modflow.load( namfile, model_ws=mpth, verbose=False, load_only=["dis", "bas6"] ) - filenametocheck = "freyberg_CBC_KPER1_KSTP1.vtu" + filenametocheck = "freyberg_CBC_000000.vtu" # export and check with point scalar otfolder = os.path.join(cpth, "freyberg_CBCTEST") - vtk.export_cbc( - m, - cbcfile, - otfolder, - kstpkper=[(0, 0), (0, 1), (0, 2)], - point_scalars=True, - ) + + vtkobj = Vtk(m, binary=False, xml=True, pvd=True, point_scalars=True) + vtkobj.add_cell_budget(cbc, kstpkper=[(0, 0), (0, 1), (0, 2)]) + vtkobj.write(os.path.join(otfolder, "freyberg_CBC")) + filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==2630875) nlines = count_lines_in_file(filetocheck) - assert nlines == 19093 + assert nlines == 39243 # with point scalars and binary otfolder = os.path.join(cpth, "freyberg_CBCTEST_bin") - vtk.export_cbc( - m, - cbcfile, - otfolder, - kstpkper=[(0, 0), (0, 1), (0, 2)], - point_scalars=True, - binary=True, - ) + vtkobj = Vtk(m, xml=True, pvd=True, point_scalars=True) + vtkobj.add_cell_budget(cbc, kstpkper=[(0, 0), (0, 1), (0, 2)]) + vtkobj.write(os.path.join(otfolder, "freyberg_CBC")) filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==1205818) - # nlines1 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines1==3088) assert os.path.exists(filetocheck) - # with point scalars and binary, only one budget component - otfolder = os.path.join(cpth, "freyberg_CBCTEST_bin2") - vtk.export_cbc( - m, - cbcfile, - otfolder, - kstpkper=(0, 0), - text="CONSTANT HEAD", - point_scalars=True, - binary=True, - ) - filetocheck = os.path.join(otfolder, filenametocheck) - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==10142) - # nlines2 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines2==66) - assert os.path.exists(filetocheck) - - return - def test_vtk_vector(): + try: + import vtk + except ImportError: + return + from flopy.utils import postprocessing as pp from flopy.utils import HeadFile, CellBudgetFile @@ -452,519 +357,261 @@ def test_vtk_vector(): filenametocheck = "discharge.vtu" # export and check with point scalar - vtk.export_vector(m, q, output_dir, "discharge", point_scalars=True) + vtkobj = Vtk(m, xml=True, binary=False, point_scalars=True) + vtkobj.add_vector(q, "discharge") + vtkobj.write(os.path.join(output_dir, filenametocheck)) + filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==2247857) nlines = count_lines_in_file(filetocheck) - assert nlines == 10605 + assert nlines == 36045 # with point scalars and binary - vtk.export_vector( - m, q, f"{output_dir}_bin", "discharge", point_scalars=True, binary=True - ) + vtkobj = Vtk(m, point_scalars=True) + vtkobj.add_vector(q, "discharge") + vtkobj.write(os.path.join(f"{output_dir}_bin", filenametocheck)) filetocheck = os.path.join(f"{output_dir}_bin", filenametocheck) - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==942413) - # nlines1 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines1==3824) assert os.path.exists( filetocheck ), f"file (0) does not exist: {filetocheck}" - # with values directly given at vertices - q = pp.get_specific_discharge(vectors, m, head, position="vertices") - nancount = np.count_nonzero(np.isnan(q[0])) - assert nancount == 472, f"nancount != 472 ({nancount})" - overall = np.nansum(q[0]) + np.nansum(q[1]) + np.nansum(q[2]) - assert np.allclose( - overall, -15.849639024891047 - ), f"vertices overall = {overall}" + # test at cell centers + q = pp.get_specific_discharge(vectors, m, head) + output_dir = os.path.join(cpth, "freyberg_vector") filenametocheck = "discharge_verts.vtu" - vtk.export_vector(m, q, output_dir, "discharge_verts") + vtkobj = Vtk(m, xml=True, binary=False) + vtkobj.add_vector(q, "discharge") + vtkobj.write(os.path.join(output_dir, filenametocheck)) + filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==1990047) nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 10598, f"nlines != 10598 ({nlines2})" + assert nlines2 == 27645, f"nlines != 10598 ({nlines2})" # with values directly given at vertices and binary - vtk.export_vector( - m, q, f"{output_dir}_bin", "discharge_verts", binary=True - ) + vtkobj = Vtk(m, xml=True, binary=False) + vtkobj.add_vector(q, "discharge") + vtkobj.write(os.path.join(f"{output_dir}_bin", filenametocheck)) + filetocheck = os.path.join(f"{output_dir}_bin", filenametocheck) - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==891486) - # nlines3 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines3==3012) assert os.path.exists( filetocheck ), f"file (1) does not exist: {filetocheck}" - return - - -def test_vtk_vti(): - # create model with regular and equal grid spacing in x, y and z directions - name = "test_vti" - m = flopy.modflow.Modflow(name) - nlay, nrow, ncol = 2, 3, 4 - delr = np.ones(ncol) - delc = np.ones(nrow) - top = 2.0 * np.ones((nrow, ncol)) - botm1 = np.ones((1, nrow, ncol)) - botm2 = np.zeros((1, nrow, ncol)) - botm = np.concatenate((botm1, botm2)) - dis = flopy.modflow.ModflowDis( - m, nlay, nrow, ncol, delr=delr, delc=delc, top=top, botm=botm - ) - output_dir = os.path.join(cpth, m.name) - filenametocheck = "DIS.vti" - - # export and check - dis.export(output_dir, fmt="vtk") - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==1075) - nlines = count_lines_in_file(filetocheck) - assert nlines == 17, f"nlines ({nlines}) not equal to 17" - - # with point scalar - dis.export(f"{output_dir}_points", fmt="vtk", point_scalars=True) - filetocheck = os.path.join(f"{output_dir}_points", filenametocheck) - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==2474) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 29, f"nlines1 ({nlines1}) not equal to 29" - - # with binary - dis.export(f"{output_dir}_bin", fmt="vtk", binary=True) - filetocheck = os.path.join(f"{output_dir}_bin", filenametocheck) - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==1144) - # nlines2 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines2==18) - assert os.path.exists(filetocheck) - - # force .vtr - filenametocheck = "DIS.vtr" - dis.export(output_dir, fmt="vtk", vtk_grid_type="RectilinearGrid") - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==1606) - nlines3 = count_lines_in_file(filetocheck) - assert nlines3 == 37, f"nlines3 ({nlines3}) not equal to 37" - - # force .vtu - filenametocheck = "DIS.vtu" - dis.export(output_dir, fmt="vtk", vtk_grid_type="UnstructuredGrid") - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes4 = os.path.getsize(filetocheck) - # assert(totalbytes4==5723) - nlines4 = count_lines_in_file(filetocheck) - assert nlines4 == 125, f"nlines4 ({nlines4}) not equal to 125" - - # vector - filenametocheck = "vect.vti" - ones_array = np.ones(m.modelgrid.shape) - v = (ones_array, 2.0 * ones_array, 3.0 * ones_array) - vtk.export_vector(m, v, output_dir, "vect", point_scalars=True) - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes5 = os.path.getsize(filetocheck) - # assert(totalbytes5==1578) - nlines5 = count_lines_in_file(filetocheck) - assert nlines5 == 20 - - # vector with point scalars and binary - vtk.export_vector( - m, v, f"{output_dir}_bin", "vect", point_scalars=True, binary=True - ) - filetocheck = os.path.join(f"{output_dir}_bin", filenametocheck) - # totalbytes6 = os.path.getsize(filetocheck) - # assert(totalbytes6==2666) - # nlines6 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines6==18) - assert os.path.exists(filetocheck) - - return - - -def test_vtk_vtr(): - # test mf 2005 l1a2k - mpth = os.path.join("..", "examples", "data", "mf2005_test") - namfile = "l1a2k.nam" - m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False) - output_dir = os.path.join(cpth, m.name) - filenametocheck = "EVT_01.vtr" - - # export and check - m.export(output_dir, fmt="vtk") - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==79953) - nlines = count_lines_in_file(filetocheck) - assert nlines == 87 - - # with point scalar - m.export(f"{output_dir}_points", fmt="vtk", point_scalars=True) - filetocheck = os.path.join(f"{output_dir}_points", filenametocheck) - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==182168) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 121 - - # with binary - m.export(f"{output_dir}_bin", fmt="vtk", binary=True) - filetocheck = os.path.join(f"{output_dir}_bin", filenametocheck) - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==47874) - # nlines2 = count_lines_in_file(filetocheck, binary=True) - # assert(nlines2==28) - assert os.path.exists(filetocheck) - - # force .vtu - filenametocheck = "EVT_01.vtu" - m.export(output_dir, fmt="vtk", vtk_grid_type="UnstructuredGrid") - filetocheck = os.path.join(output_dir, filenametocheck) - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==78762) - nlines3 = count_lines_in_file(filetocheck) - assert nlines3 == 1105 - - return - - -def test_vtk_export_true2d_regular(): - mpath = os.path.join("..", "examples", "data", "mf2005_test") - output_dir = os.path.join(cpth, "true2d_regular") - - # test mf 2005 test1ss, which has one layer with non-constant elevations - namfile = "test1ss.nam" - m = flopy.modflow.Modflow.load( - namfile, model_ws=mpath, verbose=False, load_only=["dis", "bas6"] - ) - - # export and check (.vti, with point scalars) - m.dis.botm.export( - output_dir, - name="test1ss_botm", - fmt="vtk", - point_scalars=True, - true2d=True, - ) - filetocheck = os.path.join(output_dir, "test1ss_botm.vti") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==3371) - nlines = count_lines_in_file(filetocheck) - assert nlines == 32 - - # vector (.vti, with point scalars) - vect = (m.dis.botm.array, m.dis.botm.array) - vtk.export_vector( - m, - vect, - output_dir, - "test1ss_botm_vect", - point_scalars=True, - true2d=True, - ) - filetocheck = os.path.join(output_dir, "test1ss_botm_vect.vti") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==6022) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 32 - - # vector directly at vertices (.vti) - vect = [m.dis.botm.array, m.dis.botm.array] - for i, vcomp in enumerate(vect): - vect[i] = m.modelgrid.array_at_verts(vcomp) - vtk.export_vector(m, vect, output_dir, "test1ss_botm_vectv", true2d=True) - filetocheck = os.path.join(output_dir, "test1ss_botm_vectv.vti") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==3496) - nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 27 - - # export and check (force .vtu, with point scalars) - m.dis.botm.export( - output_dir, - name="test1ss_botm", - fmt="vtk", - point_scalars=True, - vtk_grid_type="UnstructuredGrid", - true2d=True, - ) - filetocheck = os.path.join(output_dir, "test1ss_botm.vtu") - # totalbytes3 = os.path.getsize(filetocheck) - # assert(totalbytes3==23827) - nlines3 = count_lines_in_file(filetocheck) - assert nlines3 == 608 - - # test mf 2005 swiex3, which has one row - namfile = "swiex3.nam" - m = flopy.modflow.Modflow.load( - namfile, model_ws=mpath, verbose=False, load_only=["dis", "bas6"] - ) - - # export and check (.vtr) - m.dis.botm.export(output_dir, name="swiex3_botm", fmt="vtk", true2d=True) - filetocheck = os.path.join(output_dir, "swiex3_botm.vtr") - # totalbytes4 = os.path.getsize(filetocheck) - # assert(totalbytes4==8022) - nlines4 = count_lines_in_file(filetocheck) - assert nlines4 == 229 - - # export and check (force .vtu) - m.dis.botm.export( - output_dir, - name="swiex3_botm", - fmt="vtk", - vtk_grid_type="UnstructuredGrid", - true2d=True, - ) - filetocheck = os.path.join(output_dir, "swiex3_botm.vtu") - # totalbytes5 = os.path.getsize(filetocheck) - # assert(totalbytes5==85446) - nlines5 = count_lines_in_file(filetocheck) - assert nlines5 == 2426 - - return - - -def test_vtk_export_true2d_nonregxy(): - import flopy.utils.binaryfile as bf - from flopy.utils import postprocessing as pp - output_dir = os.path.join(cpth, "true2d_nonregxy") - cbc_unit_nb = 53 - - # model with one layer, non-regular grid in x and y - name = "nonregxy" - m = flopy.modflow.Modflow(name, model_ws=output_dir, exe_name="mf2005") - nlay, nrow, ncol = 1, 10, 10 - delr = np.concatenate((np.ones((5,)), 2.0 * np.ones((5,)))) - delc = delr - top = 50.0 - botm = 0.0 - dis = flopy.modflow.ModflowDis( - m, nlay, nrow, ncol, delr=delr, delc=delc, top=top, botm=botm - ) - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - ibound[:, :, 0] = -1 - ibound[:, :, -1] = -1 - strt = np.linspace(1.0, 0.0, ncol).reshape(1, 1, ncol) - strt = strt * np.ones((nlay, nrow, ncol)) - bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt) - lpf = flopy.modflow.ModflowLpf(m, hk=1.0, vka=1.0, ipakcb=cbc_unit_nb) - spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} - oc = flopy.modflow.ModflowOc(m, stress_period_data=spd, compact=True) - pcg = flopy.modflow.ModflowPcg(m) - m.write_input() - m.run_model(silent=True) - - # export and check head with point scalar - hdsfile = os.path.join(output_dir, f"{name}.hds") - hds = bf.HeadFile(hdsfile) - head = hds.get_data() - vtk.export_array( - m, head, output_dir, f"{name}_head", point_scalars=True, true2d=True - ) - filetocheck = os.path.join(output_dir, f"{name}_head.vtr") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==4997) - nlines = count_lines_in_file(filetocheck) - assert nlines == 59 - - # export and check specific discharge given at vertices - cbcfile = os.path.join(output_dir, f"{name}.cbc") - cbc = bf.CellBudgetFile(cbcfile) - keys = ["FLOW RIGHT FACE", "FLOW FRONT FACE"] - vectors = [cbc.get_data(text=t)[0] for t in keys] - q = pp.get_specific_discharge(vectors, m, position="vertices") - vtk.export_vector( - m, q, output_dir, f"{name}_q", point_scalars=True, true2d=True - ) - filetocheck = os.path.join(output_dir, f"{name}_q.vtr") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==5772) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 54 - return - - -def test_vtk_export_true2d_nonregxz(): - import flopy.utils.binaryfile as bf - from flopy.utils import postprocessing as pp - - output_dir = os.path.join(cpth, "true2d_nonregxz") - cbc_unit_nb = 53 - - # model with one row, non-regular grid in x and stepwise z - name = "nonregxz" - m = flopy.modflow.Modflow(name, model_ws=output_dir, exe_name="mf2005") - nlay, nrow, ncol = 2, 1, 10 - delr = np.concatenate((np.ones((5,)), 2.0 * np.ones((5,)))) - delc = 1.0 - top = np.linspace(2.0, 3.0, ncol).reshape((1, 1, ncol)) - botm1 = np.linspace(1.0, 2.5, ncol).reshape((1, 1, ncol)) - botm2 = np.linspace(0.0, 0.5, ncol).reshape((1, 1, ncol)) - botm = np.concatenate((botm1, botm2)) - dis = flopy.modflow.ModflowDis( - m, nlay, nrow, ncol, delr=delr, delc=delc, top=top, botm=botm - ) - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - ibound[:, :, 0] = -1 - ibound[:, :, -1] = -1 - strt = np.linspace(0.0, 1.0, ncol).reshape(1, 1, ncol) - strt = strt * np.ones((nlay, nrow, ncol)) - bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt) - lpf = flopy.modflow.ModflowLpf(m, hk=1.0, vka=1.0, ipakcb=cbc_unit_nb) - spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} - oc = flopy.modflow.ModflowOc(m, stress_period_data=spd, compact=True) - pcg = flopy.modflow.ModflowPcg(m) - m.write_input() - m.run_model(silent=True) - - # export and check head - hdsfile = os.path.join(output_dir, f"{name}.hds") - hds = bf.HeadFile(hdsfile) - head = hds.get_data() - vtk.export_array(m, head, output_dir, f"{name}_head", true2d=True) - filetocheck = os.path.join(output_dir, f"{name}_head.vtu") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==4217) - nlines = count_lines_in_file(filetocheck) - assert nlines == 105 - - # export and check head with point scalar - hdsfile = os.path.join(output_dir, f"{name}.hds") - hds = bf.HeadFile(hdsfile) - head = hds.get_data() - vtk.export_array( - m, - head, - output_dir, - f"{name}_head_points", - point_scalars=True, - true2d=True, - ) - filetocheck = os.path.join(output_dir, f"{name}_head_points.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==6155) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 129 +def test_vtk_unstructured(): + try: + import vtk + from vtk.util import numpy_support + except ImportError: + return + + def load_verts(fname): + verts = np.genfromtxt( + fname, dtype=[int, float, float], names=["iv", "x", "y"] + ) + verts["iv"] -= 1 # zero based + return verts + + def load_iverts(fname): + f = open(fname, "r") + iverts = [] + xc = [] + yc = [] + for line in f: + ll = line.strip().split() + iverts.append([int(i) - 1 for i in ll[4:]]) + xc.append(float(ll[1])) + yc.append(float(ll[2])) + return iverts, np.array(xc), np.array(yc) + + u_data_ws = os.path.join("..", "examples", "data", "unstructured") + outfile = os.path.join("temp", "t050", "vtk_disu_model", "disu_grid.vtu") + + # load vertices + fname = os.path.join(u_data_ws, "ugrid_verts.dat") + verts = load_verts(fname) + + # load the index list into iverts, xc, and yc + fname = os.path.join(u_data_ws, "ugrid_iverts.dat") + iverts, xc, yc = load_iverts(fname) + + # create a 3 layer model grid + ncpl = np.array(3 * [len(iverts)]) + nnodes = np.sum(ncpl) + + top = np.ones( + (nnodes), + ) + botm = np.ones( + (nnodes), + ) + + # set top and botm elevations + i0 = 0 + i1 = ncpl[0] + elevs = [100, 0, -100, -200] + for ix, cpl in enumerate(ncpl): + top[i0:i1] *= elevs[ix] + botm[i0:i1] *= elevs[ix + 1] + i0 += cpl + i1 += cpl + + # create the modelgrid + modelgrid = flopy.discretization.UnstructuredGrid( + vertices=verts, + iverts=iverts, + xcenters=xc, + ycenters=yc, + top=top, + botm=botm, + ncpl=ncpl, + ) + + vtkobj = Vtk( + modelgrid=modelgrid, vertical_exageration=2, binary=True, smooth=False + ) + vtkobj.add_array(modelgrid.top, "top") + vtkobj.add_array(modelgrid.botm, "botm") + vtkobj.write(outfile) + + if not os.path.exists(outfile): + raise FileNotFoundError("VTK DISU test file not written") + + reader = vtk.vtkUnstructuredGridReader() + reader.SetFileName(outfile) + reader.ReadAllFieldsOn() + reader.Update() + + data = reader.GetOutput() + + top2 = numpy_support.vtk_to_numpy(data.GetCellData().GetArray("top")) + + if not np.allclose(np.ravel(top), top2): + raise AssertionError("Field data not properly written") + + +def test_vtk_vertex(): + try: + import vtk + from vtk.util import numpy_support + except ImportError: + return + + # disv test + workspace = os.path.join("..", "examples", "data", "mf6", + "test003_gwfs_disv") + # outfile = os.path.join("vtk_transient_test", "vtk_pacakages") + outfile = os.path.join("temp", "t050", "vtk_disv", "disv.vtk") + sim = flopy.mf6.MFSimulation.load(sim_ws=workspace) + gwf = sim.get_model("gwf_1") + sim.run_simulation() + + vtkobj = Vtk(model=gwf, binary=True, smooth=False) + vtkobj.add_model(gwf) + vtkobj.write(outfile) + + outfile = outfile.split(".")[0] + "_000000.vtk" + if not os.path.exists(outfile): + raise FileNotFoundError("Vertex VTK File was not written") + + reader = vtk.vtkUnstructuredGridReader() + reader.SetFileName(outfile) + reader.ReadAllFieldsOn() + reader.Update() + + data = reader.GetOutput() + + hk2 = numpy_support.vtk_to_numpy(data.GetCellData().GetArray("k")) + hk = gwf.npf.k.array + hk[gwf.modelgrid.idomain == 0] = np.nan + + if not np.allclose(np.ravel(hk), hk2, equal_nan=True): + raise AssertionError("Field data not properly written") + + +def test_vtk_pathline(): + try: + import vtk + from vtk.util import numpy_support + except ImportError: + return + # pathline test for vtk + ws = os.path.join("..", "examples", "data", "freyberg") + modelpth = os.path.join("temp", "t050") + outfile = os.path.join(modelpth, "pathline_test", "pathline.vtk") + ml = flopy.modflow.Modflow.load("freyberg.nam", model_ws=ws, + exe_name="mf2005") + ml.change_model_ws(new_pth=modelpth) + ml.write_input() + ml.run_model() + + mpp = flopy.modpath.Modpath6( + "freybergmpp", modflowmodel=ml, model_ws=modelpth, exe_name="mp6" + ) + mpbas = flopy.modpath.Modpath6Bas( + mpp, + hnoflo=ml.bas6.hnoflo, + hdry=ml.lpf.hdry, + ibound=ml.bas6.ibound.array, + prsity=0.2, + prsityCB=0.2, + ) + sim = mpp.create_mpsim( + trackdir="backward", simtype="pathline", packages="WEL" + ) + mpp.write_input() + mpp.run_model() - # export and check specific discharge given at vertices - cbcfile = os.path.join(output_dir, f"{name}.cbc") - cbc = bf.CellBudgetFile(cbcfile) - keys = ["FLOW RIGHT FACE", "FLOW LOWER FACE"] - vectors = [cbc.get_data(text=t)[0] for t in keys] - vectors.insert(1, None) - q = pp.get_specific_discharge(vectors, m, position="vertices") - vtk.export_vector( - m, q, output_dir, f"{name}_q", point_scalars=True, true2d=True - ) - filetocheck = os.path.join(output_dir, f"{name}_q.vtu") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==7036) - nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 123 - return + pthfile = os.path.join(modelpth, mpp.sim.pathline_file) + pthobj = flopy.utils.PathlineFile(pthfile) + travel_time_max = 200.0 * 365.25 * 24.0 * 60.0 * 60.0 + plines = pthobj.get_alldata(totim=travel_time_max, ge=False) + vtkobj = Vtk(model=ml, binary=True, vertical_exageration=50, smooth=False) + vtkobj.add_model(ml) + vtkobj.add_pathline_points(plines) + vtkobj.write(outfile) + + outfile = outfile.split(".")[0] + "_pathline.vtk" + if not os.path.exists(outfile): + raise FileNotFoundError("Pathline VTK file not properly written") + + reader = vtk.vtkUnstructuredGridReader() + reader.SetFileName(outfile) + reader.ReadAllFieldsOn() + reader.Update() -def test_vtk_export_true2d_nonregyz(): - import flopy.utils.binaryfile as bf - from flopy.utils import postprocessing as pp - - output_dir = os.path.join(cpth, "true2d_nonregyz") - cbc_unit_nb = 53 - - # model with one col, non-regular grid in y and stepwise z - name = "nonregyz" - m = flopy.modflow.Modflow(name, model_ws=output_dir, exe_name="mf2005") - nlay, nrow, ncol = 2, 10, 1 - delr = 1.0 - delc = np.concatenate((2.0 * np.ones((5,)), np.ones((5,)))) - top = np.linspace(3.0, 2.0, nrow).reshape((1, nrow, 1)) - botm1 = np.linspace(2.5, 1.0, nrow).reshape((1, nrow, 1)) - botm2 = np.linspace(0.5, 0.0, nrow).reshape((1, nrow, 1)) - botm = np.concatenate((botm1, botm2)) - dis = flopy.modflow.ModflowDis( - m, nlay, nrow, ncol, delr=delr, delc=delc, top=top, botm=botm - ) - ibound = np.ones((nlay, nrow, ncol), dtype=np.int32) - ibound[:, 0, :] = -1 - ibound[:, -1, :] = -1 - strt = np.linspace(1.0, 0.0, nrow).reshape(1, nrow, 1) - strt = strt * np.ones((nlay, nrow, ncol)) - bas = flopy.modflow.ModflowBas(m, ibound=ibound, strt=strt) - lpf = flopy.modflow.ModflowLpf(m, hk=1.0, vka=1.0, ipakcb=cbc_unit_nb) - spd = {(0, 0): ["print head", "print budget", "save head", "save budget"]} - oc = flopy.modflow.ModflowOc(m, stress_period_data=spd, compact=True) - pcg = flopy.modflow.ModflowPcg(m) - m.write_input() - m.run_model(silent=True) - - # export and check head - hdsfile = os.path.join(output_dir, f"{name}.hds") - hds = bf.HeadFile(hdsfile) - head = hds.get_data() - vtk.export_array(m, head, output_dir, f"{name}_head", true2d=True) - filetocheck = os.path.join(output_dir, f"{name}_head.vtu") - # totalbytes = os.path.getsize(filetocheck) - # assert(totalbytes==4217) - nlines = count_lines_in_file(filetocheck) - assert nlines == 105 - - # export and check head with point scalar - hdsfile = os.path.join(output_dir, f"{name}.hds") - hds = bf.HeadFile(hdsfile) - head = hds.get_data() - vtk.export_array( - m, - head, - output_dir, - f"{name}_head_points", - point_scalars=True, - true2d=True, - ) - filetocheck = os.path.join(output_dir, f"{name}_head_points.vtu") - # totalbytes1 = os.path.getsize(filetocheck) - # assert(totalbytes1==6155) - nlines1 = count_lines_in_file(filetocheck) - assert nlines1 == 129 - - # export and check specific discharge given at vertices - cbcfile = os.path.join(output_dir, f"{name}.cbc") - cbc = bf.CellBudgetFile(cbcfile) - keys = ["FLOW FRONT FACE", "FLOW LOWER FACE"] - vectors = [cbc.get_data(text=t)[0] for t in keys] - vectors.insert(0, None) - q = pp.get_specific_discharge(vectors, m, position="vertices") - vtk.export_vector( - m, q, output_dir, f"{name}_q", point_scalars=True, true2d=True - ) - filetocheck = os.path.join(output_dir, f"{name}_q.vtu") - # totalbytes2 = os.path.getsize(filetocheck) - # assert(totalbytes2==7032) - nlines2 = count_lines_in_file(filetocheck) - assert nlines2 == 123 - return + data = reader.GetOutput() + + totim = numpy_support.vtk_to_numpy(data.GetCellData().GetArray("time")) + pid = numpy_support.vtk_to_numpy(data.GetCellData().GetArray("particleid")) + + maxtime = 0 + for p in plines: + if np.max(p['time']) > maxtime: + maxtime = np.max(p['time']) + + if not len(totim) == 12054: + raise AssertionError("Array size is incorrect for modpath VTK") + + if not np.abs(np.max(totim) - maxtime) < 100: + raise AssertionError("time values are incorrect for modpath VTK") + + if not np.max(pid) == len(plines): + raise AssertionError( + "number of particles are incorrect for modpath VTK" + ) if __name__ == "__main__": - # test_vtk_export_array2d() - # test_vtk_export_array3d() - # test_vtk_transient_array_2d() + test_vtk_export_array2d() + test_vtk_export_array3d() + test_vtk_transient_array_2d() test_vtk_export_packages() - # test_vtk_mf6() - # test_vtk_binary_head_export() - # test_vtk_cbc() - # test_vtk_vector() - test_vtk_vti() - # test_vtk_vtr() - # test_vtk_export_true2d_regular() - # test_vtk_export_true2d_nonregxy() - # test_vtk_export_true2d_nonregxz() - # test_vtk_export_true2d_nonregyz() + test_vtk_mf6() + test_vtk_binary_head_export() + test_vtk_cbc() + test_vtk_vector() + test_vtk_unstructured() + test_vtk_vertex() + test_vtk_pathline() \ No newline at end of file diff --git a/etc/environment.yml b/etc/environment.yml index 9e7b494c6..e78f3169b 100644 --- a/etc/environment.yml +++ b/etc/environment.yml @@ -26,3 +26,4 @@ dependencies: - shapely - geos - geojson + - vtk \ No newline at end of file diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 0ded8e005..150a2a918 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -622,25 +622,30 @@ def model_export(f, ml, fmt=None, **kwargs): elif fmt == "vtk": # call vtk model export - nanval = kwargs.get("nanval", -1e20) + name = kwargs.get("name", ml.name) + xml = kwargs.get("xml", False) + masked_values = kwargs.get("masked_values", None) + pvd = kwargs.get("pvd", False) smooth = kwargs.get("smooth", False) point_scalars = kwargs.get("point_scalars", False) - vtk_grid_type = kwargs.get("vtk_grid_type", "auto") - true2d = kwargs.get("true2d", False) - binary = kwargs.get("binary", False) + binary = kwargs.get("binary", True) + vertical_exageration = kwargs.get("vertical_exageration", 1) kpers = kwargs.get("kpers", None) - vtk.export_model( + package_names = kwargs.get("package_names", None) + + vtkobj = vtk.Vtk( ml, - f, - package_names=package_names, - nanval=nanval, + vertical_exageration=vertical_exageration, + binary=binary, + xml=xml, + pvd=pvd, smooth=smooth, point_scalars=point_scalars, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - kpers=kpers, ) + vtkobj.add_model( + ml, masked_values=masked_values, selpaklist=package_names + ) + vtkobj.write(os.path.join(f, name), kpers) else: raise NotImplementedError(f"unrecognized export argument:{f}") @@ -714,13 +719,30 @@ def package_export(f, pak, fmt=None, **kwargs): elif fmt == "vtk": # call vtk array export to folder - nanval = kwargs.get("nanval", -1e20) + name = kwargs.get("name", pak.name[0]) + xml = kwargs.get("xml", False) + masked_values = kwargs.get("masked_values", None) + pvd = kwargs.get("pvd", False) smooth = kwargs.get("smooth", False) point_scalars = kwargs.get("point_scalars", False) - vtk_grid_type = kwargs.get("vtk_grid_type", "auto") - true2d = kwargs.get("true2d", False) - binary = kwargs.get("binary", False) + binary = kwargs.get("binary", True) + vertical_exageration = kwargs.get("vertical_exageration", 1) kpers = kwargs.get("kpers", None) + + vtkobj = vtk.Vtk( + pak.parent, + vertical_exageration=vertical_exageration, + binary=binary, + xml=xml, + pvd=pvd, + smooth=smooth, + point_scalars=point_scalars, + ) + + vtkobj.add_package(pak, masked_values=masked_values) + vtkobj.write(os.path.join(f, name), kper=kpers) + + """ vtk.export_package( pak.parent, pak.name, @@ -733,6 +755,7 @@ def package_export(f, pak, fmt=None, **kwargs): binary=binary, kpers=kpers, ) + """ else: raise NotImplementedError(f"unrecognized export argument:{f}") @@ -1083,27 +1106,35 @@ def transient2d_export(f, t2d, fmt=None, **kwargs): elif fmt == "vtk": name = kwargs.get("name", t2d.name) - nanval = kwargs.get("nanval", -1e20) + xml = kwargs.get("xml", False) + masked_values = kwargs.get("masked_values", None) + pvd = kwargs.get("pvd", False) smooth = kwargs.get("smooth", False) point_scalars = kwargs.get("point_scalars", False) - vtk_grid_type = kwargs.get("vtk_grid_type", "auto") - true2d = kwargs.get("true2d", False) - binary = kwargs.get("binary", False) + binary = kwargs.get("binary", True) + vertical_exageration = kwargs.get("vertical_exageration", 1) kpers = kwargs.get("kpers", None) - vtk.export_transient( + + if t2d.array is not None: + if hasattr(t2d, "transient_2ds"): + d = t2d.transient_2ds + else: + d = {ix: i for ix, i in enumerate(t2d.array)} + else: + raise AssertionError("No data available to export") + + vtkobj = vtk.Vtk( t2d.model, - t2d.array, - f, - name, - nanval=nanval, + vertical_exageration=vertical_exageration, + binary=binary, + xml=xml, + pvd=pvd, smooth=smooth, point_scalars=point_scalars, - array2d=True, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - kpers=kpers, ) + vtkobj.add_transient_array(d, name=name, masked_values=masked_values) + vtkobj.write(os.path.join(f, name), kper=kpers) + else: raise NotImplementedError(f"unrecognized export argument:{f}") @@ -1253,27 +1284,26 @@ def array3d_export(f, u3d, fmt=None, **kwargs): elif fmt == "vtk": # call vtk array export to folder name = kwargs.get("name", u3d.name) - nanval = kwargs.get("nanval", -1e20) + xml = kwargs.get("xml", False) + masked_values = kwargs.get("masked_values", None) smooth = kwargs.get("smooth", False) point_scalars = kwargs.get("point_scalars", False) - vtk_grid_type = kwargs.get("vtk_grid_type", "auto") - true2d = kwargs.get("true2d", False) - binary = kwargs.get("binary", False) + binary = kwargs.get("binary", True) + vertical_exageration = kwargs.get("vertical_exageration", 1) + if isinstance(name, list) or isinstance(name, tuple): name = name[0] - vtk.export_array( + vtkobj = vtk.Vtk( u3d.model, - u3d.array, - f, - name, - nanval=nanval, smooth=smooth, point_scalars=point_scalars, - vtk_grid_type=vtk_grid_type, - true2d=true2d, binary=binary, + xml=xml, + vertical_exageration=vertical_exageration, ) + vtkobj.add_array(u3d.array, name, masked_values=masked_values) + vtkobj.write(os.path.join(f, name)) else: raise NotImplementedError(f"unrecognized export argument:{f}") @@ -1395,28 +1425,29 @@ def array2d_export(f, u2d, fmt=None, **kwargs): elif fmt == "vtk": - # call vtk array export to folder name = kwargs.get("name", u2d.name) - nanval = kwargs.get("nanval", -1e20) + xml = kwargs.get("xml", False) + masked_values = kwargs.get("masked_values", None) smooth = kwargs.get("smooth", False) point_scalars = kwargs.get("point_scalars", False) - vtk_grid_type = kwargs.get("vtk_grid_type", "auto") - true2d = kwargs.get("true2d", False) - binary = kwargs.get("binary", False) - vtk.export_array( + binary = kwargs.get("binary", True) + vertical_exageration = kwargs.get("vertical_exageration", 1) + + if isinstance(name, list) or isinstance(name, tuple): + name = name[0] + + vtkobj = vtk.Vtk( u2d.model, - u2d.array, - f, - name, - nanval=nanval, smooth=smooth, point_scalars=point_scalars, - array2d=True, - vtk_grid_type=vtk_grid_type, - true2d=true2d, binary=binary, + xml=xml, + vertical_exageration=vertical_exageration, ) - + array = np.ones((modelgrid.nnodes,)) * np.nan + array[0 : u2d.array.size] = np.ravel(u2d.array) + vtkobj.add_array(array, name, masked_values=masked_values) + vtkobj.write(os.path.join(f, name)) else: raise NotImplementedError(f"unrecognized export argument:{f}") diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 898fcd572..1b7392c5e 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -1,2142 +1,1937 @@ -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 - -np_to_vtk_type = { - "int8": "Int8", - "uint8": "UInt8", - "int16": "Int16", - "uint16": "UInt16", - "int32": "Int32", - "uint32": "UInt32", - "int64": "Int64", - "uint64": "UInt64", - "float32": "Float32", - "float64": "Float64", -} - -np_to_struct = { - "int8": "b", - "uint8": "B", - "int16": "h", - "uint16": "H", - "int32": "i", - "uint32": "I", - "int64": "q", - "uint64": "Q", - "float32": "f", - "float64": "d", -} - - -class XmlWriterInterface: - """ - Helps writing vtk files. - - Parameters - ---------- - - file_path : str - output file path - """ - - def __init__(self, file_path): - # class attributes - self.open_tag = False - self.current = [] - self.indent_level = 0 - self.indent_char = " " - - # open file and initialize - self.f = self._open_file(file_path) - self.write_string('') - - # open VTKFile element - self.open_element("VTKFile").add_attributes(version="0.1") - - def _open_file(self, file_path): - """ - Open the file for writing. - - Return - ------ - File object. - """ - raise NotImplementedError("must define _open_file in child class") - - def write_string(self, string): - """ - Write a string to the file. - """ - raise NotImplementedError("must define write_string in child class") - - def open_element(self, tag): - if self.open_tag: - self.write_string(">") - indent = self.indent_level * self.indent_char - self.indent_level += 1 - tag_string = f"\n{indent}<{tag}" - self.write_string(tag_string) - self.open_tag = True - self.current.append(tag) - return self - - def close_element(self, tag=None): - self.indent_level -= 1 - if tag: - assert self.current.pop() == tag - if self.open_tag: - self.write_string(">") - self.open_tag = False - indent = self.indent_level * self.indent_char - tag_string = f"\n{indent}" - self.write_string(tag_string) - else: - self.write_string("/>") - self.open_tag = False - self.current.pop() - return self - - def add_attributes(self, **kwargs): - assert self.open_tag - for key in kwargs: - st = f' {key}="{kwargs[key]}"' - self.write_string(st) - return self - - def write_line(self, text): - if self.open_tag: - self.write_string(">") - self.open_tag = False - self.write_string("\n") - indent = self.indent_level * self.indent_char - self.write_string(indent) - self.write_string(text) - return self - - def write_array(self, array, actwcells=None, **kwargs): - """ - Write an array to the file. - - Parameters - ---------- - array : ndarray - the data array being output - actwcells : array - array of the active cells - kwargs : dictionary - Attributes to be added to the DataArray element - """ - raise NotImplementedError("must define write_array in child class") - - def final(self): - """ - Finalize the file. Must be called. - """ - self.close_element("VTKFile") - assert not self.open_tag - self.f.close() - - -class XmlWriterAscii(XmlWriterInterface): - """ - Helps writing ascii vtk files. - - Parameters - ---------- - - file_path : str - output file path - """ - - def __init__(self, file_path): - super().__init__(file_path) - - def _open_file(self, file_path): - """ - Open the file for writing. - - Return - ------ - File object. - """ - return open(file_path, "w") - - def write_string(self, string): - """ - Write a string to the file. - """ - self.f.write(string) - - def write_array(self, array, actwcells=None, **kwargs): - """ - Write an array to the file. - - Parameters - ---------- - array : ndarray - the data array being output - actwcells : array - array of the active cells - kwargs : dictionary - Attributes to be added to the DataArray element - """ - # open DataArray element with relevant attributes - self.open_element("DataArray") - vtk_type = np_to_vtk_type[array.dtype.name] - self.add_attributes(type=vtk_type) - self.add_attributes(**kwargs) - self.add_attributes(format="ascii") - - # write the data - nlay = array.shape[0] - for lay in range(nlay): - if actwcells is not None: - idx = actwcells[lay] != 0 - array_lay_flat = array[lay][idx].flatten() - else: - array_lay_flat = array[lay].flatten() - # replace NaN values by -1e9 as there is a bug is Paraview when - # reading NaN in ASCII mode - # https://gitlab.kitware.com/paraview/paraview/issues/19042 - # this may be removed in the future if they fix the bug - array_lay_flat[np.isnan(array_lay_flat)] = -1e9 - self.write_line(" ".join(str(val) for val in array_lay_flat)) - - # close DataArray element - self.close_element("DataArray") - return - - -class XmlWriterBinary(XmlWriterInterface): - """ - Helps writing binary vtk files. - - Parameters - ---------- - - file_path : str - output file path - - """ - - def __init__(self, file_path): - super().__init__(file_path) - - if sys.byteorder == "little": - self.byte_order = "<" - self.add_attributes(byte_order="LittleEndian") - else: - self.byte_order = ">" - self.add_attributes(byte_order="BigEndian") - self.add_attributes(header_type="UInt64") - - # class attributes - self.offset = 0 - self.byte_count_size = 8 - self.processed_arrays = [] - - def _open_file(self, file_path): - """ - Open the file for writing. - - Return - ------ - File object. - """ - return open(file_path, "wb") - - def write_string(self, string): - """ - Write a string to the file. - """ - self.f.write(str.encode(string)) - - def write_array(self, array, actwcells=None, **kwargs): - """ - Write an array to file. - - Parameters - ---------- - array : ndarray - the data array being output - actwcells : array - array of the active cells - kwargs : dictionary - Attributes to be added to the DataArray element - """ - # open DataArray element with relevant attributes - self.open_element("DataArray") - vtk_type = np_to_vtk_type[array.dtype.name] - self.add_attributes(type=vtk_type) - self.add_attributes(**kwargs) - self.add_attributes(format="appended", offset=self.offset) - - # store array for later writing (appended data section) - if actwcells is not None: - array = array[actwcells != 0] - a = np.ascontiguousarray(array.ravel()) - array_size = array.size * array[0].dtype.itemsize - self.processed_arrays.append([a, array_size]) - - # calculate the offset of the start of the next piece of data - # offset is calculated from beginning of data section - self.offset += array_size + self.byte_count_size - - # close DataArray element - self.close_element("DataArray") - return - - def _write_size(self, block_size): - # size is a 64 bit unsigned integer - byte_order = f"{self.byte_order}Q" - block_size = struct.pack(byte_order, block_size) - self.f.write(block_size) - - def _append_array_binary(self, data): - # see vtk documentation and more details here: - # https://vtk.org/Wiki/VTK_XML_Formats#Appended_Data_Section - assert data.flags["C_CONTIGUOUS"] or data.flags["F_CONTIGUOUS"] - assert data.ndim == 1 - data_format = ( - self.byte_order + str(data.size) + np_to_struct[data.dtype.name] - ) - binary_data = struct.pack(data_format, *data) - self.f.write(binary_data) - - def final(self): - """ - Finalize the file. Must be called. - """ - # build data section - self.open_element("AppendedData") - self.add_attributes(encoding="raw") - self.write_line("_") - for a, block_size in self.processed_arrays: - self._write_size(block_size) - self._append_array_binary(a) - self.close_element("AppendedData") - - # call super final - super().final() - - -class _Array: - # 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: - """ - 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 layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output file will be binary, default is False - - Attributes - ---------- - - arrays : dict - Stores data arrays added to VTK object - """ - - def __init__( - self, - model, - verbose=None, - nanval=-1e20, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=False, - binary=False, - ): - - if point_scalars: - smooth = True - - if verbose is None: - verbose = model.verbose - self.verbose = verbose - - # set up variables - self.model = model - self.modelgrid = model.modelgrid - self.nlay = self.modelgrid.nlay - if hasattr(self.model, "dis") and hasattr(self.model.dis, "laycbd"): - self.nlay = self.nlay + np.sum(self.model.dis.laycbd.array > 0) - self.nrow = self.modelgrid.nrow - self.ncol = self.modelgrid.ncol - self.shape = (self.nlay, self.nrow, self.ncol) - self.shape2d = (self.shape[1], self.shape[2]) - self.shape_verts = ( - self.shape[0] + 1, - self.shape[1] + 1, - self.shape[2] + 1, - ) - self.shape_verts2d = (self.shape_verts[1], self.shape_verts[2]) - self.nanval = nanval - - self.arrays = {} - self.vectors = {} - - self.smooth = smooth - self.point_scalars = point_scalars - self.has_cell_data = False - self.has_point_data = False - - # 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 - - self.true2d = true2d - self.nx = self.modelgrid.ncol - self.ny = self.modelgrid.nrow - self.nz = self.modelgrid.nlay - if self.true2d: - if self.nz == 1: - self.nz = 0 - elif self.ny == 1: - self.ny = 0 - elif self.nx == 1: - self.nx = 0 - else: - raise ValueError( - "The option true2d was used but the model is not 2d." - ) - self.cell_type = 8 - else: - self.cell_type = 11 - - self.vtk_grid_type, self.file_extension = self._vtk_grid_type( - vtk_grid_type - ) - - self.binary = binary - - return - - def _vtk_grid_type(self, vtk_grid_type="auto"): - """ - Determines the vtk grid type and corresponding file extension. - - Parameters - ---------- - vtk_grid_type : str - Specific vtk_grid_type or 'auto'. Possible specific values are - 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - - Returns - ---------- - (vtk_grid_type, file_extension) : tuple of two strings - """ - # if 'auto', determine the vtk grid type automatically - if vtk_grid_type == "auto": - if self.modelgrid.grid_type == "structured": - if ( - self.modelgrid.is_regular - or (self.modelgrid.is_regular_xy and self.nz == 0) - or (self.modelgrid.is_regular_xz and self.ny == 0) - or (self.modelgrid.is_regular_yz and self.nx == 0) - ): - vtk_grid_type = "ImageData" - elif self.modelgrid.is_rectilinear or self.nz == 0: - vtk_grid_type = "RectilinearGrid" - else: - vtk_grid_type = "UnstructuredGrid" - else: - vtk_grid_type = "UnstructuredGrid" - # otherwise, check the validity of the passed vtk_grid_type - else: - allowable_types = [ - "ImageData", - "RectilinearGrid", - "UnstructuredGrid", - ] - if not any(vtk_grid_type in s for s in allowable_types): - raise ValueError( - f'"{vtk_grid_type}" is not a correct vtk_grid_type.' - ) - if ( - vtk_grid_type == "ImageData" - or vtk_grid_type == "RectilinearGrid" - ) and not self.modelgrid.grid_type == "structured": - raise NotImplementedError( - 'vtk_grid_type cannot be "' - + vtk_grid_type - + '" for a grid ' - "that is not structured" - ) - if ( - vtk_grid_type == "ImageData" - and not self.modelgrid.is_regular - and not (self.modelgrid.is_regular_xy and self.nz == 0) - and not (self.modelgrid.is_regular_xz and self.ny == 0) - and not (self.modelgrid.is_regular_yz and self.nx == 0) - ): - raise ValueError( - 'vtk_grid_type cannot be "ImageData" for a ' - "non-regular grid spacing" - ) - if ( - vtk_grid_type == "RectilinearGrid" - and not self.modelgrid.is_rectilinear - and not self.nz == 0 - ): - raise ValueError( - 'vtk_grid_type cannot be "RectilinearGrid" ' - "for a non-rectilinear grid spacing" - ) - - # determine the file extension - if vtk_grid_type == "ImageData": - file_extension = ".vti" - elif vtk_grid_type == "RectilinearGrid": - file_extension = ".vtr" - # else vtk_grid_type == 'UnstructuredGrid' - else: - file_extension = ".vtu" - - # return vtk grid type and file extension - return (vtk_grid_type, file_extension) - - def _format_array(self, a, array2d=False): - """ - Formats array for vtk output. - - 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 - - Return - ------ - Formatted array (note a copy is made) - """ - # if array is 2d reformat to 3d array - if array2d: - if a.shape == self.shape2d: - array = np.full(self.shape, self.nanval) - elif a.shape == self.shape_verts2d: - array = np.full(self.shape_verts, self.nanval) - else: - raise ValueError("Incompatible array size") - array[0, :, :] = a - a = array - - # deal with inactive cells - inactive3d = self.ibound == 0 - if a.shape == self.shape: - # set to nan where nanval or where ibound==0 - where_to_nan = np.logical_or(a == self.nanval, inactive3d) - self.has_cell_data = True - elif a.shape == self.shape_verts: - # set to nan where ibound==0 at all 8 neighbors - where_to_nan = np.full(self.shape_verts, True) - where_to_nan[:-1, :-1, :-1] = inactive3d - where_to_nan[:-1, :-1, 1:] = np.logical_and( - where_to_nan[:-1, :-1, 1:], inactive3d - ) - where_to_nan[:-1, 1:, :-1] = np.logical_and( - where_to_nan[:-1, 1:, :-1], inactive3d - ) - where_to_nan[:-1, 1:, 1:] = np.logical_and( - where_to_nan[:-1, 1:, 1:], inactive3d - ) - where_to_nan[1:, :-1, :-1] = np.logical_and( - where_to_nan[1:, :-1, :-1], inactive3d - ) - where_to_nan[1:, :-1, 1:] = np.logical_and( - where_to_nan[1:, :-1, 1:], inactive3d - ) - where_to_nan[1:, 1:, :-1] = np.logical_and( - where_to_nan[1:, 1:, :-1], inactive3d - ) - where_to_nan[1:, 1:, 1:] = np.logical_and( - where_to_nan[1:, 1:, 1:], inactive3d - ) - self.has_point_data = True - self.smooth = True - else: - # incompatible size, skip this array - return None - a = np.where(where_to_nan, np.nan, a) - - return a - - 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. - The shape should match either grid cells or grid vertices. - array2d : bool - true if the array is 2d and represents the first layer, - default is False - """ - # format array - a = self._format_array(a, array2d) - - # add to self.arrays - if a is not None: - self.arrays[name] = a - - return - - def add_vector(self, name, v, array2d=False): - """ - Adds a vector (i.e., a tuple of arrays) to the vtk object. - - Parameters - ---------- - - name : str - Name of the vector. - v : tuple of arrays - The vector to be added to the vtk object. The shape of each - component should match either grid cells or grid vertices. - array2d : bool - true if the vector components are 2d arrays and represent the first - layer, default is False - - Notes - ----- - If the grid is rotated, the vector will be rotated too, assuming that - the first and second components are along x and y directions, - respectively. - """ - # format each component of the vector - vf = () - for vcomp in v: - vcomp = self._format_array(vcomp, array2d=array2d) - if vcomp is None: - return - vf = vf + (vcomp,) - - # rotate the vector according to grid - if self.modelgrid.angrot_radians != 0.0: - from ..utils import geometry - - vf = list(vf) - vf[0], vf[1] = geometry.rotate( - vf[0], vf[1], 0.0, 0.0, self.modelgrid.angrot_radians - ) - vf = tuple(vf) - - # add to self.vectors - self.vectors[name] = vf - - return - - def write(self, output_file, timeval=None): - """ - Writes the stored arrays to vtk file in XML format. - - Parameters - ---------- - - output_file : str - output file name without extension (extension is determined - automatically) - timeval : scalar - model time value to be stored in the time section of the vtk - file, default is None - """ - # output file - output_file = output_file + self.file_extension - if self.verbose: - print(f"Writing vtk file: {output_file}") - - # initialize xml file - if self.binary: - xml = XmlWriterBinary(output_file) - else: - xml = XmlWriterAscii(output_file) - xml.add_attributes(type=self.vtk_grid_type) - - # grid type - xml.open_element(self.vtk_grid_type) - - # if time value write time section - if timeval: - xml.open_element("FieldData") - xml.write_array( - np.array([timeval]), - Name="TimeValue", - NumberOfTuples="1", - RangeMin="{0}", - RangeMax="{0}", - ) - xml.close_element("FieldData") - - if self.vtk_grid_type == "UnstructuredGrid": - # get the active data cells based on the data arrays and ibound - actwcells3d = self._configure_data_arrays() - - # get the verts and iverts to be output - verts, iverts, _ = self._get_3d_vertex_connectivity( - actwcells=actwcells3d - ) - - # check if there is data to be written out - if len(verts) == 0: - # if nothing, cannot write file - return - - # get the total number of cells and vertices - ncells = len(iverts) - if self.true2d: - npoints = ncells * 4 - else: - npoints = ncells * 8 - if self.verbose: - print( - f"Number of point is {npoints}, Number of cells is {ncells}\n" - ) - - # piece - xml.open_element("Piece") - xml.add_attributes(NumberOfPoints=npoints, NumberOfCells=ncells) - - # points - xml.open_element("Points") - verts = np.array(list(verts.values())) - verts.reshape(npoints, 3) - xml.write_array(verts, Name="points", NumberOfComponents="3") - xml.close_element("Points") - - # cells - xml.open_element("Cells") - - # connectivity - iverts = np.array(list(iverts.values())) - xml.write_array( - iverts, Name="connectivity", NumberOfComponents="1" - ) - - # offsets - offsets = np.empty((iverts.shape[0]), np.int32) - icount = 0 - for index, row in enumerate(iverts): - icount += len(row) - offsets[index] = icount - xml.write_array(offsets, Name="offsets", NumberOfComponents="1") - - # types - types = np.full((iverts.shape[0]), self.cell_type, dtype=np.uint8) - xml.write_array(types, Name="types", NumberOfComponents="1") - - # end cells - xml.close_element("Cells") - - elif self.vtk_grid_type == "ImageData": - # note: in vtk, "extent" actually means indices of grid lines - vtk_extent_str = f"0 {self.nx} 0 {self.ny} 0 {self.nz}" - xml.add_attributes(WholeExtent=vtk_extent_str) - grid_extent = self.modelgrid.xyzextent - vtk_origin_str = ( - f"{grid_extent[0]} {grid_extent[2]} {grid_extent[4]}" - ) - xml.add_attributes(Origin=vtk_origin_str) - vtk_spacing_str = "{} {} {}".format( - self.modelgrid.delr[0], - self.modelgrid.delc[0], - self.modelgrid.top[0, 0] - self.modelgrid.botm[0, 0, 0], - ) - xml.add_attributes(Spacing=vtk_spacing_str) - - # piece - xml.open_element("Piece").add_attributes(Extent=vtk_extent_str) - - elif self.vtk_grid_type == "RectilinearGrid": - # note: in vtk, "extent" actually means indices of grid lines - vtk_extent_str = f"0 {self.nx} 0 {self.ny} 0 {self.nz}" - xml.add_attributes(WholeExtent=vtk_extent_str) - - # piece - xml.open_element("Piece").add_attributes(Extent=vtk_extent_str) - - # grid coordinates - xml.open_element("Coordinates") - - # along x - xedges = self.modelgrid.xyedges[0] - xml.write_array(xedges, Name="coord_x", NumberOfComponents="1") - - # along y - yedges = np.flip(self.modelgrid.xyedges[1]) - xml.write_array(yedges, Name="coord_y", NumberOfComponents="1") - - # along z - zedges = np.flip(self.modelgrid.zedges) - xml.write_array(zedges, Name="coord_z", NumberOfComponents="1") - - # end coordinates - xml.close_element("Coordinates") - - if self.has_cell_data: - # cell data - xml.open_element("CellData") - - # loop through stored arrays - for name, a in self.arrays.items(): - if a.shape == self.shape_verts: - # these are dealt with later - continue - if self.vtk_grid_type == "UnstructuredGrid": - xml.write_array( - a, - actwcells=actwcells3d, - Name=name, - NumberOfComponents="1", - ) - else: - # flip "a" so coordinates increase along with indices as in - # vtk - a = np.flip(a, axis=[0, 1]) - xml.write_array(a, Name=name, NumberOfComponents="1") - - # loop through stored vectors - for name, v in self.vectors.items(): - if v[0].shape == self.shape_verts: - # these are dealt with later - continue - ncomp = len(v) - v_as_array = np.moveaxis(np.array(v), 0, -1) - if self.vtk_grid_type == "UnstructuredGrid": - shape4d = actwcells3d.shape + (ncomp,) - actwcells4d = actwcells3d.reshape(actwcells3d.shape + (1,)) - actwcells4d = np.broadcast_to(actwcells4d, shape4d) - xml.write_array( - v_as_array, - actwcells=actwcells4d, - Name=name, - NumberOfComponents=ncomp, - ) - else: - # flip "v" so coordinates increase along with indices as in - # vtk - v_as_array = np.flip(v_as_array, axis=[0, 1]) - xml.write_array( - v_as_array, Name=name, NumberOfComponents=ncomp - ) - - # end cell data - xml.close_element("CellData") - - if self.point_scalars or self.has_point_data: - # point data (i.e., values at vertices) - xml.open_element("PointData") - - # loop through stored arrays - for name, a in self.arrays.items(): - if a.shape == self.shape: - if not self.point_scalars: - continue - # get the array values onto vertices - if self.vtk_grid_type == "UnstructuredGrid": - _, _, averts = self._get_3d_vertex_connectivity( - actwcells=actwcells3d, zvalues=a - ) - a = np.array(list(averts.values())) - else: - a = self.modelgrid.array_at_verts(a) - a = np.flip(a, axis=[0, 1]) - # deal with true2d - if self.true2d: - if self.nz == 0: - a = a[0, :, :] - elif self.ny == 0: - a = a[:, 0, :] - elif self.nz == 0: - a = a[:, :, 0] - else: - if self.vtk_grid_type == "UnstructuredGrid": - # still need to do this to be consistent with - # connectivity (i.e. 8 points for every cell) - _, _, averts = self._get_3d_vertex_connectivity( - actwcells=actwcells3d, zvalues=a - ) - a = np.array(list(averts.values())) - else: - # flip "a" so coordinates increase along with indices - # as in vtk - a = np.flip(a, axis=[0, 1]) - # deal with true2d - if self.true2d: - if self.nz == 0: - a = a[0, :, :] - elif self.ny == 0: - a = a[:, 0, :] - elif self.nz == 0: - a = a[:, :, 0] - xml.write_array(a, Name=name, NumberOfComponents="1") - - # loop through stored vectors - for name, v in self.vectors.items(): - if v[0].shape == self.shape: - if not self.point_scalars: - continue - # get the vector values onto vertices - v_verts = () - for vcomp in v: - if self.vtk_grid_type == "UnstructuredGrid": - _, _, averts = self._get_3d_vertex_connectivity( - actwcells=actwcells3d, zvalues=vcomp - ) - vcomp = np.array(list(averts.values())) - else: - vcomp = self.modelgrid.array_at_verts(vcomp) - vcomp = np.flip(vcomp, axis=[0, 1]) - # deal with true2d - if self.true2d: - if self.nz == 0: - vcomp = vcomp[0, :, :] - elif self.ny == 0: - vcomp = vcomp[:, 0, :] - elif self.nz == 0: - vcomp = vcomp[:, :, 0] - v_verts = v_verts + (vcomp,) - v = v_verts - else: - v_verts = () - for vcomp in v: - if self.vtk_grid_type == "UnstructuredGrid": - # still need to do this to be consistent with - # connectivity (i.e. 8 points for every cell) - _, _, averts = self._get_3d_vertex_connectivity( - actwcells=actwcells3d, zvalues=vcomp - ) - vcomp = np.array(list(averts.values())) - else: - vcomp = np.flip(vcomp, axis=[0, 1]) - # deal with true2d - if self.true2d: - if self.nz == 0: - vcomp = vcomp[0, :, :] - elif self.ny == 0: - vcomp = vcomp[:, 0, :] - elif self.nz == 0: - vcomp = vcomp[:, :, 0] - v_verts = v_verts + (vcomp,) - v = v_verts - # write to file - ncomp = len(v) - v_as_array = np.moveaxis(np.array(v), 0, -1) - xml.write_array( - v_as_array, Name=name, NumberOfComponents=ncomp - ) - - # end point data - xml.close_element("PointData") - - # end piece - xml.close_element("Piece") - - # end vtk_grid_type - xml.close_element(self.vtk_grid_type) - - # finalize and close xml file - xml.final() - - # clear arrays - self.arrays.clear() - self.vectors.clear() - - def _configure_data_arrays(self): - """ - Compares arrays and active cells to find where active data - exists, and what cells to output. - """ - # build 1d index array - shape1d = self.shape[0] * self.shape[1] * self.shape[2] - actwcells1d = np.zeros(shape1d, dtype=int) - if self.has_point_data: - shape1d_verts = ( - self.shape_verts[0] * self.shape_verts[1] * self.shape_verts[2] - ) - actwcells1d_verts = np.zeros(shape1d_verts, dtype=int) - - # loop through arrays - for a in self.arrays.values(): - # make array 1d - a1d = a.ravel() - # get the indexes where there is data - idxs = np.argwhere(np.logical_not(np.isnan(a1d))) - # set the active array to 1 - if a.shape == self.shape: - actwcells1d[idxs] = 1 - elif self.has_point_data: - actwcells1d_verts[idxs] = 1 - - # loop through vectors - for v in self.vectors.values(): - for vcomp in v: - # make array 1d - vcomp1d = vcomp.ravel() - # get the indexes where there is data - idxs = np.argwhere(np.logical_not(np.isnan(vcomp1d))) - # set the active array to 1 - if vcomp.shape == self.shape: - actwcells1d[idxs] = 1 - elif self.has_point_data: - actwcells1d_verts[idxs] = 1 - - # reshape to 3D array - actwcells3d = actwcells1d.reshape(self.shape) - if self.has_point_data: - actwcells3d_verts = actwcells1d_verts.reshape(self.shape_verts) - # activate cells that are neighbor of 8 active vertices - activate = np.full(self.shape, True) - activate[actwcells3d_verts[:-1, :-1, :-1] == 0] = False - activate[actwcells3d_verts[:-1, :-1, 1:] == 0] = False - activate[actwcells3d_verts[:-1, 1:, :-1] == 0] = False - activate[actwcells3d_verts[:-1, 1:, 1:] == 0] = False - activate[actwcells3d_verts[1:, :-1, :-1] == 0] = False - activate[actwcells3d_verts[1:, :-1, 1:] == 0] = False - activate[actwcells3d_verts[1:, 1:, :-1] == 0] = False - activate[actwcells3d_verts[1:, 1:, 1:] == 0] = False - activate[self.ibound == 0] = False - actwcells3d[activate] = 1 - - return actwcells3d - - 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: - if zvalues.shape == self.shape: - # interpolate using the given values - zVertices = self.modelgrid.array_at_verts(zvalues) - else: - # in this case the given values are already at vertices - zVertices = zvalues - else: - zVertices = self.modelgrid.zverts_smooth - 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 - - # determine z values - if self.nz == 0 and zvalues is None: - elev = np.nanmin( - self.modelgrid.top_botm_withnan[k + 1, :, :] - ) - zvals = [ - [elev, elev, elev, elev], - [elev, elev, elev, elev], - ] - elif not self.smooth: - zbot = self.modelgrid.top_botm[k + 1, i, j] - ztop = self.modelgrid.top_botm[k, i, j] - zvals = [ - [zbot, zbot, zbot, zbot], - [ztop, ztop, ztop, ztop], - ] - else: - zvals = [ - [ - zVertices[k + 1, i + 1, j], - zVertices[k + 1, i + 1, j + 1], - zVertices[k + 1, i, j], - zVertices[k + 1, i, j + 1], - ], - [ - zVertices[k, i + 1, j], - zVertices[k, i + 1, j + 1], - zVertices[k, i, j], - zVertices[k, i, j + 1], - ], - ] - - # fill in the output lists - if self.nz == 0: - verts.append([pt1[0], pt1[1], zvals[0][0]]) - verts.append([pt2[0], pt2[1], zvals[0][1]]) - verts.append([pt0[0], pt0[1], zvals[0][2]]) - verts.append([pt3[0], pt3[1], zvals[0][3]]) - ivert.extend( - [ipoint, ipoint + 1, ipoint + 2, ipoint + 3] - ) - zverts.extend( - [ - zvals[0][0], - zvals[0][1], - zvals[0][2], - zvals[0][3], - ] - ) - ipoint += 4 - elif self.ny == 0: - verts.append([pt1[0], pt1[1], zvals[0][0]]) - verts.append([pt2[0], pt2[1], zvals[0][1]]) - verts.append([pt1[0], pt1[1], zvals[1][0]]) - verts.append([pt2[0], pt2[1], zvals[1][1]]) - ivert.extend( - [ipoint, ipoint + 1, ipoint + 2, ipoint + 3] - ) - zverts.extend( - [ - zvals[0][0], - zvals[0][1], - zvals[1][0], - zvals[1][1], - ] - ) - ipoint += 4 - elif self.nx == 0: - verts.append([pt1[0], pt1[1], zvals[0][0]]) - verts.append([pt0[0], pt0[1], zvals[0][2]]) - verts.append([pt1[0], pt1[1], zvals[1][0]]) - verts.append([pt0[0], pt0[1], zvals[1][2]]) - ivert.extend( - [ipoint, ipoint + 1, ipoint + 2, ipoint + 3] - ) - zverts.extend( - [ - zvals[0][0], - zvals[0][2], - zvals[1][0], - zvals[1][2], - ] - ) - ipoint += 4 - else: - for zvals_l in zvals: - verts.append([pt1[0], pt1[1], zvals_l[0]]) - verts.append([pt2[0], pt2[1], zvals_l[1]]) - verts.append([pt0[0], pt0[1], zvals_l[2]]) - verts.append([pt3[0], pt3[1], zvals_l[3]]) - ivert.extend( - [ipoint, ipoint + 1, ipoint + 2, ipoint + 3] - ) - zverts.extend(zvals_l) - ipoint += 4 - vertsdict[cellid] = verts - ivertsdict[cellid] = ivert - zvertsdict[cellid] = zverts - return vertsdict, ivertsdict, zvertsdict - - -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 - - -def export_cbc( - model, - cbcfile, - otfolder, - precision="single", - verbose=False, - nanval=-1e20, - kstpkper=None, - text=None, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=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 - Precision of data in the cell by cell file: 'single' or 'double'. - Default is 'single'. - verbose : bool - If True, write information to the screen. Default is False. - 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, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output 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 - pvdfilename = f"{model.name}_CBC.pvd" - pvdfile = open(os.path.join(otfolder, pvdfilename), "w") - - pvdfile.write( - """ - - \n""" - ) - - # load cbc - cbb = bf.CellBudgetFile(cbcfile, precision=precision, verbose=verbose) - - # 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 - - # get the export times - if kstpkper is not None: - if isinstance(kstpkper, tuple): - kstpkper = [kstpkper] - elif not isinstance(kstpkper, list) or not isinstance( - kstpkper[0], tuple - ): - raise Exception( - "kstpkper must be a tuple (kstp, kper) or a list of tuples" - ) - else: - kstpkper = cbb.get_kstpkper() - - # get model name - model_name = model.name - - vtk = Vtk( - model, - nanval=nanval, - smooth=smooth, - point_scalars=point_scalars, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - - # export data - addarray = False - count = 1 - for kstpkper_i in kstpkper: - ot_base = ( - f"{model_name}_CBC_KPER{kstpkper_i[1] + 1}_KSTP{kstpkper_i[0] + 1}" - ) - otfile = os.path.join(otfolder, ot_base) - pvdfile.write( - f""" -""" - ) - for name in keylist: - - try: - rec = cbb.get_data(kstpkper=kstpkper_i, 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=kstpkper_i, 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 currently supported for cbc output" - ) - # print('Data type not currently 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, - text="head", - precision="auto", - verbose=False, - nanval=-1e20, - kstpkper=None, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=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 - text : string - Name of the text string in the head file. Default is 'head'. - precision : str - Precision of data in the head file: 'auto', 'single' or 'double'. - Default is 'auto'. - verbose : bool - If True, write information to the screen. Default is False. - 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, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output 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 - pvdfilename = f"{model.name}_{text}.pvd" - pvdfile = open(os.path.join(otfolder, pvdfilename), "w") - - pvdfile.write( - """ - - \n""" - ) - - # get the heads - hds = HeadFile(hdsfile, text=text, precision=precision, verbose=verbose) - - # get the export times - if kstpkper is not None: - if isinstance(kstpkper, tuple): - kstpkper = [kstpkper] - elif not isinstance(kstpkper, list) or not isinstance( - kstpkper[0], tuple - ): - raise Exception( - "kstpkper must be a tuple (kstp, kper) or a list of tuples" - ) - else: - kstpkper = hds.get_kstpkper() - - # set upt the vtk - vtk = Vtk( - model, - smooth=smooth, - point_scalars=point_scalars, - nanval=nanval, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - - # output data - count = 0 - for kstpkper_i in kstpkper: - hdarr = hds.get_data(kstpkper_i) - vtk.add_array(text, hdarr) - ot_base = ( - f"{model.name}_{text}_" - f"KPER{kstpkper_i[1] + 1}_KSTP{kstpkper_i[0] + 1}" - ) - otfile = os.path.join(otfolder, ot_base) - # vtk.write(otfile, timeval=totim_dict[(kstp, kper)]) - vtk.write(otfile) - pvdfile.write( - f""" -""" - ) - count += 1 - - pvdfile.write( - """ -""" - ) - - pvdfile.close() - - -def export_array( - model, - array, - output_folder, - name, - nanval=-1e20, - array2d=False, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=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 the array is 2d and represents the first layer, default is - False - smooth : bool - if True, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output 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_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - vtk.add_array(name, array, array2d=array2d) - otfile = os.path.join(output_folder, name) - vtk.write(otfile) - - return - - -def export_vector( - model, - vector, - output_folder, - name, - nanval=-1e20, - array2d=False, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=False, - binary=False, -): - - """ - - Export vector (i.e., a tuple of arrays) to vtk - - Parameters - ---------- - - model : flopy model instance - the flopy model instance - vector : tuple of arrays - vector to be exported - output_folder : str - output folder to write the data - name : str - name of vector - nanval : scalar - no data value, default value is -1e20 - array2d : bool - true if the vector components are 2d arrays and represent the first - layer, default is False - smooth : bool - if True, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output 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_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - vtk.add_vector(name, vector, array2d=array2d) - otfile = os.path.join(output_folder, name) - vtk.write(otfile) - - return - - -def export_transient( - model, - array, - output_folder, - name, - nanval=-1e20, - array2d=False, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=False, - binary=False, - kpers=None, -): - """ - 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, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output file will be binary, default is False - kpers : iterable of int - Stress periods to export. If None (default), all stress periods will be - exported. - """ - - 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, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - - if name.endswith("_"): - separator = "" - else: - separator = "_" - - if kpers is None: - kpers = range(array.shape[0]) - else: - assert isinstance(kpers, list) or isinstance(kpers, np.ndarray) - - if array2d: - for kper in kpers: - - 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) - - otname = f"{name}{separator}0{kper + 1}" - otfile = os.path.join(output_folder, otname) - vtk.write(otfile, timeval=to_tim[kper]) - - else: - for kper in kpers: - vtk.add_array(name, array[kper]) - - otname = f"{name}{separator}0{kper + 1}" - otfile = os.path.join(output_folder, otname) - vtk.write(otfile, 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, - otfolder, - vtkobj=None, - nanval=-1e20, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=False, - binary=False, - kpers=None, -): - """ - Exports package to vtk - - Parameters - ---------- - - pak_model : flopy model instance - the model of the package - pak_name : str - the name of the package - otfolder : 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, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output file will be binary, default is False - kpers : iterable of int - Stress periods to export. If None (default), all stress periods will be - exported. - """ - - # 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, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - else: - # otherwise use the vtk object that was supplied - vtk = vtkobj - - if not os.path.exists(otfolder): - os.mkdir(otfolder) - - # 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: - otfile = os.path.join(otfolder, pak_name) - vtk.write(otfile) - - # write transient data - if vtk_trans_dict: - - # only retain requested stress periods - if kpers is not None: - assert isinstance(kpers, list) or isinstance(kpers, np.ndarray) - vtk_trans_dict = {kper: vtk_trans_dict[kper] for kper in kpers} - - # 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 - otfile = os.path.join(otfolder, f"{pak_name}_0{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(otfile, timeval=time) - vtk.write(otfile) - return - - -def export_model( - model, - otfolder, - package_names=None, - nanval=-1e20, - smooth=False, - point_scalars=False, - vtk_grid_type="auto", - true2d=False, - binary=False, - kpers=None, -): - """ - Exports model to vtk - - Parameters - ---------- - - 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, will create smooth layer elevations, default is False - point_scalars : bool - if True, will also output array values at cell vertices, default is - False; note this automatically sets smooth to True - vtk_grid_type : str - Specific vtk_grid_type or 'auto' (default). Possible specific values - are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. - If 'auto', the grid type is automatically determined. Namely: - * A regular grid (in all three directions) will be saved as an - 'ImageData'. - * A rectilinear (in all three directions), non-regular grid - will be saved as a 'RectilinearGrid'. - * Other grids will be saved as 'UnstructuredGrid'. - true2d : bool - If True, the model is expected to be 2d (1 layer, 1 row or 1 column) - and the data will be exported as true 2d data, default is False. - binary : bool - if True the output file will be binary, default is False - kpers : iterable of int - Stress periods to export. If None (default), all stress periods will be - exported. - """ - vtk = Vtk( - model, - nanval=nanval, - smooth=smooth, - point_scalars=point_scalars, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - ) - - 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 model.packagelist] - - if not os.path.exists(otfolder): - os.mkdir(otfolder) - - for pak_name in package_names: - export_package( - model, - pak_name, - otfolder, - vtkobj=vtk, - nanval=nanval, - smooth=smooth, - point_scalars=point_scalars, - vtk_grid_type=vtk_grid_type, - true2d=true2d, - binary=binary, - kpers=kpers, - ) +import numpy as np +import os +from flopy.datbase import DataType, DataInterface +from flopy.utils import Util3d +import warnings + +warnings.simplefilter("always", DeprecationWarning) + + +VTKIGNORE = ("vkcb", "perlen", "steady", "tsmult", "nstp") + + +class Pvd: + """ + Simple class to build a Paraview Data File (PVD) + """ + + def __init__(self): + self.__data = [ + '\n', + '\n', + "\n", + ] + + def add_timevalue(self, file, timevalue): + """ + Method to add a Dataset record to the pvd file + + Parameters + ---------- + file : str + vtu file name + timevalue : float + time step value in model time + """ + if not file.endswith(".vtu"): + raise AssertionError("File name must end with .vtu ") + + file = os.path.split(file)[-1] + + record = ( + f'\n' + ) + self.__data.append(record) + + def write(self, f): + """ + Method to write a pvd file from the PVD object + + Paramaters + ---------- + f : str + PVD file name + + """ + if not f.endswith(".pvd"): + f += ".pvd" + + with open(f, "w") as foo: + foo.writelines(self.__data) + foo.write("\n") + foo.write("") + + +class Vtk: + """ + Class that builds VTK objects and exports models to VTK files + + + Parameters: + ---------- + model : flopy.ModelInterface object + any flopy model object, example flopy.modflow.Modflow() object + modelgrid : flopy.discretization.Grid object + any flopy modelgrid object, example. VertexGrid + vertical_exageration : float + floating point value to scale vertical exageration of the vtk points + default is 1. + binary : bool + flag that indicates if Vtk will write a binary or text file. Binary + is prefered as paraview has a bug (8/4/2021) where is cannot represent + NaN values from ASCII (non xml) files. In this case no-data values + are set to 1e+30. + xml : bool + flag to write xml based VTK files + pvd : bool + boolean flag to write a paraview pvd file for transient data. This + file maps vtk files to a model time. + shared_points : bool + boolean flag to share points in grid polyhedron construction. Default + is False, as paraview has a bug (8/4/2021) where some polyhedrons will + not properly close when using shared points. If shared_points is True + file size will be slightly smaller. + smooth : bool + boolean flag to interpolate vertex elevations based on shared cell + elevations. Default is False. + point_scalars : bool + boolen flag to write interpolated data at each point based "shared + vertices". + + """ + + def __init__( + self, + model=None, + modelgrid=None, + vertical_exageration=1, + binary=True, + xml=False, + pvd=False, + shared_points=False, + smooth=False, + point_scalars=False, + ): + + try: + import vtk + except ImportError: + err = "vtk not installed, use pip install vtk" + raise ImportError(err) + + if model is None and modelgrid is None: + raise AssertionError( + "A model or modelgrid must be provided to use Vtk" + ) + + elif model is not None: + self.modelgrid = model.modelgrid + self.modeltime = model.modeltime + else: + self.modelgrid = modelgrid + self.modeltime = None + + self.binary = binary + self.xml = xml + self.pvd = pvd + + if self.pvd and not self.xml: + print( + "Switching to xml, ASCII and standard binary are not " + "supported by Paraview's PVD reader" + ) + self.xml = True + + self.vertical_exageration = vertical_exageration + self.shared_points = shared_points + self.smooth = smooth + self.point_scalars = point_scalars + + self.verts = self.modelgrid.verts + self.iverts = self.modelgrid.iverts + self.nnodes = self.modelgrid.nnodes + + nvpl = 0 + for iv in self.iverts: + nvpl += len(iv) - 1 + + self.nvpl = nvpl + + # method to accomodate DISU grids, do not use modelgrid.ncpl! + self.ncpl = len(self.iverts) + if self.nnodes == len(self.iverts): + self.nlay = 1 + else: + self.nlay = self.modelgrid.nlay + + self._laycbd = self.modelgrid.laycbd + if self._laycbd is None: + self._laycbd = np.zeros((self.nlay,), dtype=int) + + self._active = [] + self._ncbd = 0 + for i in range(self.nlay): + self._active.append(1) + if self._laycbd[i] != 0: + self._active.append(0) + self._ncbd += 1 + + # USG trap + if self.modelgrid.top.size == self.nnodes: + self.top = self.modelgrid.top.reshape(self.nnodes) + else: + self.top = self.modelgrid.top.reshape(self.ncpl) + self.botm = self.modelgrid.botm.reshape(-1, self.ncpl) + + if self.modelgrid.idomain is not None: + self.idomain = self.modelgrid.idomain.reshape(self.nnodes) + else: + self.idomain = np.ones((self.nnodes,), dtype=int) + + if self.modeltime is not None: + perlen = self.modeltime.perlen + self._totim = np.add.accumulate(perlen) + + self.points = [] + self.faces = [] + self.vtk_grid = None + self.vtk_polygons = None + self.vtk_pathlines = None + self._pathline_points = [] + self._point_scalar_graph = None + self._point_scalar_numpy_graph = None + self._idw_weight_graph = None + self._idw_total_weight_graph = None + + self._vtk_geometry_set = False + self.__transient_output_data = False + self.__transient_data = {} + self.__transient_vector = {} + self.__pathline_transient_data = {} + self.__vtk = vtk + + if self.point_scalars: + self._create_point_scalar_graph() + + def _create_smoothed_elevation_graph(self, adjk, top_elev=True): + """ + Method to create a dictionary of shared point elevations + + Parameters + ---------- + adjk : int + confining bed adjusted layer + + Returns + ------- + dict : {vertex number: elevation} + + """ + elevations = {} + for i, iv in enumerate(self.iverts): + iv = iv[1:] + for v in iv: + if v is None: + continue + if not top_elev: + zv = self.botm[adjk][i] * self.vertical_exageration + elif adjk == 0: + zv = self.top[i] * self.vertical_exageration + else: + if self.top.size == self.nnodes: + adji = (adjk * self.ncpl) + i + zv = self.top[adji] * self.vertical_exageration + else: + zv = self.botm[adjk - 1][i] * self.vertical_exageration + + if v in elevations: + elevations[v].append(zv) + else: + elevations[v] = [zv] + + for key in elevations: + elevations[key] = np.mean(elevations[key]) + + return elevations + + def _create_point_scalar_graph(self): + """ + Method to create a point scalar graph to map cells to points + + """ + graph = {} + v0 = 0 + v1 = 0 + nvert = len(self.verts) + shared_points = self.shared_points + if len(self._active) != self.nlay: + shared_points = False + + for k in range(self.nlay): + for i, iv in enumerate(self.iverts): + iv = iv[1:] + adji = (k * self.ncpl) + i + for v in iv: + if v is None: + continue + xvert = self.verts[v, 0] + yvert = self.verts[v, 1] + adjv = v + v1 + if adjv not in graph: + graph[adjv] = { + "vtk_points": [v0], + "idx": [adji], + "xv": [xvert], + "yv": [yvert], + } + else: + graph[adjv]["vtk_points"].append(v0) + graph[adjv]["idx"].append(adji) + graph[adjv]["xv"].append(xvert) + graph[adjv]["yv"].append(yvert) + v0 += 1 + + v1 += nvert + + if k == self.nlay - 1 or not shared_points: + for i, iv in enumerate(self.iverts): + iv = iv[1:] + adji = (k * self.ncpl) + i + for v in iv: + if v is None: + continue + xvert = self.verts[v, 0] + yvert = self.verts[v, 1] + adjv = v + v1 + if adjv not in graph: + graph[adjv] = { + "vtk_points": [v0], + "idx": [adji], + "xv": [xvert], + "yv": [yvert], + } + else: + graph[adjv]["vtk_points"].append(v0) + graph[adjv]["idx"].append(adji) + graph[adjv]["xv"].append(xvert) + graph[adjv]["yv"].append(yvert) + v0 += 1 + v1 += nvert + + # convert this to a numpy representation + max_shared = 0 + num_points = v0 + for k, d in graph.items(): + if len(d["vtk_points"]) > max_shared: + max_shared = len(d["vtk_points"]) + + numpy_graph = np.ones((max_shared, num_points), dtype=int) * -1 + xvert = np.ones((max_shared, num_points)) * np.nan + yvert = np.ones((max_shared, num_points)) * np.nan + for k, d in graph.items(): + for ix, pt in enumerate(d["vtk_points"]): + for ixx, value in enumerate(d["idx"]): + numpy_graph[ixx, pt] = value + xvert[ixx, pt] = d["xv"][ixx] + yvert[ixx, pt] = d["yv"][ixx] + + # now create the IDW weights for point scalars + xc = np.ravel(self.modelgrid.xcellcenters) + yc = np.ravel(self.modelgrid.ycellcenters) + + if xc.size != self.nnodes: + xc = np.tile(xc, self.nlay) + yc = np.tile(yc, self.nlay) + + xc = np.where(numpy_graph != -1, xc[numpy_graph], np.nan) + yc = np.where(numpy_graph != -1, yc[numpy_graph], np.nan) + + asq = (xvert - xc) ** 2 + bsq = (yvert - yc) ** 2 + + weights = np.sqrt(asq + bsq) + tot_weights = np.nansum(weights, axis=0) + + self._point_scalar_graph = graph + self._point_scalar_numpy_graph = numpy_graph + self._idw_weight_graph = weights + self._idw_total_weight_graph = tot_weights + + def _build_grid_geometry(self): + """ + Method that creates lists of vertex points and cell faces + + """ + points = [] + faces = [] + v0 = 0 + v1 = 0 + ncb = 0 + shared_points = self.shared_points + if len(self._active) != self.nlay: + shared_points = False + + for k in range(self.nlay): + adjk = k + ncb + if k != self.nlay - 1: + if self._active[adjk + 1] == 0: + ncb += 1 + + if self.smooth: + elevations = self._create_smoothed_elevation_graph(adjk) + + for i, iv in enumerate(self.iverts): + iv = iv[1:] + for v in iv: + if v is None: + continue + xv = self.verts[v, 0] + yv = self.verts[v, 1] + if self.smooth: + zv = elevations[v] + elif k == 0: + zv = self.top[i] * self.vertical_exageration + else: + if self.top.size == self.nnodes: + adji = (adjk * self.ncpl) + i + zv = self.top[adji] * self.vertical_exageration + else: + zv = ( + self.botm[adjk - 1][i] + * self.vertical_exageration + ) + + points.append([xv, yv, zv]) + v1 += 1 + + cell_faces = [ + [v for v in range(v0, v1)], + [v + self.nvpl for v in range(v0, v1)], + ] + + for v in range(v0, v1): + if v != v1 - 1: + cell_faces.append( + [v + 1, v, v + self.nvpl, v + self.nvpl + 1] + ) + else: + cell_faces.append( + [v0, v, v + self.nvpl, v0 + self.nvpl] + ) + + v0 = v1 + faces.append(cell_faces) + + if k == self.nlay - 1 or not shared_points: + if self.smooth: + elevations = self._create_smoothed_elevation_graph( + adjk, top_elev=False + ) + + for i, iv in enumerate(self.iverts): + iv = iv[1:] + for v in iv: + if v is None: + continue + xv = self.verts[v, 0] + yv = self.verts[v, 1] + if self.smooth: + zv = elevations[v] + else: + zv = self.botm[adjk][i] * self.vertical_exageration + + points.append([xv, yv, zv]) + v1 += 1 + + v0 = v1 + + self.points = points + self.faces = faces + + def _set_vtk_grid_geometry(self): + """ + Method to set vtk's geometry and add it to the vtk grid object + + """ + if self._vtk_geometry_set: + return + + if not self.faces: + self._build_grid_geometry() + + self.vtk_grid = self.__vtk.vtkUnstructuredGrid() + + points = self.__vtk.vtkPoints() + for point in self.points: + points.InsertNextPoint(point) + + self.vtk_grid.SetPoints(points) + + for node in range(self.nnodes): + cell_faces = self.faces[node] + nface = len(cell_faces) + fid_list = self.__vtk.vtkIdList() + fid_list.InsertNextId(nface) + for face in cell_faces: + fid_list.InsertNextId(len(face)) + [fid_list.InsertNextId(i) for i in face] + + self.vtk_grid.InsertNextCell(self.__vtk.VTK_POLYHEDRON, fid_list) + + self._vtk_geometry_set = True + + def _build_hfbs(self, pkg): + """ + Method to add hfb planes to the vtk object + + Parameters + ---------- + pkg : object + flopy hfb object + + """ + from vtk.util import numpy_support + + # check if modflow 6 or modflow 2005 + if hasattr(pkg, "hfb_data"): + mf6 = False + hfb_data = pkg.hfb_data + else: + # asssume that there is no transient hfb data for now + hfb_data = pkg.stress_period_data.array[0] + mf6 = True + + points = [] + faces = [] + array = [] + cnt = 0 + verts = self.modelgrid.cross_section_vertices + verts = np.dstack((verts[0], verts[1])) + botm = np.ravel(self.botm) + for hfb in hfb_data: + if self.modelgrid.grid_type == "structured": + if not mf6: + k = hfb["k"] + cell1 = list(hfb[["k", "irow1", "icol1"]]) + cell2 = list(hfb[["k", "irow2", "icol2"]]) + else: + cell1 = hfb["cellid1"] + cell2 = hfb["cellid2"] + k = cell1[0] + n0, n1 = self.modelgrid.get_node([cell1, cell2]) + + else: + cell1 = hfb["cellid1"] + cell2 = hfb["cellid2"] + if len(cell1) == 2: + k, n0 = cell1 + _, n1 = cell2 + else: + n0 = cell1[0] + n1 = cell1[1] + k = 0 + + array.append(hfb["hydchr"]) + adjn0 = n0 - (k * self.ncpl) + adjn1 = n1 - (k * self.ncpl) + v1 = verts[adjn0] + v2 = verts[adjn1] + + # get top and botm elevations, use max and min: + if k == 0 or self.top.size == self.nnodes: + tv = np.max([self.top[n0], self.top[n1]]) + bv = np.min([botm[n0], botm[n1]]) + else: + tv = np.max([botm[n0 - self.ncpl], botm[n1 - self.ncpl]]) + bv = np.min([botm[n0], botm[n1]]) + + tv *= self.vertical_exageration + bv *= self.vertical_exageration + + pts = [] + for v in v1: + # ix = np.where(v2 == v) + ix = np.where((v2.T[0] == v[0]) & (v2.T[1] == v[1])) + if len(ix[0]) > 0 and len(pts) < 2: + pts.append(v2[ix[0][0]]) + + pts = np.sort(pts)[::-1] + + # create plane, counter-clockwise order + pts = [ + (pts[0, 0], pts[0, 1], tv), + (pts[1, 0], pts[1, 1], tv), + (pts[1, 0], pts[1, 1], bv), + (pts[0, 0], pts[0, 1], bv), + ] + + # add to points and faces + plane_face = [] + for pt in pts: + points.append(pt) + plane_face.append(cnt) + cnt += 1 + + faces.append(plane_face) + + # now create the vtk geometry + vtk_points = self.__vtk.vtkPoints() + for point in points: + vtk_points.InsertNextPoint(point) + + # now create an UnstructuredGrid object + polydata = self.__vtk.vtkUnstructuredGrid() + polydata.SetPoints(vtk_points) + + # create the vtk polygons + for face in faces: + polygon = self.__vtk.vtkPolygon() + polygon.GetPointIds().SetNumberOfIds(4) + for ix, iv in enumerate(face): + polygon.GetPointIds().SetId(ix, iv) + polydata.InsertNextCell( + polygon.GetCellType(), polygon.GetPointIds() + ) + + # and then set the hydchr data + vtk_arr = numpy_support.numpy_to_vtk( + num_array=np.array(array), array_type=self.__vtk.VTK_FLOAT + ) + vtk_arr.SetName("hydchr") + polydata.GetCellData().SetScalars(vtk_arr) + self.vtk_polygons = polydata + + def _build_point_scalar_array(self, array): + """ + Method to build a point scalar array using an inverse distance + weighting scheme + + Parameters + ---------- + array : np.ndarray + ndarray of shape nnodes + + Returns + ------- + np.ndarray of shape n vtk points + """ + isint = False + if np.issubdtype(array[0], np.dtype(int)): + isint = True + + if isint: + # maybe think of a way to do ints with the numpy graph, + # looping through the dict is very inefficient + ps_array = np.zeros((len(self.points),), dtype=array.dtype) + for _, value in self._point_scalar_graph.items(): + for ix, pt in enumerate(value["vtk_points"]): + ps_array[pt] = array[value["idx"][ix]] + else: + ps_graph = self._point_scalar_numpy_graph + ps_array = np.where(ps_graph >= 0, array[ps_graph], np.nan) + + # do inverse distance weighting and apply mask to retain + # nan valued cells because numpy returns 0 when all vals are nan + weighted_vals = self._idw_weight_graph * ps_array + mask = np.isnan(weighted_vals).all(axis=0) + weighted_vals = np.nansum(weighted_vals, axis=0) + weighted_vals[mask] = np.nan + ps_array = weighted_vals / self._idw_total_weight_graph + + return ps_array + + def _add_timevalue(self, index, fname): + """ + Method to add a TimeValue to a vtk object, used with transient arrays + + Parameters + ---------- + index : int, tuple + integer representing kper or a tuple of (kstp, kper) + fname : str + path to the vtu file + + """ + if not self.pvd: + return + + try: + timeval = self._totim[index] + except (IndexError, KeyError): + return + + self.pvd.add_timevalue(fname, timeval) + + def _mask_values(self, array, masked_values): + """ + Method to mask values with nan + + Parameters + ---------- + array : np.ndarray + numpy array of values + masked_values : list + values to convert to nan + + Returns + ------- + numpy array + """ + if masked_values is not None: + try: + for mv in masked_values: + array[array == mv] = np.nan + except ValueError: + pass + + return array + + def add_array(self, array, name, masked_values=None, dtype=None): + """ + Method to set an array to the vtk grid + + Parameters + ---------- + array : np.ndarray, list + list of array values to set to the vtk_array + name : str + array name for vtk + masked_values : list, None + list of values to set equal to nan + dtype : vtk datatype object + method to supply and force a vtk datatype + + """ + from vtk.util import numpy_support + + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + array = np.ravel(array) + + if array.size != self.nnodes: + raise AssertionError("array must be the size as the modelgrid") + + if self.idomain is not None: + try: + array[self.idomain == 0] = np.nan + except ValueError: + pass + + array = self._mask_values(array, masked_values) + + if not self.binary and not self.xml: + # ascii does not properly render nan values + array = np.nan_to_num(array, nan=1e30) + + if self.point_scalars: + array = self._build_point_scalar_array(array) + + if dtype is None: + dtype = self.__vtk.VTK_FLOAT + if np.issubdtype(array[0], np.dtype(int)): + dtype = self.__vtk.VTK_INT + + vtk_arr = numpy_support.numpy_to_vtk(num_array=array, array_type=dtype) + vtk_arr.SetName(name) + + if self.point_scalars: + self.vtk_grid.GetPointData().AddArray(vtk_arr) + else: + self.vtk_grid.GetCellData().AddArray(vtk_arr) + + def add_transient_array(self, d, name=None, masked_values=None): + """ + Method to add transient array data to the vtk object + + Parameters + ---------- + d: dict + dictionary of array2d, arry3d data or numpy array data + name : str, None + parameter name, required when user provides a dictionary + of numpy arrays + masked_values : list, None + list of values to set equal to nan + + Returns + ------- + + """ + if self.__transient_output_data: + raise AssertionError( + "Transient arrays cannot be mixed with transient output, " + "Please create a seperate vtk object for transient package " + "data" + ) + + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + k = list(d.keys())[0] + transient = dict() + if isinstance(d[k], DataInterface): + if d[k].data_type in (DataType.array2d, DataType.array3d): + if name is None: + name = d[k].name + if isinstance(name, list): + name = name[0] + + for kper, value in d.items(): + if value.array.size != self.nnodes: + array = np.zeros(self.nnodes) * np.nan + array[: value.array.size] = np.ravel(value.array) + else: + array = value.array + + array = self._mask_values(array, masked_values) + transient[kper] = array + else: + if name is None: + raise ValueError( + "name must be specified when providing numpy arrays" + ) + for kper, trarray in d.items(): + if trarray.size != self.nnodes: + array = np.zeros(self.nnodes) * np.nan + array[: trarray.size] = np.ravel(trarray) + else: + array = trarray + + array = self._mask_values(array, masked_values) + transient[kper] = array + + for k, v in transient.items(): + if k not in self.__transient_data: + self.__transient_data[k] = {name: v} + else: + self.__transient_data[k][name] = v + + def add_transient_list(self, mflist, masked_values=None): + """ + Method to add transient list data to a vtk object + + Parameters + ---------- + mflist : flopy.utils.MfList object + masked_values : list, None + list of values to set equal to nan + + """ + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + pkg_name = mflist.package.name[0] + mfl = mflist.array + if isinstance(mfl, dict): + for arr_name, arr4d in mflist.array.items(): + d = {kper: array for kper, array in enumerate(arr4d)} + name = f"{pkg_name}_{arr_name}" + self.add_transient_array(d, name) + else: + export = {} + for kper in range(mflist.package.parent.nper): + try: + arr_dict = mflist.to_array(kper, mask=True) + except ValueError: + continue + + if arr_dict is None: + continue + + for arr_name, array in arr_dict.items(): + if arr_name not in export: + export[arr_name] = {kper: array} + else: + export[arr_name][kper] = array + + for arr_name, d in export.items(): + name = f"{pkg_name}_{arr_name}" + self.add_transient_array(d, name, masked_values=masked_values) + + def add_vector(self, vector, name, masked_values=None): + """ + Method to add vector data to vtk + + Parameters + ---------- + vector : array + array of dimension (3, nnodes) + name : str + name of the vector to be displayed in vtk + masked_values : list, None + list of values to set equal to nan + + """ + from vtk.util import numpy_support + + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + if isinstance(vector, (tuple, list)): + vector = np.array(vector) + + if vector.size != 3 * self.nnodes: + if vector.size == 3 * self.ncpl: + vector = np.reshape(vector, (3, self.ncpl)) + tv = np.full((3, self.nnodes), np.nan) + for ix, q in enumerate(vector): + tv[ix, : self.ncpl] = q + vector = tv + else: + raise AssertionError( + "Size of vector must be 3 * nnodes or 3 * ncpl" + ) + else: + vector = np.reshape(vector, (3, self.nnodes)) + + if self.point_scalars: + tmp = [] + for v in vector: + tmp.append(self._build_point_scalar_array(v)) + vector = np.array(tmp) + + vector = self._mask_values(vector, masked_values) + + vtk_arr = numpy_support.numpy_to_vtk( + num_array=vector, array_type=self.__vtk.VTK_FLOAT + ) + vtk_arr.SetName(name) + vtk_arr.SetNumberOfComponents(3) + + if self.point_scalars: + self.vtk_grid.GetPointData().SetVectors(vtk_arr) + else: + self.vtk_grid.GetCellData().SetVectors(vtk_arr) + + def add_transient_vector(self, d, name, masked_values=None): + """ + Method to add transient vector data to vtk + + Paramters + --------- + d : dict + dictionary of vectors + name : str + name of vector to be displayed in vtk + masked_values : list, None + list of values to set equal to nan + + """ + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + if self.__transient_data: + k = list(self.__transient_data.keys())[0] + if len(d) != len(self.__transient_data[k]): + print( + "Transient vector not same size as transient arrays time " + "stamp may be unreliable for vector data in VTK file" + ) + + if isinstance(d, dict): + cnt = 0 + for key, value in d.items(): + if not isinstance(value, np.ndarray): + value = np.array(value) + + if ( + value.size != 3 * self.ncpl + or value.size != 3 * self.nnodes + ): + raise AssertionError( + "Size of vector must be 3 * nnodes or 3 * ncpl" + ) + + value = self._mask_values(value, masked_values) + self.__transient_vector[cnt] = {name: value} + cnt += 1 + + def add_package(self, pkg, masked_values=None): + """ + Method to set all arrays within a package to VTK arrays + + Parameters + ---------- + pkg : flopy.pakbase.Package object + flopy package object, example ModflowWel + masked_values : list, None + list of values to set equal to nan + + """ + if not self._vtk_geometry_set: + self._set_vtk_grid_geometry() + + if "hfb" in pkg.name[0].lower(): + self._build_hfbs(pkg) + return + + for item, value in pkg.__dict__.items(): + if item in VTKIGNORE: + continue + + if isinstance(value, list): + for v in value: + if isinstance(v, Util3d): + if value.array.size != self.nnodes: + continue + self.add_array(v.array, item, masked_values) + + if isinstance(value, DataInterface): + if value.data_type in (DataType.array2d, DataType.array3d): + if value.array is not None: + if value.array.size < self.nnodes: + if value.array.size < self.ncpl: + continue + + array = np.zeros(self.nnodes) * np.nan + array[: value.array.size] = np.ravel(value.array) + + elif value.array.size > self.nnodes and self._ncbd > 0: + # deal with confining beds + array = np.array( + [ + value.array[ix] + for ix, i in enumerate(self._active) + if i != 0 + ] + ) + + else: + array = value.array + + self.add_array(array, item, masked_values) + + elif value.data_type == DataType.transient2d: + if value.array is not None: + if hasattr(value, "transient_2ds"): + self.add_transient_array( + value.transient_2ds, item, masked_values + ) + else: + d = {ix: i for ix, i in enumerate(value.array)} + self.add_transient_array(d, item, masked_values) + + elif value.data_type == DataType.transient3d: + if value.array is not None: + self.add_transient_array( + value.transient_3ds, item, masked_values + ) + + elif value.data_type == DataType.transientlist: + self.add_transient_list(value, masked_values) + + else: + pass + + def add_model(self, model, selpaklist=None, masked_values=None): + """ + Method to add all array and transient list data from a modflow model + to a timeseries of vtk files + + Parameters + ---------- + model : fp.modflow.ModelInterface + any flopy model object + selpaklist : list, None + optional parameter where the user can supply a list of packages + to export. + masked_values : list, None + list of values to set equal to nan + + """ + for package in model.packagelist: + if selpaklist is not None: + if package.name[0] not in selpaklist: + continue + + self.add_package(package, masked_values) + + def add_pathline_points(self, pathlines, timeseries=False): + """ + Method to add Modpath output from a pathline or timeseries file + to the grid. Colors will be representative of totim. + + Parameters: + ---------- + pathlines : np.recarray or list + pathlines accepts a numpy recarray of a particle pathline or + a list of numpy reccarrays associated with pathlines + + timeseries : bool + method to plot data as a series of vtk timeseries files for + animation or as a single static vtk file. Default is false + + """ + if isinstance(pathlines, (np.recarray, np.ndarray)): + pathlines = [pathlines] + + keys = ["particleid", "time"] + if not timeseries: + arrays = {key: [] for key in keys} + points = [] + for recarray in pathlines: + recarray["z"] *= self.vertical_exageration + for rec in recarray: + points.append(tuple(rec[["x", "y", "z"]])) + for key in keys: + arrays[key].append(rec[key]) + + self._set_modpath_point_data(points, arrays) + + else: + self.vtk_pathlines = self.__vtk.vtkUnstructuredGrid() + timeseries_data = {} + points = {} + for recarray in pathlines: + recarray["z"] *= self.vertical_exageration + for rec in recarray: + time = rec["time"] + if time not in points: + points[time] = [tuple(rec[["x", "y", "z"]])] + t = {key: [] for key in keys} + timeseries_data[time] = t + + else: + points[time].append(tuple(rec[["x", "y", "z"]])) + + for key in keys: + timeseries_data[time][key].append(rec[key]) + + self.__pathline_transient_data = timeseries_data + self._pathline_points = points + + def add_heads(self, hds, kstpkper=None, masked_values=None): + """ + Method to add head data to a vtk file + + Parameters + ---------- + hds : flopy.utils.LayerFile object + Binary or Formatted HeadFile type object + kstpkper : tuple, list of tuples, None + tuple or list of tuples of kstpkper, if kstpkper=None all + records are selected + masked_values : list, None + list of values to set equal to nan + + """ + if not self.__transient_output_data and self.__transient_data: + raise AssertionError( + "Head data cannot be mixed with transient package data, " + "Please create a seperate vtk object for transient head data" + ) + + if kstpkper is None: + kstpkper = hds.get_kstpkper() + elif isinstance(kstpkper, (list, tuple)): + if not isinstance(kstpkper[0], (list, tuple)): + kstpkper = [kstpkper] + else: + pass + + # reset totim based on values read from head file + times = hds.get_times() + kstpkpers = hds.get_kstpkper() + self._totim = {ki: time for (ki, time) in zip(kstpkpers, times)} + + text = hds.text.decode() + + d = dict() + for ki in kstpkper: + d[ki] = hds.get_data(ki) + + self.__transient_output_data = False + self.add_transient_array(d, name=text, masked_values=masked_values) + self.__transient_output_data = True + + def add_cell_budget( + self, cbc, text=None, kstpkper=None, masked_values=None + ): + """ + Method to add cell budget data to vtk + + Parameters + ---------- + cbc : flopy.utils.CellBudget object + flopy binary CellBudget object + text : str or None + The text identifier for the record. Examples include + 'RIVER LEAKAGE', 'STORAGE', 'FLOW RIGHT FACE', etc. If text=None + all compatible records are exported + kstpkper : tuple, list of tuples, None + tuple or list of tuples of kstpkper, if kstpkper=None all + records are selected + masked_values : list, None + list of values to set equal to nan + + """ + if not self.__transient_output_data and self.__transient_data: + raise AssertionError( + "Binary data cannot be mixed with transient package data, " + "Please create a seperate vtk object for transient head data" + ) + + records = cbc.get_unique_record_names(decode=True) + imeth_dict = { + record: imeth for (record, imeth) in zip(records, cbc.imethlist) + } + if text is None: + keylist = records + else: + if not isinstance(text, list): + keylist = [text] + else: + keylist = text + + if kstpkper is None: + kstpkper = cbc.get_kstpkper() + elif isinstance(kstpkper, tuple): + if not isinstance(kstpkper[0], (list, tuple)): + kstpkper = [kstpkper] + else: + pass + + # reset totim based on values read from budget file + times = cbc.get_times() + kstpkpers = cbc.get_kstpkper() + self._totim = {ki: time for (ki, time) in zip(kstpkpers, times)} + + for name in keylist: + d = {} + for i, ki in enumerate(kstpkper): + try: + array = cbc.get_data(kstpkper=ki, text=name, full3D=True) + if len(array) == 0: + continue + + array = np.ma.filled(array, np.nan) + if array.size < self.nnodes: + if array.size < self.ncpl: + raise AssertionError( + "Array size must be equal to " + "either ncpl or nnodes" + ) + + array = np.zeros(self.nnodes) * np.nan + array[: array.size] = np.ravel(array) + + except ValueError: + if imeth_dict[name] == 6: + array = np.full((self.nnodes,), np.nan) + rec = cbc.get_data(kstpkper=ki, text=name)[0] + for [node, q] in zip(rec["node"], rec["q"]): + array[node] = q + else: + continue + + d[ki] = array + + self.__transient_output_data = False + self.add_transient_array(d, name, masked_values) + self.__transient_output_data = True + + def _set_modpath_point_data(self, points, d): + """ + Method to build the vtk point geometry and set arrays for + modpath pathlines + + Parameters + ---------- + points : list + list of (x, y, z) points + d : dict + dictionary of numpy arrays to add to vtk + + """ + from vtk.util import numpy_support + + nverts = len(points) + + self.vtk_pathlines = self.__vtk.vtkUnstructuredGrid() + + vtk_points = self.__vtk.vtkPoints() + for point in points: + vtk_points.InsertNextPoint(point) + + self.vtk_pathlines.SetPoints(vtk_points) + + # create a Vertex instance for each point data add to grid + for i in range(nverts): + vertex = self.__vtk.vtkPolyVertex() + vertex.GetPointIds().SetNumberOfIds(1) + vertex.GetPointIds().SetId(0, i) + + # set data to the pathline grid + self.vtk_pathlines.InsertNextCell( + vertex.GetCellType(), vertex.GetPointIds() + ) + + # process arrays and add arrays to grid. + for name, array in d.items(): + array = np.array(array) + vtk_array = numpy_support.numpy_to_vtk( + num_array=array, array_type=self.__vtk.VTK_FLOAT + ) + vtk_array.SetName(name) + self.vtk_pathlines.GetCellData().AddArray(vtk_array) + + def write(self, f, kper=None): + """ + Method to write a vtk file from the VTK object + + Parameters + ---------- + f : str + vtk file name + kpers : int, list, tuple + stress period or list of stress periods to write to vtk. This + parameter only applies to transient package data. + + """ + grids = [self.vtk_grid, self.vtk_polygons, self.vtk_pathlines] + suffix = ["", "_hfb", "_pathline"] + + extension = ".vtk" + if self.pvd: + self.pvd = Pvd() + extension = ".vtu" + + fpth, _ = os.path.split(f) + if not os.path.exists(fpth): + os.mkdir(fpth) + + if kper is not None: + if isinstance(kper, (int, float)): + kper = [int(kper)] + + for ix, grid in enumerate(grids): + if grid is None: + continue + + if not f.endswith(".vtk") and not f.endswith(".vtu"): + foo = f"{f}{suffix[ix]}{extension}" + else: + foo = f"{f[:-4]}{suffix[ix]}{f[-4:]}" + + if not self.xml: + w = self.__vtk.vtkUnstructuredGridWriter() + if self.binary: + w.SetFileTypeToBinary() + else: + w = self.__vtk.vtkXMLUnstructuredGridWriter() + if not self.binary: + w.SetDataModeToAscii() + + if self.__pathline_transient_data and ix == 2: + stp = 0 + for time, d in self.__pathline_transient_data.items(): + tf = self.__create_transient_vtk_path(foo, stp) + points = self._pathline_points[time] + self._set_modpath_point_data(points, d) + + w.SetInputData(self.vtk_pathlines) + w.SetFileName(tf) + w.Update() + stp += 1 + + else: + w.SetInputData(grid) + + if ( + self.__transient_data or self.__transient_vector + ) and ix == 0: + if self.__transient_data: + cnt = 0 + for per, d in self.__transient_data.items(): + if kper is not None: + if per not in kper: + continue + + if self.__transient_output_data: + tf = self.__create_transient_vtk_path(foo, cnt) + else: + tf = self.__create_transient_vtk_path(foo, per) + self._add_timevalue(per, tf) + for name, array in d.items(): + self.add_array(array, name) + + if per in self.__transient_vector: + d = self.__transient_vector[d] + for name, vector in d.items(): + self.add_vector(vector, name) + + w.SetFileName(tf) + w.Update() + cnt += 1 + else: + cnt = 0 + for per, d in self.__transient_vector.items(): + if kper is not None: + if per not in kper: + continue + + if self.__transient_output_data: + tf = self.__create_transient_vtk_path(foo, cnt) + else: + tf = self.__create_transient_vtk_path(foo, per) + self._add_timevalue(per) + for name, vector in d.items(): + self.add_vector(vector, name) + + w.SetFileName(tf) + w.update() + cnt += 1 + else: + w.SetFileName(foo) + w.Update() + + if not type(self.pvd) == bool: + if not f.endswith(".vtu") or f.endswith(".vtk"): + pvdfile = f"{f}.pvd" + else: + pvdfile = f"{f[:-4]}.pvd" + + self.pvd.write(pvdfile) + + def __create_transient_vtk_path(self, path, kper): + """ + Method to set naming convention for transient vtk file series + + Parameters + ---------- + path : str + vtk file path + kper : int + zero based stress period number + + Returns + ------- + updated vtk file path of format _{:06d}.vtk where + {:06d} represents the six zero padded stress period time + """ + pth = ".".join(path.split(".")[:-1]) + if pth.endswith("_"): + pth = pth[:-1] + extension = path.split(".")[-1] + return f"{pth}_{kper :06d}.{extension}" + + +def export_model( + model, + otfolder, + package_names=None, + nanval=-1e20, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, + kpers=None, +): + """ + DEPRECATED method to export model to vtk + + Parameters + ---------- + + model : flopy model instance + flopy model + otfolder : 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, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + kpers : iterable of int + Stress periods to export. If None (default), all stress periods will be + exported. + """ + warnings.warn("export_model is deprecated, please use Vtk.add_model()") + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + + vtk = Vtk( + model, + vertical_exageration=1, + binary=binary, + smooth=smooth, + point_scalars=point_scalars, + ) + + vtk.add_model(model, package_names) + + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + name = model.name + vtk.write(os.path.join(otfolder, name)) + + +def export_package( + pak_model, + pak_name, + otfolder, + vtkobj=None, + nanval=-1e20, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, + kpers=None, +): + """ + DEPRECATED method to export package to vtk + + Parameters + ---------- + + pak_model : flopy model instance + the model of the package + pak_name : str + the name of the package + otfolder : 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, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + kpers : iterable of int + Stress periods to export. If None (default), all stress periods will be + exported. + """ + warnings.warn("export_package is deprecated, use Vtk.add_package()") + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + + if not vtkobj: + vtk = Vtk( + pak_model, + binary=binary, + smooth=smooth, + point_scalars=point_scalars, + ) + else: + vtk = vtkobj + + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + if isinstance(pak_name, list): + pak_name = pak_name[0] + + p = pak_model.get_package(pak_name) + vtk.add_package(p) + + vtk.write(os.path.join(otfolder, pak_name)) + + +def export_transient( + model, + array, + output_folder, + name, + nanval=-1e20, + array2d=False, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, + kpers=None, +): + """ + DEPRECATED method to export transient arrays and lists 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, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + kpers : iterable of int + Stress periods to export. If None (default), all stress periods will be + exported. + """ + warnings.warn( + "export_transient is deprecated, use Vtk.add_transient_array() or " + "Vtk.add_transient_list()" + ) + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + if array2d: + warnings.warn( + "array2d parameter is deprecated, 2d arrays are " + "handled automatically" + ) + + if not os.path.exists(output_folder): + os.mkdir(output_folder) + + vtk = Vtk(model, binary=binary, smooth=smooth, point_scalars=point_scalars) + + if array.data_type == DataType.transient2d: + if array.array is not None: + if hasattr(array, "transient_2ds"): + vtk.add_transient_array(array.transient_2ds, name) + else: + d = {ix: i for ix, i in enumerate(array.array)} + vtk.add_transient_array(d, name) + + elif array.data_type == DataType.transient3d: + if array.array is None: + vtk.add_transient_array(array.transient_3ds, name) + + elif array.data_type == DataType.transientlist: + vtk.add_transient_list(array) + + else: + raise TypeError(f"type {type(array)} not valid for export_transient") + + +def export_array( + model, + array, + output_folder, + name, + nanval=-1e20, + array2d=False, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, +): + """ + DEPRECATED method to 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 the array is 2d and represents the first layer, default is + False + smooth : bool + if True, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + """ + warnings.warn("export_array is deprecated, please use Vtk.add_array()") + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + if array2d: + warnings.warn( + "array2d parameter is deprecated, 2d arrays are " + "handled automatically" + ) + + if not os.path.exists(output_folder): + os.mkdir(output_folder) + + if array.size < model.modelgrid.nnodes: + if array.size < model.modelgrid.ncpl: + raise AssertionError( + "Array size must be equal to either ncpl or nnodes" + ) + + array = np.zeros(model.modelgrid.nnodes) * np.nan + array[: array.size] = np.ravel(array) + + vtk = Vtk(model, binary=binary, smooth=smooth, point_scalars=point_scalars) + + vtk.add_array(array, name) + vtk.write(os.path.join(output_folder, name)) + + +def export_heads( + model, + hdsfile, + otfolder, + text="head", + precision="auto", + verbose=False, + nanval=-1e20, + kstpkper=None, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, +): + """ + DEPRECATED method + + Exports binary head file to vtk + + Parameters + ---------- + + model : MFModel + the flopy model instance + hdsfile : str, HeadFile object + binary head file path or object + otfolder : str + output folder to write the data + text : string + Name of the text string in the head file. Default is 'head'. + precision : str + Precision of data in the head file: 'auto', 'single' or 'double'. + Default is 'auto'. + verbose : bool + If True, write information to the screen. Default is False. + 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, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + """ + warnings.warn("export_heads is deprecated, use Vtk.add_heads()") + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + + from flopy.utils import HeadFile + + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + if not isinstance(hdsfile, HeadFile): + hds = HeadFile( + hdsfile, text=text, precision=precision, verbose=verbose + ) + else: + hds = hdsfile + + vtk = Vtk(model, binary=binary, smooth=smooth, point_scalars=point_scalars) + + vtk.add_heads(hds, kstpkper) + name = f"{model.name}_{text}" + vtk.write(os.path.join(otfolder, name)) + + +def export_cbc( + model, + cbcfile, + otfolder, + precision="single", + verbose=False, + nanval=-1e20, + kstpkper=None, + text=None, + smooth=False, + point_scalars=False, + vtk_grid_type="auto", + true2d=False, + binary=True, +): + """ + DEPRECATED method + + 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 + Precision of data in the cell by cell file: 'single' or 'double'. + Default is 'single'. + verbose : bool + If True, write information to the screen. Default is False. + 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, will create smooth layer elevations, default is False + point_scalars : bool + if True, will also output array values at cell vertices, default is + False; note this automatically sets smooth to True + vtk_grid_type : str + Specific vtk_grid_type or 'auto' (default). Possible specific values + are 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'. + If 'auto', the grid type is automatically determined. Namely: + * A regular grid (in all three directions) will be saved as an + 'ImageData'. + * A rectilinear (in all three directions), non-regular grid + will be saved as a 'RectilinearGrid'. + * Other grids will be saved as 'UnstructuredGrid'. + true2d : bool + If True, the model is expected to be 2d (1 layer, 1 row or 1 column) + and the data will be exported as true 2d data, default is False. + binary : bool + if True the output file will be binary, default is False + """ + warnings.warn("export_cbc is deprecated, use Vtk.add_cell_budget()") + + if nanval != -1e20: + warnings.warn("nanval is deprecated, setting to np.nan") + if true2d: + warnings.warn("true2d is no longer supported, setting to False") + if vtk_grid_type != "auto": + warnings.warn("vtk_grid_type is deprecated, setting to binary") + + from flopy.utils import CellBudgetFile + + if not os.path.exists(otfolder): + os.mkdir(otfolder) + + if not isinstance(cbcfile, CellBudgetFile): + cbc = CellBudgetFile(cbcfile, precision=precision, verbose=verbose) + else: + cbc = cbcfile + + vtk = Vtk(model, binary=binary, smooth=smooth, point_scalars=point_scalars) + + vtk.add_cell_budget(cbc, text, kstpkper) + fname = f"{model.name}_CBC" + vtk.write(os.path.join(otfolder, fname)) From 329f08a986aa7b9392e6cbc684377c8cc338196e Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 28 Sep 2021 14:02:45 -0700 Subject: [PATCH 3/5] minor update to deprecated methods --- flopy/export/vtk.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 1b7392c5e..0c5ca8de9 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -1493,7 +1493,7 @@ def export_model( os.mkdir(otfolder) name = model.name - vtk.write(os.path.join(otfolder, name)) + vtk.write(os.path.join(otfolder, name), kper=kpers) def export_package( @@ -1577,7 +1577,7 @@ def export_package( p = pak_model.get_package(pak_name) vtk.add_package(p) - vtk.write(os.path.join(otfolder, pak_name)) + vtk.write(os.path.join(otfolder, pak_name), kper=kpers) def export_transient( @@ -1675,6 +1675,8 @@ def export_transient( else: raise TypeError(f"type {type(array)} not valid for export_transient") + vtk.write(os.path.join(output_folder, name), kper=kpers) + def export_array( model, From 1d72214b311c88a0cf01e5d3c1c14973b15ea842 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 28 Sep 2021 15:35:46 -0700 Subject: [PATCH 4/5] updates to test_vtk_vector() and test_vtk_pathline() --- autotest/t050_test.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 9d9501063..5064f3fd1 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -504,7 +504,6 @@ def test_vtk_vertex(): outfile = os.path.join("temp", "t050", "vtk_disv", "disv.vtk") sim = flopy.mf6.MFSimulation.load(sim_ws=workspace) gwf = sim.get_model("gwf_1") - sim.run_simulation() vtkobj = Vtk(model=gwf, binary=True, smooth=False) vtkobj.add_model(gwf) @@ -597,7 +596,7 @@ def test_vtk_pathline(): if not np.abs(np.max(totim) - maxtime) < 100: raise AssertionError("time values are incorrect for modpath VTK") - if not np.max(pid) == len(plines): + if not len(np.unique(pid)) == len(plines): raise AssertionError( "number of particles are incorrect for modpath VTK" ) From 73e5b63052844ea7292e826e53b2f07c8594c885 Mon Sep 17 00:00:00 2001 From: Joshua Larsen Date: Tue, 28 Sep 2021 17:14:19 -0700 Subject: [PATCH 5/5] Minor codacy changes --- flopy/export/vtk.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py index 0c5ca8de9..12f44f07e 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -1,8 +1,13 @@ -import numpy as np +""" +The vtk module provides functionality for exporting model inputs and +outputs to VTK. +""" + import os +import warnings +import numpy as np from flopy.datbase import DataType, DataInterface from flopy.utils import Util3d -import warnings warnings.simplefilter("always", DeprecationWarning) @@ -13,6 +18,7 @@ class Pvd: """ Simple class to build a Paraview Data File (PVD) + """ def __init__(self): @@ -331,7 +337,7 @@ def _create_point_scalar_graph(self): xvert = np.ones((max_shared, num_points)) * np.nan yvert = np.ones((max_shared, num_points)) * np.nan for k, d in graph.items(): - for ix, pt in enumerate(d["vtk_points"]): + for _, pt in enumerate(d["vtk_points"]): for ixx, value in enumerate(d["idx"]): numpy_graph[ixx, pt] = value xvert[ixx, pt] = d["xv"][ixx]