How to generate a vector grid of polygons in QGIS using Python

In this task we will use python scripting to generate a vector grid of polygons (each one of them will have specific attributes) and save it as a memory layer. The Processing algorithm Vector grid does the same thing, but we want to do it using Python.

We need some input parameters:

  1. the extent of the grid that will be generated;
  2. the horizontal spacing, i.e. the horizontal side length for the features in the grid;
  3. the vertical spacing, i.e. the vertical side length for the features in the grid.

We can use the following code (if you don’t know how to create a new script from Processing, please see here):

##Grid_extent=extent
##Horizontal_spacing=number 10
##Vertical_spacing=number 10

from qgis.core import *
from qgis.PyQt.QtCore import QVariant
from qgis.utils import iface

extent = Grid_extent.split(',')
(xmin, xmax, ymin, ymax) = (float(extent[0]), float(extent[1]), float(extent[2]), float(extent[3]))
hspacing = Horizontal_spacing
vspacing = Vertical_spacing

# Create the grid layer
crs = iface.mapCanvas().mapSettings().destinationCrs().toWkt()
vector_grid = QgsVectorLayer('Polygon?crs='+ crs, 'vector_grid' , 'memory')
prov = vector_grid.dataProvider()

# Add ids and coordinates fields
fields = QgsFields()
fields.append(QgsField('ID', QVariant.Int, '', 10, 0))
fields.append(QgsField('XMIN', QVariant.Double, '', 24, 6))
fields.append(QgsField('XMAX', QVariant.Double, '', 24, 6))
fields.append(QgsField('YMIN', QVariant.Double, '', 24, 6))
fields.append(QgsField('YMAX', QVariant.Double, '', 24, 6))
prov.addAttributes(fields)

# Generate the features for the vector grid
id = 0
y = ymax
while y >= ymin:
    x = xmin
    while x <= xmax:
        point1 = QgsPoint(x, y)
        point2 = QgsPoint(x + hspacing, y)
        point3 = QgsPoint(x + hspacing, y - vspacing)
        point4 = QgsPoint(x, y - vspacing)
        vertices = [point1, point2, point3, point4] # Vertices of the polygon for the current id
        inAttr = [id, x, x + hspacing, y - vspacing, y]
        feat = QgsFeature()
        feat.setGeometry(QgsGeometry().fromPolygon([vertices])) # Set geometry for the current id
        feat.setAttributes(inAttr) # Set attributes for the current id
        prov.addFeatures([feat])
        x = x + hspacing
        id += 1
    y = y - vspacing

# Update fields for the vector grid
vector_grid.updateFields()

# Add the layer to the Layers panel
QgsMapLayerRegistry.instance().addMapLayers([vector_grid])

In addition to the creation of the geometries, the script will automatically assign an identifier and the coordinates of the vertices for each feature in the attribute table.

As an example, we test the script on a dataset available from DIVA-GIS. Before running the script, please do these steps:

  1. first of all, download the administrative areas shapefile for France;
  2. load the shapefile in the Layers Panel of QGIS (in particular, load the shapefile with “adm0” as final part of the filename which represents the boundaries of the country itself);
  3. save the shapefile using a proper Projected Coordinate System Reference (in this task, we will use the EPSG:32631). This is a necessary operation because the original shapefile was in a Geographic Coordinate System Reference (i.e. the EPSG:4326).

Starting from this:

input

and setting the spacing of 50000 m (= 50km) for both sides of the polygons, this will be the result:

grid

If we change the opacity of the vector grid to show the map beneath, we obtain:

final

with this attribute table:

attributes

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