bug: GridIntersect's intersect method may throw IndexError when LineString starts outside model grid

[vincentpost]

trafficstars

This error can (not will) occur when a LineString starts outside the model area for model grids that have a rotated grid or an x,y offset. I ran into it when using FloPy version 3.3.5.

It is hard to reproduce but as far as I can see the bug is caused after the LineString gets clipped (in my case in _intersect_linestring_structured, line 821 of grid_intersect.py: lineclip = shp.intersection(pl)). When _get_nodes_intersecting_linestring is called next, the first point of the clipped line is transformed to real-world coordinates using transform(...) and then intersect is called on line 1009: (i, j) = self.intersect(shapely_geo.Point(x0[0], y0[0])).cellids[0] to determine the cellid of the point. The IndexError results when the intersect method returns an empty recarray. This can happen because when the coordinates of the first point of the line are back-transformed to local coordinates (within _intersect_point_structured line 725 in gridintersect.py: rx, ry = transform(... ), the local x coordinate rx is not 0 but (in my case) -1.4551915228366852e-11. As a result, the point falls outside the bounds of the grid in local coordinates and the returned recarray is empty, causing the IndexError.

A possible solution might be to use <A href="https://docs.python.org/3/whatsnew/3.5.html#pep-485-a-function-for-testing-approximate-equality">math.isclose()</A> within find_position_in_array in gridintersect.py. But before changing the code I thought I'd throw it out here for discussion.

[dbrakenhoff]

Hi @vincentpost,

Thanks for posting this issue and getting to the bottom of the cause :)! I think your suggestion sounds reasonable, but using numpy's np.isclose instead. If you want to submit a PR, I'd be happy to review it.

Out of curiosity, does the bug occur when you use method="vertex"?.

Cheers,

Davíd

[vincentpost]

Hi Davíd,

Thanks! With method='vertex' it does work. Nevertheless I'll see if I can fix it with isclose (NumPy's version, of course) and submit a PR.

Merry Christmas, Vincent

[dbrakenhoff]

@vincentpost, did you ever get around to fixing this :innocent: ?

[wpbonelli]

We are trying to shrink the backlog before releasing 3.7.0 later this month so I took a shot at this in #2173 with @vincentpost's proposed solution. I didn't add any tests, so it would be good to confirm this fixes the original case. If anyone has a simple script reproducing the issue it could be added as a test case, feel free to push directly to the PR.

Do we need a way to adjust comparison tolerance, or are the defaults used by np.isclose() OK?