TransWikia.com

Moving rotated polygons horizontally to avoid overlap

Geographic Information Systems Asked on May 1, 2021

To rotate rectangular polygons, I use the following (courtesy from this answer):

rotation = 50
vlayer = iface.activeLayer()
provider = vlayer.dataProvider()

couples_id_geom = []
for feature in vlayer.getFeatures():
    geom = feature.geometry()
    centroid = feature.geometry().centroid().asPoint()
    geom.rotate(rotation, centroid)
    # accumulate args to avoid rotation feature by feature
    couples_id_geom.append([feature.id(), geom])

# Change the layer features rotation in one go
provider.changeGeometryValues({
  couple_id_geom[0]: couple_id_geom[1] for couple_id_geom in couples_id_geom
})

# Refresh to see the changes
vlayer.triggerRepaint()

However, when rotating adjacent polygons, overlaps occur. For example at 50 degrees:

Rectangles

How can I set it like the last pair of rectangles shown in the image so that when the polygons are rotated, they are moved at a minimum horizontal distance to avoid overlap and print the distance required?

EDIT

I asked a question on Math SE about calculating the minimum horizontal distance but I am unsure how to program this into Python.

One Answer

I tried to use derived formula in Math SE answer but it was not possible. I think there is an undetected error in its demonstration. However, you have an obvious answer in your image:

enter image description here

Distance for moving blue rectangle (assuming height > width), independent of any rotation angle, is always X. It is easily calculated as intersect length of green dotted line (by using line between polygons centroid) with intersection area of two polygons. Following code can be used.

rotation = 50
vlayer = iface.activeLayer()

feats = [ feat for feat in vlayer.getFeatures() ]

centroids = [ feat.geometry().centroid().asPoint() for feat in feats ]

new_geom = []

for feat in feats:
    geom = feat.geometry()
    geom.rotate(rotation, feat.geometry().centroid().asPoint())
    new_geom.append(geom.asWkt())

intersect1 = QgsGeometry.fromWkt(new_geom[0]).intersection(QgsGeometry.fromWkt(new_geom[1]))

line = QgsGeometry.fromPolylineXY(centroids)

distance = intersect1.intersection(line).length()

new_layer = [QgsGeometry.fromWkt(new_geom[1])]

geom = QgsGeometry.fromWkt(new_geom[0])

geom.translate(distance, 0, 0, 0)

new_layer.append(geom)

epsg = vlayer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon',
                           'memory')

prov = mem_layer.dataProvider()

new_feats = [ QgsFeature() for i in range(len(feats)) ]

for i, feat in enumerate(new_feats):
    feat.setAttributes([i])
    feat.setGeometry(new_layer[i])

prov.addFeatures(new_feats)

QgsProject.instance().addMapLayer(mem_layer)

After creating two adjacent polygons as follows:

enter image description here

above code was running in Python Console of QGIS 3 with following result (as memory layer):

enter image description here

When rotation angle was changed for 90º it was obtained below result. I tried it out with different angle values and it also worked as expected.

enter image description here

When height < width, you have to move blue rectangle (in your above image) in opposite sense (sign of distance would be negative).

Editing Note:

By using a complete grid as in following image:

enter image description here

the code, with a different algorithm, is as follows:

rotation = 50
vlayer = iface.activeLayer()

feats = [ feat for feat in vlayer.getFeatures() ]

centroids = [ feat.geometry().centroid().asPoint() for feat in feats ]

new_geom = []

for i in range(2): #for calculating translation distance
    geom = feats[i].geometry()
    geom.rotate(rotation, feats[i].geometry().centroid().asPoint())
    new_geom.append(geom.asWkt())

intersect1 = QgsGeometry.fromWkt(new_geom[0]).intersection(QgsGeometry.fromWkt(new_geom[1]))

line = QgsGeometry.fromPolylineXY(centroids)

distance = intersect1.intersection(line).length()

geom_rot = []

for i in range(len(feats)):
    geom = feats[i].geometry()
    geom.rotate(rotation, feats[i].geometry().centroid().asPoint())
    geom_rot.append(geom.asWkt())

cols = 4  #you need in this case column number

geom_trans = []

for i, item in enumerate(geom_rot):
    mod = i%cols
    geom = QgsGeometry.fromWkt(geom_rot[i])
    if mod == 0:
        k = 1
        geom_trans.append(geom)
    else:
        g = geom
        g.translate(k*distance, 0, 0, 0)
        geom_trans.append(g)
        k += 1

epsg = vlayer.crs().postgisSrid()

uri = "Polygon?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'polygon',
                           'memory')

prov = mem_layer.dataProvider()

new_feats = [ QgsFeature() for i in range(len(geom_trans)) ]

for i, feat in enumerate(new_feats):
    feat.setAttributes([i])
    feat.setGeometry(geom_trans[i])

prov.addFeatures(new_feats)

QgsProject.instance().addMapLayer(mem_layer)

After running above code, obtained result (for rotation = 50) is as expected:

enter image description here

However, it appears a Y displacement, as it was commented previously, and it must be corrected for complete matching but, it would be another question.

Correct answer by xunilk on May 1, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP