3

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?

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

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

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)

enter image description here

  • 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 at 6:28
  • 1
    You are too nice :). Im sure it can be done with less code! – BERA Sep 6 at 6:56
1

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 ta AS ( WITH at AS (SELECT id, (ST_ExteriorRing(((ST_Dump(geom)).geom))) as 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)) as geom FROM at, intervals GROUP BY id, intervals.steps, geom), tb as (SELECT (ST_Dump(ST_VoronoiPolygons(ST_Collect(geom)))).geom bgeom FROM ta), tc as (SELECT ST_Intersection (a.bgeom, b.geom) cgeom FROM tb a JOIN adm_polygons b ON ST_Intersects (a.bgeom,b.geom)) SELECT ST_PointonSurface(cgeom) geom FROM tc;

See the figure below for the result enter image description here

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

I count on your success,

With respect,

OS ...

P.S. I will name this tool - ST_EqualNumberPointsInPolygons.

  • Is it only applicable in PostGIS? – Taras Sep 6 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 Sep 6 at 20:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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