bug: Grid intersect taking excessively long on simple voronoi grid

[christianlangevin]

Describe the bug
The grid intersect routine is timing out on a relatively simple voronoi grid

To Reproduce

This script demonstrates the problem. It creates a simple voronoi grid and plots it. The next block, however, takes 2 mins on my computer, and doesn't finish in a reasonable amount of time if the grid if more finely discretized.

import pathlib as pl
import matplotlib.pyplot as plt
import numpy as np
from shapely.geometry import Polygon, LineString
import flopy
import flopy.utils.triangle
import flopy.utils.voronoi

ws = pl.Path('.')
tempdir = ws / "temp"
tempdir.mkdir(parents=True, exist_ok=True)

theta = np.arange(0.0, 2 * np.pi, 0.1)
radius = 80000.0
xoff = radius
yoff = radius
x = xoff + radius * np.cos(theta)
y = yoff + radius * np.sin(theta)
circle_poly = [(x, y) for x, y in zip(x, y)]

tri = flopy.utils.triangle.Triangle(maximum_area=2500**2, angle=30, model_ws=tempdir)
tri.add_polygon(circle_poly)
tri.build(verbose=False)

# And now for the voronoi grid
voronoi_grid = flopy.utils.voronoi.VoronoiGrid(tri)

fig = plt.figure(figsize=(5, 5))
ax = plt.subplot(1, 1, 1, aspect="equal")
voronoi_grid.plot(ax=ax, facecolor="none")

Now do some intersecting...

modelgrid = flopy.discretization.VertexGrid(**voronoi_grid.get_gridprops_vertexgrid(), nlay=1)
gi = flopy.utils.GridIntersect(modelgrid)

ls = LineString([p for p in circle_poly])
gi.intersects(ls)

Expected behavior
Seems like this intersection should only take a few seconds. Note, that it would be nice to be able to change the maximum_area down to 1000**2.

[dbrakenhoff]

Hi @langevin-usgs, thanks for posting this issue. I can reproduce it locally. It seems like GridIntersect is taking forever to parse the cell2d list and construct the polygons. I'm a bit confused why I haven't noticed it before, but it looks like the code to build the geometries could use an update.

Creating the geometries using the following code is near instantaneous on my pc:

geoms = [shapely.polygons(modelgrid.get_cell_vertices(node)) for node in range(modelgrid.nnodes)]

Modifying _vtx_grid_to_geoms_cellids into the following fixes the issue:

def _vtx_grid_to_geoms_cellids(self):
        shapely_geo = import_optional_dependency("shapely.geometry")
        geoms = [shapely_geo.polygons(self.mfgrid.get_cell_vertices(node)) for node in range(self.mfgrid.nnodes)]
        return np.array(geoms), np.arange(self.mfgrid.nnodes)

Changing the max area to 1000**2 runs in about 0.5s on my PC now. I'll submit a PR.

[christianlangevin]

Excellent, @dbrakenhoff. Thanks for looking at this so quickly. I was hoping it was something simple. Looking forward to the fix as it will help me right away with some dev work.

[dbrakenhoff]

I'm testing the fix, but running into some issues regarding rotated/offset grids. Do you happen to know how to efficiently obtain local x,y coordinates from a VertexGrid? Or would it be simpler to transform the resulting geometries using affine or something along those lines?

[christianlangevin]

grid.py has a get_local_coords() method and I think there are other local conversions available in the discretization types. Is there anything in there that would help?

[jlarsen-usgs]

There are a couple of options. 1) get_local_coords() as Chris mentioned, 2) you could also use flopy.utils.geometry.transform(xcoords, ycoords, xoff, yoff, angrot, inverse=True)

[dbrakenhoff]

Thanks @langevin-usgs, @jlarsen-usgs, I've submitted a PR.