4

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

  • Where have you searched this far? Not found anything useful, or? – user30184 Oct 18 '17 at 15:42
  • Did you read my answer? Did it help you? – mgri Nov 15 '17 at 20:53
6

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):

enter image description here

7

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.

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.