# 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+(x*spacing) for x in range(int(nspacesx)+1)]
ylist=[randomstart+(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

`````` • Holy shit...@Bera, I could not expect that you are that much skilled in PyQGIS as in ArcPy. Tremendous appreciation! Thank you – Taras Sep 6 '19 at 6:28
• 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