[GIS] How to generate a parameter defined rotated grid in QGIS, GDAL/OGR and/or Python

ogrpyqgisqgisvector-grid

Following this question which answers partially to my question, I'd like to go further and find out how I could generate a parameter-directed grid in QGIS, with origin X,Y coordinates, Direction/Angle, number of cells, size of cells.

My problem initially came from a need to add grid squares easily to an existing rotated grid. I could do it by hand, but using a script, I thought might later be able to pass on some attributes while creating the additional grid.

All resources I find explain how to create an horizontal grid with square or rectangle like cells, but I didn't find yet how to generate an inclined grid (with an angle and a Origin X/Y), ideally in the form of a QGIS python script with parameters in the processing toolbox for easy sharing.

I think it's just a matter of adaptation of the following script found out by @Gene to include some angle consideration. I'll try to answer my own question when I find the right formulas but if any of you have already some resource I could use or if you would like to take up the challenge… 😉

#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])

#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

Best Answer

With an already created polygon grid, next code works perfectly for my case [to rotate a grid, with 50 features, 70 degrees around QgsPoint(416047.568418, 4443411.85665)]. For particular cases, it's only necessary to change rotation axis and rotation angle. 

from math import sin, cos, pi

def rotatePoints(axis, r_angle, points):

    points_rot = []

    for point in points:
        tr = QgsPoint( point.x() - axis.x(), point.y() - axis.y() )
        xrot = tr.x()*cos(r_angle*(pi/180)) - tr.y()*sin(r_angle*(pi/180))  
        yrot = tr.x()*sin(r_angle*(pi/180)) + tr.y()*cos(r_angle*(pi/180)) 
        xnew = xrot + axis.x()
        ynew = yrot + axis.y()
        points_rot.append(QgsPoint(xnew,ynew))

    return points_rot

def createPolygonMemoryLayer(epsg, pol_rot):

    uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

    mem_layer = QgsVectorLayer(uri,
                               "rotated_grid",
                               "memory")

    QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

    prov = mem_layer.dataProvider()

    n = len(pol_rot)

    #Set geometry
    for i in range(n):
        feature = QgsFeature()
        feature.setGeometry(QgsGeometry.fromPolygon(pol_rot[i]))
        #set attributes values
        feature.setAttributes([i])
        prov.addFeatures([feature])

#Code starts here
layer = iface.activeLayer()

crs = layer.crs()
epsg = crs.postgisSrid()

geom = [ feature.geometry().asPolygon() for feature in layer.getFeatures() ]

points = geom

n = layer.featureCount()

axis = QgsPoint(416047.568418, 4443411.85665)

r_angle = 70

pol_rot = [ [ rotatePoints(axis, r_angle, points[i][0]) ]  for i in range(n) ]

createPolygonMemoryLayer(epsg, pol_rot)

After running the code at the Python Console of QGIS, I got:

enter image description here

Related Question