diff --git a/autotest/t007_test.py b/autotest/t007_test.py
index 9d8fdc8de6..7c4751d32c 100644
--- a/autotest/t007_test.py
+++ b/autotest/t007_test.py
@@ -262,43 +262,72 @@ def test_export_output():
 def test_write_shapefile():
     from flopy.discretization import StructuredGrid
     from flopy.export.shapefile_utils import shp2recarray
-    from flopy.export.shapefile_utils import write_grid_shapefile, write_grid_shapefile2
+    from flopy.export.shapefile_utils import write_grid_shapefile
 
     sg = StructuredGrid(delr=np.ones(10) *1.1,  # cell spacing along model rows
                         delc=np.ones(10) *1.1,  # cell spacing along model columns
                         epsg=26715)
-    outshp1 = os.path.join(tpth, 'junk.shp')
-    outshp2 = os.path.join(tpth, 'junk2.shp')
-    write_grid_shapefile(outshp1, sg, array_dict={})
-    write_grid_shapefile2(outshp2, sg, array_dict={})
+    outshp = os.path.join(tpth, 'junk.shp')
+    write_grid_shapefile(outshp, sg, array_dict={})
 
     # test that vertices aren't getting altered by writing shapefile
-    for outshp in [outshp1, outshp2]:
-        # check that pyshp reads integers
-        # this only check that row/column were recorded as "N"
-        # not how they will be cast by python or numpy
-        import shapefile as sf
-        sfobj = sf.Reader(outshp)
-        for f in sfobj.fields:
-            if f[0] == 'row' or f[0] == 'column':
-                assert f[1] == 'N'
-        recs = list(sfobj.records())
-        for r in recs[0]:
-            assert isinstance(r, int)
-
-        # check that row and column appear as integers in recarray
-        ra = shp2recarray(outshp)
-        assert np.issubdtype(ra.dtype['row'], np.integer)
-        assert np.issubdtype(ra.dtype['column'], np.integer)
-
-        try: # check that fiona reads integers
-            import fiona
-            with fiona.open(outshp) as src:
-                meta = src.meta
-                assert 'int' in meta['schema']['properties']['row']
-                assert 'int' in meta['schema']['properties']['column']
-        except:
-            pass
+    # check that pyshp reads integers
+    # this only check that row/column were recorded as "N"
+    # not how they will be cast by python or numpy
+    import shapefile as sf
+    sfobj = sf.Reader(outshp)
+    for f in sfobj.fields:
+        if f[0] == 'row' or f[0] == 'column':
+            assert f[1] == 'N'
+    recs = list(sfobj.records())
+    for r in recs[0]:
+        assert isinstance(r, int)
+
+    # check that row and column appear as integers in recarray
+    ra = shp2recarray(outshp)
+    assert np.issubdtype(ra.dtype['row'], np.integer)
+    assert np.issubdtype(ra.dtype['column'], np.integer)
+
+    try: # check that fiona reads integers
+        import fiona
+        with fiona.open(outshp) as src:
+            meta = src.meta
+            assert 'int' in meta['schema']['properties']['row']
+            assert 'int' in meta['schema']['properties']['column']
+    except:
+        pass
+
+
+def test_shapefile_polygon_closed():
+    import os
+    import flopy
+    try:
+        import shapefile
+    except:
+        return
+
+    xll, yll = 468970, 3478635
+    xur, yur = 681010, 3716462
+
+    spacing = 2000
+
+    ncol = int((xur - xll) / spacing)
+    nrow = int((yur - yll) / spacing)
+    print(nrow, ncol)
+
+    m = flopy.modflow.Modflow("test.nam", proj4_str="EPSG:32614", xll=xll,
+                              yll=yll)
+
+    flopy.modflow.ModflowDis(m, delr=spacing, delc=spacing, nrow=nrow,
+                             ncol=ncol)
+
+    shp_file = os.path.join(spth, "test_polygon.shp")
+    m.dis.export(shp_file)
+
+    shp = shapefile.Reader(shp_file)
+    for shape in shp.iterShapes():
+        if len(shape.points) != 5:
+            raise AssertionError("Shapefile polygon is not closed!")
 
 
 def test_export_array():
@@ -1331,4 +1360,5 @@ def test_export_contourf():
     #test_tricontour_NaN()
     #test_export_contourf()
     #test_sr()
+    # test_shapefile_polygon_closed()
     pass
