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

  • Could you generate a regular grid for the extent of your polygon layer and then clip it using the polygons, or is that not suitable?
    – Matt
    Commented Feb 19, 2022 at 0:04
  • 1
    What do you mean by "as well as from the edge of the polygon"? Commented Feb 19, 2022 at 0:23
  • 1
    I guess you need to provide an example image of what you're aiming for.
    – Erik
    Commented Feb 19, 2022 at 8:33
  • Welcome to Geographic Information Systems! Welcome to GIS SE! We're a little different from other sites; this isn't a discussion forum but a Q&A site. Your questions should as much as possible describe not just what you want to do, but precisely what you have tried and where you are stuck trying that. Please check out our short tour for more about how the site works
    – Ian Turton
    Commented Feb 19, 2022 at 10:22

2 Answers 2


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
        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]
            success = 1


enter image description here


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):
            f = QgsFeature()

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

enter image description here

enter image description here

Your Answer

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

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