Add point, polyline, and polygon intersection capabilities

[jdhughes-dev]

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']

[jdhughes-dev]

A notebook provided by @dbrakenhoff with some intersection ideas (@rubencalje)

Intersecting Polygons and Polylines with model grids using Shapely

• Rectangular regular grid
  • Polygon with regular grid
  • Polyline with regular grid
• Triangular grid
  • Polygon with triangular grid
  • Polyline with triangular grid

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()

[dbrakenhoff]

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...

[christianlangevin]

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

mfgridutil.zip

[dbrakenhoff]

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

[jdhughes-dev]

@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?

[dbrakenhoff]

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.

[vincentpost]

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.

[dbrakenhoff]

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.

Example_GridIntersection.zip

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?

[christianlangevin]

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!

[christianlangevin]

@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!

[christianlangevin]

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)

[dbrakenhoff]

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!

flopy_grid_intersect.zip

ps. just remembered your comment about the plotting of rotated grids, but haven't yet fixed my code to take that into account.

[dbrakenhoff]

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.

flopy_grid_intersect.zip