Making polygon for every point in set using Voronoi diagram

Question

I'm trying to make polygons from points, where each point will be represented by one polygon. I tried using Voronois, but the results are incomplete. While I more or less understand where this incompleteness comes from, I don't know how to solve such problem.

I also have another, bigger polygon that is an area containing all the points. At the end, I want all the small Voronois to cover 100% of this big polygon area. In other words:

I want all the empty white spaces filled with voronois/polygons and I want each voronoi/polygon to contain just one point.

My code:

from scipy.spatial import Voronoi
import numpy as np
import shapely
import geopandas as gpd

#pointset is geodataframe with points, region is geodataframe with big 
#polygon I want to fill with voronois

x = pointset.geometry.x.values
y = pointset.geometry.y.values
coords = np.vstack((x, y)).T
vor = Voronoi(coords)

lines = [shapely.geometry.LineString(vor.vertices[line]) for line in 
    vor.ridge_vertices if -1 not in line]
polys = shapely.ops.polygonize(lines)
voronois = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polys))

voronois = gpd.overlay(voronois, region)

Basically, I just use Voronoi and then intersect the result with my region polygon. Any advice about how to make result complete?

Answer 1

Create a large rectangle surrounding the polygon. Then add many points along its boundary, concat these and the points inside and Voronoi.

import geopandas as gpd
import shapely
import numpy as np
from scipy.spatial import Voronoi
import matplotlib.pyplot as plt

#Create a polygon
polygon = shapely.wkt.loads('Polygon ((526309 6885773, 526806 6884710, 526001 6883351, 527696 6884442, 527550 6886041, 526309 6885773)))')
bound = polygon.buffer(2000).envelope.boundary #Create a large rectangle surrounding it

#Create many points along the rectangle boundary. I create one every 100 m.
boundarypoints = [bound.interpolate(distance=d) for d in range(0, np.ceil(bound.length).astype(int), 100)]
boundarycoords = np.array([[p.x, p.y] for p in boundarypoints])

#Create the points inside the polygon
coords = np.array(
          [[527014.0, 6884067.0], [527219.0, 6884190.0], [526884.0, 6884718.0], [526445.0, 6885680.0], [527024.0, 6885868.0], 
          [526962.0, 6885625.0], [527062.0, 6885553.0], [527048.0, 6885540.0], [526891.0, 6885440.0], [526897.0, 6885368.0], 
          [527034.0, 6885351.0], [527226.0, 6885252.0], [526976.0, 6885207.0], [527007.0, 6885122.0], [527171.0, 6885087.0], 
          [527041.0, 6884961.0], [527082.0, 6884851.0], [527134.0, 6884817.0], [527154.0, 6884793.0], [527144.0, 6884591.0], 
          [527117.0, 6884533.0], [527058.0, 6884457.0], [527004.0, 6884334.0], [526997.0, 6884190.0], [526918.0, 6884348.0], 
          [526935.0, 6884498.0], [527014.0, 6884591.0], [526853.0, 6884385.0], [526808.0, 6884228.0], [526777.0, 6884156.0], 
          [526647.0, 6884022.0], [526593.0, 6883854.0], [526469.0, 6883748.0], [526295.0, 6883652.0], [526284.0, 6883621.0], 
          [527528.0, 6885488.0], [527517.0, 6885406.0], [527531.0, 6885327.0], [527541.0, 6885231.0], [527534.0, 6885183.0], 
          [527545.0, 6885087.0], [527569.0, 6884992.0], [527576.0, 6884899.0], [527576.0, 6884824.0], [527552.0, 6884690.0], 
          [527524.0, 6884632.0], [527439.0, 6884635.0], [527329.0, 6884865.0], [527257.0, 6884690.0], [527291.0, 6884444.0]])

all_coords = np.concatenate((boundarycoords, coords)) #Create an array of all points on the boundary and inside the polygon

vor = Voronoi(points=all_coords)
lines = [shapely.geometry.LineString(vor.vertices[line]) for line in 
    vor.ridge_vertices if -1 not in line]

polys = shapely.ops.polygonize(lines)
voronois = gpd.GeoDataFrame(geometry=gpd.GeoSeries(polys), crs="epsg:3006")

polydf = gpd.GeoDataFrame(geometry=[polygon], crs="epsg:3006")
points = gpd.GeoDataFrame(geometry=gpd.points_from_xy(x=coords[:,0], y=coords[:,1], crs="epsg:3006"))

result = gpd.overlay(df1=voronois, df2=polydf, how="intersection")

fig, ax = plt.subplots(figsize=(15, 15))
polydf.boundary.plot(ax=ax, edgecolor="blue", linewidth=6)
result.plot(ax=ax, color="red", alpha=0.3, edgecolor="black")
points.plot(ax=ax, color="maroon")