bug in gridintersect when intersecting linestrings starting in vertex with method="structured"

[dbrakenhoff]

I came across a bug in the structured GridIntersect routines, when linestrings start in a vertex (not on the model boundary) and continue towards the southeast. The initial intersect calculation to find the initial grid cell uses the starting point of the linestring. This returns the cell with the lowest cellid, which lies to the northwest of the vertex. Then its neighbors (north, south, east, west) are checked for intersection results, which misses the cell to the southeast of the vertex. This means the line cannot be followed through the grid and the intersection result is incomplete.

I have a fix that works, but doesn't feel very elegant. If the "check neighbors" method does not yield any intersections with length > 0, the southeastern neighboring cell is also checked for an intersection result. This increases the number of intersection calculations slightly for some cases, but does deal with the issue described above.

I will submit as part of the gridintersect PR coming up.

Code to reproduce:

import matplotlib.pyplot as plt
import numpy as np
from shapely.geometry import LineString, MultiLineString

import flopy.discretization as fgrid
from flopy.utils.gridintersect import GridIntersect

delc = 10 * np.ones(10, dtype=float)
delr = 10 * np.ones(10, dtype=float)

xoff = 0.0
yoff = 0.0
angrot = 0.0
sgr = fgrid.StructuredGrid(
    delc, delr, top=None, botm=None, xoff=xoff, yoff=yoff, angrot=angrot
)

ls1 = LineString([(95, 105), (30, 50)])
ls2 = LineString([(30, 50), (90, 22)])
ls3 = LineString([(90, 22), (0, 0)])
mls = MultiLineString(lines=[ls1, ls2, ls3])

ix = GridIntersect(sgr, method="structured")

r = ix.intersect(mls)

# necessary because ixshapes contains lists of geoms
r["isxhapes"] = [ishp for ishp in r["ixshapes"]]

# note missing portion of intersection result
fig, ax = plt.subplots(1, 1, figsize=(8, 8))
sgr.plot(ax=ax)
ix.plot_linestring(r, ax=ax, cmap="viridis")
plt.show()