As posted on flopy's temporary replacement repo by ebake310:
I am trying to intersect a stream shapefile with my model grid in FloPy but when I try to intersect my longer stream it throws an error. Below is the code I am using to load the shapefile, and the error I receive. The code works if rshp.geometry[0] but not for rshp.geometry[1]. Do you know what the problem may be? Is rshp.geometry[1] too large? Even if I modify the shapefile so there is only one linestring, it gives me an error if I try to intersect that same data that was originially in rshp.geometry[1].
Load river shapefile using geopandas
spth = '/Users/river.shp'
rshp = gpd.read_file(spth)
Plot the stream
stream = rshp.plot()
gwf.modelgrid.plot(ax=ax)
rshp.geometry
which produces:
0 LINESTRING (485833.667 5029980.118, 485778.491...
1 LINESTRING (485385.184 5032309.695, 485463.660...
Name: geometry, dtype: geometry
Create a intersection object from the modelgrid
ixobj = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid, method='structured')
Intersect the geometry object and the river shapefile goemetry
rint = ixobj.intersect_linestring(rshp.geometry[1])
Then I get this error:
NotImplementedError Traceback (most recent call last)
in
----> 1 rint = ixobj.intersect_linestring(rshp.geometry[1])
/opt/miniconda3/envs/pyclass/lib/python3.8/site-packages/flopy/utils/gridintersect.py in _intersect_linestring_structured(self, shp, keepzerolengths)
598 else: # linestring is fully within grid
599 nodelist, lengths, vertices, ixshapes =
--> 600 self._get_nodes_intersecting_linestring(
601 lineclip)
602 # if necessary, transform coordinates back to real
/opt/miniconda3/envs/pyclass/lib/python3.8/site-packages/flopy/utils/gridintersect.py in _get_nodes_intersecting_linestring(self, linestring)
744 (i, j) = nodelist[n]
745 node, length, verts, ixshape =
--> 746 self._check_adjacent_cells_intersecting_line(
747 linestring, (i, j), nodelist)
748
/opt/miniconda3/envs/pyclass/lib/python3.8/site-packages/flopy/utils/gridintersect.py in _check_adjacent_cells_intersecting_line(self, linestring, i_j, nodelist)
853 y = intersect.xy[1]
854 verts.append([(ixy[0], ixy[1])
--> 855 for ixy in zip(*intersect.xy)])
856 node.append((ii, jj))
857
/opt/miniconda3/envs/pyclass/lib/python3.8/site-packages/shapely/geometry/base.py in xy(self)
339 def xy(self):
340 """Separate arrays of X and Y coordinate values"""
--> 341 raise NotImplementedError
342
343 # Python feature protocol
NotImplementedError:
dbrakenhoff:
Thanks for submitting the issue. Can you upload a zip with the shape you're trying to intersect? Then I'll take a look at why this is throwing an error.
Meanwhile, have you tried intersecting using a different method? I recently made some changes so if you're on the latest develop branch it would be:
ixobj = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid, method='vertex')
or if you have flopy v3.3.0 installed:
ixobj = flopy.utils.gridintersect.GridIntersect(gwf.modelgrid, method='strtree')
Then try the intersection like you did before.
As posted on flopy's temporary replacement repo by ebake310:
dbrakenhoff: