feat(vtk): improve export at vertices

* Use bilinear interpolation to calculate the z of cell vertices (this is more appropriate than simple averaging when the grid is not regular)
* Ignore inactive cells when calculating the z of cell vertices (until now, they were not ignored, potentially yielding wrong values at the boundary of inactive cells)
* Use trilinear interpolation to calculate array values at cell vertices (this is more appropriate than simple averaging when the grid is not regular; furthermore, until now, everything was shifted downward with an awkward repetition of values at the top and bottom of the first layer)
* Improve performance by using scipy functions and avoiding loops (until now, interpolation was done by looping through all the cells in Python)
* Allow for directly passing array values at vertices to functions add_array() and add_vector() of the Vtk class
* New option "position" in the function postprocessing.get_specific_discharge() for calculating the values at "centers" (default), "faces" or "vertices"
* Cast delc, delr, top and botm of the grid to float64 (float32 was yielding significant round-off errors when rotating grid coordinates)

README.md

@@ -209,6 +209,8 @@ Additional dependencies to use optional FloPy helper methods are listed below.
 | `Raster()` in `flopy.utils.Raster`                                                   | **rasterio**, **affine**, and **scipy**            |
 | `Raster().sample_polygon()` in `flopy.utils.Raster`                                  | **shapely**                                        |
 | `Raster().crop()` in `flopy.utils.Raster`                                            | **shapely**                                        |
+| `._zverts_smooth()` in `flopy.discretization.structuredgrid` `StructuredGrid` class  | **scipy.interpolate**                              |
+| `.array_at_verts()` in `flopy.discretization.structuredgrid` `StructuredGrid` class  | **scipy.interpolate**                              |
 
 How to Cite
 -----------------------------------------------

autotest/t050_test.py

@@ -39,15 +39,15 @@ def test_vtk_export_array2d():
     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)
+    # assert(totalbytes==351846)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==2846)
 
     # with smoothing
     m.dis.top.export(output_dir, fmt='vtk', name='top_smooth', smooth=True)
     filetocheck = os.path.join(output_dir, 'top_smooth.vtu')
     totalbytes1 = os.path.getsize(filetocheck)
-    # assert(totalbytes1==351829)
+    assert(totalbytes1==357587)
     nlines1 = count_lines_in_file(filetocheck)
     assert(nlines1==2846)
 
@@ -66,7 +66,7 @@ def test_vtk_export_array3d():
     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)
+    # assert(totalbytes==992036)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==8486)
 
@@ -75,7 +75,7 @@ def test_vtk_export_array3d():
                     point_scalars=True)
     filetocheck = os.path.join(output_dir, 'hk_points.vtu')
     totalbytes1 = os.path.getsize(filetocheck)
-    # assert(totalbytes1==1321502)
+    assert(totalbytes1==1354204)
     nlines1 = count_lines_in_file(filetocheck)
     assert(nlines1==10605)
 
@@ -84,9 +84,9 @@ def test_vtk_export_array3d():
                     point_scalars=True, binary=True)
     filetocheck = os.path.join(output_dir, 'hk_points_bin.vtu')
     totalbytes2 = os.path.getsize(filetocheck)
-    # assert(totalbytes2==629401)
+    assert(totalbytes2==629401)
     nlines2 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines2==1869)
+    assert(nlines2==2257)
 
     return
 
@@ -104,27 +104,27 @@ def test_vtk_transient_array_2d():
     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)
+    # assert(totalbytes==355144)
     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)
+    assert(totalbytes1==354442)
     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)
+    assert(totalbytes2==168339)
     nlines2 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines2==762)
+    assert(nlines2==846)
     filetocheck = os.path.join(output_dir_bin, 'rech_01097.vtu')
     totalbytes3 = os.path.getsize(filetocheck)
-    # assert(totalbytes3==168339)
+    assert(totalbytes3==168339)
     nlines3 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines3==762)
+    assert(nlines3==846)
 
     return
 
@@ -141,7 +141,7 @@ def test_vtk_export_packages():
     m.dis.export(output_dir, fmt='vtk')
     filetocheck = os.path.join(output_dir, 'DIS.vtu')
     totalbytes = os.path.getsize(filetocheck)
-    # assert(totalbytes==1020397)
+    # assert(totalbytes==1019857)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==8496)
 
@@ -150,7 +150,7 @@ def test_vtk_export_packages():
     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==2485991)
+    assert(totalbytes1==2531309)
     nlines1 = count_lines_in_file(filetocheck)
     assert(nlines1==21215)
 
@@ -159,7 +159,7 @@ def test_vtk_export_packages():
     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)
+    assert(totalbytes2==1035056)
     nlines2 = count_lines_in_file(filetocheck)
     assert(nlines2==8491)
 
@@ -168,12 +168,12 @@ def test_vtk_export_packages():
     m.drn.export(output_dir, fmt='vtk')
     filetocheck = os.path.join(output_dir, 'DRN_01.vtu')
     totalbytes3 = os.path.getsize(filetocheck)
-    # assert(totalbytes3==20702)
+    assert(totalbytes3==20670)
     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)
