2

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])
2
  • it generates two outputs, but the algorithm doesn't really work as intended. Bummer.
    – willsbit
    Commented Sep 23, 2020 at 20:37
  • The same thing from this link gis.stackexchange.com/questions/212003/… But that was made for qgis 2, so I was trying to adapt it
    – willsbit
    Commented Sep 24, 2020 at 12:56

2 Answers 2

3

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] < points[i][1]:
            idx = i
        if points[i][1] == ymax and points[i-1][1] < 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:

enter image description here

Reference: Modifying polygons to be more rectangular using PyQGIS

0

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

0

Not the answer you're looking for? Browse other questions tagged or ask your own question.