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

@@ -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()

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(

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

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))