+    assert(totalbytes4==20670)
     nlines4 = count_lines_in_file(filetocheck)
     assert(nlines4==191)
 
@@ -182,18 +182,18 @@ def test_vtk_export_packages():
     m.dis.export(output_dir, fmt='vtk', binary=True)
     filetocheck = os.path.join(output_dir, 'DIS.vtu')
     totalbytes5 = os.path.getsize(filetocheck)
-    # assert(totalbytes5==519516)
+    assert(totalbytes5==519516)
     nlines5 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines5==1545)
+    assert(nlines5==1797)
 
     # upw with point scalars and binary
     output_dir = os.path.join(cpth, 'UPW_bin')
     m.upw.export(output_dir, fmt='vtk', point_scalars=True, binary=True)
     filetocheck = os.path.join(output_dir, 'UPW.vtu')
     totalbytes6 = os.path.getsize(filetocheck)
-    # assert(totalbytes6==1349801)
+    assert(totalbytes6==1349801)
     nlines6 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines6==4004)
+    assert(nlines6==4392)
 
     return
 
@@ -245,7 +245,7 @@ def test_vtk_binary_head_export():
                                                                       1089)])
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes = os.path.getsize(filetocheck)
-    # assert(totalbytes==993755)
+    # assert(totalbytes==993215)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==8486)
 
@@ -257,7 +257,7 @@ def test_vtk_binary_head_export():
                      point_scalars=True, nanval=-999.99)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes1 = os.path.getsize(filetocheck)
-    # assert(totalbytes1==1332153)
+    assert(totalbytes1==1365320)
     nlines1 = count_lines_in_file(filetocheck)
     assert(nlines1==10605)
 
@@ -269,7 +269,7 @@ def test_vtk_binary_head_export():
                      smooth=True, nanval=-999.99)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes2 = os.path.getsize(filetocheck)
-    # assert(totalbytes2==993551)
+    assert(totalbytes2==1026553)
     nlines2 = count_lines_in_file(filetocheck)
     assert(nlines2==8486)
 
@@ -281,9 +281,9 @@ def test_vtk_binary_head_export():
                      smooth=True, binary=True, nanval=-999.99)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes3 = os.path.getsize(filetocheck)
-    # assert(totalbytes3==493853)
+    assert(totalbytes3==493853)
     nlines3 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines3==1529)
+    assert(nlines3==1825)
 
     # with smoothing and binary, single time
     otfolder = os.path.join(cpth, 'heads_test_4')
@@ -292,9 +292,9 @@ def test_vtk_binary_head_export():
                      nanval=-999.99)
     filetocheck = os.path.join(otfolder, 'freyberg_Heads_KPER1_KSTP1.vtu')
     totalbytes4 = os.path.getsize(filetocheck)
-    # assert(totalbytes4==493853)
+    assert(totalbytes4==493853)
     nlines4 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines4==1535)
+    assert(nlines4==1831)
 
     return
 
@@ -314,7 +314,7 @@ def test_vtk_cbc():
                    kstpkper=[(0, 0), (0, 1), (0, 2)], point_scalars=True)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes = os.path.getsize(filetocheck)
-    # assert(totalbytes==2626880)
+    # assert(totalbytes==2664009)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==19093)
 
@@ -325,9 +325,9 @@ def test_vtk_cbc():
                    binary=True)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes1 = os.path.getsize(filetocheck)
-    # assert(totalbytes1==1205818)
+    assert(totalbytes1==1205818)
     nlines1 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines1==2514)
+    assert(nlines1==3310)
 
     # with point scalars and binary, only one budget component
     otfolder = os.path.join(cpth, 'freyberg_CBCTEST_bin2')
@@ -336,41 +336,68 @@ def test_vtk_cbc():
                    point_scalars=True,  binary=True)
     filetocheck = os.path.join(otfolder, filenametocheck)
     totalbytes2 = os.path.getsize(filetocheck)
-    # assert(totalbytes2==10142)
+    assert(totalbytes2==10142)
     nlines2 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines2==62)
+    assert(nlines2==66)
 
     return
 
 def test_vtk_vector():
     from flopy.utils import postprocessing as pp
+    from flopy.utils import geometry
     # test mf 2005 freyberg
     mpth = os.path.join('..', 'examples', 'data',
                         'freyberg_multilayer_transient')
     namfile = 'freyberg.nam'
     cbcfile = os.path.join(mpth, 'freyberg.cbc')
+    hdsfile = os.path.join(mpth, 'freyberg.hds')
     m = flopy.modflow.Modflow.load(namfile, model_ws=mpth, verbose=False,
                                    load_only=['dis', 'bas6'])
-    q = pp.get_specific_discharge(m, cbcfile=cbcfile)
+    q = list(pp.get_specific_discharge(m, cbcfile=cbcfile))
+    q[0], q[1] = geometry.rotate(q[0], q[1], 0., 0.,
+                                 m.modelgrid.angrot_radians)
     output_dir = os.path.join(cpth, 'freyberg_vector')
     filenametocheck = 'discharge.vtu'
 
     # export and check with point scalar
     vtk.export_vector(m, q, output_dir, 'discharge', point_scalars=True)
     filetocheck = os.path.join(output_dir, filenametocheck)
     totalbytes = os.path.getsize(filetocheck)
