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