Modifying polygons to be more rectangular using PyQGIS

Question

I have some buildings that were quickly digitized from an aerial image. As you can see most of the polygons are not rectangular.

How can I modify the polygons to make them more rectangular?

Is there some equation that can solve a rectangle from a group of points?

I do not want to re-digitize the polygons. The polygons are supposed to represent buildings. I am using QGIS.

I tried code from @Matte but resulted in

Answer 1

Next code uses analytical geometry to change each polygon quasi rectangular in rectangular and it could be used:

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(QgsPoint(x1,y1))

    #x,y coordinates of second point
    x2 = points[idx+1][0] 
    y2 = points[idx+1][1]

    rectangle.append(QgsPoint(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(QgsPoint(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([QgsPoint(x5, y5),QgsPoint(x1, y1)])

    polygon = []

    polygon.append(rectangle)

    geom = QgsGeometry.fromPolygon(polygon)

    feat = QgsFeature()

    feat.setAttributes([i])
    feat.setGeometry(geom)
    prov.addFeatures( [feat] )

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

I tried it out with shapefile of next image:

After running the code at the Python Console of QGIS I got:

Answer 2

You can try to use the following approach that changes the geometry based on the bounding box and the angle of the first digitized edge. You can of course alter the angle with the one from the longest edge or something. Only works well when your polygons are near rectangular already (as it looks like).

Input in the console in Qgis and the layer needs to be selected and editable:

import shapely

from shapely import affinity

from shapely.wkb import loads

layer = qgis.utils.iface.activeLayer()

for feature in layer.getFeatures():

    azimuth = feature.geometry().vertexAt(0).azimuth(feature.geometry().vertexAt(1))

    bbox = QgsGeometry.fromRect(feature.geometry().boundingBox())

    input = loads(bbox.asWkb())

    shape = shapely.geometry.asShape(input)

    rotated = affinity.rotate(shape, azimuth-90.0)

    new_geom = QgsGeometry.fromWkt(rotated.wkt)

    layer.changeGeometry(feature.id(),new_geom)

Answer 3

For each 4-vertex polygon:

Compute an estimate of the centre of the rectangle, C, as the mean of the 4 vertices.

Compute an estimate of the diagonal length of the rectangle, R, as the mean of the four distances from the estimated centre to each point times two.

Compute an estimate of the rotation angle of the rectangle by taking the mean of the angle of the line between points 1 and 3 and between points 2 and 4. Or possibly the mean of (the mean of the angle between (1 and C) and (C and 3)) and (the mean of the angle between (2 and C) and (4 and C). Basically an estimate of the angle midway between the two diagonals.

Compute the rectangle aspect ratio by computing the average of the unsigned angles of a line from each point to the centre point and the rotation angle. So if the rectangle is square this is 45 degrees.

Compute all the four points from the knowledge of the centre, the diagonal length, the aspect ratio, and the rotation angle.

I think this is robust to any convex arrangement of four points - I don't know what it might do if given a concave polygon...

Any chance you can dump your sample data and I might try and code this up....