Create rotated point grid inside polygon

Question

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 08 '19 at 00: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 08 '19 at 06: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 08 '19 at 08: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 08 '19 at 09:18
@Techie_Gus: yes you are right...i've entered sin(30) into google...that was a mistake...so ignore my comment — eurojam, Nov 08 '19 at 09:51
If you want rotated points you also need to use rotate QgsGeometry method. I hope my answer helps. — xunilk, Nov 08 '19 at 14:37

score 2 · Answer 1 · answered Nov 08 '19 at 14:35

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.

Create rotated point grid inside polygon

1 Answers1