Add point, polyline, and polygon intersection capabilities
jdhughes-usgs opened this issue · 18 comments
Here are some thoughts related to March 2019 developer meeting at Delft.
Possible Dependencies
geopandas, shapely, others
Projections, sampling, etc
Initial thoughts for syntax of calls
Point shape to grid for boundary conditions
>>> mg = modelgrid
>>> int_dict = mg.get_cellids(shp_xy, type=POINT)
>>> list(int_dict.keys())
['cellids', 'point_id']
>>> wel_spd = [((cellid), q) for cellid, q in zip(int_dict['cellids'], qs)]
or
>>> wel_spd = [((cellid), q) for cellid, q in zip(int_dict.cellids, qs)]
Polyline shape to grid for boundary conditions
>>> mg = modelgrid
>>> int_dict = mg.get_cellids(shp_poly, type=POLYLINE)
>>> list(int_dict.keys())
['cellids', 'poly_id', 'length']
Polygon shape to grid for boundary conditions
>>> mg = modelgrid
>>> int_dict = mg.get_cellids(shp_polygon, type=POLYGON)
>>> list(int_dict.keys())
['cellids', 'polyg_id', 'area']
>>> mg = modelgrid
>>> int_dict = mg.get_cellids(shp_polygon, type=POLYGON, vertices=True)
>>> list(int_dict.keys())
['cellids', 'polyg_id', 'area', 'vert']
A notebook provided by @dbrakenhoff with some intersection ideas (@rubencalje)
Intersecting Polygons and Polylines with model grids using Shapely
import numpy as np
import shapely.geometry as geom
from shapely.geometry import Polygon, LineString, MultiLineString, MultiPolygon
from shapely.strtree import STRtree
from descartes import PolygonPatch
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
General intersect function based on list of Shapely geometries that represent a grid, a Shapely Polygon or Polyline to intersect with and a declaration of shptype (could probably also be inferred from shp).
def intersect(mglist, shp, shptype="POLYGON"):
intersect_dict = {}
s = STRtree(mglist)
result = s.query(shp)
isectshp = []
cellids = []
vertices = []
areas = []
lengths = []
for i, r in enumerate(result):
intersect = shp.intersection(r)
if shptype == "POLYGON":
if intersect.area > 0.0:
cellids.append(mglist.index(r))
isectshp.append(intersect)
areas.append(intersect.area)
vertices.append(intersect.__geo_interface__["coordinates"])
elif shptype == "POLYLINE":
try:
isect_iter = iter(intersect)
except TypeError:
isect_iter = [intersect]
for isect in isect_iter:
if isect.length > 0.0:
cellids.append(mglist.index(r))
isectshp.append(isect)
lengths.append(isect.length)
vertices.append(isect.__geo_interface__["coordinates"])
else:
raise NotImplementedError("shptype '{}' is not supported!".format(shptype))
intersect_dict["intersect"] = isectshp
intersect_dict["cellids"] = cellids
intersect_dict["vertices"] = vertices
if shptype == "POLYGON":
intersect_dict["areas"] = areas
elif shptype == "POLYLINE":
intersect_dict["lengths"] = lengths
return intersect_dict
Rectangular regular grid
grid_cells = []
nx = 50
ny = 50
x0, y0 = 0, 0
dx, dy = 1000./nx, 1000./ny
for i in range(ny):
for j in range(nx):
x = x0+j*dx
y = y0+i*dy
xy = [[x, y], [x+dx, y], [x+dx, y+dy], [x, y+dy], [x, y]]
grid_cells.append(Polygon(xy))
Shapes to intersect with:
poly_intersect = Polygon(shell=[(150, 150), (200, 500), (350, 750), (810, 440), (400, 50), (150, 120)],
holes=[[(250, 250), (250, 450), (450, 450), (450, 250)]])
ls1 = LineString([(50, 50), (350, 475)])
ls2 = LineString([(350, 475), (875, 225)])
mls = MultiLineString(lines=[ls1, ls2])
poly_intersect
mls
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.set_aspect("equal")
polys = []
for g in grid_cells:
pp = PolygonPatch(g, facecolor="C0", alpha=0.5)
ax.add_patch(pp)
polys.append(g)
pp2 = PolygonPatch(poly_intersect, alpha=0.5, facecolor="red")
ax.add_patch(pp2)
for ig in mls.geoms:
ax.plot(ig.xy[0], ig.xy[1])
ax.set_xlim(-dx, nx*1000./nx+dx)
ax.set_ylim(-dy, ny*1000./ny+dy)
(-20.0, 1020.0)
Polygon with regular grid
%%time
result = intersect(polys, poly_intersect, shptype="POLYGON")
print(result.keys())
dict_keys(['intersect', 'cellids', 'vertices', 'areas'])
Wall time: 22.6 s
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.set_aspect("equal")
for g in grid_cells:
pp = PolygonPatch(g, edgecolor="k", alpha=1.0, facecolor="none")
ax.add_patch(pp)
for i, ishp in enumerate(result["intersect"]):
ppi = PolygonPatch(ishp, facecolor="C{}".format(i%10), alpha=0.5)
ax.add_patch(ppi)
for cid in result["cellids"]:
c = polys[cid].centroid
ax.plot(c.x, c.y, "r.")
ax.set_xlim(-dx, (nx)*1000/nx+dx)
ax.set_ylim(-dy, (ny)*1000/ny+dy)
plt.show()
Polyline with regular grid
%%time
result = intersect(polys, mls, shptype="POLYLINE")
print(result.keys())
dict_keys(['intersect', 'cellids', 'vertices', 'lengths'])
Wall time: 2.17 s
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.set_aspect("equal")
for g in grid_cells:
pp = PolygonPatch(g, edgecolor="k", alpha=1.0, facecolor="none")
ax.add_patch(pp)
for i, ishp in enumerate(result["intersect"]):
ax.plot(ishp.xy[0], ishp.xy[1], ls="-", c="C{}".format(i%10))
for cid in result["cellids"]:
c = polys[cid].centroid
ax.plot(c.x, c.y, "r.")
ax.set_xlim(-dx, (nx)*1000./nx+dx)
ax.set_ylim(-dy, (ny)*1000./ny+dy)
plt.show()
Triangular grid
from flopy.utils.triangle import Triangle as Triangle
import config
flopy is installed in C:\Users\dbrak\Anaconda3\lib\site-packages\flopy
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mf2005.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mfnwt.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mp7.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mt3dms.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mt3dusgs.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mf6.exe
Executable file could not be found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\mf6beta.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\gridgen.exe
Executable file found: C:\GitHub\mf6flopy2019_classrepo\exercises\bin\win64\triangle.exe
maximum_area = 500.
domainpoly = [(x0, y0), (x0, y0+ny*dy), (x0+nx*dx, y0+ny*dy), (x0+nx*dx, y0)]
tri = Triangle(maximum_area=maximum_area, angle=30, model_ws=".", exe_name=config.triangleexe)
tri.add_polygon(domainpoly)
tri.build(verbose=False)
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(1, 1, 1, aspect='equal')
pc = tri.plot(ax=ax)
tri.ncpl
3224
iverts, verts = tri.iverts, tri.verts
%%time
ptchs = []
tridict = {}
for icell, ivertlist in enumerate(iverts):
points = []
for iv in ivertlist:
points.append((verts[iv, 0], verts[iv, 1]))
# close the polygon, if necessary
if ivertlist[0] != ivertlist[-1]:
iv = ivertlist[0]
points.append((verts[iv, 0], verts[iv, 1]))
ptchs.append(Polygon(points))
Wall time: 28 ms
Polygon with triangular grid
%%time
result = intersect(ptchs, poly_intersect, shptype="POLYGON")
print(result.keys())
dict_keys(['intersect', 'cellids', 'vertices', 'areas'])
Wall time: 36.9 s
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(1, 1, 1, aspect='equal')
pc = tri.plot(ax=ax)
for i, ishp in enumerate(result["intersect"]):
ppi = PolygonPatch(ishp, facecolor="C{}".format(i%10), alpha=0.5)
ax.add_patch(ppi)
for cid in result["cellids"]:
c = ptchs[cid].centroid
ax.plot(c.x, c.y, "r.")
ax.set_xlim(-dx, (nx)*1000./nx+dx)
ax.set_ylim(-dy, (ny)*1000./ny+dy)
plt.show()
Polyline with triangular grid
%%time
result = intersect(ptchs, mls, shptype="POLYLINE")
print(result.keys())
dict_keys(['intersect', 'cellids', 'vertices', 'lengths'])
Wall time: 4.85 s
fig = plt.figure(figsize=(10, 10))
ax = plt.subplot(1, 1, 1, aspect='equal')
pc = tri.plot(ax=ax)
for i, ishp in enumerate(result["intersect"]):
ax.plot(ishp.xy[0], ishp.xy[1], ls="-", c="C{}".format(i%10))
for cid in result["cellids"]:
c = ptchs[cid].centroid
ax.plot(c.x, c.y, "r.")
ax.set_xlim(-dx, (nx)*1000./ny+dx)
ax.set_ylim(-dy, (ny)*1000./nx+dy)
plt.show()
As @langevin-usgs mentioned in an email, the timings on intersecting with the unstructured grid are quite slow. The times shown in the notebook are what i get on my laptop. When I was writing the notebook I hadn't really noticed initially, but I think i was working with a much smaller grid then.
I'll look into why this is the case... could be that STRtree isn't very efficient for unstructured grids, or shapely just isn't as fast as I hoped it would be...
This is something that we wrote a while back that has some accelerated intersection capabilities, at least for structured grids. It still uses shapely, but it limits the number of cells that it has to check by following the line through the grid.
'''
Classes and methods for performing spatial analyses with a ModflowGrid
object or a ModflowGrid2D object.
'''
class PointGridIntersection(object):
'''
'''
def __init__(self, grid, point, **kwargs):
'''
Find the point in the grid. Return None if not found in range,
otherwise return (row, column) in zero-based indexing. If x, y, fall
directly on a cell edge, return the smaller cell indices.
Arguments:
*point*: A tuple containing the x and y coordinates of the point,
(x, y) or a tuple containing the x, y, and z coordinates of
the point (x, y, z).
Store the node (ipos, jpos) or (kpos, ipos, jpos) in self.nodelist.
Set self.nodelist = None if point is not in grid.
'''
from mfgridutil.gridutil import ModflowGridIndices
self.nodelist = None
#two dimensional point
Xe = grid.Xe
Ye = grid.Ye
x = point[0]
jpos = ModflowGridIndices.find_position_in_array(Xe, x)
y = point[1]
ipos = ModflowGridIndices.find_position_in_array(Ye, y)
self.nodelist = (ipos, jpos)
#three dimensional point
if len(point) == 3:
#find k
z = point[2]
kpos = ModflowGridIndices.find_position_in_array(
grid.botm[:, ipos, jpos], z)
if kpos is None:
self.nodelist = None
else:
self.nodelist = (kpos, ipos, jpos)
return
class LineGridIntersection(object):
'''
This class contains the methods for intersecting a polyline with a
MODFLOW grid object. The nodes intersecting the line and the lengths of
each line are stored in the nodes and lengths properties of the class
object.
'''
def __init__(self, grid, line, keepzerolengths=False):
'''
Create the line intersection object. Store the list of nodes that
intersect the line in self.nodelist. Store the corresponding lengths in
self.lengths.
Nodes are represented as (row, col).
Arguments:
*grid*: A ModflowGrid or ModflowGrid2D object.
*line*: A tuple or list of points defining a line.
*keepzerolengths*: A true or false flag indicating whether line
segments that have zero lengths should be included in the
list of nodes and list of lengths. Zero length line segments
can occur when a line touches a cell edge.
'''
from shapely.geometry import LineString, Polygon, MultiLineString, box
self.grid = grid
self.line = line
(xmin, ymin), (xmax, ymax) = self.grid.get_extent()
pl = box(xmin, ymin, xmax, ymax)
ls = LineString(self.line)
lineclip = ls.intersection(pl)
self.nodelist = []
self.lengths = []
if lineclip.length == 0.:
return
if lineclip.geom_type is 'MultiLineString': #there are multiple lines
for l in lineclip:
self.get_nodes_intersecting_linestring(l)
else:
self.get_nodes_intersecting_linestring(lineclip)
#eliminate any nodes that have a zero length
if not keepzerolengths:
tempnodes = []
templengths = []
for i in range(len(self.nodelist)):
if self.lengths[i] > 0:
tempnodes.append(self.nodelist[i])
templengths.append(self.lengths[i])
self.nodelist = tempnodes
self.lengths = templengths
return
def get_nodes_intersecting_linestring(self, linestring):
'''
Intersect the linestring with the model grid and return a list of
node indices and the length of the line in that node.
'''
from shapely.geometry import LineString, Polygon, MultiLineString, box
#start at the beginning of the line
x, y = linestring.xy
x0 = x[0]
y0 = y[0]
(i, j) = self.grid.intersection((x0, y0), 'point')
xmin = self.grid.Xe[j]
xmax = self.grid.Xe[j + 1]
ymax = self.grid.Ye[i]
ymin = self.grid.Ye[i + 1]
pl = box(xmin, ymin, xmax, ymax)
length = linestring.intersection(pl).length
self.lengths.append(length)
self.nodelist.append( (i, j) )
n = 0
while True:
(i, j) = self.nodelist[n]
self.check_adjacent_cells_intersecting_line(linestring, (i, j))
if n == len(self.nodelist) - 1:
break
n += 1
return
def check_adjacent_cells_intersecting_line(self, linestring, i_j):
i,j = i_j
from shapely.geometry import LineString, Polygon, box
#check to left
if j > 0:
ii = i
jj = j - 1
if (ii, jj) not in self.nodelist:
xmin = self.grid.Xe[jj]
xmax = self.grid.Xe[jj + 1]
ymax = self.grid.Ye[ii]
ymin = self.grid.Ye[ii + 1]
pl = box(xmin, ymin, xmax, ymax)
if linestring.intersects(pl):
length = linestring.intersection(pl).length
self.nodelist.append( (ii,jj) )
self.lengths.append(length)
#check to right
if j < self.grid.ncol - 1:
ii = i
jj = j + 1
if (ii, jj) not in self.nodelist:
xmin = self.grid.Xe[jj]
xmax = self.grid.Xe[jj + 1]
ymax = self.grid.Ye[ii]
ymin = self.grid.Ye[ii + 1]
pl = box(xmin, ymin, xmax, ymax)
if linestring.intersects(pl):
length = linestring.intersection(pl).length
self.nodelist.append( (ii,jj) )
self.lengths.append(length)
#check to back
if i > 0:
ii = i - 1
jj = j
if (ii, jj) not in self.nodelist:
xmin = self.grid.Xe[jj]
xmax = self.grid.Xe[jj + 1]
ymax = self.grid.Ye[ii]
ymin = self.grid.Ye[ii + 1]
pl = box(xmin, ymin, xmax, ymax)
if linestring.intersects(pl):
length = linestring.intersection(pl).length
self.nodelist.append( (ii,jj) )
self.lengths.append(length)
#check to front
if i < self.grid.nrow - 1:
ii = i + 1
jj = j
if (ii, jj) not in self.nodelist:
xmin = self.grid.Xe[jj]
xmax = self.grid.Xe[jj + 1]
ymax = self.grid.Ye[ii]
ymin = self.grid.Ye[ii + 1]
pl = box(xmin, ymin, xmax, ymax)
if linestring.intersects(pl):
length = linestring.intersection(pl).length
self.nodelist.append( (ii,jj) )
self.lengths.append(length)
return
class RectangleGridIntersection(object):
'''
'''
def __init__(self, grid, rectangle, **kwargs):
'''
Given a rectangle defined as [(xmin, ymin), (xmax, ymax)]
return the cells (k, i, j) that are within or touching
the rectangle. This is faster than using the more generic
PolygonGridIntersect approach.
Arguments:
*grid*: The ModflowGrid or ModflowGrid2D object.
*rectangle*: A tuple containing ((xmin, ymin), (xmax, ymax))
'''
from shapely.geometry import Point, Polygon, box
from mfgridutil.gridutil import ModflowGridIndices
self.nodelist = []
#return if rectangle does not contain any cells
(xmin, ymin), (xmax, ymax) = grid.get_extent()
bgrid = box(xmin, ymin, xmax, ymax)
(xmin, ymin), (xmax, ymax) = rectangle
b = box(xmin, ymin, xmax, ymax)
if not b.intersects(bgrid):
#return with nodelist as an empty list
return
Xe = grid.Xe
Ye = grid.Ye
jmin = ModflowGridIndices.find_position_in_array(Xe, xmin)
if jmin is None:
if xmin <= Xe[0]:
jmin = 0
elif xmin >= Xe[-1]:
jmin = grid.ncol - 1
jmax = ModflowGridIndices.find_position_in_array(Xe, xmax)
if jmax is None:
if xmax <= Xe[0]:
jmax = 0
elif xmax >= Xe[-1]:
jmax = grid.ncol - 1
imin = ModflowGridIndices.find_position_in_array(Ye, ymax)
if imin is None:
if ymax >= Ye[0]:
imin = 0
elif ymax <= Ye[-1]:
imin = grid.nrow - 1
imax = ModflowGridIndices.find_position_in_array(Ye, ymin)
if imax is None:
if ymin >= Ye[0]:
imax = 0
elif ymin <= Ye[-1]:
imax = grid.nrow - 1
for i in range(imin, imax + 1):
for j in range(jmin, jmax + 1):
self.nodelist.append( (i, j) )
return
class PolygonGridIntersection(object):
'''
'''
def __init__(self, grid, polygon, **kwargs):
'''
Find the nodes within the polygon.
'''
from mfgridutil.gridutil import ModflowGridIndices
from shapely.geometry import Point, Polygon
holes = None
if 'holes' in kwargs.keys(): holes = kwargs['holes']
#initialize the result arrays
self.nodelist = []
self.areas = []
self.containscentroid = []
pg = Polygon(polygon, holes=holes)
#set X and Y
X = grid.X
Y = grid.Y
#use the bounds of the polygon to restrict the cell search
minx, miny, maxx, maxy = pg.bounds
rectangle = ((minx, miny), (maxx, maxy))
nodelist = RectangleGridIntersection(grid, rectangle).nodelist
for (i, j) in nodelist:
nodenumber = ModflowGridIndices.nn0_from_kij(0, i, j,
grid.nrow, grid.ncol)
#node = grid.get_nodeobj(nodenumber)
#node_polygon = Polygon(node.vertices)
node_polygon = Polygon(grid.get_vertices(nodenumber))
if pg.intersects(node_polygon):
area = pg.intersection(node_polygon).area
if area > 0.:
self.nodelist.append( (i, j) )
self.areas.append(area)
#pt = Point(node.position)
#if pg.contains(pt):
# self.containscentroid = True
#else:
# self.containscentroid = False
return
I figured it out. The time for intersecting the polygon with the unstructured grid is now 150 ms instead of 35 seconds. The lookup in the list of modelgrid shapes to get the cellid from the modelgrid is extremely slow.
If you comment out those lookups in the intersect function it's a lot faster.
def intersect(mglist, shp, shptype="POLYGON"):
intersect_dict = {}
s = STRtree(mglist)
result = s.query(shp)
isectshp = []
cellids = []
vertices = []
areas = []
lengths = []
for i, r in enumerate(result):
intersect = shp.intersection(r)
if shptype == "POLYGON":
if intersect.area > 0.0:
# cellids.append(mglist.index(r)) # this is extremely slow
isectshp.append(intersect)
areas.append(intersect.area)
vertices.append(intersect.__geo_interface__["coordinates"])
elif shptype == "POLYLINE":
try:
isect_iter = iter(intersect)
except TypeError:
isect_iter = [intersect]
for isect in isect_iter:
if isect.length > 0.0:
# cellids.append(mglist.index(r)) # this is extremely slow
isectshp.append(isect)
lengths.append(isect.length)
vertices.append(isect.__geo_interface__["coordinates"])
else:
raise NotImplementedError("shptype '{}' is not supported!".format(shptype))
intersect_dict["intersect"] = isectshp
intersect_dict["cellids"] = cellids
intersect_dict["vertices"] = vertices
if shptype == "POLYGON":
intersect_dict["areas"] = areas
elif shptype == "POLYLINE":
intersect_dict["lengths"] = lengths
return intersect_dict
@dbrakenhoff and @rubencalje, @jdhughes-usgs and @langevin-usgs discussed using a np.recarray
instead of the dictionary for the return. This would make it as extensible as a dictionary and may make it easier to iterate over.
Thoughts?
I am working on combining the code @langevin-usgs posted and the notebook I made. I'm fine with the np.recarray
as a return, so I'll implement it and see how that works.
On the speed issue related to the line
cellids.append(mglist.index(r))
in intersect()
A possible fix is to replace it with (see: shapely/shapely#618)
cellids.append(r.name)
where name
contains the id of the cell that one assigns when creating the polygons. For example for the rectangular regular grid example above:
k = 0
for i in range(ny):
for j in range(nx):
x = x0+j*dx
y = y0+i*dy
xy = [[x, y], [x+dx, y], [x+dx, y+dy], [x, y+dy], [x, y]]
p = Polygon(xy)
p.name = k
grid_cells.append(p)
k+=1
This worked for me so I thought I'd share it here.
Thanks for the tip @vincentpost !
It's taken a while, but I've attached a notebook that contains a GridIntersect class which structures the intersect functionality from the first notebook. The way it works is you create an GridIntersect object from your flopy model grid. This converts the gridcells to Shapely geometries and builds an STRTree.
Intersects can then be performed using intersect_point()
, intersect_polyline()
and intersect_polygon()
. The result is a numpy.rec.array containing cellids, vertices (coordinates of the shape within that gridcell), the intersection Shapely shape, and areas or lengths depending on what shape you're intersecting with your grid.
I've added some tests and some timing at the bottom of the notebook.
Curious to hear what you think!
If this is something worth including, do you think adding it to the existing Grid objects like the example below would be a good idea? Or would you prefer adding as something separate?
class StructuredGrid:
def intersect_polylines(self, list_of_lines):
ix = GridIntersect(self)
result = []
for line in list_of_lines:
result.append(ix.intersect_polyline(line))
return result
This example does mean the conversion to Shapely and building the STRTree are performed for each intersect call (but you can do multiple intersects in one go). It does mean that changing the grid after creation won't result in unexpected differences.
ps. @langevin-usgs I wanted to compare with your gridops.py file, but I couldn't find mfgridutil.gridutil
. Any idea where I can find those functions?
There is a version of the mfgridutil here: https://github.com/usgs/gw-general-models/tree/master/general-models/MFGrid/mfgrid.
Excited to take a look at your intersection capabilities. Will try to do so over the next few days. Thanks for putting it together!
@dbrakenhoff, I'm starting to play with your nice intersect class. I ran into one quick thing early on, that is just related to plotting. Trying passing in a rotation angle like:
sgr = fgrid.StructuredGrid(delc, delr, top=None, botm=None, xoff=0.0, yoff=0.0, angrot=45)
In order for things to show up right, you need to replace [irow, 0] and [0, icol] with [irow, icol], but this is trivial of course. We are at the MODFLOW conference next week, so maybe we can talk it over with some people and see what they think. Will keep you posted!
One thing to think about: what is a reasonable amount of time for some common operations? For example, the following search for 100 random points in a million cell model takes a bit of time (a couple minutes I think). Do we need to have an optimized search for structured grids in this case?
nrnc = 1000
delc = np.ones(nrnc, dtype=np.float)
delr = np.ones(nrnc, dtype=np.float)
sgr = fgrid.StructuredGrid(delc, delr, top=None, botm=None)
ix = GridIntersect(sgr)
npoints = 100
points = np.random.random((npoints, 2)) * nrnc
mp = MultiPoint(points=[Point(x, y) for x, y in points])
ix.intersect_point(mp, return_all_ix=True)
I agree we should make use of the optimized search if at all possible.
I've added the methods you sent me earlier to the GridIntersect class. The methods used now depend on the method keyword argument passed when you instantiate the class. By setting method="structured"
it uses the functions you wrote for structured grids.
I've added more tests, and included some timing comparisons in the notebook in the attached zipfile.
Curious to hear what you think!
ps. just remembered your comment about the plotting of rotated grids, but haven't yet fixed my code to take that into account.
Just submitted a PR #610 .
I also have a notebook containing some of the functionality, but not sure if I can add that to flopy notebooks, because it cannot be run unless shapely is available.
Here's a zipfile containing the gridintersect code, the tests and the example notebook which can be run seperately.
Can you try to update the install block in the .travis.yml
file to install shapely and test the notebook on travis?
- pip install shapely[vectorized]
- pip install nbconvert
- pip install nose-timer
- pip install coveralls
- pip install pylint
You should use the same initial code block from the other notebooks so that flopy loads correctly.
According to PyPi the [vectorized]
option adds a few extra speedups that depend on Numpy. If the install fails you can try it without the option[vectorized]
(i.e., - pip install shapely
).
That seems to work. There are some failures on Python 2.7 for grid intersection and there is some unrelated Notebook failing.
Does the failing notebook run for you? Also is your fork current with develop on upstream?
Added by #610 . Issue can be closed.
Functionality has been added with #610. Closing...