Reconciling point intersection methods

[jlarsen-usgs]

I'm starting this issue/discussion as part of my review of pull requests related to point intersection within the modflow grid. Currently there are two approaches that tackle this problem. My personal opinion is to select one of the two methods for intersection to ease the code maintenance burden.

#2654 which is an optimization for the grid.intersect() method

and

#2657 which extends the GridIntersect class and adds a points_to_cellids() method.

Both methods accomplish similar goals but take different approaches. I've also done some profiling of these methods as well as a geopandas based method to get cellids and here are the results (StructuredGrid, Python 3.12, Windows 11 OS):

Results here show that convex hull method #2654 is slightly slower than the grid intersect points_to_cellids method #2657 for up to 10,000 points. Around 100,000 points these two methods perform similarly. At 1,000,000 pts the convex hull method is about 6x faster. Profiling revealed about half of the execution time of the grid intersect method was due to recarray construction (~1.7 seconds for 1 million pts).

The geopandas method that I tested for comparison used spatial joins between geodataframes of the grid and points. This method was the fastest for up to 100,000 points of intersection, at 1,000,000 points it was slightly slower than the convex hull method. It was also only a few lines of code, but I didn't do any data prep to formalize the method as it was a test (so it would be slightly slower and more than a few lines of code if it was updated for "production").

Profiling results for Vertex/UnstructuredGrid's are very different, however.

Performance of #2657 (grid intersect) applied to a voronoi grid (3 layer, 2317 nodes per layer) is consistent to the performance it had on a structured grid. Performance of #2654 (convex hull) is poor on a voronoi grid and profiling was halted at 5000 points of intersection, due to slow performance. The difference in performance between structured and voronoi grid for #2654 is due to very different approaches taken for different grid types. For structured grids, a vectorized "sweep-line" type search is performed which is fast. For unstructured and vertex grids the convex hull method is applied and performance is slow, likely due to a lot of looping in python to search for the grid cell in an indexed neighborhood.

Consideration should also be given to the maintenance burden of the code. #2657 reduces the total amount of code in the repository by about 200 lines. #2654 adds about 1,100 lines of code to the repository, with the core feature being around 964 lines of code.

Given the results outlined here, #2657 (grid intersect improvements) seems like the best possible option moving forward. It is consistently fast (although not the fastest in all regards) across grid types. With the addition of .points_to_cellids() on GridIntersect I think we could use the .intersection() functionality on the grid as a thin wrapper that calls GridIntersect() instead of maintaining a separate intersection routine.

My recommendation is to close #2654 and not merge it in its current state as it does not perform well for vertex and unstructured grids. Any additional code review and development work should be focused on #2657 and improving the existing GridIntersect utility.

Comments and thoughts are welcome.

[aacovski]

I'll review my kdtree convexhull, because it should be performing quickly and I think it's also worth testing # of cells as well.

About half my code is adding more test edge cases.

I'll put it under draft for now as I haven't had a chance to work on it recently and compare myself against the other methods, but if your profiling is correct, 1 second for 1000000 points is fast enough and you are right that it should be closed.

[wpbonelli]

@aacovski even if your PR doesn't merge as is, I still think there is a lot of value in it — maybe with some tweaks it could move into the experimental 4.x repo.

The approach seems to shine for huge grids, where not testing every cell as in the vectorized version is a win. And the DISV/DISU performance could be improved by leaning on geopandas for containment checking? Seems worth continuing to prototype/benchmark.

Thanks a lot for your work. Anything we learn from it will inform ongoing development.

[jlarsen-usgs]

@aacovski

I've seen this happen before with ray casting algorithms and how they are vectorized where they perform very fast under certain conditions and degrade under other conditions. A lot of point-in-polygon algorithms can be optimized for either large grids or a large number of points; the challenge with that is finding a way set up the equations that balance those two constraints and be fast for large grids and a large number of points.

As @wpbonelli stated, there is a lot of value in this work and some of it could integrated into the FloPy 4.x work that's ongoing.

[aacovski]

Okay, my methods should be faster now with points greater than 1000. 2D convex hull was not as effective for 2D as other methods, which I have since replaced with a faster vectorized edge method. As for 3D (varies_by_layer), convex hull still most effective (I assume). You are welcome to checkout and test again. Furthermore, worth noting that I don't use a convexhull approach to the structuredgrid. I still need to test the geopandas optimizations