7

existing buildings

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 enter image description here

  • 1
    Have a search for 'orthogonalize' on this site - there are some other questions relating to this, none of which are really duplicates IMHO. – Simbamangu Sep 26 '16 at 16:25
  • How are your skills concerning python or sql? You could calculate the bounding box, the angle of your longest edge and then rotate the bounding box according to it. Also with PostGIS it would be relative easy. – Matte Sep 26 '16 at 16:38
  • Most of my data is in shapefile format so I would like to stay away from SQL. I can work with Python. – GreyHippo Sep 26 '16 at 17:01
  • OpenJump has a tool to create Minimum Bounding Rectangles, gis.stackexchange.com/questions/160259/…. – klewis Sep 26 '16 at 20:07
  • 2
    The OpenStreetMap editor JOSM also has a nice rectify tool. Unfortunately, it does not read shapefiles, so you have to convert from shp to osm and later backwards. Or look into the open source code how they did it. – AndreJ Sep 27 '16 at 17:46
4

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:

enter image description here

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

enter image description here

  • Does that method always end up with rectangles that have one side the same as one side of the original quadrilateral? – Spacedman Sep 28 '16 at 21:01
  • No. If you comment these lines in the code, if points[i][1] == ymax and points[i-1][1] < points[i][1] and idx = i-1, it is obtained this result: postimg.org/image/ao0zrelsr . This block is the "smart" part of code and it could be refined by using other bounding box parameters. – xunilk Sep 29 '16 at 8:33
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)
  • I get an error. name 'shapely' is not defined. – GreyHippo Sep 27 '16 at 17:22
  • I have installed shapely but not sure if qgis can see it. – GreyHippo Sep 27 '16 at 17:28
  • A sorry, i forgot an import (added in the post). Shapely is installed alongside Qgis nowadays i think. – Matte Sep 27 '16 at 17:29
  • See me edit above – GreyHippo Sep 27 '16 at 17:37
  • Well thats probably why i took the first edge for the angle. And i´m not sure what your first edge is. I assumed the one from lower left to lower right for the example code. The idea from Spacedman is therfore another good idea to get the angle for rotating and also not to use the bounding box as it is always horizontaly aligned an therefore creates to large polygons when the original polygon was further away from horizontal. – Matte Sep 27 '16 at 18:04
0

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....

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.