# Insert n points equally distributed into many polygons

I know the "Random Points inside Polygons" tool in QGIS. It generates the below output.

What I want to achieve however are 8 equally distributed points for each of the polygons. I thought of dividing the polygons Voronoi-style by the random points and then generating new (equally distributed) points, but I don't seem to be able to do so. Using "Voronoi Polygons" in the toolbox, I only get a Voronoi diagram that's bounded by a rectangle. No way to divide the whole polygon by it:

How can I get 8 equally distributed points inside each of the polygons?

• You are very close, just for a start you need to create its borders for each polygon, use the <polygon-in-line> tool, then divide each line evenly into 8 parts, and only then run the <Voronoi Polygons> tool for each polygon and then split your polygons with the result <Voronoi Polygons>, and then extract the centroids ... If a programmer think how to automate this process, for example on a python ... Sep 2, 2019 at 17:09
• Grid pattern of points would be ok, too! :) Prerequisite is that there are 8 points in each polygon. Sep 3, 2019 at 8:48

I implemented Cyril's comment. This is the input test polygon:

1. Polygon to Lines

1. Split lines by maximum length using \$length / 7 in the expression builder. (8 vertices)

1. Extract specific vertices (0 stands for start point)

1. Voronoi Polygons with 100% buffer to make sure that parts of the polygon do not lie outside the generated Voronoi (in theory they may still lie outside, but 100% is quite big)

1. Intersect Voronoi and test polygon

1. Point on surface (centroids could fall into holes)

The resulting 8 points look like this inside the test polygon:

I will now try to automate this procedure and edit my answer when I have results.

Update 5.12.2019: I haven't been able to automate this task.
Update 4.1.2022: Here is an implementation in Shapely:

from itertools import count, islice

from shapely.geometry import MultiPoint, Polygon
from shapely.ops import substring, voronoi_diagram

def _get_voronoi_starting_points(polygon: Polygon, point_count: int) -> MultiPoint:
perimeter = polygon.exterior
segment_length = perimeter.length / point_count
segment_starts = islice(count(0, segment_length), 0, point_count)

return MultiPoint([substring(perimeter, start, start) for start in segment_starts])

polygon =  # shapely.geometry.Polygon
point_count = 8

voronoi = voronoi_diagram(_get_voronoi_starting_points(polygon, point_count))
equally_distributed_points = MultiPoint(
[part.intersection(polygon).representative_point() for part in voronoi.geoms]
)
• I suppose that QGIS's graphical modeler might be a suitable solution Sep 6, 2019 at 6:26
• Can I feed the graphical modeler with features one by one? Sep 6, 2019 at 6:39
• Good question. By so far I would say IDK. Because 'iterate over this layer' is not available in modeler Iterative execution of algorithms but what if "Split vector layer"will be used, which is available in QGIS graphical modeler. I have not tried it yet Sep 6, 2019 at 6:48

I use this code to ~randomly place a grid of sample plots for forest inventory:

import numpy as np
from itertools import product

polylayer = iface.activeLayer() #Highlight polygon layer in layer tree
npoints = 8 #Change

pointlayer = QgsVectorLayer('Point?crs=epsg:3006', 'point' , 'memory') #Change epsg
prov = pointlayer.dataProvider()

for poly in polylayer.getFeatures():
geom = poly.geometry()
bbox = poly.geometry().boundingBox()
xmin, xmax, ymin, ymax = bbox.xMinimum(),bbox.xMaximum(),bbox.yMinimum(),bbox.yMaximum()
success = 0
while success == 0:
featlist = []
spacing = ((geom.area()/npoints)**0.5)*np.random.uniform(0.7,1.3) #Adjust random range. If code runs forever increase range.
nspacesx = np.ceil((xmax-xmin)/spacing)
nspacesy = np.ceil((ymax-ymin)/spacing)
randomstart = [xmin-spacing*np.random.random(),ymin-spacing*np.random.random()]
xlist=[randomstart[0]+(x*spacing) for x in range(int(nspacesx)+1)]
ylist=[randomstart[1]+(y*spacing) for y in range(int(nspacesy)+1)]
for x,y in product(xlist,ylist):
feat = QgsFeature()
feat.setGeometry(QgsPoint(x,y))
featlist.append(feat)
points_inside = [1 if f.geometry().intersects(geom.buffer(-5,10)) else 0 for f in featlist] #I use negative buffer to prevent points to end up near polygon edge. You need to adjust -5 m or remove the buffering.
#print(sum(points_inside))
if sum(points_inside)==npoints:
featlist = [p for p,i in zip(featlist, points_inside) if i==1]
success = 1

I’m moving from theory to practice, the initial data is a layer (table) with the name adm_polygons,

run the script:

WITH tbla AS (
WITH atbl AS (SELECT id, (ST_ExteriorRing(((ST_Dump(geom)).geom))) geom FROM adm_polygons),
intervals AS (SELECT generate_series (0, 8) as steps)
SELECT steps AS stp, ST_LineInterpolatePoint(geom, steps/(SELECT count(steps)::float-1 FROM intervals)) geom FROM atbl, intervals GROUP BY id, intervals.steps, geom),
tblb AS (SELECT (ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom geom FROM tbla),
tblc AS (SELECT ST_Intersection (a.geom, b.geom) geom FROM tblb a JOIN adm_polygons b ON ST_Intersects (a.geom, b.geom))
SELECT ST_PointonSurface(geom) geom FROM tblc;

See the figure below for the result

Warning: additional points may appear near the borders of the polygons!

This script is called - ST_EqualNumberPointsInPolygons...

• Is it only applicable in PostGIS? Sep 6, 2019 at 18:22
• @Taras, I think this is possible both in python and using SQL for QGIS, but it is advisable to check this ... Sep 6, 2019 at 20:43