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

Screenshot of agricultural parcels with some points

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:

Screenshot of agricultural parcels with some points and Voronoi

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


3 Answers 3


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

test polygon

  1. Polygon to Lines

test polygon to lines

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

test polygon split at max length

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

test extract specific vertices

  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)

test Voronoi from vertices

  1. Intersect Voronoi and test polygon

test intersect Voronoi and polygon

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

test point on surface

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

result: test polygon with points

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
    – Taras
    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
    – Taras
    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()
        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.
        if sum(points_inside)==npoints:
            featlist = [p for p,i in zip(featlist, points_inside) if i==1]
            success = 1


enter image description here


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 enter image description here

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

This script is called - ST_EqualNumberPointsInPolygons...

  • Is it only applicable in PostGIS?
    – Taras
    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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.