Voronoi Grid

[jonathanqv]

Hi everyone, I was not able to run the VoronoiGrid class/function until I installed the last GitHub available version. I found an error in my grid that may be useful for improving this. Will show it with images. Thanks!

I am able to remove that triangle by changing grid sizes, the general grid was 100 meters and refinement 20. With a grid size of 50 meters and refinement of 10 works nice. Thanks for the efforts!

[langevin-usgs]

Hi @jonathanqv, we are glad you reported this. Can you post a short script to this issue that would allow us to regenerate this strange case? Thank you.

[jonathanqv]

Hi @langevin-usgs , thanks for your reply. I am a big fan of your work.

I have the shapefiles here:

https://github.com/jonathanqv/VoronoiError

The script goes as follows:

`import os, sys, flopy import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt import geopandas as gpd from shapely.geometry import LineString, Point from flopy.utils.triangle import Triangle as Triangle from flopy.utils.voronoi import VoronoiGrid from flopy.discretization import VertexGrid

area = gpd.read_file('../Data/shp/area.shp')

x,y = area['geometry'].exterior[0].coords.xy area_coords = np.dstack((x,y))[0]

area_centroid = area.geometry[0].centroid.coords.xy

xmin,ymin,xmax,ymax = area.geometry.total_bounds

ref = gpd.read_file('../Data/shp/ref.shp')

x,y = ref['geometry'].exterior[0].coords.xy ref_coords = np.dstack((x,y))[0]

ref_centroid = ref.geometry[0].centroid.coords.xy

angle_min = 30 area_max = 10000.0 delr = area_max * 0.5 ncol = xmax / delr nrow = ymax / delr nodes = ncol nrow

tri = Triangle(angle=angle_min, model_ws='../Model/Grid',exe_name='../Exe/triangle.exe')#maximum_area=area_max, tri.add_polygon(area_coords) tri.add_polygon(ref_coords) tri.add_region((xmax-1, ymax-1), 0, maximum_area=100) tri.add_region((ref_centroid[0][0], ref_centroid[1][0]), 1, maximum_area=20) tri.build(verbose=False)

fig = plt.figure(figsize=(10, 10)) ax = plt.subplot(1, 1, 1, aspect="equal") pc = tri.plot(ax=ax);

voronoi_grid = VoronoiGrid(tri) fig = plt.figure(figsize=(10, 10)) ax = plt.subplot(1, 1, 1, aspect="equal") voronoi_grid.plot(ax=ax, facecolor="none");`

---

Moreover, when I try to generate the "intersect" object with this grid I get the following error:

I tried to follow your other voronoi grids examples in the "voronoi.ipynb" file and I get the same problem with the last (more complex) voronoi grid. I can avoid the error in "flopy.utils.GridIntersect" by specifying {method='vertex',rtree = False}, but when I try to intersect an object later I get the same error.

kind regards

[langevin-usgs]

So I took a quick look and was able to reproduce the issue using the data files and code that you sent. I think the problem is in some of the low-level voronoi calculations that require a certain level of precision. In your case, you have very large x and y values, so the easy solution is to just translate your gpd shapes prior to creating and working with the grid. So for example, I modified your code to be

xoff = -400000.
yoff = -4.e6

area = gpd.read_file(f'{datapath}/area.shp')
area = area.translate(xoff=xoff, yoff=yoff)
x,y = area.exterior[0].coords.xy
area_coords = np.dstack((x,y))[0]
area_centroid = area.geometry[0].centroid.coords.xy
xmin,ymin,xmax,ymax = area.geometry.total_bounds

and

ref = gpd.read_file(f'{datapath}/ref.shp')
ref = ref.translate(xoff=xoff, yoff=yoff)
x,y = ref.exterior[0].coords.xy
ref_coords = np.dstack((x,y))[0]
ref_centroid = ref.geometry[0].centroid.coords.xy

When we translate these values, the grid seems to turn out okay

This is probably good practice in general as you send these vertices directly into MODFLOW (with DISV) and may be unnecessarily losing quite a bit of precision without translating closer to a (0,0) cartesian origin.

Hope this helps. I will close for now. Feel free to reopen if there are still issues.

[jonathanqv]

@langevin-usgs Oh, that's a nice trick. Thanks! The last error I have is when I try to create the intersect object "flopy.utils.GridIntersect". I get an "list index out of range" error. I was thinking that it may be a shapely library error as I am using a windows computer with the ".whl" files for Shapely.

[dbrakenhoff]

I'll check today if I can figure out what's going wrong in GridIntersect. It seems to be tripping over an empty list of points while converting voronoi cells to shapely Polygons.

[dbrakenhoff]

There is one entry in VertexGrid._cell2d that is empty it seems. This means GridIntersect can't convert that cell to a Polygon which is the error you're seeing.

vgrid._cell2d._cell2d[14]

For me this results in:

# [icpl, xc, yc, nv, ... nv vertices numbers]
[14, 79895.47744963976, 400346.5514225885, 0]

Not sure if this is a bug in VoronoiGrid? This the location of the xc, yc (listed above) in the grid.

I can probably add a check in GridIntersect to avoid trying to create polygons for empty grid cells, but it seems cleaner to just add a grid that doesnt contain an empty cell...? Also this is probably a separate issue?

[langevin-usgs]

@dbrakenhoff, thank you for looking into this. This must be an error in the voronoi mesh generator. I'll take this on.

[jonathanqv]

If it helps, I was getting the same error with this geometry from the flopy3_voronoi.ipynb example file (the last one in the examples)

The other examples work okay.