diff --git a/autotest/t050_test.py b/autotest/t050_test.py
index 02eb9fc684..bbcd77428b 100644
--- a/autotest/t050_test.py
+++ b/autotest/t050_test.py
@@ -1,116 +1,207 @@
import shutil
import os
-import numpy as np
import flopy
from flopy.export 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 now.
+
# create output directory
cpth = os.path.join('temp', 't050')
if os.path.isdir(cpth):
shutil.rmtree(cpth)
os.makedirs(cpth)
-# binary output directory
-binot = os.path.join(cpth, 'bin')
-if os.path.isdir(binot):
- shutil.rmtree(binot)
-os.makedirs(binot)
+def count_lines_in_file(filepath):
+ f = open(filepath, 'r')
+ n = len(f.readlines())
+ return n
+def count_lines_in_file_bin(filepath):
+ f = open(filepath, 'rb')
+ n = len(f.readlines())
+ return n
def test_vtk_export_array2d():
- """Export 2d array"""
+ # test mf 2005 freyberg
mpath = os.path.join('..', 'examples', 'data',
'freyberg_multilayer_transient')
namfile = 'freyberg.nam'
- m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False)
- m.dis.top.export(os.path.join(cpth, 'array_2d_test'), fmt='vtk')
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False,
+ load_only=['dis', 'bas6'])
+ 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==352026)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==2846)
+
# with smoothing
- m.dis.top.export(os.path.join(cpth, 'array_2d_test'), fmt='vtk',
- name='top_smooth', smooth=True)
+ 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==351829)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==2846)
+ return
def test_vtk_export_array3d():
- """Vtk export 3d array"""
+ # test mf 2005 freyberg
mpath = os.path.join('..', 'examples', 'data',
'freyberg_multilayer_transient')
namfile = 'freyberg.nam'
- m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False)
- m.upw.hk.export(os.path.join(cpth, 'array_3d_test'), fmt='vtk')
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False,
+ load_only=['dis', 'bas6', 'upw'])
+ 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==992576)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==8486)
+
# with point scalars
- m.upw.hk.export(os.path.join(cpth, 'array_3d_test'), fmt='vtk',
- name='hk_points', point_scalars=True)
- # binary test
- m.upw.hk.export(os.path.join(binot, 'array_3d_test'), fmt='vtk',
- name='hk_points', point_scalars=True, binary=True)
+ m.upw.hk.export(output_dir, fmt='vtk', name='hk_points',
+ point_scalars=True)
+ filetocheck = os.path.join(output_dir, 'hk_points.vtu')
+ # totalbytes1 = os.path.getsize(filetocheck)
+ # assert(totalbytes1==1321292)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==10605)
+
+ # with point scalars and binary
+ m.upw.hk.export(output_dir, 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==637861)
+ nlines2 = count_lines_in_file_bin(filetocheck)
+ assert(nlines2==1848)
return
-
def test_vtk_transient_array_2d():
- """VTK export transient 2d array"""
+ # test mf 2005 freyberg
mpath = os.path.join('..', 'examples', 'data',
'freyberg_multilayer_transient')
namfile = 'freyberg.nam'
- m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False)
- m.rch.rech.export(os.path.join(cpth, 'transient_2d_test'), fmt='vtk')
-
- # binary test
- m.rch.rech.export(os.path.join(binot, 'transient_2d_test'), fmt='vtk',
- binary=True)
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False,
+ load_only=['dis', 'bas6', 'rch'])
+ output_dir = os.path.join(cpth, 'transient_2d_test')
+ output_dir_bin = os.path.join(cpth, 'transient_2d_test_bin')
+
+ # export and check
+ r = m.rch.rech
+ m.rch.rech.export(output_dir, fmt='vtk')
+ filetocheck = os.path.join(output_dir, 'rech_01.vtu')
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==355324)
+ 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==354622)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==2851)
+
+ # with binary
+ m.rch.rech.export(output_dir_bin, fmt='vtk', binary=True)
+ filetocheck = os.path.join(output_dir_bin, 'rech_01.vtu')
+ # totalbytes2 = os.path.getsize(filetocheck)
+ # assert(totalbytes2==168339)
+ nlines2 = count_lines_in_file_bin(filetocheck)
+ assert(nlines2==762)
+ filetocheck = os.path.join(output_dir_bin, 'rech_01097.vtu')
+ # totalbytes3 = os.path.getsize(filetocheck)
+ # assert(totalbytes3==168339)
+ nlines3 = count_lines_in_file_bin(filetocheck)
+ assert(nlines3==762)
return
-
def test_vtk_export_packages():
- """testing vtk package export"""
+ # test mf 2005 freyberg
mpath = os.path.join('..', 'examples', 'data',
'freyberg_multilayer_transient')
namfile = 'freyberg.nam'
- m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False)
- # test dis export
- m.dis.export(os.path.join(cpth, 'DIS'), fmt='vtk')
- # upw with point scalar output
- m.upw.export(os.path.join(cpth, 'UPW'), fmt='vtk', point_scalars=True)
- # bas with smoothing on
- m.bas6.export(os.path.join(cpth, 'BAS'), fmt='vtk', smooth=True)
- # transient package drain
- m.drn.export(os.path.join(cpth, 'DRN'), fmt='vtk')
- # binary test
- m.dis.export(os.path.join(binot, 'DIS'), fmt='vtk', binary=True)
- # upw with point scalar output
- m.upw.export(os.path.join(binot, 'UPW'), fmt='vtk', point_scalars=True,
- binary=True)
-
- return
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpath, verbose=False,
+ load_only=['dis', 'bas6', 'upw', 'DRN'])
+
+ # 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')
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==1010527)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==8496)
+ # 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==2485295)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==21215)
-# add mf2005 model exports
-def test_export_mf2005_vtk():
- """test vtk model export mf2005"""
- pth = os.path.join('..', 'examples', 'data', 'mf2005_test')
- namfiles = [namfile for namfile in os.listdir(pth) if
- namfile.endswith('.nam')]
- skip = ['bcf2ss.nam']
- for namfile in namfiles:
- if namfile in skip:
- continue
- print('testing namefile', namfile)
- m = flopy.modflow.Modflow.load(namfile, model_ws=pth, verbose=False)
- m.export(os.path.join(cpth, m.name), fmt='vtk')
+ # 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==1002054)
+ nlines2 = count_lines_in_file(filetocheck)
+ assert(nlines2==8491)
- # binary test
- m.export(os.path.join(binot, m.name), fmt='vtk', binary=True)
+ # transient package drain
+ output_dir = os.path.join(cpth, 'DRN')
+ m.drn.export(output_dir, fmt='vtk')
+ filetocheck = os.path.join(output_dir, 'DRN_01.vtu')
+ # totalbytes3 = os.path.getsize(filetocheck)
+ # assert(totalbytes3==20702)
+ 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==20702)
+ nlines4 = count_lines_in_file(filetocheck)
+ assert(nlines4==191)
+
+ # 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==536436)
+ nlines5 = count_lines_in_file_bin(filetocheck)
+ assert(nlines5==1545)
+
+ # 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==1400561)
+ nlines6 = count_lines_in_file_bin(filetocheck)
+ assert(nlines6==1868)
return
-
def test_vtk_mf6():
+ # test mf6
mf6expth = os.path.join('..', 'examples', 'data', 'mf6')
- # test vtk mf6 export
- mf6sims = ['test045_lake1ss_table',
- 'test036_twrihfb', 'test045_lake2tr', 'test006_2models_mvr']
- # mf6sims = ['test005_advgw_tidal']
- # mf6sims = ['test036_twrihfb']
+ mf6sims = ['test045_lake1ss_table', 'test036_twrihfb', 'test045_lake2tr',
+ 'test006_2models_mvr']
for simnm in mf6sims:
print(simnm)
@@ -124,91 +215,234 @@ def test_vtk_mf6():
m = loaded_sim.get_model(mname)
m.export(os.path.join(cpth, m.name), fmt='vtk')
+ # check one
+ filetocheck = os.path.join(cpth, 'twrihfb2015', 'npf.vtr')
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==21609)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==76)
+
return
def test_vtk_binary_head_export():
+ # test mf 2005 freyberg
+ mpth = os.path.join('..', 'examples', 'data',
+ 'freyberg_multilayer_transient')
+ namfile = 'freyberg.nam'
+ hdsfile = os.path.join(mpth, 'freyberg.hds')
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False,
+ load_only=['dis', 'bas6'])
+ filenametocheck = 'freyberg_Heads_KPER455_KSTP1.vtu'
- """test vet export of heads"""
-
- freyberg_pth = os.path.join('..', 'examples', 'data',
- 'freyberg_multilayer_transient')
-
- hdsfile = os.path.join(freyberg_pth, 'freyberg.hds')
-
- m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_pth,
- verbose=False)
+ # 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)])
- # test with points
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==993755)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==8486)
+
+ # 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)
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes1 = os.path.getsize(filetocheck)
+ # assert(totalbytes1==1332346)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==10605)
- # test vtk export heads with smoothing and no point scalars
+ # 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)],
- point_scalars=False, smooth=True, nanval=-999.99)
-
- # test binary output
+ smooth=True, nanval=-999.99)
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes2 = os.path.getsize(filetocheck)
+ # assert(totalbytes2==993551)
+ 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)],
- point_scalars=False, smooth=True, binary=True,
- nanval=-999.99)
+ smooth=True, binary=True, nanval=-999.99)
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes3 = os.path.getsize(filetocheck)
+ # assert(totalbytes3==502313)
+ nlines3 = count_lines_in_file_bin(filetocheck)
+ assert(nlines3==1538)
+ # 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)
+ point_scalars=False, smooth=True, binary=True,
+ nanval=-999.99)
+ filetocheck = os.path.join(otfolder, 'freyberg_Heads_KPER1_KSTP1.vtu')
+ # totalbytes4 = os.path.getsize(filetocheck)
+ # assert(totalbytes4==502313)
+ nlines4 = count_lines_in_file_bin(filetocheck)
+ assert(nlines4==1528)
return
-
def test_vtk_cbc():
# test mf 2005 freyberg
- freyberg_cbc = os.path.join('..', 'examples', 'data',
- 'freyberg_multilayer_transient',
- 'freyberg.cbc')
-
- freyberg_mpth = os.path.join('..', 'examples', 'data',
- 'freyberg_multilayer_transient')
-
- m = flopy.modflow.Modflow.load('freyberg.nam', model_ws=freyberg_mpth,
- verbose=False)
+ mpth = os.path.join('..', 'examples', 'data',
+ 'freyberg_multilayer_transient')
+ namfile = 'freyberg.nam'
+ cbcfile = os.path.join(mpth, 'freyberg.cbc')
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False)
+ filenametocheck = 'freyberg_CBC_KPER1_KSTP1.vtu'
- vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST'),
+ # 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)
-
- vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth, 'freyberg_CBCTEST_bin'),
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==2496132)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==19093)
+
+ # 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)
-
- vtk.export_cbc(m, freyberg_cbc, os.path.join(cpth,
- 'freyberg_CBCTEST_bin2'),
+ filetocheck = os.path.join(otfolder, filenametocheck)
+ # totalbytes1 = os.path.getsize(filetocheck)
+ # assert(totalbytes1==1248118)
+ nlines1 = count_lines_in_file_bin(filetocheck)
+ assert(nlines1==2609)
+
+ # 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==10262)
+ nlines2 = count_lines_in_file_bin(filetocheck)
+ assert(nlines2==61)
+
+ return
+
+def test_vtk_vti():
+ # test mf 2005 ibs2k
+ mpth = os.path.join('..', 'examples', 'data', 'mf2005_test')
+ namfile = 'ibs2k.nam'
+ m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False)
+ output_dir = os.path.join(cpth, m.name)
+ filenametocheck = 'DIS.vti'
+
+ # export and check
+ m.export(output_dir, fmt='vtk')
+ filetocheck = os.path.join(output_dir, filenametocheck)
+ # totalbytes = os.path.getsize(filetocheck)
+ # assert(totalbytes==6322)
+ nlines = count_lines_in_file(filetocheck)
+ assert(nlines==21)
+
+ # with point scalar
+ m.export(output_dir + '_points', fmt='vtk', point_scalars=True)
+ filetocheck = os.path.join(output_dir + '_points', filenametocheck)
+ # totalbytes1 = os.path.getsize(filetocheck)
+ # assert(totalbytes1==16382)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==38)
+
+ # with binary
+ m.export(output_dir + '_bin', fmt='vtk', binary=True)
+ filetocheck = os.path.join(output_dir + '_bin', filenametocheck)
+ # totalbytes2 = os.path.getsize(filetocheck)
+ # assert(totalbytes2==6537)
+ nlines2 = count_lines_in_file_bin(filetocheck)
+ assert(nlines2==18)
+
+ # force .vtr
+ filenametocheck = 'DIS.vtr'
+ m.export(output_dir, fmt='vtk', vtk_grid_type='RectilinearGrid')
+ filetocheck = os.path.join(output_dir, filenametocheck)
+ # totalbytes3 = os.path.getsize(filetocheck)
+ # assert(totalbytes3==7146)
+ nlines3 = count_lines_in_file(filetocheck)
+ assert(nlines3==56)
+
+ # force .vtu
+ filenametocheck = 'DIS.vtu'
+ m.export(output_dir, fmt='vtk', vtk_grid_type='UnstructuredGrid')
+ filetocheck = os.path.join(output_dir, filenametocheck)
+ # totalbytes4 = os.path.getsize(filetocheck)
+ # assert(totalbytes4==67065)
+ nlines4 = count_lines_in_file(filetocheck)
+ assert(nlines4==993)
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(output_dir + '_points', fmt='vtk', point_scalars=True)
+ filetocheck = os.path.join(output_dir + '_points', filenametocheck)
+ # totalbytes1 = os.path.getsize(filetocheck)
+ # assert(totalbytes1==182472)
+ nlines1 = count_lines_in_file(filetocheck)
+ assert(nlines1==121)
+
+ # with binary
+ m.export(output_dir + '_bin', fmt='vtk', binary=True)
+ filetocheck = os.path.join(output_dir + '_bin', filenametocheck)
+ # totalbytes2 = os.path.getsize(filetocheck)
+ # assert(totalbytes2==47778)
+ nlines2 = count_lines_in_file_bin(filetocheck)
+ assert(nlines2==28)
+
+ # 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
if __name__ == '__main__':
test_vtk_export_array2d()
test_vtk_export_array3d()
test_vtk_transient_array_2d()
test_vtk_export_packages()
- test_export_mf2005_vtk()
test_vtk_mf6()
test_vtk_binary_head_export()
test_vtk_cbc()
+ test_vtk_vti()
+ test_vtk_vtr()
diff --git a/flopy/discretization/grid.py b/flopy/discretization/grid.py
index dcab127bed..9853a95738 100644
--- a/flopy/discretization/grid.py
+++ b/flopy/discretization/grid.py
@@ -284,6 +284,12 @@ def extent(self):
raise NotImplementedError(
'must define extent in child class')
+ @property
+ def xyzextent(self):
+ return (np.min(self.xyzvertices[0]), np.max(self.xyzvertices[0]),
+ np.min(self.xyzvertices[1]), np.max(self.xyzvertices[1]),
+ np.min(self.xyzvertices[2]), np.max(self.xyzvertices[2]))
+
@property
def grid_lines(self):
raise NotImplementedError(
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index 23222da45d..f8a77b45e1 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -177,6 +177,23 @@ def xyedges(self):
else:
return self._cache_dict[cache_index].data_nocopy
+ @property
+ def zedges(self):
+ """
+ Return zedges for (column, row)==(0, 0).
+ """
+ cache_index = 'zedges'
+ if cache_index not in self._cache_dict or \
+ self._cache_dict[cache_index].out_of_date:
+ zedge = np.concatenate((np.array([self.top[0, 0]]),
+ self.botm[:, 0, 0]))
+ self._cache_dict[cache_index] = \
+ CachedData(zedge)
+ if self._copy_cache:
+ return self._cache_dict[cache_index].data
+ else:
+ return self._cache_dict[cache_index].data_nocopy
+
@property
def xyzcellcenters(self):
"""
@@ -421,6 +438,57 @@ def write_shapefile(self, filename='grid.shp', epsg=None, prj=None):
write_grid_shapefile(filename, self, array_dict={}, nan_val=-1.0e9,
epsg=epsg, prj=prj)
+ def is_regular(self):
+ """
+ Test whether the grid spacing is regular or not (including in the
+ vertical direction).
+ """
+ # Relative tolerance to use in test
+ rel_tol = 1.e-5
+
+ # Regularity test in x direction
+ rel_diff_x = (self.delr - self.delr[0]) / self.delr[0]
+ is_regular_x = np.count_nonzero(rel_diff_x > rel_tol) == 0
+
+ # Regularity test in y direction
+ rel_diff_y = (self.delc - self.delc[0]) / self.delc[0]
+ is_regular_y = np.count_nonzero(rel_diff_y > rel_tol) == 0
+
+ # Regularity test in z direction
+ thickness = (self.top[0, 0] - self.botm[0, 0, 0])
+ rel_diff_z1 = (self.top - self.botm[0, :, :] - thickness) / thickness
+ failed = np.abs(rel_diff_z1) > rel_tol
+ is_regular_z = np.count_nonzero(failed) == 0
+ for k in range(self.nlay - 1):
+ rel_diff_zk = (self.botm[k, :, :] - self.botm[k + 1, :, :] -
+ thickness) / thickness
+ failed = np.abs(rel_diff_zk) > rel_tol
+ is_regular_z = is_regular_z and np.count_nonzero(failed) == 0
+
+ return is_regular_x and is_regular_y and is_regular_z
+
+ def is_rectilinear(self):
+ """
+ Test whether the grid is rectilinear (it is always so in the x and
+ y directions, but not necessarily in the z direction).
+ """
+ # Relative tolerance to use in test
+ rel_tol = 1.e-5
+
+ # Rectilinearity test in z direction
+ thickness = (self.top[0, 0] - self.botm[0, 0, 0])
+ rel_diff_z1 = (self.top - self.botm[0, :, :] - thickness) / thickness
+ failed = np.abs(rel_diff_z1) > rel_tol
+ is_rectilinear_z = np.count_nonzero(failed) == 0
+ for k in range(self.nlay - 1):
+ thickness_k = (self.botm[k, 0, 0] - self.botm[k + 1, 0, 0])
+ rel_diff_zk = (self.botm[k, :, :] - self.botm[k + 1, :, :] -
+ thickness_k) / thickness_k
+ failed = np.abs(rel_diff_zk) > rel_tol
+ is_rectilinear_z = is_rectilinear_z and \
+ np.count_nonzero(failed) == 0
+
+ return is_rectilinear_z
if __name__ == "__main__":
import matplotlib.pyplot as plt
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index 4b4a9eba89..1e23712374 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -512,32 +512,17 @@ def model_export(f, ml, fmt=None, **kwargs):
or dictionary of ....
ml : flopy.modflow.mbase.ModelInterface object
flopy model object
- fmt: str
- output format flag. set to 'vtk' to export to vtk .vtu file
-
+ fmt : str
+ output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
- modelgrid: flopy.discretization.Grid
+ modelgrid: flopy.discretization.Grid
user supplied modelgrid object which will supercede the built
in modelgrid object
epsg : int
epsg projection code
prj : str
prj file name
-
- smooth: bool
- For vtk, if set to True will output smooth surface
-
- point_scalars: bool
- For vtk, if set to True will create point scalar values along
- with cell values.
-
- binary: bool
- For vtk, if set to Ture will output binary .vtu files. Default
- is False which exports standard text xml .vtu files.
-
- nanval: For vtk, no data value
-
-
+ if fmt is set to 'vtk', parameters of vtk.export_model
"""
assert isinstance(ml, ModelInterface)
@@ -567,13 +552,14 @@ def model_export(f, ml, fmt=None, **kwargs):
elif fmt == 'vtk':
# call vtk model export
+ nanval = kwargs.get('nanval', -1e20)
smooth = kwargs.get('smooth', False)
point_scalars = kwargs.get('point_scalars', False)
+ vtk_grid_type = kwargs.get('vtk_grid_type', 'auto')
binary = kwargs.get('binary', False)
- nanval = kwargs.get('nanval', -1e20)
- vtk.export_model(ml, f, package_names=package_names, smooth=smooth,
- point_scalars=point_scalars, binary=binary,
- nanval=nanval)
+ vtk.export_model(ml, f, package_names=package_names, nanval=nanval,
+ smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
else:
raise NotImplementedError("unrecognized export argument:{0}".format(f))
@@ -591,8 +577,8 @@ def package_export(f, pak, fmt=None, **kwargs):
output file name (ends in .shp for shapefile or .nc for netcdf)
pak : flopy.pakbase.Package object
package to export
- fmt: str
- output format flag. set to 'vtk' to export to vtk .vtu file
+ fmt : str
+ output format flag. 'vtk' will export to vtk
** kwargs : keword arguments
modelgrid: flopy.discretization.Grid
user supplied modelgrid object which will supercede the built
@@ -601,19 +587,7 @@ def package_export(f, pak, fmt=None, **kwargs):
epsg projection code
prj : str
prj file name
-
- smooth: bool
- For vtk, if set to True will output smooth surface
-
- point_scalars: bool
- For vtk, if set to True will create point scalar values along
- with cell values.
-
- binary: bool
- For vtk, if set to Ture will output binary .vtu files. Default
- is False which exports standard text xml .vtu files.
-
- nanval: For vtk, no data value
+ if fmt is set to 'vtk', parameters of vtk.export_package
Returns
-------
@@ -656,14 +630,14 @@ def package_export(f, pak, fmt=None, **kwargs):
elif fmt == 'vtk':
# call vtk array export to folder
+ nanval = kwargs.get('nanval', -1e20)
smooth = kwargs.get('smooth', False)
point_scalars = kwargs.get('point_scalars', False)
+ vtk_grid_type = kwargs.get('vtk_grid_type', 'auto')
binary = kwargs.get('binary', False)
- nanval = kwargs.get('nanval', -1e20)
- vtk.export_package(pak.parent, pak.name, f, smooth=smooth,
- point_scalars=point_scalars, binary=binary,
- nanval=nanval)
-
+ vtk.export_package(pak.parent, pak.name, f, nanval=nanval,
+ smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
else:
raise NotImplementedError("unrecognized export argument:{0}".format(f))
@@ -874,18 +848,14 @@ def transient2d_export(f, t2d, fmt=None, **kwargs):
f : str
filename or existing export instance type (NetCdf only for now)
u2d : Transient2d instance
- fmt: set to 'vtk' to export a .vtu file
+ fmt : str
+ output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
min_valid : minimum valid value
max_valid : maximum valid value
modelgrid : flopy.discretization.Grid
model grid instance which will supercede the flopy.model.modelgrid
- smooth: for vtk to output a smooth represenation of the model
- point_scalars: for vtk to output point value scalars as well as cell
- name: for vtk to set a specific name for array and output file
- binary: bool
- For vtk, if set to True will output binary .vtu files. Default
- is False which exports standard text xml .vtu files.
+ if fmt is set to 'vtk', parameters of vtk.export_transient
"""
@@ -986,13 +956,16 @@ def transient2d_export(f, t2d, fmt=None, **kwargs):
return f
elif fmt == 'vtk':
- smooth = kwargs.get('smooth', False)
- point_scalars = kwargs.get('point_scalars', False)
name = kwargs.get('name', t2d.name)
nanval = kwargs.get('nanval', -1e20)
- vtk.export_transient(t2d.model, t2d.array, f, name, smooth=smooth,
- point_scalars=point_scalars, array2d=True,
- nanval=nanval)
+ smooth = kwargs.get('smooth', False)
+ point_scalars = kwargs.get('point_scalars', False)
+ vtk_grid_type = kwargs.get('vtk_grid_type', 'auto')
+ binary = kwargs.get('binary', False)
+ vtk.export_transient(t2d.model, t2d.array, f, name, nanval=nanval,
+ smooth=smooth, point_scalars=point_scalars,
+ array2d=True, vtk_grid_type=vtk_grid_type,
+ binary=binary)
else:
raise NotImplementedError("unrecognized export argument:{0}".format(f))
@@ -1006,18 +979,14 @@ def array3d_export(f, u3d, fmt=None, **kwargs):
f : str
filename or existing export instance type (NetCdf only for now)
u2d : Util3d instance
- fmt: set to 'vtk' to export a .vtu file
+ fmt : str
+ output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
min_valid : minimum valid value
max_valid : maximum valid value
modelgrid : flopy.discretization.Grid
model grid instance which will supercede the flopy.model.modelgrid
- smooth: for vtk to output a smooth represenation of the model
- point_scalars: for vtk to output point value scalars as well as cell
- name: for vtk to set a specific name for array and output file
- binary: bool
- For vtk, if set to Ture will output binary .vtu files. Default
- is False which exports standard text xml .vtu files.
+ if fmt is set to 'vtk', parameters of vtk.export_array
"""
@@ -1137,17 +1106,18 @@ 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)
smooth = kwargs.get('smooth', False)
point_scalars = kwargs.get('point_scalars', False)
- name = kwargs.get('name', u3d.name)
+ vtk_grid_type = kwargs.get('vtk_grid_type', 'auto')
binary = kwargs.get('binary', False)
- nanval = kwargs.get('nanval', -1e20)
if isinstance(name, list) or isinstance(name, tuple):
name = name[0]
- vtk.export_array(u3d.model, u3d.array, f, name, smooth=smooth,
- point_scalars=point_scalars, binary=binary,
- nanval=nanval)
+ vtk.export_array(u3d.model, u3d.array, f, name, nanval=nanval,
+ smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
else:
raise NotImplementedError("unrecognized export argument:{0}".format(f))
@@ -1162,18 +1132,14 @@ def array2d_export(f, u2d, fmt=None, **kwargs):
f : str
filename or existing export instance type (NetCdf only for now)
u2d : Util2d instance
- fmt: set to 'vtk' to export a .vtu file
+ fmt : str
+ output format flag. 'vtk' will export to vtk
**kwargs : keyword arguments
min_valid : minimum valid value
max_valid : maximum valid value
modelgrid : flopy.discretization.Grid
model grid instance which will supercede the flopy.model.modelgrid
- smooth: for vtk to output a smooth represenation of the model
- point_scalars: for vtk to output point value scalars as well as cell
- name: for vtk to set a specific name for array and output file
- binary: bool
- For vtk, if set to Ture will output binary .vtu files. Default
- is False which exports standard text xml .vtu files.
+ if fmt is set to 'vtk', parameters of vtk.export_array
"""
assert isinstance(u2d, DataInterface), "util2d_helper only helps " \
@@ -1269,14 +1235,16 @@ 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)
smooth = kwargs.get('smooth', False)
point_scalars = kwargs.get('point_scalars', False)
- name = kwargs.get('name', u2d.name)
+ vtk_grid_type = kwargs.get('vtk_grid_type', 'auto')
binary = kwargs.get('binary', False)
- nanval = kwargs.get('nanval', -1e20)
- vtk.export_array(u2d.model, u2d.array, f, name, smooth=smooth,
- point_scalars=point_scalars, array2d=True,
- binary=binary, nanval=nanval)
+ vtk.export_array(u2d.model, u2d.array, f, name, nanval=nanval,
+ smooth=smooth, point_scalars=point_scalars,
+ array2d=True, vtk_grid_type=vtk_grid_type,
+ binary=binary)
else:
raise NotImplementedError("unrecognized export argument:{0}".format(f))
diff --git a/flopy/export/vtk.py b/flopy/export/vtk.py
index dee14c1327..1fed49a8a7 100644
--- a/flopy/export/vtk.py
+++ b/flopy/export/vtk.py
@@ -11,7 +11,17 @@
# Module for exporting vtk from flopy
-# BINARY *********************************************
+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',
@@ -24,127 +34,293 @@
'float64': 'd'}
-class BinaryXml:
+class XmlWriterInterface:
"""
- Helps write binary vtk files
+ Helps writing vtk files.
Parameters
----------
file_path : str
output file path
-
"""
def __init__(self, file_path):
- self.stream = open(file_path, "wb")
+ # class attributes
self.open_tag = False
self.current = []
- self.stream.write(b'')
- if sys.byteorder == "little":
- self.byte_order = '<'
- else:
- self.byte_order = '>'
+ self.indent_level = 0
+ self.indent_char=' '
- def write_size(self, block_size):
- # size is a 64 bit unsigned integer
- byte_order = self.byte_order + 'Q'
- block_size = struct.pack(byte_order, block_size)
- self.stream.write(block_size)
+ # open file and initialize
+ self.f = self._open_file(file_path)
+ self.write_string('')
- def write_array(self, data):
- assert (data.flags['C_CONTIGUOUS'] or data.flags['F_CONTIGUOUS'])
+ # open VTKFile element
+ self.open_element('VTKFile').add_attributes(version='0.1')
- # ravel in fortran order
- dd = np.ravel(data, order='F')
-
- data_format = self.byte_order + str(data.size) + np_to_struct[
- data.dtype.name]
- binary_data = struct.pack(data_format, *dd)
- self.stream.write(binary_data)
-
- def write_coord_arrays(self, x, y, z):
- # check that arrays are the same shape and data type
- assert (x.size == y.size == z.size)
- assert (x.dtype.itemsize == y.dtype.itemsize == z.dtype.itemsize)
-
- # check if arrays are contiguous
- assert (x.flags['C_CONTIGUOUS'] or x.flags['F_CONTIGUOUS'])
- assert (y.flags['C_CONTIGUOUS'] or y.flags['F_CONTIGUOUS'])
- assert (z.flags['C_CONTIGUOUS'] or z.flags['F_CONTIGUOUS'])
-
- data_format = self.byte_order + str(1) + \
- np_to_struct[x.dtype.name]
-
- xrav = np.ravel(x, order='F')
- yrav = np.ravel(y, order='F')
- zrav = np.ravel(z, order='F')
+ def _open_file(self, file_path):
+ """
+ Open the file for writing.
- for idx in range(x.size):
- bx = struct.pack(data_format, xrav[idx])
- by = struct.pack(data_format, yrav[idx])
- bz = struct.pack(data_format, zrav[idx])
- self.stream.write(bx)
- self.stream.write(by)
- self.stream.write(bz)
+ Return
+ ------
+ File object.
+ """
+ raise NotImplementedError('must define _open_file in child class')
- def close(self):
- assert (not self.open_tag)
- self.stream.close()
+ def write_string(self, string):
+ """
+ Write a string into the file.
+ """
+ raise NotImplementedError('must define write_string in child class')
def open_element(self, tag):
if self.open_tag:
- self.stream.write(b">")
- tag_string = "\n<%s" % tag
- self.stream.write(str.encode(tag_string))
+ self.write_string(">")
+ indent = self.indent_level * self.indent_char
+ self.indent_level += 1
+ tag_string = "\n" + indent + "<%s" % 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.stream.write(b">")
+ self.write_string(">")
self.open_tag = False
- string = "\n" % tag
- self.stream.write(str.encode(string))
+ indent = self.indent_level * self.indent_char
+ tag_string = "\n" + indent + "" % tag
+ self.write_string(tag_string)
else:
- self.stream.write(b"/>")
+ self.write_string("/>")
self.open_tag = False
self.current.pop()
return self
- def add_text(self, text):
- if self.open_tag:
- self.stream.write(b">\n")
- self.open_tag = False
- self.stream.write(str.encode(text))
- return self
-
def add_attributes(self, **kwargs):
assert self.open_tag
for key in kwargs:
st = ' %s="%s"' % (key, kwargs[key])
- self.stream.write(str.encode(st))
+ 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
-# END BINARY *********************************************
+ def write_array(self, array, actwcells=None, **kwargs):
+ """
+ Writes an array to the output vtk 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()
-def start_tag(f, tag, indent_level, indent_char=' '):
- # starts xml tag
- s = indent_level * indent_char + tag
- indent_level += 1
- f.write(s + '\n')
- return indent_level
+class XmlWriterAscii(XmlWriterInterface):
+ """
+ Helps writing ascii vtk files.
+ Parameters
+ ----------
-def end_tag(f, tag, indent_level, indent_char=' '):
- # ends xml tag
- indent_level -= 1
- s = indent_level * indent_char + tag
- f.write(s + '\n')
- return indent_level
+ file_path : str
+ output file path
+ """
+ def __init__(self, file_path):
+ super(XmlWriterAscii, self).__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 into the file.
+ """
+ self.f.write(string)
+
+ def write_array(self, array, actwcells=None, **kwargs):
+ """
+ Writes an array to the output vtk 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 -1.e9 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)] = -1.e9
+ s = ' '.join(['{}'.format(val) for val in array_lay_flat])
+ self.write_line(s)
+
+ # 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(XmlWriterBinary, self).__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 into the file.
+ """
+ self.f.write(str.encode(string))
+
+ def write_array(self, array, actwcells=None, **kwargs):
+ """
+ Writes an array to the output vtk 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 = 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(XmlWriterBinary, self).final()
class _Array(object):
@@ -173,24 +349,34 @@ class Vtk(object):
model : MFModel
flopy model instance
verbose : bool
- if True, stdout is verbose
+ If True, stdout is verbose
nanval : float
no data value, default is -1e20
smooth : bool
- If True will create smooth output surface
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- if True will output point scalar values, this will set smooth to True.
+ 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'. Possible specific values are
+ 'ImageData', 'RectilinearGrid', and 'UnstructuredGrid'.
+ If 'auto', the grid type is automatically determined. Namely:
+ * A regular grid (i.e., with cubic cells) will be saved as an
+ 'ImageData'.
+ * A rectilinear, non-regular grid will be saved as a
+ 'RectilinearGrid'.
+ * Other grids will be saved as 'UnstructuredGrid'.
+ 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=-1e+20, smooth=False,
- point_scalars=False):
+ point_scalars=False, vtk_grid_type='auto', binary=False):
if point_scalars:
smooth = True
@@ -202,7 +388,6 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False,
# set up variables
self.model = model
self.modelgrid = model.modelgrid
- self.arrays = {}
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)
@@ -246,12 +431,80 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False,
self.verts, self.iverts, self.zverts = \
self.get_3d_vertex_connectivity()
+ self.vtk_grid_type, self.file_extension = \
+ self._vtk_grid_type(vtk_grid_type)
+
+ self.binary = binary
+
return
- def add_array(self, name, a, array2d=False):
+ def _vtk_grid_type(self, vtk_grid_type='auto'):
+ """
+ Determine 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():
+ vtk_grid_type = 'ImageData'
+ elif self.modelgrid.is_rectilinear():
+ 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('"' + 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():
+ 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():
+ 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 add_array(self, name, a, array2d=False):
+ """
Adds an array to the vtk object
Parameters
@@ -263,12 +516,7 @@ def add_array(self, name, a, array2d=False):
the array to be added to the vtk object
array2d : bool
True if the array is 2d
-
"""
-
- if name == 'ibound':
- return
-
# if array is 2d reformat to 3d array
if array2d:
assert a.shape == self.shape2d
@@ -277,392 +525,205 @@ def add_array(self, name, a, array2d=False):
a = array
try:
-
assert a.shape == self.shape
except AssertionError:
return
- a = np.where(a == self.nanval, np.nan, a)
- a = a.astype(float)
- # idxs = np.argwhere(a == self.nanval)
+ # assign nan where nanval or ibound==0
+ where_to_nan = np.logical_or(a==self.nanval, self.ibound==0)
+ a = np.where(where_to_nan, np.nan, a)
- # add array to self.arrays
+ # add a copy of the array to self.arrays
+ a = a.astype(float)
self.arrays[name] = a
return
def write(self, output_file, timeval=None):
"""
-
- writes the stored arrays to vtk file
+ Writes the stored arrays to vtk file in XML format.
Parameters
----------
output_file : str
- output file name to write the vtk data
-
+ 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
"""
-
- # make sure file ends with vtu
- assert output_file.lower().endswith(".vtu")
-
- # get the active data cells based on the data arrays and ibound
- actwcells3d = self._configure_data_arrays()
- actwcells = actwcells3d.ravel()
-
- # get the indexes of the active cells
- idxs = np.argwhere(actwcells != 0).ravel()
-
- # get the verts and iverts to be output
- verts = [self.verts[idx] for idx in idxs]
- iverts = self._build_iverts(verts)
-
- # get the total number of cells and vertices
- ncells = len(iverts)
- npoints = ncells * 8
-
+ # output file
+ output_file = output_file + self.file_extension
if self.verbose:
print('Writing vtk file: ' + output_file)
- print('Number of point is {}, Number of cells is {}\n'.format(
- npoints, ncells))
-
- # open output file for writing
- f = open(output_file, 'w')
- # write xml
- indent_level = 0
- s = ''
- f.write(s + '\n')
- indent_level = start_tag(f, '',
- indent_level)
+ # initialize xml file
+ if self.binary:
+ xml = XmlWriterBinary(output_file)
+ else:
+ xml = XmlWriterAscii(output_file)
+ xml.add_attributes(type=self.vtk_grid_type)
- indent_level = start_tag(f, '', indent_level)
+ # grid type
+ xml.open_element(self.vtk_grid_type)
# if time value write time section
if timeval:
-
- indent_level = start_tag(f, '', indent_level)
-
- s = ''
- indent_level = start_tag(f, s, indent_level)
-
- f.write(indent_level * ' ' + '{}\n'.format(timeval))
-
- indent_level = end_tag(f, '', indent_level)
-
- indent_level = end_tag(f, '', indent_level)
-
- # piece
- s = ''.format(npoints, ncells)
- indent_level = start_tag(f, s, indent_level)
-
- # points
- s = ''
- indent_level = start_tag(f, s, indent_level)
-
- s = ''
- indent_level = start_tag(f, s, indent_level)
- assert (isinstance(self.modelgrid, StructuredGrid))
- for cell in verts:
- for row in cell:
- s = indent_level * ' ' + '{} {} {} \n'.format(*row)
- f.write(s)
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- # cells
- s = ''
- indent_level = start_tag(f, s, indent_level)
-
- s = ''
- indent_level = start_tag(f, s, indent_level)
- for row in iverts:
- s = indent_level * ' ' + ' '.join([str(i) for i in row]) + '\n'
- f.write(s)
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- s = ''
- indent_level = start_tag(f, s, indent_level)
- icount = 0
- for row in iverts:
- icount += len(row)
- s = indent_level * ' ' + '{} \n'.format(icount)
- f.write(s)
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- s = ''
- indent_level = start_tag(f, s, indent_level)
- for row in iverts:
- s = indent_level * ' ' + '{} \n'.format(self.cell_type)
- f.write(s)
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- # add cell data
- s = ''
- indent_level = start_tag(f, s, indent_level)
-
- # write data arrays to file
- for arrayName, arrayValues in self.arrays.items():
- self.write_data_array(f, indent_level, arrayName,
- arrayValues, actwcells3d)
-
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- # write arrays as point data if point scalars is set to True
- if self.point_scalars:
- s = ''
- indent_level = start_tag(f, s, indent_level)
- for array_name, array_values in self.arrays.items():
- self.write_point_value(f, indent_level, array_values,
- array_name, actwcells3d)
-
- s = ''
- indent_level = end_tag(f, s, indent_level)
-
- else:
- pass
-
- # end piece
- indent_level = end_tag(f, '', indent_level)
-
- # end unstructured grid
- indent_level = end_tag(f, '', indent_level)
-
- # end xml
- indent_level = end_tag(f, '', indent_level)
-
- # end file
- f.close()
- self.arrays.clear()
- return
-
- def write_binary(self, output_file):
-
- """
-
- outputs binary .vtu file
-
- Parameters
- ----------
-
- output_file : str
- vtk output file
-
- """
-
- # make sure file ends with vtu
- assert output_file.lower().endswith(".vtu")
-
- if self.verbose:
- print('writing binary vtk file')
-
- xml = BinaryXml(output_file)
- offset = 0
- grid_type = 'UnstructuredGrid'
-
- # get the active data cells based on the data arrays and ibound
- actwcells3d = self._configure_data_arrays()
- actwcells = actwcells3d.ravel()
-
- # get the indexes of the active cells
- idxs = np.argwhere(actwcells != 0).ravel()
-
- # get the verts and iverts to be output
- verts = [self.verts[idx] for idx in idxs]
- iverts = self._build_iverts(verts)
-
- # check if there is data to be written out
- if len(verts) == 0:
- # if not cannot write binary .vtu file
- return
-
- # get the total number of cells and vertices
- ncells = len(iverts)
- npoints = ncells * 8
-
- if self.verbose:
- print('Writing vtk file: ' + output_file)
- print('Number of point is {}, Number of cells is {}\n'.format(
- npoints, ncells))
-
- # format verts and iverts
- verts = np.array(verts)
- verts.reshape(npoints, 3)
- iverts = np.ascontiguousarray(iverts, np.float64)
-
- # write xml file info
- xml.open_element("VTKFile"). \
- add_attributes(type=grid_type, version="1.0",
- byte_order=self._get_byte_order(),
- header_type="UInt64")
- # unstructured grid
- xml.open_element(grid_type)
-
- # piece
- xml.open_element('Piece')
- xml.add_attributes(NumberOfPoints=npoints, NumberOfCells=ncells)
-
- # points
- xml.open_element('Points')
-
- xml.open_element('DataArray')
- xml.add_attributes(Name='points', NumberOfComponents='3',
- type='Float64',
- format='appended', offset=offset)
-
- # calculate the offset of the start of the next piece of data
- # offset is calculated from beginning of data section
- points_size = verts.size * verts[0].dtype.itemsize
- offset += points_size + 8
-
- xml.close_element('DataArray')
-
- xml.close_element('Points')
-
- # cells
- xml.open_element('Cells')
-
- # connectivity
- xml.open_element('DataArray')
- xml.add_attributes(Name='connectivity', NumberOfComponents='1',
- type='Float64',
- format='appended', offset=offset)
- conn_size = iverts.size * iverts[0].dtype.itemsize
- offset += conn_size + 8
-
- xml.close_element('DataArray')
-
- xml.open_element('DataArray')
- xml.add_attributes(Name='offsets', NumberOfComponents='1',
- type='Float64',
- format='appended', offset=offset)
- offsets_size = iverts.shape[0] * iverts[0].dtype.itemsize
- offset += offsets_size + 8
-
- xml.close_element('DataArray')
-
- xml.open_element('DataArray')
- xml.add_attributes(Name='types', NumberOfComponents='1',
- type='Float64',
- format='appended', offset=offset)
- types_size = iverts.shape[0] * iverts[0].dtype.itemsize
- offset += types_size + 8
-
- xml.close_element('DataArray')
-
- xml.close_element('Cells')
-
+ xml.open_element('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)
+ npoints = ncells * 8
+ if self.verbose:
+ print('Number of point is {}, Number of cells is {}\n'.format(
+ npoints, ncells))
+
+ # 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 = '0' + ' ' + str(self.modelgrid.ncol) + ' ' + \
+ '0' + ' ' + str(self.modelgrid.nrow) + ' ' + \
+ '0' + ' ' + str(self.modelgrid.nlay)
+ xml.add_attributes(WholeExtent=vtk_extent_str)
+ grid_extent = self.modelgrid.xyzextent
+ vtk_origin_str = str(grid_extent[0]) + ' ' + \
+ str(grid_extent[2]) + ' ' + \
+ str(grid_extent[4])
+ xml.add_attributes(Origin=vtk_origin_str)
+ vtk_spacing_str = str(self.modelgrid.delr[0]) + ' ' + \
+ str(self.modelgrid.delc[0]) + ' ' + \
+ str(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 = '0' + ' ' + str(self.modelgrid.ncol) + ' ' + \
+ '0' + ' ' + str(self.modelgrid.nrow) + ' ' + \
+ '0' + ' ' + str(self.modelgrid.nlay)
+ 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')
+
+ # cell data
xml.open_element('CellData')
- xml.add_attributes(Scalars='scalars')
- # format data arrays and store for later output
- processed_arrays = []
+ # loop through stored arrays
for name, a in self.arrays.items():
- a = a.ravel()[idxs]
- xml.open_element('DataArray')
- xml.add_attributes(Name=name, NumberOfComponents='1',
- type='Float64',
- format='appended', offset=offset)
- a = np.ascontiguousarray(a, np.float64)
- processed_arrays.append([a, a.size * a[0].dtype.itemsize])
- offset += processed_arrays[-1][-1] + 8
- xml.close_element('DataArray')
+ 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)
+ a = np.flip(a, axis=1)
+ xml.write_array(a, Name=name, NumberOfComponents='1')
+ # end cell data
xml.close_element('CellData')
- # for data array point scalars
if self.point_scalars:
-
+ # point data (i.e., values at vertices)
xml.open_element('PointData')
- xml.add_attributes(Scalars='scalars')
- # get output point arrays
# loop through stored arrays
for name, a in self.arrays.items():
# get the array values onto vertices
- verts_info = self.get_3d_vertex_connectivity(
- actwcells=actwcells3d, zvalues=a)
- # get values
- point_values_dict = verts_info[2]
- a = np.array([point_values_dict[cellid] for cellid in
- sorted(point_values_dict.keys())]).ravel()
-
- xml.open_element('DataArray')
- xml.add_attributes(Name=name, NumberOfComponents='1',
- type='Float64',
- format='appended', offset=offset)
- a = np.ascontiguousarray(a, np.float64)
- processed_arrays.append([a, a.size * a[0].dtype.itemsize])
- offset += processed_arrays[-1][-1] + 8
-
- xml.close_element('DataArray')
+ if self.vtk_grid_type == 'UnstructuredGrid':
+ _, _, zverts = self.get_3d_vertex_connectivity(
+ actwcells=actwcells3d, zvalues=a)
+ a = np.array(list(zverts.values()))
+ else:
+ a = self.extendedDataArray(a)
+ # flip "a" so coordinates increase along with indices as in
+ # vtk
+ a = np.flip(a, axis=0)
+ a = np.flip(a, axis=1)
+ xml.write_array(a, Name=name, NumberOfComponents='1')
+
+ # end point data
xml.close_element('PointData')
# end piece
xml.close_element('Piece')
- # end unstructured grid
- xml.close_element('UnstructuredGrid')
+ # end vtk_grid_type
+ xml.close_element(self.vtk_grid_type)
+
+ # finalize and close xml file
+ xml.final()
- # build data section
- xml.open_element("AppendedData").add_attributes(
- encoding="raw").add_text("_")
-
- xml.write_size(points_size)
- # format verts for output
- verts_x = np.ascontiguousarray(np.ravel(verts[:, :, 0]),
- np.float64)
- verts_y = np.ascontiguousarray(np.ravel(verts[:, :, 1]),
- np.float64)
- verts_z = np.ascontiguousarray(np.ravel(verts[:, :, 2]),
- np.float64)
- # write coordinates
- xml.write_coord_arrays(verts_x, verts_y, verts_z)
-
- # write iverts
- xml.write_size(conn_size)
- rav_iverts = np.ascontiguousarray(np.ravel(iverts), np.float64)
- xml.write_array(rav_iverts)
-
- xml.write_size(offsets_size)
- data = np.empty((iverts.shape[0]), np.float64)
- icount = 0
- for index, row in enumerate(iverts):
- icount += len(row)
- data[index] = icount
- xml.write_array(data)
-
- # write cell types (11)
- xml.write_size(types_size)
- data = np.empty((iverts.shape[0]), np.float64)
- data.fill(self.cell_type)
- xml.write_array(data)
-
- # write out the array scalars and array point scalars
- for a, block_size in processed_arrays:
- xml.write_size(block_size)
- xml.write_array(a)
-
- # end xml
- xml.close_element("AppendedData")
- xml.close_element('VTKFile')
- xml.close()
# clear arrays
self.arrays.clear()
@@ -671,12 +732,12 @@ def _configure_data_arrays(self):
Compares arrays and active cells to find where active data
exists, and what cells to output.
"""
-
# get 1d shape
shape1d = self.shape[0] * self.shape[1] * self.shape[2]
# build index array
ot_idx_array = np.zeros(shape1d, dtype=np.int)
+
# loop through arrays
for name in self.arrays:
array = self.arrays[name]
@@ -693,15 +754,11 @@ def _configure_data_arrays(self):
# reset the shape of the active data array
ot_idx_array = ot_idx_array.reshape(self.shape)
- # where the ibound is 0 set the active array to 0
- ot_idx_array[self.ibound == 0] = 0
return ot_idx_array
def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None):
-
"""
-
Builds x,y,z vertices
Parameters
@@ -720,7 +777,6 @@ def get_3d_vertex_connectivity(self, actwcells=None, zvalues=None):
dictionary of iverts
zvertsdict : dict
dictionary of zverts
-
"""
# set up active cells
if actwcells is None:
@@ -827,125 +883,12 @@ def extendedDataArray(self, dataArray):
index[1]])
neighList = np.array(neighList)
if neighList[neighList != self.nanval].shape[0] > 0:
- headMean = neighList[neighList != self.nanval].mean()
+ valMean = neighList[neighList != self.nanval].mean()
else:
- headMean = self.nanval
- matrix[lay, row, col] = headMean
+ valMean = self.nanval
+ matrix[lay, row, col] = valMean
return matrix
- @staticmethod
- def _get_byte_order():
- if sys.byteorder == "little":
- return "LittleEndian"
- else:
- return "BigEndian"
-
- @staticmethod
- def write_data_array(f, indent_level, arrayName, arrayValues,
- actWCells):
- """
-
- Writes the data array to the output vtk file
-
- Parameters
- ----------
- f : file object
- output vtk file
- indent_level : int
- current indent of the xml
- arrayName : str
- name of the output array
- arrayValues : array
- the data array being output
- actWCells : array
- array of the active cells
-
- """
-
- s = ''.format(
- arrayName)
- indent_level = start_tag(f, s, indent_level)
-
- # data
- nlay = arrayValues.shape[0]
-
- for lay in range(nlay):
- s = indent_level * ' '
- f.write(s)
- idx = (actWCells[lay] != 0)
- arrayValuesLay = arrayValues[lay][idx].flatten()
- for layValues in arrayValuesLay:
- s = ' {}'.format(layValues)
- f.write(s)
- f.write('\n')
-
- s = ''
- indent_level = end_tag(f, s, indent_level)
- return
-
- def write_point_value(self, f, indent_level, data_array, array_name,
- actwcells):
- """
- Writes the data array to the output vtk file as point scalars
- """
- # header tag
- s = ''.format(
- array_name)
- indent_level = start_tag(f, s, indent_level)
-
- # data
- verts_info = self.get_3d_vertex_connectivity(
- actwcells=actwcells, zvalues=data_array)
-
- zverts = verts_info[2]
-
- for cellid in sorted(zverts):
- for z in zverts[cellid]:
- s = indent_level * ' '
- f.write(s)
- s = ' {}'.format(z)
- f.write(s)
- f.write('\n')
-
- # ending tag
- s = ''
- indent_level = end_tag(f, s, indent_level)
- return
-
- @staticmethod
- def _build_iverts(verts):
- """
-
- Builds the iverts based on the vertices being output
-
- Parameters
- ----------
- verts : array
- vertices being output
-
- Returns
- -------
-
- iverts : array
- array of ivert values
-
- """
- ncells = len(verts)
- npoints = ncells * 8
- iverts = []
- ivert = []
- count = 1
- for i in range(npoints):
- ivert.append(i)
- if count == 8:
- iverts.append(ivert)
- ivert = []
- count = 0
- count += 1
- iverts = np.array(iverts)
-
- return iverts
-
def _get_names(in_list):
ot_list = []
@@ -958,10 +901,9 @@ def _get_names(in_list):
def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20,
- kstpkper=None, text=None, smooth=False,
- point_scalars=False, binary=False):
+ kstpkper=None, text=None, smooth=False, point_scalars=False,
+ vtk_grid_type='auto', binary=False):
"""
-
Exports cell by cell file to vtk
Parameters
@@ -984,13 +926,21 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20,
The text identifier for the record. Examples include
'RIVER LEAKAGE', 'STORAGE', 'FLOW RIGHT FACE', etc.
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
-
+ if True the output file will be binary, default is False
"""
mg = model.modelgrid
@@ -1054,7 +1004,8 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20,
# get model name
model_name = model.name
- vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars)
+ vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
# export data
addarray = False
@@ -1062,7 +1013,7 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20,
for kper in kperlist:
for kstp in kstplist:
- ot_base = '{}_CBC_KPER{}_KSTP{}.vtu'.format(
+ ot_base = '{}_CBC_KPER{}_KSTP{}'.format(
model_name, kper + 1, kstp + 1)
otfile = os.path.join(otfolder, ot_base)
pvdfile.write("""
@@ -1120,10 +1068,9 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20,
def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None,
- smooth=False, point_scalars=False,
+ smooth=False, point_scalars=False, vtk_grid_type='auto',
binary=False):
"""
-
Exports binary head file to vtk
Parameters
@@ -1141,13 +1088,21 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None,
A tuple containing the time step and stress period (kstp, kper).
The kstp and kper values are zero based.
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
-
+ if True the output file will be binary, default is False
"""
# setup output folder
@@ -1185,7 +1140,8 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None,
kstplist = list(set([x[0] for x in hds.get_kstpkper() if x[0] > -1]))
# set upt the vtk
- vtk = Vtk(model, smooth=smooth, point_scalars=point_scalars, nanval=nanval)
+ vtk = Vtk(model, smooth=smooth, point_scalars=point_scalars, nanval=nanval,
+ vtk_grid_type=vtk_grid_type, binary=binary)
# output data
count = 0
@@ -1193,14 +1149,11 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None,
for kstp in kstplist:
hdarr = hds.get_data((kstp, kper))
vtk.add_array('head', hdarr)
- ot_base = '{}_Heads_KPER{}_KSTP{}.vtu'.format(
+ ot_base = '{}_Heads_KPER{}_KSTP{}'.format(
model.name, kper + 1, kstp + 1)
otfile = os.path.join(otfolder, ot_base)
# vtk.write(otfile, timeval=totim_dict[(kstp, kper)])
- if binary:
- vtk.write_binary(otfile)
- else:
- vtk.write(otfile)
+ vtk.write(otfile)
pvdfile.write("""\n""".format(count, ot_base))
count += 1
@@ -1213,10 +1166,8 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None,
def export_array(model, array, output_folder, name, nanval=-1e+20,
array2d=False, smooth=False, point_scalars=False,
- binary=False):
-
+ vtk_grid_type='auto', binary=False):
"""
-
Export array to vtk
Parameters
@@ -1235,34 +1186,39 @@ def export_array(model, array, output_folder, name, nanval=-1e+20,
array2d : bool
True if array is 2d, default is False
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
-
+ 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 = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
vtk.add_array(name, array, array2d=array2d)
- otfile = os.path.join(output_folder, '{}.vtu'.format(name))
- if binary:
- vtk.write_binary(otfile)
- else:
- vtk.write(otfile)
+ otfile = os.path.join(output_folder, '{}'.format(name))
+ vtk.write(otfile)
return
def export_transient(model, array, output_folder, name, nanval=-1e+20,
array2d=False, smooth=False, point_scalars=False,
- binary=False):
+ vtk_grid_type='auto', binary=False):
"""
-
Export transient 2d array to vtk
Parameters
@@ -1281,13 +1237,21 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20,
array2d : bool
True if array is 2d, default is False
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
-
+ if True the output file will be binary, default is False
"""
if not os.path.exists(output_folder):
@@ -1295,10 +1259,15 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20,
to_tim = model.dis.get_totim()
- vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars)
+ vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
- if array2d:
+ if name.endswith('_'):
+ separator = ''
+ else:
+ separator = '_'
+ if array2d:
for kper in range(array.shape[0]):
t2d_array_kper = array[kper]
@@ -1308,21 +1277,17 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20,
vtk.add_array(name, t2d_array_input, array2d=True)
- ot_name = '{}_0{}'.format(name, kper + 1)
- ot_file = os.path.join(output_folder, '{}.vtu'.format(ot_name))
- vtk.write(ot_file, timeval=to_tim[kper])
+ otname = '{}'.format(name) + separator + '0{}'.format(kper + 1)
+ otfile = os.path.join(output_folder, '{}'.format(otname))
+ vtk.write(otfile, timeval=to_tim[kper])
else:
-
for kper in range(array.shape[0]):
vtk.add_array(name, array[kper])
- ot_name = '{}_0{}'.format(name, kper + 1)
- ot_file = os.path.join(output_folder, '{}.vtu'.format(ot_name))
- if binary:
- vtk.write_binary(ot_file)
- else:
- vtk.write(ot_file, timeval=to_tim[kper])
+ otname = '{}'.format(name) + separator + '0{}'.format(kper + 1)
+ otfile = os.path.join(output_folder, '{}'.format(otname))
+ vtk.write(otfile, timeval=to_tim[kper])
return
@@ -1340,10 +1305,9 @@ def trans_dict(in_dict, name, trans_array, array2d=False):
return in_dict
-def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
+def export_package(pak_model, pak_name, otfolder, vtkobj=None,
nanval=-1e+20, smooth=False, point_scalars=False,
- binary=False):
-
+ vtk_grid_type='auto', binary=False):
"""
Exports package to vtk
@@ -1354,7 +1318,7 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
the model of the package
pak_name : str
the name of the package
- ot_folder : str
+ otfolder : str
output folder to write the data
vtkobj : VTK instance
a vtk object (allows export_package to be called from
@@ -1362,12 +1326,21 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
nanval : scalar
no data value, default value is -1e20
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
+ if True the output file will be binary, default is False
"""
@@ -1375,13 +1348,14 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
if not vtkobj:
# if not build one
vtk = Vtk(pak_model, nanval=nanval, smooth=smooth,
- point_scalars=point_scalars)
+ point_scalars=point_scalars, vtk_grid_type=vtk_grid_type,
+ binary=binary)
else:
# otherwise use the vtk object that was supplied
vtk = vtkobj
- if not os.path.exists(ot_folder):
- os.mkdir(ot_folder)
+ if not os.path.exists(otfolder):
+ os.mkdir(otfolder)
# is there output data
has_output = False
@@ -1478,11 +1452,8 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
# write out data
# write array data
if len(vtk.arrays) > 0:
- ot_file = os.path.join(ot_folder, '{}.vtu'.format(pak_name))
- if binary:
- vtk.write_binary(ot_file)
- else:
- vtk.write(ot_file)
+ otfile = os.path.join(otfolder, '{}'.format(pak_name))
+ vtk.write(otfile)
# write transient data
if vtk_trans_dict:
@@ -1497,7 +1468,7 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
# else:
# time = None
# set up output file
- ot_file = os.path.join(ot_folder, '{} _0{}.vtu'.format(
+ otfile = os.path.join(otfolder, '{}_0{}'.format(
pak_name, kper + 1))
for name, array in sorted(array_dict.items()):
if array.array2d:
@@ -1506,17 +1477,19 @@ def export_package(pak_model, pak_name, ot_folder, vtkobj=None,
else:
a = array.array
vtk.add_array(name, a, array.array2d)
- # vtk.write(ot_file, timeval=time)
- if binary:
- vtk.write_binary(ot_file)
- else:
- vtk.write(ot_file)
+ # vtk.write(otfile, timeval=time)
+ vtk.write(otfile)
return
-def export_model(model, ot_folder, package_names=None, nanval=-1e+20,
- smooth=False, point_scalars=False, binary=False):
+def export_model(model, otfolder, package_names=None, nanval=-1e+20,
+ smooth=False, point_scalars=False, vtk_grid_type='auto',
+ binary=False):
"""
+ Exports model to vtk
+
+ Parameters
+ ----------
model : flopy model instance
flopy model
@@ -1529,15 +1502,24 @@ def export_model(model, ot_folder, package_names=None, nanval=-1e+20,
array2d : bool
True if array is 2d, default is False
smooth : bool
- If true a smooth surface will be output, default is False
+ if True, will create smooth layer elevations, default is False
point_scalars : bool
- If True point scalar values will be written, default is False
+ 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'. 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'.
binary : bool
- if True the output .vtu file will be binary, default is
- False.
-
+ if True the output file will be binary, default is False
"""
- vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars)
+ vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars,
+ vtk_grid_type=vtk_grid_type, binary=binary)
if package_names is not None:
if not isinstance(package_names, list):
@@ -1545,10 +1527,10 @@ def export_model(model, ot_folder, package_names=None, nanval=-1e+20,
else:
package_names = [pak.name[0] for pak in model.packagelist]
- if not os.path.exists(ot_folder):
- os.mkdir(ot_folder)
+ if not os.path.exists(otfolder):
+ os.mkdir(otfolder)
for pak_name in package_names:
- export_package(model, pak_name, ot_folder, vtkobj=vtk, nanval=nanval,
+ export_package(model, pak_name, otfolder, vtkobj=vtk, nanval=nanval,
smooth=smooth, point_scalars=point_scalars,
- binary=binary)
+ vtk_grid_type=vtk_grid_type, binary=binary)