From c97c6c83419eb5b0fe7f33f4edf3fcbc4f6c85de Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Fri, 21 Jun 2019 13:09:29 -0500 Subject: [PATCH 1/2] feat(export_contourf): add capability to export filled contours to a shapefile Required changes to flopy.utils.geometry so that the Polygon.pyshp_prts handled holes correctly for multiple shapefile versions --- autotest/t007_test.py | 20 ++++++++- flopy/export/utils.py | 95 ++++++++++++++++++++++++++++++++++++++++- flopy/utils/geometry.py | 27 +++++++++++- 3 files changed, 137 insertions(+), 5 deletions(-) diff --git a/autotest/t007_test.py b/autotest/t007_test.py index 10e6caaec0..84537f9988 100644 --- a/autotest/t007_test.py +++ b/autotest/t007_test.py @@ -1217,6 +1217,21 @@ def test_export_array_contours(): return +def test_export_contourf(): + try: + import shapely + except: + return + import matplotlib.pyplot as plt + from flopy.export.utils import export_contourf + filename = os.path.join(spth, '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' + return + + if __name__ == '__main__': #test_shapefile() # test_shapefile_ibound() @@ -1248,6 +1263,7 @@ def test_export_array_contours(): #test_wkt_parse() #test_get_rc_from_node_coordinates() # test_export_array() - test_export_array_contours() - test_tricontour_NaN() + #test_export_array_contours() + #test_tricontour_NaN() + test_export_contourf() pass diff --git a/flopy/export/utils.py b/flopy/export/utils.py index 80af2aa469..ccc4bb86ba 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -1080,8 +1080,8 @@ def export_contours(modelgrid, filename, contours, df : dataframe of shapefile contents """ - from flopy.utils.geometry import LineString - from flopy.export.shapefile_utils import recarray2shp + from ..utils.geometry import LineString + from .shapefile_utils import recarray2shp if not isinstance(contours, list): contours = [contours] @@ -1105,6 +1105,97 @@ def export_contours(modelgrid, filename, contours, dtype=[(fieldname, float)]).view(np.recarray) recarray2shp(ra, geoms, filename, epsg=epsg, prj=prj, **kwargs) + return + + +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. + + Parameters + ---------- + filename : str + name of output shapefile (e.g. myshp.shp) + contours : matplotlib.contour.QuadContourSet or list of them + (object returned by matplotlib.pyplot.contourf) + fieldname : str + Name of shapefile attribute field to contain the contour level. The + fieldname column in the attribute table will contain the lower end of + the range represented by the polygon. Default is 'level'. + epsg : int + EPSG code. See https://www.epsg-registry.org/ or spatialreference.org + prj : str + Existing projection file to be used with new shapefile. + + **kwargs : keyword arguments to flopy.export.shapefile_utils.recarray2shp + + Returns + ------- + None + + Examples + -------- + >>> import flopy + >>> import matplotlib.pyplot as plt + >>> from flopy.export.utils import export_contourf + >>> a = np.random.random((10, 10)) + >>> cs = plt.contourf(a) + >>> export_contourf('myfilledcontours.shp', cs) + + """ + + try: + from shapely import geometry + except: + raise ImportError('export_contourf requires python shapely package') + + from ..utils.geometry import Polygon + from .shapefile_utils import recarray2shp + + shapelygeoms = [] + level = [] + + if not isinstance(contours, list): + contours = [contours] + + for c in contours: + levels = c.levels + for i, 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" + 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)]) + 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[i]) + + 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] + + print('Writing {} polygons'.format(len(level))) + + # Create recarray + ra = np.array(level, dtype=[(fieldname, float)]).view(np.recarray) + + recarray2shp(ra, geoms, filename, epsg=epsg, prj=prj, **kwargs) + return def export_array_contours(modelgrid, filename, a, diff --git a/flopy/utils/geometry.py b/flopy/utils/geometry.py index 366e6c99ec..a4e50e0eb5 100644 --- a/flopy/utils/geometry.py +++ b/flopy/utils/geometry.py @@ -85,7 +85,32 @@ def geojson(self): @property def pyshp_parts(self): - return [list(self.exterior) + [list(i) for i in self.interiors]] + from ..export.shapefile_utils import import_shapefile, shapefile_version + + # 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: + exterior.reverse() + + interiors = [] + for i in self.interiors: + il = list(i) + if shapefile.signed_area(il) < 0: + 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) + return result @property def patch(self): From b99c872f9df4cd6c2510564a5efad5e0146a241d Mon Sep 17 00:00:00 2001 From: Christian Langevin Date: Fri, 21 Jun 2019 13:43:25 -0500 Subject: [PATCH 2/2] codacy fixes --- flopy/export/utils.py | 3 ++- flopy/utils/geometry.py | 3 ++- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/flopy/export/utils.py b/flopy/export/utils.py index ccc4bb86ba..030be51561 100644 --- a/flopy/export/utils.py +++ b/flopy/export/utils.py @@ -1171,7 +1171,8 @@ def export_contourf(filename, contours, fieldname='level', epsg=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)]) + new_shape = geometry.Polygon([(i[0], i[1]) + for i in zip(x, y)]) if ncp == 0: poly = new_shape else: diff --git a/flopy/utils/geometry.py b/flopy/utils/geometry.py index a4e50e0eb5..84e1c570c3 100644 --- a/flopy/utils/geometry.py +++ b/flopy/utils/geometry.py @@ -85,7 +85,8 @@ def geojson(self): @property def pyshp_parts(self): - from ..export.shapefile_utils import import_shapefile, shapefile_version + from ..export.shapefile_utils import (import_shapefile, + shapefile_version) # exterior ring must be clockwise (negative area) # interiors rings must be counter-clockwise (positive area)