3

I want to create a rotated grid inside a polygon using a starting point and a rotation angle. I have started with the python script found here.

I have altered the original script to accommodate changes to QGIS 3.0 and introduced changes for the rotation and the loops. Obviously, something is amiss since the script is getting stuck in an infinite loop and I have to crash QGIS to stop.

Any idea where I went wrong?

Here is the script:

from qgis.core import *
from qgis.utils import iface
import math

# Set the spacing
spacing = 2.5
rot = math.radians(35.392)
# Set the inset
inset_x = 0
inset_y = 0

# Get the Coordinate Reference System and the extent from the loaded layer
layer = iface.activeLayer() # load the layer as you want
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()

# Create the coordinates of the points in the grid
points = []

for ft in layer.getFeatures():
    feat_geom = ft.geometry()
    # Set the extent of the new layer
    #xmin = ext.xMinimum() + inset_x
    xmin = 755358.66
    xmax = 755687.28
    ymin = ext.yMinimum()
    ymax = 3726718.02 # - inset_y
    #y = ymax
    #y = ymax
    xstart = xmax
    while xstart >= xmin:
        y = ymax
        x = xstart
        while y >= ymin:
            geom = QgsGeometry().fromPointXY(QgsPointXY(x, y))
            feat = QgsFeature()
            point = QgsPointXY(x,y)
            feat.setGeometry(QgsGeometry.fromPointXY(point))
            if feat_geom.contains(feat.geometry()):
                points.append(feat)
            x = x + (spacing * math.cos(rot))
            y = y - (spacing * math.sin(rot))
        xstart = xstart - spacing * math.cos(rot)

prov.addFeatures(points)
points_layer.updateExtents()

# Add the layer to the map
QgsProject.instance().addMapLayer(points_layer)
  • The line y = y - (spacing * math.sin(rot)) is inside of second while. I think that it should be out. – xunilk Nov 8 at 0:17
  • Moving that line out of the second loop to the first one will create lines parallel to the Y axis, which is not what I want. Say I want the lines rotated at a 45 degrees angle from the Y axis. – Techie_Gus Nov 8 at 6:43
  • since math.sin(rot) is a negative value your y value will never be less then ymin. Change y = y - (spacing * math.sin(rot)) to y = y + (spacing * math.sin(rot)) – eurojam Nov 8 at 8:52
  • @eurojam I'm a bit confused as to why math.sin(rot) is a negative value if rot is less than 90 degrees. – Techie_Gus Nov 8 at 9:18
  • 1
    @Techie_Gus: yes you are right...i've entered sin(30) into google...that was a mistake...so ignore my comment – eurojam Nov 8 at 9:51
2

If you want rotated points you also need to use rotate QgsGeometry method. Following code uses it.

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures()]

points = []

spacing_y = 5000
spacing_x = 5000

for feat in feats:
    centroid = feat.geometry().centroid().asPoint()
    extent = feat.geometry().boundingBox()
    xmin, ymin, xmax, ymax = extent.toRectF().getCoords()

    rows = int((ymax - ymin)/spacing_y)
    cols = int((xmax - xmin)/spacing_x)

    x = xmin
    y = ymax

    geom_feat = feat.geometry()

    for i in range(rows+1):
        for j in range(cols+1):
            pt = QgsPointXY(x,y)
            tmp_pt = QgsGeometry.fromPointXY(pt)
            tmp_pt.rotate(-45, centroid)
            if tmp_pt.within(geom_feat):
                points.append(tmp_pt.asPoint())
            x += spacing_x
        x = xmin
        y -= spacing_y

epsg = layer.crs().postgisSrid()

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

mem_layer = QgsVectorLayer(uri,
                           'point',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(points)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPointXY(points[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

After running above code with layer of following image, I got points rotated 45 degrees around centroid of feature layer.

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.