diff --git a/flopy/discretization/structuredgrid.py b/flopy/discretization/structuredgrid.py
index fdab4e9551..524c98e057 100644
--- a/flopy/discretization/structuredgrid.py
+++ b/flopy/discretization/structuredgrid.py
@@ -386,11 +386,11 @@ def from_gridspec(cls, gridspec_file, lenuni=0):
     # Exporting
     def write_shapefile(self, filename='grid.shp', epsg=None, prj=None):
         """Write a shapefile of the grid with just the row and column attributes"""
-        from ..export.shapefile_utils import write_grid_shapefile2
+        from ..export.shapefile_utils import write_grid_shapefile
         if epsg is None and prj is None:
             epsg = self.epsg
-        write_grid_shapefile2(filename, self, array_dict={}, nan_val=-1.0e9,
-                              epsg=epsg, prj=prj)
+        write_grid_shapefile(filename, self, array_dict={}, nan_val=-1.0e9,
+                             epsg=epsg, prj=prj)
 
 
 if __name__ == "__main__":
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
index 9decf5a6b3..529f3b1f79 100755
--- a/flopy/export/shapefile_utils.py
+++ b/flopy/export/shapefile_utils.py
@@ -25,20 +25,6 @@ def import_shapefile():
                         "importing shapefile - try pip install pyshp")
 
 
-def shapefile_version(sf):
-    """
-    Return the shapefile major version number
-    Parameters
-    ----------
-    sf : shapefile package
-
-    Returns
-    -------
-    int
-    """
-    return int(sf.__version__.split('.')[0])
-
-
 def write_gridlines_shapefile(filename, mg):
     """
     Write a polyline shapefile of the grid lines - a lightweight alternative
@@ -56,11 +42,7 @@ def write_gridlines_shapefile(filename, mg):
 
     """
     shapefile = import_shapefile()
-    sfv = shapefile_version(shapefile)
-    if sfv < 2:
-        wr = shapefile.Writer(shapeType=shapefile.POLYLINE)
-    else:
-        wr = shapefile.Writer(filename, shapeType=shapefile.POLYLINE)
+    wr = shapefile.Writer(filename, shapeType=shapefile.POLYLINE)
     wr.field("number", "N", 18, 0)
     if isinstance(mg, SpatialReference):
         grid_lines = mg.get_grid_lines()
@@ -73,116 +55,38 @@ def write_gridlines_shapefile(filename, mg):
     for i, line in enumerate(grid_lines):
         wr.poly([line])
         wr.record(i)
-    if sfv < 2:
-        wr.save(filename)
-    else:
-        wr.close()
+
+    wr.close()
     return
 
 
-def write_grid_shapefile(filename, mg, array_dict, nan_val=None):  # -1.0e9):
+def write_grid_shapefile(filename, mg, array_dict, nan_val=np.nan,  # -1.0e9,
+                         epsg=None, prj=None):
     """
-    Write a grid shapefile array_dict attributes.
+    Method to write a shapefile of gridded input data
 
     Parameters
     ----------
-    filename : string
-        name of the shapefile to write
-    mg : model grid instance
-        object for model grid
+    filename : str
+        shapefile file name path
+    mg : flopy.discretization.Grid object
+        flopy model grid
     array_dict : dict
-       Dictionary of name and 2D array pairs.  Additional 2D arrays to add as
-       attributes to the grid shapefile.
+        dictionary of model input arrays
+    nan_val : float
+        value to fill nans
+    epsg : str, int
+        epsg code
+    prj : str
+        projection file name path
 
     Returns
     -------
     None
 
     """
-
     shapefile = import_shapefile()
