diff --git a/autotest/t007_test.py b/autotest/t007_test.py index edcd02f85..8f2b572de 100644 --- a/autotest/t007_test.py +++ b/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(): diff --git a/flopy/export/utils.py b/flopy/export/utils.py index af6d560d3..2421ba178 100644 --- a/flopy/export/utils.py +++ b/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") diff --git a/flopy/plot/crosssection.py b/flopy/plot/crosssection.py index 404968494..73f6e4189 100644 --- a/flopy/plot/crosssection.py +++ b/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) diff --git a/flopy/plot/map.py b/flopy/plot/map.py index fafaf6eb4..158b1f996 100644 --- a/flopy/plot/map.py +++ b/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)