Creating n randomly distributed points within given (multi)polygon with PyQGIS?

Question

Is there a way to create n randomly distributed point features within a given polygon using QGIS, PyQGIS and/or PostGIS?

My first approach was to implement something like (fantasy language)

n = 0
repeat until n = 20
    create point feature within polygons boundingbox
    if point within polgon
        increment n
        keep point

It's easy to create randomly distributed points within a rectangle, but depending on the shape of the polygon, this algorithm may create a huge amount of points that are not within the polygon and may run a while, if it terminates at all...

Answer 1

If there isn't the need to use PyQGIS, maybe the Random points inside polygons (fixed) algorithm from the Processing Toolbox should return what you are looking for: it allows to specify the number of the points or the density as a parameter (in addition to a minimum distance among them).

Instead, using PyQGIS, you may run this simple code from the Python Console:

import random
from qgis.PyQt.QtCore import QVariant

layer=iface.activeLayer()
crs = layer.crs().toWkt()

# Create the output layer
outLayer = QgsVectorLayer('Point?crs='+ crs, 'random_points' , 'memory')
prov = outLayer.dataProvider()
prov.addAttributes([QgsField('ID', QVariant.Int, '', 10, 0)])
outLayer.updateFields()


ext=layer.extent()
xmin = ext.xMinimum()
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum()

points = 50 # set as you want

first = True
for feat in layer.getFeatures():
    if first:
        outFeat = QgsFeature()
        outGeom = QgsGeometry(feat.geometry())
        first = False
    else:
        outGeom = outGeom.combine(feat.geometry())
outFeat.setGeometry(outGeom)

id = 0
p = 1
while p <= points:
    x_coord = random.uniform(xmin, xmax)
    y_coord = random.uniform(ymin, ymax)
    pt = QgsPoint(x_coord, y_coord)
    tmp_geom = QgsGeometry.fromPoint(pt)
    if tmp_geom.intersects(outFeat.geometry()):
        outGeom = QgsFeature()
        outGeom.setGeometry(tmp_geom)
        outGeom.setAttributes([id])
        prov.addFeatures([outGeom])
        id += 1
        p += 1

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayer(outLayer)

It will create a new memory layer that stores as many points as specified (there isn't any check for the minimum distance between points or a control about the density, but this is not necessary for answering the question).

For example (some points seem to be outside the polygons but, actually, it depends on the width used when rendering them):

Answer 2

The ST_GeneratePoints function in PostGIS does this quite well.

You simply provide it with a geometry and the desired number of points, for example:

SELECT ST_GeneratePoints(geom, 100);

It's implemented following the same approach as your pseudocode, but it prepares (internally indexes) the polygon so that the if point within polygon tests are fast.

Answer 3

Since QGIS 3.10 QgsGeometry has the function randomPointsInPolygon (https://qgis.org/pyqgis/3.10/core/QgsGeometry.html#qgis.core.QgsGeometry.randomPointsInPolygon). Parameters available in pyqgis are a) Number of points to generate and b) Random seed.