-    sfv = shapefile_version(shapefile)
-    if sfv < 2:
-        wr = shapefile.Writer(shapeType=shapefile.POLYGON)
-    else:
-        wr = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
-    if isinstance(mg, SpatialReference):
-        warnings.warn(
-            "SpatialReference has been deprecated. Use StructuredGrid"
-            " instead.",
-            category=DeprecationWarning)
-        wr.field("row", "N", 10, 0)
-        wr.field("column", "N", 10, 0)
-    elif mg.grid_type == 'structured':
-        wr.field("row", "N", 10, 0)
-        wr.field("column", "N", 10, 0)
-    elif mg.grid_type == 'vertex':
-        wr.field("node", "N", 10, 0)
-    else:
-        raise Exception('Grid type {} not supported.'.format(mg.grid_type))
-
-    arrays = []
-    names = list(array_dict.keys())
-    names.sort()
-    # for name,array in array_dict.items():
-    for name in names:
-        array = array_dict[name]
-        if array.ndim == 3:
-            assert array.shape[0] == 1
-            array = array[0, :, :]
-        assert array.shape == (mg.nrow, mg.ncol)
-        if array.dtype in [np.float, np.float32, np.float64]:
-            array[np.where(np.isnan(array))] = nan_val
-        else:
-            j = 2
-        # if array.dtype in [np.int, np.int32, np.int64]:
-        #    wr.field(name, "N", 18, 0)
-        # else:
-        #    wr.field(name, "F", 18, 12)
-        wr.field(name, *get_pyshp_field_info(array.dtype.name))
-        arrays.append(array)
-
-    if isinstance(mg, SpatialReference) or mg.grid_type == 'structured':
-        for i in range(mg.nrow):
-            for j in range(mg.ncol):
-                try:
-                    pts = mg.get_cell_vertices(i, j)
-                except AttributeError:
-                    # support old style SR object
-                    pts = mg.get_vertices(i, j)
-
-                wr.poly([pts])
-                rec = [i + 1, j + 1]
-                for array in arrays:
-                    rec.append(array[i, j])
-                wr.record(*rec)
-    elif mg.grid_type == 'vertex':
-        for i in range(mg.ncpl):
-            pts = mg.get_cell_vertices(i)
-
-            wr.poly([pts])
-            rec = [i + 1]
-            for array in arrays:
-                rec.append(array[i])
-            wr.record(*rec)
-
-    # close or write the file
-    if sfv < 2:
-        wr.save(filename)
-    else:
-        wr.close()
-    print('wrote {}'.format(filename))
-    return
-
-
-def write_grid_shapefile2(filename, mg, array_dict, nan_val=np.nan,  # -1.0e9,
-                          epsg=None, prj=None):
-    shapefile = import_shapefile()
-    sfv = shapefile_version(shapefile)
-    if sfv < 2:
-        w = shapefile.Writer(shapeType=shapefile.POLYGON)
-    else:
-        w = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
+    w = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
     w.autoBalance = 1
 
     if isinstance(mg, SpatialReference):
@@ -201,55 +105,55 @@ def write_grid_shapefile2(filename, mg, array_dict, nan_val=np.nan,  # -1.0e9,
     else:
         raise Exception('Grid type {} not supported.'.format(mg.grid_type))
 
-    # set up the attribute fields
+    # set up the attribute fields and arrays of attributes
     if isinstance(mg, SpatialReference) or mg.grid_type == 'structured':
         names = ['node', 'row', 'column'] + list(array_dict.keys())
-        names = enforce_10ch_limit(names)
         dtypes = [('node', np.dtype('int')),
                   ('row', np.dtype('int')),
                   ('column', np.dtype('int'))] + \
-                 [(enforce_10ch_limit([name])[0], arr.dtype)
-                  for name, arr in array_dict.items()]
-    elif mg.grid_type == 'vertex':
-        names = ['node'] + list(array_dict.keys())
-        names = enforce_10ch_limit(names)
-        dtypes = [('node', np.dtype('int'))] + \
-                 [(enforce_10ch_limit([name])[0], arr.dtype)
-                  for name, arr in array_dict.items()]
-
-    fieldinfo = {name: get_pyshp_field_info(dtype.name) for name, dtype in
-                 dtypes}
-    for n in names:
-        w.field(n, *fieldinfo[n])
-
-    if isinstance(mg, SpatialReference) or mg.grid_type == 'structured':
-        # set-up array of attributes of shape ncells x nattributes
+                 [(enforce_10ch_limit([name])[0], array_dict[name].dtype)
+                  for name in names[3:]]
         node = list(range(1, mg.ncol * mg.nrow + 1))
         col = list(range(1, mg.ncol + 1)) * mg.nrow
         row = sorted(list(range(1, mg.nrow + 1)) * mg.ncol)
         at = np.vstack(
             [node, row, col] +
-            [arr.ravel() for arr in array_dict.values()]).transpose()
-        if at.dtype in [np.float, np.float32, np.float64]:
-            at[np.isnan(at)] = nan_val
+            [array_dict[name].ravel() for name in names[3:]]).transpose()
+
+        names = enforce_10ch_limit(names)
+
     elif mg.grid_type == 'vertex':
