feat(GridIntersect): add optional kwarg to .intersect() to keep all c…

[mkennard-aquaveo]

…ells with identical vertices

Add "keep_stacked_cells" kwarg to GridIntersect.intersect() which is used to prevent removal of
cells from the intersection results which have identical vertices as a cell already included in
the results. This kwarg is only used when intersecting using shapely, and is needed when dealing
with a DISU grid which doesn't have layers. With DIS and DISV, the intersection results include
cellids which have the layer number (either LAY, ROW, COL, or LAY, CELL2D) so cells above or
below the intersected cell can be easily found just be changing the layer number. But DISU
cellids do not include the layer so it is not at all easy to find cells above or below the
intersected cell. Thus having this kwarg which keeps all the cells in the results is needed.

[codecov]

Codecov Report

Merging #1214 (db121d2) into develop (8ec8d09) will increase coverage by 1.164%.
The diff coverage is 100.000%.

❗ Current head db121d2 differs from pull request most recent head 18b87e7. Consider uploading reports for the commit 18b87e7 to get more accurate results

@@              Coverage Diff              @@
##           develop     #1214       +/-   ##
=============================================
+ Coverage   74.226%   75.391%   +1.164%     
=============================================
  Files          230       226        -4     
  Lines        50192     49418      -774     
=============================================
+ Hits         37256     37257        +1     
+ Misses       12936     12161      -775     

Impacted Files	Coverage Δ
flopy/utils/gridintersect.py	74.967% <100.000%> (+0.033%)	⬆️
flopy/mf6/utils/reference.py
flopy/utils/voronoi.py
flopy/mt3d/mtcts.py
flopy/mf6/utils/testutils.py

[dbrakenhoff]

Thanks for submitting your proposed changes. I added a few comments for you to take a look at.

I think we need to think of a better name for the option, because keep_stacked_cells doesn't quite capture the full behavior. If you set the option to True, points or linestrings on the boundary between two cells will result in 2 intersection results (even if there are no stacked cells). So maybe keep_all_intersections, return_all_intersections, ...? I do think it's a good idea to briefly describe both stacked cells case and the boundary case in the docstring for this option as examples.

Also it would be good to maybe adjust or add one of the tests to ensure that what I describe up here works as expected. There are tests for intersections on boundaries for points and linestrings, maybe we can extend those to check whether the intersection result with the option set to True does yield 2 intersection results instead of 1.

[dbrakenhoff]

keep_stacked_cells doesn't quite capture the full behavior of this option. Even without stacked cells, setting this option to True will e.g. return 2 intersection results for a point or linestring on a boundary between two cells.

also ensure <80 characters per line in docsstring

[jdhughes-dev]

I like keep_all_intersections and return_all_intersections. If I had to pick I would chose return_all_intersections.

[dbrakenhoff]

Do we want to keep this syntax with **kwargs, or should we add the keyword arguments to the intersect method explicitly? (Even if they're sometimes not used, depending on the GridIntersect settings).

[jdhughes-dev]

I would vote for an actual argument rather than a **kwargs. I would leave **kwargs that are passed on to lower level functions like matplotlib.

[dbrakenhoff]

docstring <80 chars

[christianlangevin]

Good comments, @dbrakenhoff! The use case so far has focused on the DISU keeping of stacked cells, but this has a substantial impact on the regular intersection behavior that we should probably think about more closely.

Couple quick thoughts here:

1. The keyword naming is tricky, but maybe we can think about it from the other direction. By default, we are culling the results so that our model features don't appear more than we want just because they fall on a cell edge. What about something like cullReturnedCellids=True or something like that? As @dbrakenhoff mentioned, we should describe this keyword for both use cases, which are keeping of stacked cells as well as returning cellids for any cells that touch the feature (perhaps a more common case).
2. We should make some explicit tests for the touching feature case. One that comes to mind right away is intersecting a point with a structured grid, and having the point coincide with a cell vertex. In that case, we should get the lowest number cellid with cullReturnedCellids=True, and all four cells with cullReturnedCellids=False. Might even programmatically test all cell vertices for a very small structured grid to ensure all cases work as expected. There should be a similar test for a line that runs along a column or row edge of a structured grid.
3. And I also like explicitly adding this argument instead of passing it in as a **kwarg.

[mkennard-aquaveo]

I now realize that I didn't understand all the ramifications of my proposed change. I'm going to need to spend some more time on this and figure out what we really want.

[mkennard-aquaveo]

After some discussion here at Aquavaeo we think that this change isn't really going to help us with what we want to do with DISU intersections. I think I'd just like to just close this pull request.