Moving rotated polygons horizontally to avoid overlap

Question

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:

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.

@Erik It would be the same issue. If the polygons are extremely sparse this is easy enough to code, but increasing feature density quickly becomes an NP-Hard problem. — Vince, Jul 29 '20 at 12:02
@Erik - If I have a horizontal row of rectangles and rotate around corners, the row will no longer be horizontal. — user35594, Jul 29 '20 at 12:02
Isn't this just going to be cos(angle)*width/2 or maybe sin? — Ian Turton, Jul 29 '20 at 12:17
@IanTurton - The answer given in the Math link suggests it's more involved but I'm not a Math expert. — user35594, Jul 29 '20 at 12:21
It isn't really more involved. The only trick is that you need to shift each part in half the delta-X in different directions. You could solve this empirically by constructing a line from the corner point in the other feature along a constant Y and intersecting it with the polygon, the taking half the distance, but @IanTurton already gave the answer -- 0.5cos(rotation)base_width — Vince, Jul 29 '20 at 14:46

xunilk · Accepted Answer · 2020-07-31T14:23:58.923

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:

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:

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

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.

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:

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:

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.

This is exactly what I was looking for! Thank you very much!!! — user35594, Jul 30 '20 at 14:57
Only issue is that this works for 2 rectangles. If I have more, I assume I would need another for loop using range to find intersect1 between the [0] and [1] and then [1] and [2] etc? — user35594, Jul 30 '20 at 15:11
Yes, it is correct. If you have some issues for implementing that loop comment again. — xunilk, Jul 30 '20 at 16:07
By the way, for same row in a grid it is not difficult to implement the loop because distance is always the same. However, for a different row you need an additional vertical movement for matching completely. It needs to be calculated. — xunilk, Jul 30 '20 at 16:28
The translation part of the code doesn't seem to work on memory layers. The features get rotated but are not moved to the new position. Do you have any idea why? — user35594, Jul 31 '20 at 09:44
It works perfectly on memory layers as you can see in my Editing Note. I used in this case a complete grid. However, as it was commented previously, for a different row you need an additional vertical movement for matching completely. It needs to be calculated. — xunilk, Jul 31 '20 at 14:09
You have gone above and beyond with your help! Thank you very much for this! — user35594, Aug 02 '20 at 12:59

Moving rotated polygons horizontally to avoid overlap

Change the layer features rotation in one go

Refresh to see the changes

1 Answers1