Wednesday, 7 September 2016

Generating grid in PyQGIS without GUI?



Is there a QGIS Python binding for generating a grid programmatically without using the GUI?



Answer



From an existing layer:


enter image description here


1) you can use processing with the algorithm "qgis:creategrid" as in the "Hex grid from layer bounds" in Processing/Scripts/Example scripts. (look at Processing-Help / qgis / creategrid.rst)

squared grid


import processing
layer = canvas.layer(0)
cellsize = 1000
input = processing.getobject(layer.name())

centerx = (input.extent().xMinimum() + input.extent().xMaximum()) / 2
centery = (input.extent().yMinimum() + input.extent().yMaximum()) / 2
width = (input.extent().xMaximum() - input.extent().xMinimum())
height = (input.extent().yMaximum() - input.extent().yMinimum())
grid="/yourpath/grid.shp"
processing.runandload("qgis:creategrid", cellsize, cellsize, width, height, centerx, centery, 1, input.crs().authid(), grid)

enter image description here


rectangular grid


.....

cellsizex = 500
cellsizey = 1000
.....
processing.runandload("qgis:creategrid", cellsizex, cellsizey, width, height, centerx, centery, 1, input.crs().authid(), grid)

enter image description here


For other types of grids change the value in runalg (0,1,2,3 = Rectangle (line), Rectangle (polygon), Diamond (polygon), Hexagon (polygon)):


For Hexagonal grid


....
processing.runandload("qgis:creategrid", cellsize, cellsize, width, height, centerx, centery, 3, input.crs().authid(), grid)


enter image description here


2) You can adapt the script of ustroetz in Create a square grid polygon shapefile with python


part of my class to create a memory layer limited here to the geometry of Polygons:


class Crea_layer(object):
def __init__(self,name,type):
self.type=type
self.name = name
self.layer = QgsVectorLayer(self.type, self.name , "memory")
self.pr =self.layer.dataProvider()

def create_poly(self,points):
self.seg = QgsFeature()
self.seg.setGeometry(QgsGeometry.fromPolygon([points]))
self.pr.addFeatures( [self.seg] )
self.layer.updateExtents()
@property
def disp_layer(self):
QgsMapLayerRegistry.instance().addMapLayers([self.layer])

The script to create a squared grid layer from an existing layer:



from math import ceil
canvas= qgis.utils.iface.mapCanvas()
# first layer
layer = canvas.layer(0)
xmin,ymin,xmax,ymax = layer.extent().toRectF().getCoords()
gridWidth = 1000
gridHeight = 1000
rows = ceil((ymax-ymin)/gridHeight)
cols = ceil((xmax-xmin)/gridWidth)
ringXleftOrigin = xmin

ringXrightOrigin = xmin + gridWidth
ringYtopOrigin = ymax
ringYbottomOrigin = ymax-gridHeight
pol = Crea_layer("grid", "Polygon")
for i in range(int(cols)):
# reset envelope for rows
ringYtop = ringYtopOrigin
ringYbottom =ringYbottomOrigin
for j in range(int(rows)):
poly = [QgsPoint(ringXleftOrigin, ringYtop),QgsPoint(ringXrightOrigin, ringYtop),QgsPoint(ringXrightOrigin, ringYbottom),QgsPoint(ringXleftOrigin, ringYbottom),QgsPoint(ringXleftOrigin, ringYtop)]

pol.create_poly(poly)
ringYtop = ringYtop - gridHeight
ringYbottom = ringYbottom - gridHeight
ringXleftOrigin = ringXleftOrigin + gridWidth
ringXrightOrigin = ringXrightOrigin + gridWidth

pol.disp_layer

Result:


enter image description here



But you can modify the script using numpy in place of math ceil and int


import numpy as np
rows = (ymax-ymin)/gridHeight
cols = (xmax-xmin)/gridWidth
for i in np.arange(cols):
....
for j in np.arange(rows):
....

and if you don't understand xmin,ymin,xmax,ymax = layer.extent().toRectF().getCoords()



xmin = layer.extent().xMinimum()
xmax = layer.extent().xMaximum()
ymin = layer.extent().yMinimum()
ymax = layer.extent().yMaximum()
# = in one line
xmin,ymin,xmax,ymax = layer.extent().toRectF().getCoords()

and to go faster, you can even use the gridding functions of numpy


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...