Two issues related to creation of Voronoi-based grids using shapefiles

[nikobenho]

Hi,

I've encountered two situations where I was unable to create Voronoi grids from a polygon shapefile. Since I suspect this is a rather common approach for delineating a model domain, I thought I should share my experiences and how I managed to work around the issues, in case anyone is experiencing the same problems.

Issue 1:
Polygon shapefiles created using ArcGIS will have the starting AND closing vertice at the same location. e.g.

polygon_shp = np.array([
    [ 1000.0, 1000.0], #First and last points are located at the same position
    [ 1000.0, 2000.0],
    [ 2000.0, 2000.0],
    [ 2000.0, 1000.0],
    [ 1000.0, 1000.0] #First and last points are located at the same position
])

Although creating a triangular mesh using this array works fine, instantiating a Voronoi grid based on that triangular mesh will throw an error, e.g.

voronoi_grid = VoronoiGrid(tri.verts) will trigger:

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

AttributeError: 'list' object has no attribute '__array_interface__'

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
<ipython-input-15-2cc6e6ba8ade> in <module>
----> 1 voronoi_grid = VoronoiGrid(tri.verts)
      2 fig = plt.figure(figsize=(10,10))
      3 ax = plt.subplot(1, 1, 1, aspect='equal')
      4 voronoi_grid.plot(ax=ax, facecolor='none')

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\flopy\utils\voronoi.py in __init__(self, points, **kwargs)
    158 
    159     def __init__(self, points, **kwargs):
--> 160         verts, iverts = get_voronoi_grid(points, **kwargs)
    161         self.points = points
    162         self.verts = verts

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\flopy\utils\voronoi.py in get_voronoi_grid(points, **kwargs)
    116             poly = sort_vertices(poly)
    117             xc, yc = points[n]
--> 118             if not point_in_cell((xc, yc), poly):
    119                 new_vor_vertices.append((xc, yc))
    120                 iv = len(new_vor_vertices) - 1

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\flopy\utils\voronoi.py in point_in_cell(point, vertices)
     33 
     34     p = Point(point)
---> 35     poly = Polygon(vertices)
     36     if p.intersects(poly):
     37         return True

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\shapely\geometry\polygon.py in __init__(self, shell, holes)
    241 
    242         if shell is not None:
--> 243             ret = geos_polygon_from_py(shell, holes)
    244             if ret is not None:
    245                 self._geom, self._ndim = ret

c:\users\nat12nho\appdata\local\programs\python\python38\lib\site-packages\shapely\geometry\polygon.py in geos_polygon_from_py(shell, holes)
    507 
    508     if shell is not None:
--> 509         ret = geos_linearring_from_py(shell)
    510         if ret is None:
    511             return None

shapely\speedups\_speedups.pyx in shapely.speedups._speedups.geos_linearring_from_py()

ValueError: A LinearRing must have at least 3 coordinate tuples

However, omitting either the firs or last point of the polygon will work:

will_cause_error_shp = np.array([
    [ 1000.0, 1000.0], #First and last points are located at the same position
    [ 1000.0, 2000.0],
    [ 2000.0, 2000.0],
    [ 2000.0, 1000.0],
    [ 1000.0, 1000.0] #First and last points are located at the same position
])

will_work_shp = np.array([
    [ 1000.0, 1000.0], 
    [ 1000.0, 2000.0],
    [ 2000.0, 2000.0],
    [ 2000.0, 1000.0],
])

The second problem seem to arise from any situation where decimals are present. For example:

will_work_shp = np.array([
[ 1000.1, 1000.0], #The decimal will cause the creation of a malformed voronoi grid
[ 1000.0, 2000.0],
[ 2000.0, 2000.0],
[ 2000.0, 1000.0],
])

See img for example of the second problem:

Since a polygon shapefile from a real-world location is likely to contain decimals, the work-around is to round off the points to the nearest integer.

I hope these suggestions help anyone with similar problems.

[christianlangevin]

Thanks for sharing your experiences, @nikobenho. Any chance you could copy and paste an image of your voronoi grid into this issue so we can see the finished product?

[christianlangevin]

I've confirmed that there is an issue with current voronoi implemention. The issue can be encountered with some irregular domain polygons. Working to implement additional geometry processing to create closed cell polygons.