fix(gis.py): Ensure valid geometries for intersection functions

docs/source/api/index.rst

@@ -4,6 +4,7 @@ Modules
 .. toctree::
 
    Flows Module <sfrmaker.flows>
+   GIS Module <sfrmaker.gis>
    Grid Module <sfrmaker.grid>
    Lines Module <sfrmaker.lines>
    MODFLOW-2005 to 6 Module <sfrmaker.mf5to6>

docs/source/conf.py

@@ -109,9 +109,6 @@
 # a list of builtin themes.
 #
 html_theme = 'sphinx_rtd_theme'
-import sphinx_rtd_theme
-
-html_theme_path = [sphinx_rtd_theme.get_html_theme_path()]
 
 # Theme options are theme-specific and customize the look and feel of a theme
 # further.  For a list of options available for each theme, see the
@@ -206,5 +203,6 @@
     'flopy': ('https://flopy.readthedocs.io/en/latest/', None),
     'rasterstats': ('https://pythonhosted.org/rasterstats/', None),
     'shapely': ('https://shapely.readthedocs.io/en/latest/', None),
-    'pyproj': ('http://pyproj4.github.io/pyproj/stable/', None)
+    'pyproj': ('http://pyproj4.github.io/pyproj/stable/', None),
+    'rtree': ('https://rtree.readthedocs.io/en/latest/', None)
 }

sfrmaker/gis.py

@@ -10,6 +10,7 @@
 import pyproj
 from shapely.geometry import shape, Polygon, box
 from shapely.ops import unary_union
+from shapely.validation import make_valid
 import gisutils
 from gisutils import df2shp, shp2df, project, get_shapefile_crs, get_authority_crs
 import sfrmaker
@@ -49,7 +50,8 @@ def get_crs(prjfile=None, crs=None, **kwargs):
 
 
 def build_rtree_index(geom):
-    """Builds an rtree index. Useful for multiple intersections with same index.
+    """Builds an :class:`rtree.index.Index` (spatial index) object. 
+    Useful for multiple intersections with same index.
 
     Parameters
     ==========
@@ -98,33 +100,41 @@ def intersect_rtree(geom1, geom2, index=None):
     """Intersect features in geom1 with those in geom2. For each feature in geom2, return a list of
      the indices of the intersecting features in geom1.
 
-    Parameters:
+    Parameters
     ----------
     geom1 : list
         list of shapely geometry objects
     geom2 : list
-        list of shapely polygon objects to be intersected with features in geom1
-    index :
-        use an index that has already been created
+        list of shapely geometry objects to be intersected with features in geom1
+    index : rtree spatial index
+        Option to use an r-tree spatial index that has already been created for geom1.
+        The :func:`sfrmaker.gis.build_rtree_index` function can be used to
+        create a spatial index for a list of shapely geometriy objects.
+        by default, None.
 
-    Returns:
+    Returns
     -------
     A list of the same length as geom2; containing for each feature in geom2,
     a list of indicies of intersecting geometries in geom1.
     """
+
+    #make certain that the objects in geom1 and geom2 are all considered valid shapely geometries using "make_valid()"
+    for i in range(len(geom1)): geom1[i] = make_valid(geom1[i])
+    for i in range(len(geom2)): geom2[i] = make_valid(geom2[i])
+
     if index is None:
         idx = build_rtree_index(geom1)
     else:
         idx = index
     isfr = []
     print('\nIntersecting {} features...'.format(len(geom2)))
     ta = time.time()
-    for pind, poly in enumerate(geom2):
+    for pind, geom in enumerate(geom2):
         print('\r{}'.format(pind + 1), end='')
-        # test for intersection with bounding box of each polygon feature in geom2 using spatial index
-        inds = [i for i in idx.intersection(poly.bounds)]
-        # test each feature inside the bounding box for intersection with the polygon geometry
-        inds = [i for i in inds if geom1[i].intersects(poly)]
+        # test for intersection with bounding box of each feature in geom2 using spatial index
+        inds = [i for i in idx.intersection(geom.bounds)]
+        # test each feature inside the bounding box for intersection with the geometry
+        inds = [i for i in inds if geom1[i].intersects(geom)]
         isfr.append(inds)
     print("\nfinished in {:.2f}s".format(time.time() - ta))
     return isfr
@@ -134,19 +144,23 @@ def intersect(geom1, geom2):
     """Same as intersect_rtree, except without spatial indexing. Fine for smaller datasets,
     but scales by 10^4 with the side of the problem domain.
 
-    Parameters:
+    Parameters
     ----------
     geom1 : list
         list of shapely geometry objects
     geom2 : list
-        list of shapely polygon objects to be intersected with features in geom1
+        list of shapely geometry objects to be intersected with features in geom1
 
-    Returns:
+    Returns
     -------
     A list of the same length as geom2; containing for each feature in geom2,
     a list of indicies of intersecting geometries in geom1.
     """
 
+    #make certain that the objects in geom1 and geom2 are all considered valid shapely geometries using "make_valid()"
+    for i in range(len(geom1)): geom1[i] = make_valid(geom1[i])
+    for i in range(len(geom2)): geom2[i] = make_valid(geom2[i])
+
     isfr = []
     ngeom1 = len(geom1)
     print('Intersecting {} features...'.format(len(geom2)))