-        # set-up array of attributes of shape ncells x nattributes
+        names = ['node'] + list(array_dict.keys())
+        dtypes = [('node', np.dtype('int'))] + \
+                 [(enforce_10ch_limit([name])[0], array_dict[name].dtype)
+                  for name in names[1:]]
         node = list(range(1, mg.ncpl + 1))
         at = np.vstack(
             [node] +
-            [arr.ravel() for arr in array_dict.values()]).transpose()
+            [array_dict[name].ravel() for name in names[1:]]).transpose()
+
+        names = enforce_10ch_limit(names)
+
+    # flag nan values and explicitly set the dtypes
     if at.dtype in [np.float, np.float32, np.float64]:
         at[np.isnan(at)] = nan_val
+    at = np.array([tuple(i) for i in at], dtype=dtypes)
+
+    # write field information
+    fieldinfo = {name: get_pyshp_field_info(dtype.name) for name, dtype in
+                 dtypes}
+    for n in names:
+        w.field(n, *fieldinfo[n])
 
     for i, r in enumerate(at):
+        # check if polygon is closed, if not close polygon for QGIS
+        if verts[i][-1] != verts[i][0]:
+            verts[i] = verts[i] + [verts[i][0]]
         w.poly([verts[i]])
         w.record(*r)
 
     # close
-    if sfv < 2:
-        w.save(filename)
-    else:
-        w.close()
+    w.close()
     print('wrote {}'.format(filename))
     # write the projection file
     write_prj(filename, mg, epsg, prj)
