feat(GridIntersect): add array support, fast point locating and support for points with z-coordinates

[dbrakenhoff]

Major update of GridIntersect module.

@Reviewer(s), I have a few things I wonder about, especially the naming of things:

• A new method called points_to_cellids() for locating points in the grid. Is this name clear? And is cellid the correct name for (row, column) or (cellid,) output?
• I use "cellids" as a field name in the result but when return_nodenumbers=True, the intersects result can also return node numbers instead of (row, columns) for structured grids.
• Anyway, a critical look at the naming of stuff in the module would be greatly appreciated!

Support for numpy arrays of shapely geometries

Introduces support for arrays of shapely geometries in GridIntersect.intersects() method.

gi = GridIntersect(grid)
shapes = shapely.polygons(...)
gi.intersects(shapes)

Returns rec-array with a shp_id field containing the index of the original shape in shapes and the cellids of the candidate cells the polygons intersect with, for example:

rec.array([(0, (np.int64(1), np.int64(0))),
           (0, (np.int64(1), np.int64(1))),
           (1, (np.int64(0), np.int64(0))),
           (1, (np.int64(0), np.int64(1)))],
          dtype=[('shp_ids', '<i8'), ('cellids', 'O')])

So in that example output polygon 0 intersects with 2 grid cells: (1, 0) and (1, 1), and polygon 1 intersects with 2 grid cells: (0, 0), and (0, 1).

Fast point locating

GridIntersect was missing a method to obtain one cellid for each input point. The implementation of GridIntersect.intersects() drops all points outside the model grid, and potentially returns all intersection results for each point, meaning there is no 1:1 output between input points and output cellids. The GridIntersect.intersect() method allows a MultiPoint as input, but is inefficient, and groups points in the same cell as MultiPoints, also losing the 1:1 relation.

So a new method is introduced GridIntersect.points_to_cellids(). This method guarantees len(input) == len(output). Cellids are returned for each point and NaN is returned if the point is outside the model grid. For points on boundaries the lowest grid cellid is returned.

Support for points with z-coordinates

Adds support for points with z-coordinates in GridIntersect.intersect() and GridIntersect.points_to_cellids() with the handle_z-kwarg.

• The default is "ignore" which ignores the z-coordinates, similar to what GridIntersect did before.
• The second option is "drop" which removes the intersection results if the point lies above or below the model grid. This option is only included in intersect().
• The final option is "return" which returns the computed layer position for each point. This returns NaN if the point is outside the model domain, +inf if the point is above the model and -inf the point is below the model.

Deprecations

Deprecates everything related to method="structured"

Other stuff

• Cleaned up the module somewhat with some renaming since the structured methods dropped out
• Used shapely plotting routines for linestrings and points to simplify plotting calls, make static plot methods accept shapely geometries instead of rec-array results so they can also be used to plot input geometries.
• Update tests to test new features and clean up old tests
• Update example notebook.

[dbrakenhoff]

Tests failing because MF6 examples still use method kwarg. Getting tests to pass, will require a change in modflow6-examples repo:

MODFLOW-ORG/modflow6-examples#315

[aacovski]

@dbrakenhoff First off, apologies for overlapping with your work...

searchsorted (structured) and kdtree (vertex, unstructured) in my PR already provides 1:1 point-to-cell mapping - each query point returns exactly one cell (or nan if outside, can change from nan). We could add wrappers or extend grid.intersect to use this approach while maintaining the existing API... So I'm a bit confused by the timing of this PR.

The key difference from shapely: shapely handles non-convex cell geometries, whereas our approach assumes convex hulls. This assumption is valid for typical MODFLOW grids (rectangles, Voronoi cells, triangles) and enables significant simplifications and speedups. More importantly, the convex hull approach generalizes cleanly to higher-dimensional shapes - the same hyperplane equations can work for lines, polygons, and polyhedrons without needing shapely at all. I can work with you in the future once I'm done with the point to cell PR for higher-dimensional shapes.

[dbrakenhoff]

@dbrakenhoff First off, apologies for overlapping with your work...

No worries, there were already 2 places in flopy with grid intersection operations, so our simultaneous work hasn't changed that, but they have become slightly more similar now that they both support arrays of points.

searchsorted (structured) and kdtree (vertex, unstructured) in my PR already provides 1:1 point-to-cell mapping - each query point returns exactly one cell (or nan if outside, can change from nan). We could add wrappers or extend grid.intersect to use this approach while maintaining the existing API... So I'm a bit confused by the timing of this PR.

I realize there is some overlap with your work, especially between the new points_to_cellids() method and the array-based grid.intersect() methods you built for grids. The inefficient method for intersecting points in GridIntersect has bugged me for quite some time, plus as you can see in the list of issues and the TODOs that were fixed in the code itself, there was a bunch of work we still wanted to do on the module, which is all included in this PR. Also I felt it was weird for GridIntersect to have good support for linestrings and polygons but not have a decent implementation for points. Yes I could have just used self.mfgrid.intersect() but since the whole GridIntersect module was written around shapely, I wanted to stay with shapely for now, and see what I could achieve with those tools. Shapely's approach might have some limitations, but it seems very fast, and very flexible for the most common intersection operations.

So the questions now are, do we continue to let the two live side-by-side? Or do we integrate them and implement grid.intersect() in GridIntersect for point intersection calls (I would be curious to see how the performance compares?) ? But this might also be more of a question for flopy4, where a significant redesign of flopy will take place. And in the meantime we can enjoy having two intersection methods :).