Insert n points equally distributed into many polygons

Question

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?

Answer 1

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

1. Polygon to Lines

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

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

4. 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)

5. Intersect Voronoi and test polygon

6. 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]
)

Answer 2

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]
            prov.addFeatures(featlist)
            success = 1

QgsProject.instance().addMapLayer(pointlayer)

Answer 3

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...