# 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 ... – Cyril Mikhalchenko Sep 2 '19 at 17:09
• A grid pattern of Points are ok? – BERA Sep 3 '19 at 5:13
• Grid pattern of points would be ok, too! :) Prerequisite is that there are 8 points in each polygon. – Stefan Sep 3 '19 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.

• I suppose that QGIS's graphical modeler might be a suitable solution – Taras Sep 6 '19 at 6:26
• Can I feed the graphical modeler with features one by one? – Stefan Sep 6 '19 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 – Taras Sep 6 '19 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

• You are too nice :). Im sure it can be done with less code! – BERA Sep 6 '19 at 6:56

Good day,

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!

Good luck to everyone :-),

Original solutions ...

This script is called - ST_EqualNumberPointsInPolygons...

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