How to generate regularly spaced points in QGIS using Python

In this task we will see how to generate a vector grid of regularly spaced points. The Processing algorithm Regular points does the same thing, but we want to do it using Python.

We assume to start from a layer that is already loaded in QGIS (if not, please see here) because we need to define an extent for creating the points.

For example, we will work on this raster layer:

raster

To create the points, we only need two parameters:

  1. a spacing variable which defines how far the points are in map units;
  2. an inset variable to determine how close to the edge of the layer the points start, in map units.

For example, in this task we will set a spacing of 10 m and an inset of 5 m.

If the loaded layer is named”layer”, we can code in this way:

from qgis.core import * 

# Set the spacing
spacing = 10
# Set the inset
inset = 5

# Get the Coordinate Reference System and the extent from the loaded layer
crs = layer.crs().toWkt()
ext=layer.extent()

# Create a new vector point layer
points_layer = QgsVectorLayer('Point?crs=' + crs, 'grid', "memory")
prov = points_layer.dataProvider()

# Set the extent of the new layer
xmin = ext.xMinimum() + inset
xmax = ext.xMaximum()
ymin = ext.yMinimum()
ymax = ext.yMaximum() - inset

# Create the coordinates of the points in the grid
points = []
y = ymax
while y >= ymin:
    x = xmin
    while x <= xmax:
        geom = QgsGeometry().fromPoint(QgsPoint(x, y))
        feat = QgsFeature()
        point = QgsPoint(x,y)
        feat.setGeometry(QgsGeometry.fromPoint(point))
        points.append(feat)
        x += spacing
    y = y - spacing

prov.addFeatures(points)
points_layer.updateExtents()

# Add the layer to the map
QgsMapLayerRegistry.instance().addMapLayer(points_layer)

Alternatively, we can obtain the same result by using the Numpy module:

from qgis.core import *
import numpy

# Set the spacing
spacing = 10
# Set the inset
inset = 5

# Get the Coordinate Reference System and the extent from the loaded layer
crs = layer.crs().toWkt()
ext=layer.extent()

# Create a new vector point layer
points_layer = QgsVectorLayer('Point?crs=' + crs, 'grid', "memory")
prov = points_layer.dataProvider()
# Set the extent of the new layer
xmin = ext.xMinimum() + inset
xmax = ext.xMaximum() - inset
ymin = ext.yMinimum() + inset
ymax = ext.yMaximum() - inset
# Create the coordinates of the points in the grid
coords = [(x,y) for x in (i for i in numpy.arange(xmin, xmax, spacing)) for y in (j for j in numpy.arange(ymin, ymax, spacing))]

# Add the coordinates in the new point layer
points = []
for x,y in coords:
    feat = QgsFeature()
    point = QgsPoint(x,y)
    feat.setGeometry(QgsGeometry.fromPoint(point))
    points.append(feat)
prov.addFeatures(points)
points_layer.updateExtents()

# Add the layer to the map
QgsMapLayerRegistry.instance().addMapLayer(points_layer)

Regardless of the method used, this will be the result:

points

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s