1

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
2

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

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.