Creating set number of regularly spaced points inside of polygon in QGIS

Question

My application for this would be for making maps with pre-established plot points for cruising a stand of timber, but I could see this being used for other applications.

I am a little surprised that there are Random Points Inside Polygons and Regular Points under Vector > Research Tools, but no way to do Regular Points Inside Polygons.

I want to be able to input a number of points (or better yet, calculate this number based on the area of the polygon) and have the program generate a grid of points inside a polygon with the maximum possible spacing between individual points, as well as from the edge of the polygon.

I was hoping I could write an expression that would make this happen in the Geometry Generator, but I haven't been able to figure it out.

My starting point was the expression given in the answer to Halftone effect gridded points inside polygon with QGIS geometry generator

Although that expression returned NULL for me.

I'm using QGIS 3.16

Answer 1

I use this pyqgis code to create a fixed number of sample plots in forest stands/polygons. It will create a point layer in memory. Click/highlight your polygon layer, adjust the number of plots below and execute the code below. If you have polygons with very complex shapes you might want to add a counter to break the while loop or it can keep going for a long time / forever.

import numpy as np
from itertools import product

polylayer = iface.activeLayer() #Click/highlight your layer
npoints = 8 #Number of points in each polygon

pointlayer = QgsVectorLayer('Point?crs={0}'.format(polylayer.crs().authid()), 'point' , 'memory')
prov = pointlayer.dataProvider()

for poly in polylayer.getFeatures(): #For each polygon
    geom = poly.geometry()
    bbox = poly.geometry().boundingBox() #Create a bounding box
    xmin, xmax, ymin, ymax = bbox.xMinimum(),bbox.xMaximum(),bbox.yMinimum(),bbox.yMaximum() #Find max/min x/y coordinates
    success = 0
    while success == 0:
        featlist = []
        spacing = ((geom.area()/npoints)**0.5)*np.random.uniform(0.7,1.3) #Create random spacings based on polygon area
        nspacesx = np.ceil((xmax-xmin)/spacing) #Find out how many points you can fit in x direction
        nspacesy = np.ceil((ymax-ymin)/spacing)
        randomstart = [xmin-spacing*np.random.random(),ymin-spacing*np.random.random()] #Create a random start coordinate
        xlist=[randomstart[0]+(x*spacing) for x in range(int(nspacesx)+1)] #Calculate x coordinates 
        ylist=[randomstart[1]+(y*spacing) for y in range(int(nspacesy)+1)] #And y
        for x,y in product(xlist,ylist): #For each combination of x,y coordinate
            feat = QgsFeature()
            feat.setGeometry(QgsPoint(x,y)) #Create a point
            featlist.append(feat)
        points_inside = [1 if f.geometry().intersects(geom.buffer(-5,10)) else 0 for f in featlist] #Find the points inside the polygon (-5m from each polygon edge)
        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 2

I used the following script in my project. It creates a point grid (which is independent of the other point grids as seen in the first image) in the bounding box of the polygon, then adjusts the distance from the edge of the polygon (as seen in the second image).

import numpy as np

layer =  iface.activeLayer()
crs = layer.crs().authid()
plyr = QgsVectorLayer(f"Point?crs={crs}", "points", "memory")



m = 20 # interval

for feature in layer.getFeatures():
    geom = feature.geometry()
    
    bb = geom.boundingBox()
    min_x, min_y = bb.xMinimum(), bb.yMinimum()
    max_x, max_y = bb.xMaximum(), bb.yMaximum()
    dx = max_x - min_x
    dy = max_y - min_y
    start_x = min_x + (dx % m)/2
    start_y = min_y + (dy % m)/2
    end_x = max_x - (dx % m)/2 + m
    end_y = max_y - (dy % m)/2 + m
    X = np.arange(start_x, end_x, m)
    Y = np.arange(start_y, end_y, m)

    features = []
    for x in X:
        for y in Y:
            p = QgsPointXY(x, y)
            g = QgsGeometry.fromPointXY(p)
            
            if not geom.contains(p):
                continue
            
            f = QgsFeature()
            f.setGeometry(g)
            features.append(f)
            
    plyr.dataProvider().addFeatures(features)
    plyr.updateExtents()


QgsProject.instance().addMapLayer(plyr, True)