-    # assert(totalbytes==2249214)
+    # assert(totalbytes==2281411)
     nlines = count_lines_in_file(filetocheck)
     assert(nlines==10605)
 
     # with point scalars and binary
-    vtk.export_vector(m, q, output_dir + '_bin', 'discharge', point_scalars=True,
-                      binary=True)
+    vtk.export_vector(m, q, output_dir + '_bin', 'discharge',
+                      point_scalars=True, binary=True)
     filetocheck = os.path.join(output_dir + '_bin', filenametocheck)
     totalbytes1 = os.path.getsize(filetocheck)
-    # assert(totalbytes1==917033)
+    assert(totalbytes1==942413)
     nlines1 = count_lines_in_file(filetocheck, binary=True)
-    assert(nlines1==2725)
+    assert(nlines1==4014)
+
+    # with values directly given at vertices
+    q = list(pp.get_specific_discharge(m, cbcfile=cbcfile, hdsfile=hdsfile,
+                                       position='vertices'))
+    q[0], q[1] = geometry.rotate(q[0], q[1], 0., 0.,
+                                 m.modelgrid.angrot_radians)
+    output_dir = os.path.join(cpth, 'freyberg_vector')
+    filenametocheck = 'discharge_verts.vtu'
+    vtk.export_vector(m, q, output_dir, 'discharge_verts')
+    filetocheck = os.path.join(output_dir, filenametocheck)
+    totalbytes2 = os.path.getsize(filetocheck)
+    assert(totalbytes2==2039528)
+    nlines2 = count_lines_in_file(filetocheck)
+    assert(nlines2==10598)
+
+    # with values directly given at vertices and binary
+    vtk.export_vector(m, q, output_dir + '_bin', 'discharge_verts',
+                      binary=True)
+    filetocheck = os.path.join(output_dir + '_bin', filenametocheck)
+    totalbytes3 = os.path.getsize(filetocheck)
+    assert(totalbytes3==891486)
+    nlines3 = count_lines_in_file(filetocheck, binary=True)
+    assert(nlines3==3113)
 
     return
 
@@ -394,15 +421,15 @@ def test_vtk_vti():
     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)
+    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==4617)
+    assert(totalbytes2==4617)
     nlines2 = count_lines_in_file(filetocheck, binary=True)
     assert(nlines2==18)
 
@@ -411,7 +438,7 @@ def test_vtk_vti():
     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)
+    assert(totalbytes3==7146)
     nlines3 = count_lines_in_file(filetocheck)
     assert(nlines3==56)
 
@@ -420,7 +447,7 @@ def test_vtk_vti():
     m.export(output_dir, fmt='vtk', vtk_grid_type='UnstructuredGrid')
     filetocheck = os.path.join(output_dir, filenametocheck)
     totalbytes4 = os.path.getsize(filetocheck)
-    # assert(totalbytes4==67905)
+    assert(totalbytes4==67905)
     nlines4 = count_lines_in_file(filetocheck)
     assert(nlines4==993)
 
@@ -430,16 +457,16 @@ def test_vtk_vti():
     vtk.export_vector(m, T, output_dir, 'T', point_scalars=True)
     filetocheck = os.path.join(output_dir, filenametocheck)
     totalbytes5 = os.path.getsize(filetocheck)
-    # assert(totalbytes5==12621)
+    assert(totalbytes5==12621)
     nlines5 = count_lines_in_file(filetocheck)
     assert(nlines5==20)
 
-    # vector binary
+    # vector with point scalars and binary
     vtk.export_vector(m, T, output_dir + '_bin', 'T', point_scalars=True,
                       binary=True)
     filetocheck = os.path.join(output_dir + '_bin', filenametocheck)
     totalbytes6 = os.path.getsize(filetocheck)
-    # assert(totalbytes6==16716)
+    assert(totalbytes6==16716)
     nlines6 = count_lines_in_file(filetocheck, binary=True)
     assert(nlines6==18)
 
@@ -465,15 +492,15 @@ def test_vtk_vtr():
     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==182168)
+    assert(totalbytes1==182168)
     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)
+    assert(totalbytes2==47874)
     nlines2 = count_lines_in_file(filetocheck, binary=True)
     assert(nlines2==28)
 
@@ -482,7 +509,7 @@ def test_vtk_vtr():
     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)
+    assert(totalbytes3==78762)
     nlines3 = count_lines_in_file(filetocheck)
     assert(nlines3==1105)
 

flopy/discretization/grid.py

@@ -136,7 +136,11 @@ def __init__(self, grid_type=None, top=None, botm=None, idomain=None,
         LENUNI = {"u": 0, "f": 1, "m": 2, "c": 3}
         self.use_ref_coords = True
         self._grid_type = grid_type
+        if top is not None:
+            top = top.astype(float)
         self._top = top
+        if botm is not None:
+            botm = botm.astype(float)
         self._botm = botm
         self._idomain = idomain