Modifying polygons to be more rectangular using PyQGIS 3

Question

I'm trying to re-adapt an adaptation of this code but I'm getting the error stated above.

I tried reading the docs, but I honestly didn't understand it.

The input is a vector that contains many polygons and I want to make them more rectangular.

def slope(x1, y1, x2, y2):
        return (y2 - y1) / ((x2 - x1) + 0.00001) #I was getting DivisionByZero error
def intercept(y, slope, x): 
        return y - slope * x
def perpindicular(slope):
        return -1 / slope
pth = r"C:\Users\whatever...myVector.gpkg"
lyr = QgsVectorLayer(pth, "idk what i should put here", "ogr")
epsg = lyr.crs().postgisSrid()
uri = "Polygon?crs=epsg:{}&field=id:integer&index=yes".format(epsg)
rectangle = QgsVectorLayer(uri, 'Rectangles', 'memory') 
p = rectangle.dataProvider()
for f in lyr.getFeatures():
    g = f.geometry()
    xmin, ymin, xmax, ymax = g.boundingBox().toRectF().getCoords()
    pts = f.geometry().asPolygon()[0]
for i in range(len(pts)-1): 
    if pts[i][1] == ymax and pts[i+1][1] < pts[i][1]: 
        idx = i 
    if pts[i][1] == ymax and pts[i-1][1] < pts[i][1]: 
        idx = i-1 
r = [] 
x1 = pts[idx][0]

y1 = pts[idx][1] 
r.append(QgsPointXY(x1,y1)) 
x2 = pts[idx+1][0]

y2 = pts[idx+1][1] 
r.append(QgsPointXY(x2,y2))
s1 = slope(x1, y1, x2, y2) 
i1 = intercept(y1, s1, x1)
x3 = pts[idx+2][0]

y3 = pts[idx+2][1] 
i2 = intercept(y3, s1, x3) 
s3 = perpindicular(s1)
i3 = intercept(y2, s3, x2) 
x4 = (i3 - i2)/(s1 - s3)
y4 = s3 * x4 + i3 
r.append(QgsPointXY(x4, y4)) 
s4 = perpindicular(s1) 
i4 = intercept(y1, s4, x1) 
x5 = (i4 - i2)/(s1 - s4) 
y5 = s4 * x5 + i4 
r.extend([QgsPointXY(x5, y5),QgsPointXY(x1, y1)]) 
poly = [] 
poly.append(r) 
g = QgsGeometry.fromPolygonXY(poly) 
ft = QgsFeature() 
ft.setAttributes([i]) 
ft.setGeometry(g) 
p.addFeatures([ft])
QgsProject.instance().addMapLayer([rectangle, lyr])

it generates two outputs, but the algorithm doesn't really work as intended. Bummer. — willsbit, Sep 23 '20 at 20:37
The same thing from this link https://gis.stackexchange.com/questions/212003/modifying-polygons-to-be-more-rectangular-using-pyqgis/212325#212325 But that was made for qgis 2, so I was trying to adapt it — willsbit, Sep 24 '20 at 12:56

score 3 · Accepted Answer · answered Sep 24 '20 at 13:24

Use the following script. In QGIS 3, you should use QgsProject, QgsPointXY, fromPolygonXY instead of QgsMapLayerRegistry, QgsPoint, fromPolygon respectively.

layer = iface.activeLayer()
feats = [ feat for feat in layer.getFeatures() ]
n = len(feats)
crs = layer.crs()
epsg = crs.postgisSrid()
uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"
mem_layer = QgsVectorLayer(uri, 'rectangle', 'memory')
prov = mem_layer.dataProvider()
for feature in feats:
    geom = feature.geometry()
    xmin, ymin, xmax, ymax = geom.boundingBox().toRectF().getCoords()
    points = feature.geometry().asPolygon()[0]
for i in range(len(points)-1):
    if points[i][1] == ymax and points[i+1][1] &lt; points[i][1]:
        idx = i
    if points[i][1] == ymax and points[i-1][1] &lt; points[i][1]:
        idx = i-1

rectangle = []

#x,y coordinates of first point
x1 = points[idx][0] 
y1 = points[idx][1]
rectangle.append(QgsPointXY(x1,y1))

#x,y coordinates of second point
x2 = points[idx+1][0] 
y2 = points[idx+1][1]
rectangle.append(QgsPointXY(x2,y2))

#slope for first line
m1 = (y2 - y1) / (x2 - x1)
#intercept at origin for first line
int1 = y1 - m1 * x1
#slope for second line
m2 = m1
#x,y coordinates of third point
x3 = points[idx+2][0] 
y3 = points[idx+2][1]
#intercept at origin for second line
int2 = y3 - m2 * x3
#first perpendicular
m3 = -1/m1
#intercept at origin for second line
int3 = y2 - m3 * x2
#intersect point
x4 = (int3 - int2)/(m2 - m3)
y4 = m3*x4 + int3
rectangle.append(QgsPointXY(x4, y4))

#second perpendicular
m4 = -1/m1
#intercept at origin for second perpendicular
int4 = y1 - m4 * x1
#intersect point
x5 = (int4 - int2)/(m2 - m4)
y5 = m4*x5 + int4

rectangle.extend([QgsPointXY(x5, y5),QgsPointXY(x1, y1)])

polygon = []
polygon.append(rectangle)
geom = QgsGeometry.fromPolygonXY(polygon)
feat = QgsFeature()
feat.setAttributes([i])
feat.setGeometry(geom)
prov.addFeatures( [feat] )


QgsProject.instance().addMapLayer(mem_layer)

Result in QGIS 3:

Reference: Modifying polygons to be more rectangular using PyQGIS

score 0 · Answer 2 · answered Sep 24 '20 at 08:02

0

QgsProject.instance().addMapLayer([rectangle, lyr]) expects a single layer and you have added a list. You can use QgsProject.instance().addMapLayers([rectangle, lyr])

answered Sep 24 '20 at 08:02

Ian Turton

81,417
6
84
185

Modifying polygons to be more rectangular using PyQGIS 3

2 Answers2