feat(geospatial_util): geospatial consolidation utilities (#1002)

* *Added GeoSpatialUtil for geospatial conversions *Added GeoSpatialCollection for batch conversion of "Collection" objects *Updated geometry.py and added additional geometry types to geometry.py *updates to gridgen, shapefile_utils, triangle, rasters, and gridintersect modules that allow the user to supply any supported shape object to these geospatial methods. Supported shape objects are, geojson, shapely, shapefile, flopy.utils.geometry, and lists of vertices. *Added t074_test_geospatial_util.py issue #772 * remove intersect_linestring, intersect_point, and intersect_polygon from GridIntersect. *Refactored t065_test_gridintersect.py for intersect() method *Updated flopy3_grid_intersection_demo.ipynb for intersect() method
modflowpy · Oct 19, 2020 · 9ad6732 · 9ad6732
1 parent 698b4ec
commit 9ad6732
Show file tree

Hide file tree

Showing 10 changed files with 1,677 additions and 320 deletions.
diff --git a/autotest/t065_test_gridintersect.py b/autotest/t065_test_gridintersect.py
diff --git a/autotest/t074_test_geospatial_util.py b/autotest/t074_test_geospatial_util.py
diff --git a/examples/Notebooks/flopy3_grid_intersection_demo.ipynb b/examples/Notebooks/flopy3_grid_intersection_demo.ipynb
diff --git a/flopy/export/shapefile_utils.py b/flopy/export/shapefile_utils.py
@@ -90,7 +90,12 @@ def write_gridlines_shapefile(filename, mg):
 
 
 def write_grid_shapefile(
-    filename, mg, array_dict, nan_val=np.nan, epsg=None, prj=None  # -1.0e9,
+    filename,
+    mg,
+    array_dict,
+    nan_val=np.nan,
+    epsg=None,
+    prj=None,  # -1.0e9,
 ):
     """
     Method to write a shapefile of gridded input data
@@ -357,7 +362,8 @@ def model_attributes_to_shapefile(
                             for ilay in range(a.model.modelgrid.nlay):
                                 u2d = a[ilay]
                                 name = "{}_{}".format(
-                                    shape_attr_name(u2d.name), ilay + 1
+                                    shape_attr_name(u2d.name),
+                                    ilay + 1,
                                 )
                                 arr = u2d.array
                                 assert arr.shape == horz_shape
@@ -463,7 +469,12 @@ def get_pyshp_field_info(dtypename):
 
 def get_pyshp_field_dtypes(code):
     """Returns a numpy dtype for a pyshp field type."""
-    dtypes = {"N": np.int, "F": np.float, "L": np.bool, "C": np.object}
+    dtypes = {
+        "N": np.int,
+        "F": np.float,
+        "L": np.bool,
+        "C": np.object,
+    }
     return dtypes.get(code, np.object)
 
 
@@ -480,7 +491,7 @@ def shp2recarray(shpname):
     recarray : np.recarray
 
     """
-    from ..utils.geometry import shape
+    from ..utils.geospatial_utils import GeoSpatialCollection
 
     sf = import_shapefile(check_version=False)
 
@@ -489,7 +500,7 @@ def shp2recarray(shpname):
         (str(f[0]), get_pyshp_field_dtypes(f[1])) for f in sfobj.fields[1:]
     ]
 
-    geoms = [shape(s) for s in sfobj.iterShapes()]
+    geoms = GeoSpatialCollection(sfobj).flopy_geometry
     records = [
         tuple(r) + (geoms[i],) for i, r in enumerate(sfobj.iterRecords())
     ]
@@ -510,14 +521,18 @@ def recarray2shp(
 ):
     """
     Write a numpy record array to a shapefile, using a corresponding
-    list of geometries.
+    list of geometries. Method supports list of flopy geometry objects,
+    flopy Collection object, shapely Collection object, and geojson
+    Geometry Collection objects
 
     Parameters
     ----------
     recarray : np.recarray
         Numpy record array with attribute information that will go in the
         shapefile
-    geoms : list of flopy.utils.geometry objects
+    geoms : list of flopy.utils.geometry, shapely geometry collection,
+            flopy geometry collection, shapefile.Shapes,
+            list of shapefile.Shape objects, or geojson geometry collection
         The number of geometries in geoms must equal the number of records in
         recarray.
     shpname : str
@@ -536,6 +551,7 @@ def recarray2shp(
     subsequent use. See flopy.reference for more details.
 
     """
+    from ..utils.geospatial_utils import GeoSpatialCollection
 
     if len(recarray) != len(geoms):
         raise IndexError(
@@ -546,6 +562,9 @@ def recarray2shp(
         raise Exception("Recarray is empty")
 
     geomtype = None
+
+    geoms = GeoSpatialCollection(geoms).flopy_geometry
+
     for g in geoms:
         try:
             geomtype = g.shapeType
@@ -722,7 +741,10 @@ def grid_mapping_attribs(self):
         if self.wktstr is not None:
             sp = [
                 p
-                for p in [self.standard_parallel_1, self.standard_parallel_2]
+                for p in [
+                    self.standard_parallel_1,
+                    self.standard_parallel_2,
+                ]
                 if p is not None
             ]
             sp = sp if len(sp) > 0 else None
@@ -797,7 +819,12 @@ def _getvalue(self, k):
             end = s[strt:].find("]") + strt
             try:
                 return float(self.wktstr[strt:end].split(",")[1])
-            except (IndexError, TypeError, ValueError, AttributeError):
+            except (
+                IndexError,
+                TypeError,
+                ValueError,
+                AttributeError,
+            ):
                 pass
 
     def _getgcsparam(self, txt):