@@ -407,7 +311,7 @@ def model_attributes_to_shapefile(filename, ml, package_names=None,
                                 array_dict[name] = arr
 
     # write data arrays to a shapefile
-    write_grid_shapefile2(filename, grid, array_dict)
+    write_grid_shapefile(filename, grid, array_dict)
     epsg = kwargs.get('epsg', None)
     prj = kwargs.get('prj', None)
     write_prj(filename, grid, epsg, prj)
@@ -592,13 +496,9 @@ def recarray2shp(recarray, geoms, shpname='recarray.shp', mg=None,
         except AttributeError:
             continue
 
-    # set up for pyshp 1 or 2
+    # set up for pyshp 2
     shapefile = import_shapefile()
-    sfv = shapefile_version(shapefile)
-    if sfv < 2:
-        w = shapefile.Writer(shapeType=geomtype)
-    else:
-        w = shapefile.Writer(shpname, shapeType=geomtype)
+    w = shapefile.Writer(shpname, shapeType=geomtype)
     w.autoBalance = 1
 
     # set up the attribute fields
@@ -622,23 +522,12 @@ def recarray2shp(recarray, geoms, shpname='recarray.shp', mg=None,
     elif geomtype == shapefile.POINT:
         # pyshp version 2.x w.point() method can only take x and y
         # code will need to be refactored in order to write POINTZ
-        # shapes with the z attribute.  The pyshp version 1.x
-        # method w.point() took a z attribute, but only wrote it if
-        # the shapeType was shapefile.POINTZ, which it is not for
-        # flopy, even if the point is 3D.
-        if sfv < 2:
-            for i, r in enumerate(ralist):
-                w.point(*geoms[i].pyshp_parts)
-                w.record(*r)
-        else:
-            for i, r in enumerate(ralist):
-                w.point(*geoms[i].pyshp_parts[:2])
-                w.record(*r)
+        # shapes with the z attribute.
+        for i, r in enumerate(ralist):
+            w.point(*geoms[i].pyshp_parts[:2])
+            w.record(*r)
 
-    if sfv < 2:
-        w.save(shpname)
-    else:
-        w.close()
+    w.close()
     write_prj(shpname, mg, epsg, prj)
     print('wrote {}'.format(shpname))
     return
diff --git a/flopy/export/utils.py b/flopy/export/utils.py
index b1afc94b97..193abffd83 100644
--- a/flopy/export/utils.py
+++ b/flopy/export/utils.py
@@ -576,7 +576,7 @@ def mflist_export(f, mfl, **kwargs):
                         n = shapefile_utils.shape_attr_name(name, length=4)
                         aname = "{}{}{}".format(n, k + 1, int(kk) + 1)
                         array_dict[aname] = array[k]
-            shapefile_utils.write_grid_shapefile2(f, modelgrid, array_dict)
+            shapefile_utils.write_grid_shapefile(f, modelgrid, array_dict)
         else:
             from ..export.shapefile_utils import recarray2shp
             from ..utils.geometry import Polygon
@@ -693,7 +693,7 @@ def transient2d_export(f, t2d, **kwargs):
             name = '{}_{}'.format(shapefile_utils.shape_attr_name(u2d.name),
                                   kper + 1)
             array_dict[name] = u2d.array
-        shapefile_utils.write_grid_shapefile2(f, modelgrid, array_dict)
+        shapefile_utils.write_grid_shapefile(f, modelgrid, array_dict)
 
     elif isinstance(f, NetCdf) or isinstance(f, dict):
         # mask the array is defined by any row col with at lease
@@ -810,7 +810,7 @@ def array3d_export(f, u3d, **kwargs):
             name = '{}_{}'.format(
                 shapefile_utils.shape_attr_name(u2d.name), ilay + 1)
             array_dict[name] = u2d.array
-        shapefile_utils.write_grid_shapefile2(f, modelgrid, array_dict)
+        shapefile_utils.write_grid_shapefile(f, modelgrid, array_dict)
 
     elif isinstance(f, NetCdf) or isinstance(f, dict):
         var_name = u3d.name
@@ -940,7 +940,7 @@ def array2d_export(f, u2d, **kwargs):
 
     if isinstance(f, str) and f.lower().endswith(".shp"):
         name = shapefile_utils.shape_attr_name(u2d.name, keep_layer=True)
-        shapefile_utils.write_grid_shapefile2(f, modelgrid,
+        shapefile_utils.write_grid_shapefile(f, modelgrid,
                                              {name: u2d.array})
         return
 
@@ -1144,14 +1144,14 @@ def export_array(modelgrid, filename, a, nodata=-9999,
         print('wrote {}'.format(filename))
 
     elif filename.lower().endswith(".shp"):
-        from ..export.shapefile_utils import write_grid_shapefile2
+        from ..export.shapefile_utils import write_grid_shapefile
         epsg = kwargs.get('epsg', None)
         prj = kwargs.get('prj', None)
         if epsg is None and prj is None:
             epsg = modelgrid.epsg
-        write_grid_shapefile2(filename, modelgrid, array_dict={fieldname: a},
-                              nan_val=nodata,
-                              epsg=epsg, prj=prj)
+        write_grid_shapefile(filename, modelgrid, array_dict={fieldname: a},
+                             nan_val=nodata,
+                             epsg=epsg, prj=prj)
 
 
 def export_contours(modelgrid, filename, contours,
diff --git a/flopy/utils/datafile.py b/flopy/utils/datafile.py
index 3dfea109c6..a8e778ad71 100755
--- a/flopy/utils/datafile.py
+++ b/flopy/utils/datafile.py
@@ -215,8 +215,8 @@ def to_shapefile(self, filename, kstpkper=None, totim=None, mflay=None,
                 name = attrib_name + '{}'.format(k)
                 attrib_dict[name] = plotarray[k]
 
-        from ..export.shapefile_utils import write_grid_shapefile2
-        write_grid_shapefile2(filename, self.mg, attrib_dict)
+        from ..export.shapefile_utils import write_grid_shapefile
+        write_grid_shapefile(filename, self.mg, attrib_dict)
 
     def plot(self, axes=None, kstpkper=None, totim=None, mflay=None,
              filename_base=None, **kwargs):
diff --git a/flopy/utils/geometry.py b/flopy/utils/geometry.py
index 84e1c570c3..570f192ee1 100644
--- a/flopy/utils/geometry.py
+++ b/flopy/utils/geometry.py
@@ -85,14 +85,12 @@ def geojson(self):
 
     @property
     def pyshp_parts(self):
-        from ..export.shapefile_utils import (import_shapefile,
-                                              shapefile_version)
+        from ..export.shapefile_utils import import_shapefile
 
         # exterior ring must be clockwise (negative area)
         # interiors rings must be counter-clockwise (positive area)
 
         shapefile = import_shapefile()
-        sfv = shapefile_version(shapefile)
 
         exterior = list(self.exterior)
         if shapefile.signed_area(exterior) > 0:
@@ -105,12 +103,9 @@ def pyshp_parts(self):
                 il.reverse()
             interiors.append(il)
 
-        if sfv < 2:
-            result = [exterior + [i for i in interiors]]
-        else:
-            result = [exterior]
-            for i in interiors:
-                result.append(i)
+        result = [exterior]
+        for i in interiors:
+            result.append(i)
         return result
 
     @property
diff --git a/flopy/utils/gridgen.py b/flopy/utils/gridgen.py
index 2d2c986c33..40ffe0c6c4 100644
--- a/flopy/utils/gridgen.py
+++ b/flopy/utils/gridgen.py
@@ -9,11 +9,9 @@
 from .util_array import Util2d  # read1d,
 from ..export.shapefile_utils import shp2recarray
 from ..mbase import which
-from ..export.shapefile_utils import import_shapefile, shapefile_version
+from ..export.shapefile_utils import import_shapefile
 
 shapefile = import_shapefile()
-sfv = shapefile_version(shapefile)
-
 
 # todo
 # creation of line and polygon shapefiles from features (holes!)
@@ -60,39 +58,27 @@ def features_to_shapefile(features, featuretype, filename):
         raise Exception('Unrecognized feature type: {}'.format(featuretype))
 
     if featuretype.lower() == 'line':
-        if sfv < 2:
-            wr = shapefile.Writer(shapeType=shapefile.POLYLINE)
-        else:
-            wr = shapefile.Writer(filename, shapeType=shapefile.POLYLINE)
+        wr = shapefile.Writer(filename, shapeType=shapefile.POLYLINE)
         wr.field("SHAPEID", "N", 20, 0)
         for i, line in enumerate(features):
             wr.line(line)
             wr.record(i)
 
     elif featuretype.lower() == 'point':
-        if sfv < 2:
-            wr = shapefile.Writer(shapeType=shapefile.POINT)
-        else:
-            wr = shapefile.Writer(filename, shapeType=shapefile.POINT)
+        wr = shapefile.Writer(filename, shapeType=shapefile.POINT)
         wr.field("SHAPEID", "N", 20, 0)
         for i, point in enumerate(features):
             wr.point(point[0], point[1])
             wr.record(i)
 
     elif featuretype.lower() == 'polygon':
-        if sfv < 2:
-            wr = shapefile.Writer(shapeType=shapefile.POLYGON)
-        else:
-            wr = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
+        wr = shapefile.Writer(filename, shapeType=shapefile.POLYGON)
         wr.field("SHAPEID", "N", 20, 0)
         for i, polygon in enumerate(features):
             wr.poly(polygon)
             wr.record(i)
 
-    if sfv < 2:
-        wr.save(filename)
-    else:
-        wr.close()
+    wr.close()
     return
 
 
diff --git a/flopy/utils/reference.py b/flopy/utils/reference.py
index 50cd8251d7..ce10e0e090 100755
--- a/flopy/utils/reference.py
+++ b/flopy/utils/reference.py
@@ -935,11 +935,11 @@ def write_gridSpec(self, filename):
 
     def write_shapefile(self, filename='grid.shp', epsg=None, prj=None):
         """Write a shapefile of the grid with just the row and column attributes"""
-        from ..export.shapefile_utils import write_grid_shapefile2
+        from ..export.shapefile_utils import write_grid_shapefile
         if epsg is None and prj is None:
             epsg = self.epsg
-        write_grid_shapefile2(filename, self, array_dict={}, nan_val=-1.0e9,
-                              epsg=epsg, prj=prj)
+        write_grid_shapefile(filename, self, array_dict={}, nan_val=-1.0e9,
+                             epsg=epsg, prj=prj)
 
     def get_vertices(self, i, j):
         """Get vertices for a single cell or sequence if i, j locations."""
@@ -1145,14 +1145,14 @@ def export_array(self, filename, a, nodata=-9999,
             print('wrote {}'.format(filename))
 
         elif filename.lower().endswith(".shp"):
-            from ..export.shapefile_utils import write_grid_shapefile2
+            from ..export.shapefile_utils import write_grid_shapefile
             epsg = kwargs.get('epsg', None)
             prj = kwargs.get('prj', None)
             if epsg is None and prj is None:
                 epsg = self.epsg
-            write_grid_shapefile2(filename, self, array_dict={fieldname: a},
-                                  nan_val=nodata,
-                                  epsg=epsg, prj=prj)
+            write_grid_shapefile(filename, self, array_dict={fieldname: a},
+                                 nan_val=nodata,
+                                 epsg=epsg, prj=prj)
 
     def export_contours(self, filename, contours,
                         fieldname='level', epsg=None, prj=None,