From 0867b3490cb3da06b8b5dd41b8edf9745d085d22 Mon Sep 17 00:00:00 2001 From: Etienne Bresciani <53166244+etiennebresciani@users.noreply.github.com> Date: Sat, 28 Mar 2020 15:53:34 +0900 Subject: [PATCH 1/2] refactor and fix(vtk): improving code and fixing issues * Create XmlWriterInterface class and child classes (XmlWriterAscii and XmlWriterBinary), allowing to factorize lot of code. Namely, merge write() and write_binary(). binary is now a parameter of the Vtk class, consistently with other parameters. In addition, ascii and binary files are now written more consistently. Note: the written files are slightly different than before, but this only concerns minor style with no functionalities being impaired. * Use more appropriate data types for connectivity, offsets and cell types * Fix vtk.export_transient and utils.transient2d_export for binary option (it was actually saving in ascii) * Fix weird file name separators (i.e., "__", " _") in transient cases * Fix issue with inactive cells not ignored when smoothing : now ignore them (i.e., in extendedDataArray()) * Export ibound as well (it was excluded, but I see no reason why) * Improve code documentation * Simplify tests (t050) and make them more robust --- autotest/t050_test.py | 336 +++++++++----- flopy/export/utils.py | 119 ++--- flopy/export/vtk.py | 1014 ++++++++++++++++------------------------- 3 files changed, 673 insertions(+), 796 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 02eb9fc684..c7250d70e7 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) @@ -128,87 +219,126 @@ def test_vtk_mf6(): 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 - 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() diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 4b4a9eba89..ddeb410216 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,13 @@ 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) 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, + binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -591,8 +576,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 +586,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 +629,13 @@ 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) 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, + binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -874,18 +846,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 +954,14 @@ 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) + binary = kwargs.get('binary', False) + vtk.export_transient(t2d.model, t2d.array, f, name, nanval=nanval, + smooth=smooth, point_scalars=point_scalars, + array2d=True, binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -1006,18 +975,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 +1102,17 @@ 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) 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, + binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -1162,18 +1127,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 +1230,14 @@ 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) 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, 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..7878a1bf39 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() + + +class XmlWriterAscii(XmlWriterInterface): + """ + Helps writing ascii vtk files. + + Parameters + ---------- + + 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 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 + """ + 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 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 + 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,25 @@ 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 + 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, binary=False): if point_scalars: smooth = True @@ -202,7 +379,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 +422,12 @@ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, self.verts, self.iverts, self.zverts = \ self.get_3d_vertex_connectivity() + self.binary = binary + return def add_array(self, name, a, array2d=False): - """ - Adds an array to the vtk object Parameters @@ -263,12 +439,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,242 +448,73 @@ 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 + '.vtu' 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') + print('Writing vtk file: ' + output_file) - # 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='UnstructuredGrid') - indent_level = start_tag(f, '', indent_level) + # grid type + xml.open_element('UnstructuredGrid') # 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' + xml.open_element('FieldData') + xml.write_array(np.array([timeval]), Name='TimeValue', + NumberOfTuples='1', RangeMin='{0}', RangeMax='{0}') + xml.close_element('FieldData') # 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) + verts, iverts, _ = \ + self.get_3d_vertex_connectivity(actwcells=actwcells3d) # check if there is data to be written out if len(verts) == 0: - # if not cannot write binary .vtu file + # if nothing, cannot write 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) + npoints, ncells)) # piece xml.open_element('Piece') @@ -520,149 +522,69 @@ def write_binary(self, output_file): # 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') - + 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 - 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 + iverts = np.array(list(iverts.values())) + xml.write_array(iverts, Name='connectivity', + NumberOfComponents='1') - 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 + # 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') - xml.close_element('DataArray') + # 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') + # 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') + xml.write_array(a, actwcells=actwcells3d, 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') + _, _, zverts = self.get_3d_vertex_connectivity( + actwcells=actwcells3d, zvalues=a) + a = np.array(list(zverts.values())) + xml.write_array(a, Name=name, NumberOfComponents='1') + + # end point data xml.close_element('PointData') # end piece xml.close_element('Piece') - # end unstructured grid + # end grid type xml.close_element('UnstructuredGrid') - # build data section - xml.open_element("AppendedData").add_attributes( - encoding="raw").add_text("_") - - xml.write_size(points_size) - # format verts for output - verts_x = np.ascontiguousarray(np.ravel(verts[:, :, 0]), - np.float64) - verts_y = np.ascontiguousarray(np.ravel(verts[:, :, 1]), - np.float64) - verts_z = np.ascontiguousarray(np.ravel(verts[:, :, 2]), - np.float64) - # write coordinates - xml.write_coord_arrays(verts_x, verts_y, verts_z) - - # write iverts - xml.write_size(conn_size) - rav_iverts = np.ascontiguousarray(np.ravel(iverts), np.float64) - xml.write_array(rav_iverts) - - xml.write_size(offsets_size) - data = np.empty((iverts.shape[0]), np.float64) - icount = 0 - for index, row in enumerate(iverts): - icount += len(row) - data[index] = icount - xml.write_array(data) - - # write cell types (11) - xml.write_size(types_size) - data = np.empty((iverts.shape[0]), np.float64) - data.fill(self.cell_type) - xml.write_array(data) - - # write out the array scalars and array point scalars - for a, block_size in processed_arrays: - xml.write_size(block_size) - xml.write_array(a) - - # end xml - xml.close_element("AppendedData") - xml.close_element('VTKFile') - xml.close() + # finalize and close xml file + xml.final() + # clear arrays self.arrays.clear() @@ -671,12 +593,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 +615,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 +638,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 +744,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 +762,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, + binary=False): """ - Exports cell by cell file to vtk Parameters @@ -984,13 +787,12 @@ 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 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 +856,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, + binary=binary) # export data addarray = False @@ -1062,7 +865,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 +920,8 @@ 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, - binary=False): + smooth=False, point_scalars=False, binary=False): """ - Exports binary head file to vtk Parameters @@ -1141,13 +939,12 @@ 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 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 +982,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, + binary=binary) # output data count = 0 @@ -1193,14 +991,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 @@ -1214,9 +1009,7 @@ 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): - """ - Export array to vtk Parameters @@ -1235,25 +1028,22 @@ 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 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, + 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 @@ -1262,7 +1052,6 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20, array2d=False, smooth=False, point_scalars=False, binary=False): """ - Export transient 2d array to vtk Parameters @@ -1281,13 +1070,12 @@ 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 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 +1083,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, + 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 +1101,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 +1129,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): - """ Exports package to vtk @@ -1354,7 +1142,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 +1150,12 @@ 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 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 +1163,13 @@ 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, 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 +1266,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 +1282,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 +1291,18 @@ 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, +def export_model(model, otfolder, package_names=None, nanval=-1e+20, smooth=False, point_scalars=False, binary=False): """ + Exports model to vtk + + Parameters + ---------- model : flopy model instance flopy model @@ -1529,15 +1315,15 @@ 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 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, + binary=binary) if package_names is not None: if not isinstance(package_names, list): @@ -1545,10 +1331,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) From 13c1dd28b9c470c70aac624b70888f5f6badf795 Mon Sep 17 00:00:00 2001 From: Etienne Bresciani <53166244+etiennebresciani@users.noreply.github.com> Date: Sat, 28 Mar 2020 16:06:14 +0900 Subject: [PATCH 2/2] feat(vtk): export in .vti and .vtr when possible * Export data in vtk ImageData (.vti) or RectilinearGrid (.vtr) formats when possible, instead of always UnstructuredGrid (.vtu), since this can greatly reduce file volume and speed up subsequent operations (e.g. in Paraview). The user can also choose to force an export in a specific format (e.g., .vtu) by using the new parameter vtk_grid_type (by default = 'auto') * Add tests for it to t050 --- autotest/t050_test.py | 104 ++++++++ flopy/discretization/grid.py | 6 + flopy/discretization/structuredgrid.py | 68 +++++ flopy/export/utils.py | 17 +- flopy/export/vtk.py | 348 +++++++++++++++++++------ 5 files changed, 462 insertions(+), 81 deletions(-) diff --git a/autotest/t050_test.py b/autotest/t050_test.py index c7250d70e7..bbcd77428b 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -215,6 +215,13 @@ 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 @@ -334,6 +341,101 @@ def test_vtk_cbc(): 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() @@ -342,3 +444,5 @@ def test_vtk_cbc(): 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 ddeb410216..1e23712374 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -555,10 +555,11 @@ def model_export(f, ml, fmt=None, **kwargs): 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) vtk.export_model(ml, f, package_names=package_names, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -632,10 +633,11 @@ def package_export(f, pak, fmt=None, **kwargs): 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) vtk.export_package(pak.parent, pak.name, f, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -958,10 +960,12 @@ def transient2d_export(f, t2d, fmt=None, **kwargs): 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) vtk.export_transient(t2d.model, t2d.array, f, name, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - array2d=True, binary=binary) + array2d=True, vtk_grid_type=vtk_grid_type, + binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -1106,13 +1110,14 @@ def array3d_export(f, u3d, fmt=None, **kwargs): 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) if isinstance(name, list) or isinstance(name, tuple): name = name[0] vtk.export_array(u3d.model, u3d.array, f, name, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) else: raise NotImplementedError("unrecognized export argument:{0}".format(f)) @@ -1234,10 +1239,12 @@ def array2d_export(f, u2d, fmt=None, **kwargs): 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) vtk.export_array(u2d.model, u2d.array, f, name, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - array2d=True, binary=binary) + 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 7878a1bf39..1fed49a8a7 100644 --- a/flopy/export/vtk.py +++ b/flopy/export/vtk.py @@ -357,6 +357,15 @@ class Vtk(object): 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'. 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 @@ -367,7 +376,7 @@ class Vtk(object): Stores data arrays added to VTK object """ def __init__(self, model, verbose=None, nanval=-1e+20, smooth=False, - point_scalars=False, binary=False): + point_scalars=False, vtk_grid_type='auto', binary=False): if point_scalars: smooth = True @@ -422,10 +431,78 @@ 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 _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 @@ -476,19 +553,19 @@ def write(self, output_file, timeval=None): file, default is None """ # output file - output_file = output_file + '.vtu' + output_file = output_file + self.file_extension if self.verbose: - print('Writing vtk file: ' + output_file) + print('Writing vtk file: ' + output_file) # initialize xml file if self.binary: xml = XmlWriterBinary(output_file) else: xml = XmlWriterAscii(output_file) - xml.add_attributes(type='UnstructuredGrid') + xml.add_attributes(type=self.vtk_grid_type) # grid type - xml.open_element('UnstructuredGrid') + xml.open_element(self.vtk_grid_type) # if time value write time section if timeval: @@ -497,66 +574,121 @@ def write(self, output_file, timeval=None): NumberOfTuples='1', RangeMin='{0}', RangeMax='{0}') xml.close_element('FieldData') - # 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 + 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') - # 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') + # 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') # loop through stored arrays for name, a in self.arrays.items(): - xml.write_array(a, actwcells=actwcells3d, Name=name, - NumberOfComponents='1') + 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') @@ -568,9 +700,16 @@ def write(self, output_file, timeval=None): # loop through stored arrays for name, a in self.arrays.items(): # get the array values onto vertices - _, _, zverts = self.get_3d_vertex_connectivity( - actwcells=actwcells3d, zvalues=a) - a = np.array(list(zverts.values())) + 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 @@ -579,8 +718,8 @@ def write(self, output_file, timeval=None): # end piece xml.close_element('Piece') - # end grid type - xml.close_element('UnstructuredGrid') + # end vtk_grid_type + xml.close_element(self.vtk_grid_type) # finalize and close xml file xml.final() @@ -763,7 +902,7 @@ 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): + vtk_grid_type='auto', binary=False): """ Exports cell by cell file to vtk @@ -791,6 +930,15 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, 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'. 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 file will be binary, default is False """ @@ -857,7 +1005,7 @@ def export_cbc(model, cbcfile, otfolder, precision='single', nanval=-1e+20, model_name = model.name vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) # export data addarray = False @@ -920,7 +1068,8 @@ 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, binary=False): + smooth=False, point_scalars=False, vtk_grid_type='auto', + binary=False): """ Exports binary head file to vtk @@ -943,6 +1092,15 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None, 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'. 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 file will be binary, default is False """ @@ -983,7 +1141,7 @@ def export_heads(model, hdsfile, otfolder, nanval=-1e+20, kstpkper=None, # set upt the vtk vtk = Vtk(model, smooth=smooth, point_scalars=point_scalars, nanval=nanval, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) # output data count = 0 @@ -1008,7 +1166,7 @@ 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 @@ -1032,6 +1190,15 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, 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'. 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 file will be binary, default is False """ @@ -1040,7 +1207,7 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, os.mkdir(output_folder) vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) vtk.add_array(name, array, array2d=array2d) otfile = os.path.join(output_folder, '{}'.format(name)) vtk.write(otfile) @@ -1050,7 +1217,7 @@ def export_array(model, array, output_folder, name, nanval=-1e+20, def export_transient(model, array, output_folder, name, nanval=-1e+20, array2d=False, smooth=False, point_scalars=False, - binary=False): + vtk_grid_type='auto', binary=False): """ Export transient 2d array to vtk @@ -1074,6 +1241,15 @@ def export_transient(model, array, output_folder, name, nanval=-1e+20, 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'. 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 file will be binary, default is False """ @@ -1084,7 +1260,7 @@ 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, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) if name.endswith('_'): separator = '' @@ -1131,7 +1307,7 @@ def trans_dict(in_dict, name, trans_array, array2d=False): 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 @@ -1154,6 +1330,15 @@ def export_package(pak_model, pak_name, otfolder, vtkobj=None, 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'. 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 file will be binary, default is False @@ -1163,7 +1348,8 @@ def export_package(pak_model, pak_name, otfolder, vtkobj=None, if not vtkobj: # if not build one vtk = Vtk(pak_model, nanval=nanval, smooth=smooth, - point_scalars=point_scalars, binary=binary) + point_scalars=point_scalars, vtk_grid_type=vtk_grid_type, + binary=binary) else: # otherwise use the vtk object that was supplied vtk = vtkobj @@ -1297,7 +1483,8 @@ def export_package(pak_model, pak_name, otfolder, vtkobj=None, def export_model(model, otfolder, package_names=None, nanval=-1e+20, - smooth=False, point_scalars=False, binary=False): + smooth=False, point_scalars=False, vtk_grid_type='auto', + binary=False): """ Exports model to vtk @@ -1319,11 +1506,20 @@ def export_model(model, otfolder, package_names=None, nanval=-1e+20, 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'. 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 file will be binary, default is False """ vtk = Vtk(model, nanval=nanval, smooth=smooth, point_scalars=point_scalars, - binary=binary) + vtk_grid_type=vtk_grid_type, binary=binary) if package_names is not None: if not isinstance(package_names, list): @@ -1337,4 +1533,4 @@ def export_model(model, otfolder, package_names=None, nanval=-1e+20, for pak_name in package_names: 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)