update(contouring routines): updates for contour_array and export_con…

…tourf routines (#1353)

* added filled kwarg to contour_array for filled contours
* update PlotCrossSection.contour_array() to mask out edge triangles outside of cross-section boundaries
* update export_contourf, provide support for multipolygon contours and update support for holes
* update export_contourf, remove shapely dependency. Replaced with flopy geometry objects
* update t007_test.py, test_export_contourf()

autotest/t007_test.py

@@ -1731,25 +1731,40 @@ def test_export_array_contours():
 
 
 def test_export_contourf():
-    try:
-        import shapely
-    except:
-        return
     if import_shapefile() is None:
         return
-    import matplotlib.pyplot as plt
-
+    import shapefile
     from flopy.export.utils import export_contourf
+    import matplotlib.pyplot as plt
 
     ws_out = f"{base_dir}_shapefile_export_contourf"
+    filename = os.path.join(ws_out, "myfilledcontours.shp")
+
     test_setup = FlopyTestSetup(verbose=True, test_dirs=ws_out)
 
-    filename = os.path.join(ws_out, "myfilledcontours.shp")
-    a = np.random.random((10, 10))
-    cs = plt.contourf(a)
-    export_contourf(filename, cs)
-    assert os.path.isfile(filename), "did not create contourf shapefile"
+    model_ws = os.path.join("..", "examples", "data", "freyberg")
+    ml = flopy.modflow.Modflow.load("freyberg.nam", model_ws=model_ws)
+    hds_pth = os.path.join(model_ws, "freyberg.githds")
+    hds = flopy.utils.HeadFile(hds_pth)
+    head = hds.get_data()
+    levels = np.arange(10, 30, 0.5)
+
+    mapview = flopy.plot.PlotMapView(model=ml)
+    contour_set = mapview.contour_array(
+        head, masked_values=[999.0], levels=levels, filled=True
+    )
+
+    export_contourf(filename, contour_set)
     plt.close()
+    if not os.path.isfile(filename):
+        raise AssertionError("did not create contourf shapefile")
+
+    with shapefile.Reader(filename) as r:
+        shapes = r.shapes()
+        if len(shapes) != 65:
+            raise AssertionError(
+                "multipolygons were skipped in contourf routine"
+            )
 
 
 def main():

flopy/export/utils.py

@@ -1672,8 +1672,7 @@ def export_contourf(
     filename, contours, fieldname="level", epsg=None, prj=None, **kwargs
 ):
     """
-    Write matplotlib filled contours to shapefile.  This utility requires
-    that shapely is installed.
+    Write matplotlib filled contours to shapefile.
 
     Parameters
     ----------
@@ -1706,16 +1705,10 @@ def export_contourf(
     >>> export_contourf('myfilledcontours.shp', cs)
 
     """
-
-    shapely = import_optional_dependency(
-        "shapely", error_message="export_contourf requires shapely."
-    )
-    from shapely import geometry
-
-    from ..utils.geometry import Polygon
+    from ..utils.geometry import Polygon, is_clockwise
     from .shapefile_utils import recarray2shp
 
-    shapelygeoms = []
+    geoms = []
     level = []
 
     if not isinstance(contours, list):
@@ -1726,31 +1719,34 @@ def export_contourf(
         for idx, col in enumerate(c.collections):
             # Loop through all polygons that have the same intensity level
             for contour_path in col.get_paths():
-                # Create the polygon for this intensity level
-                # The first polygon in the path is the main one, the following
-                # ones are "holes"
+                # Create the polygon(s) for this intensity level
+                poly = None
                 for ncp, cp in enumerate(contour_path.to_polygons()):
                     x = cp[:, 0]
                     y = cp[:, 1]
-                    new_shape = geometry.Polygon(
-                        [(i[0], i[1]) for i in zip(x, y)]
-                    )
+                    verts = [(i[0], i[1]) for i in zip(x, y)]
+                    new_shape = Polygon(verts)
+
                     if ncp == 0:
                         poly = new_shape
                     else:
-                        # Remove the holes if there are any
-                        poly = poly.difference(new_shape)
-
-                # store shapely geometry object
-                shapelygeoms.append(poly)
-                level.append(levels[idx])
-
-    geoms = []
-    for shpgeom in shapelygeoms:
-        xa, ya = shpgeom.exterior.coords.xy
-        interiors = [s.coords for s in shpgeom.interiors]
-        pg = Polygon([(x, y) for x, y in zip(xa, ya)], interiors=interiors)
-        geoms += [pg]
+                        # check if this is a multipolygon by checking vertex
+                        # order.
+                        if is_clockwise(verts):
+                            # Clockwise is a hole, set to interiors
+                            if not poly.interiors:
+                                poly.interiors = [new_shape.exterior]
+                            else:
+                                poly.interiors.append(new_shape.exterior)
+                        else:
+                            geoms.append(poly)
+                            level.append(levels[idx])
+                            poly = new_shape
+
+                if poly is not None:
+                    # store geometry object
+                    geoms.append(poly)
+                    level.append(levels[idx])
 
     print(f"Writing {len(level)} polygons")
 

flopy/plot/crosssection.py

@@ -554,6 +554,7 @@ def contour_array(self, a, masked_values=None, head=None, **kwargs):
                     t = np.isclose(plotarray, mval)
                     ismasked += t
 
+        filled = kwargs.pop("filled", False)
         plot_triplot = kwargs.pop("plot_triplot", False)
 
         if "extent" in kwargs:
@@ -571,18 +572,33 @@ def contour_array(self, a, masked_values=None, head=None, **kwargs):
 
         if mplcontour:
             plotarray = np.ma.masked_array(plotarray, ismasked)
-            contour_set = ax.contour(xcenters, zcenters, plotarray, **kwargs)
+            if filled:
+                contour_set = ax.contourf(
+                    xcenters, zcenters, plotarray, **kwargs
+                )
+            else:
+                contour_set = ax.contour(
+                    xcenters, zcenters, plotarray, **kwargs
+                )
         else:
             triang = tri.Triangulation(xcenters, zcenters)
+            analyze = tri.TriAnalyzer(triang)
+            mask = analyze.get_flat_tri_mask(rescale=False)
 
             if ismasked is not None:
                 ismasked = ismasked.flatten()
-                mask = np.any(
+                mask2 = np.any(
                     np.where(ismasked[triang.triangles], True, False), axis=1
                 )
-                triang.set_mask(mask)
+                mask[mask2] = True
+
+            triang.set_mask(mask)
+
+            if filled:
+                contour_set = ax.tricontourf(triang, plotarray, **kwargs)
+            else:
+                contour_set = ax.tricontour(triang, plotarray, **kwargs)
 
-            contour_set = ax.tricontour(triang, plotarray, **kwargs)
             if plot_triplot:
                 ax.triplot(triang, color="black", marker="o", lw=0.75)
 

flopy/plot/map.py

@@ -205,9 +205,8 @@ def contour_array(self, a, masked_values=None, **kwargs):
             if "cmap" in kwargs.keys():
                 kwargs.pop("cmap")
 
-        plot_triplot = False
-        if "plot_triplot" in kwargs:
-            plot_triplot = kwargs.pop("plot_triplot")
+        filled = kwargs.pop("filled", False)
+        plot_triplot = kwargs.pop("plot_triplot", False)
 
         # Get vertices for the selected layer
         xcentergrid = self.mg.get_xcellcenters_for_layer(self.layer)
@@ -238,7 +237,10 @@ def contour_array(self, a, masked_values=None, **kwargs):
             )
             triang.set_mask(mask)
 
-        contour_set = ax.tricontour(triang, plotarray, **kwargs)
+        if filled:
+            contour_set = ax.tricontourf(triang, plotarray, **kwargs)
+        else:
+            contour_set = ax.tricontour(triang, plotarray, **kwargs)
 
         if plot_triplot:
             ax.triplot(triang, color="black", marker="o", lw=0.75)