Alter Minimum Bounding Box Algorithm

Question

I'm trying to create an algorithm similar to minimum bounding box (though it may end up looking nothing like it). In this case the angle will be passed as a parameter and given the angle I need the smallest rectangle covering all my points/polygons. So far my line of thought is to find the center of my points (centroid algorithm), and from there create two parallel lines with same angle as paramter angle, and two more lines perpendicular to them. Then using iteration move these lines outwards (in opposite directions) until they contain all points. Also doesn't have to be exact minimum bounding box, approximate works (I guess would depend on size of each iteration step).

Here is my code so far. I have already dissolved all my polygons into one. I then take a convex hull to reduce vertices. I then put all vertices in a list - not sure if this helps yet...

a = layer.getFeatures()
for feat in a:
 geom = feat.geometry()
a = geom.convexHull()
vertexId = QgsVertexId()
vertices = []
b = a.constGet().nextVertex(vertexId)
while b[0]:
 vertices.append(b[1])
 b = a.constGet().nextVertex(vertexId)

Notes: At some point I need to pass the angle of box. I am using QGIS 3, and need to create this in Python. Layer 'layer' has one geometry, the dissolved polygon of all other polygons - maybe iteration isn't needed to access it..

Please let me know if I should be passing on more details/info.

• FelixIP: This is straight forward task. Rotate vertices of convex hull using standard equations, stackoverflow.com/questions/20104611/… Compute minX, minY etc. Unrotate and create rectangle of 4 xy pairs.

Answer 1

Here is the complete code. It contains too many lines (much more than needed for sure) but it works. Now you can clean it if you like.

In resume the algorithm calculates the maximum distance betweeen parallel lines that have the slope defined by the rotation parameter and pass though the points. For each point there will be created a 'horizontal' and 'vertical' line. This names are just orientative as they are defined at position 0 (rotation = 0). So, for each external point there will be created this 2 posible lines and then, iteratively, the poligon will be created based on the 4 external, or said in other way, where the distance of parallel lines are maximum.

One last thing: it is made to be used in QGIS 3.8 with grass.

enter image description here

from PyQt5.QtCore import *
from qgis.core import *
from qgis.gui import *
from processing.tools import *
from qgis.utils import iface
import qgis.utils, os, glob, processing, string, time, shutil, ogr
#PARAMETERS AND LAYERS
rotation = 45 #use any value between 0 and <90 #90 would make a mess
layer1 = iface.activeLayer() # Load the layer (from active)
crs = layer1.crs().authid() #get crs
#----------------------------------------------------------------------------------------
#LINE EQUATIONS
''' 
BASIC LINE EQUATIONS
y = ax + b
a = (y2 - y1) / (x2 - x1)
b = y1 - a * x1
Distance = (| a*x1 + b*y1 + c |) / (sqrt( a*a + b*b))# Function to find straight distance betweeen line and point 
'''
# slope from angle
def sfa (a):
 return round(math.tan(math.radians(a)),12) #round to avoid problems with horizontal and vertical
# angle from slope (not used)
def afs (s):
 return (math.atan(s) / math.pi) * 180
# Function to find distance 
def shortest_distance(x1, y1, a, b, c): 
 d = round(abs((a * x1 + b * y1 + c)) / (math.sqrt(a * a + b * b)) , 12)
 return d
# Function to find interception between lines
def cross(a1,b1,a2,b2):
 x = (b2-b1) / (a1-a2)
 y = a1 * x + b1
 return (x,y)
#----------------------------------------------------------------------------------------
# GET LIST OF POINTS TO ITERATE
# Calculate convexhull to reduce the iterations between point
# This avoid calculations on 'internal' points
# process of minimum bounding geometry convexHull
MBG = processing.run("qgis:minimumboundinggeometry", {'INPUT': layer1,'FIELD':None,'TYPE':3,'OUTPUT':'TEMPORARY_OUTPUT'})
# Get vertex of MBG
MBGp = processing.run("native:extractvertices", {'INPUT':MBG['OUTPUT'],'OUTPUT':'TEMPORARY_OUTPUT'})
plist = list(MBGp['OUTPUT'].getFeatures())
lp = list()
for p in plist:
 geom = p.geometry()
 a = geom.asPoint()
 point = (a[0],a[1])
 lp.append(point)
#----------------------------------------------------------------------------------------
# PROCESS
# compare hdist and v dist betweeen each pair of point and get the most distant lines
hdist_max = 0
vdist_max = 0
index = list(range(0,len(lp))) #iteration index
bl = ['ah1','bh1','av1','bv1','ah2','bh2','av2','bv2'] #polygon lines defined by 8 parameters see below
for i in index[:-1]:
 print('i'+str(i))
 for t in index[i+1:]:
 print('t'+str(t))
 x1 = lp[i][0] #; print('x1: {}', x1)
 y1 = lp[i][1] #; print('y1: {}', y1)
 x2 = lp[t][0] #; print('x2: {}', x2)
 y2 = lp[t][1] #; print('y2: {}', y2)
 #h1 equation
 ah1 = sfa(rotation)
 bh1 = y1 - ah1 * x1
 #v1 equation
 av1 = sfa(rotation + 90) #remember that just the horizontal is the reference at 0 rotation
 bv1 = y1 - av1 * x1 
 #h2 equation
 ah2 = sfa(rotation)
 bh2 = y2 - ah2 * x2
 #v2 equation
 av2 = sfa(rotation + 90) #remember that just the horizontal is the reference
 bv2 = y2 - av2 * x2 
 # H dist
 hdist = shortest_distance(x1, y1, ah2, -1, bh2)
 vdist = shortest_distance(x1, y1, av2, -1, bv2)
 if hdist > hdist_max:
 bl[0] = ah1
 bl[1] = bh1
 bl[4] = ah2
 bl[5] = bh2
 hdist_max = hdist #update max hdist
 if vdist > vdist_max:
 bl[2] = av1
 bl[3] = bv1
 bl[6] = av2
 bl[7] = bv2
 vdist_max = vdist #update max vdist
print("Max perpendicular distance betweeen 'horizontal lines' is",hdist_max, ' m')
print("Max perpendicular distance betweeen 'verticallines' is",vdist_max, ' m')
#------------------------------------------------------------------------------------------
# GET 4 COORDS FROM BOUNDINGLINES bl
# using the slope and intercept from boundinglines can we now calculate the 4 corners of the rotated polygon
H1V1 = cross(bl[0],bl[1],bl[2],bl[3]) # H1V1
H1V2 = cross(bl[0],bl[1],bl[6],bl[7]) # H1V2
H2V1 = cross(bl[4],bl[5],bl[2],bl[3]) # H2V1
H2V2 = cross(bl[4],bl[5],bl[6],bl[7]) # H2V2
# SORT POINTS CLOCKWISE AND CREATE QgsPointXY for polygon
clist = [H1V1,H1V2,H2V1,H2V2]
points=[]
points.append(sorted(clist, key=lambda e: (e[1], e[0]))[0]); clist.remove(points[0]) #minX and minY
points.append(sorted(clist, key=lambda e: (e[0], e[1]))[0]); clist.remove(points[1]) #minY and minX
points.append(sorted(clist, key=lambda e: (e[1]), reverse=True)[0]); clist.remove(points[2]) #maxY
points.append(clist[0]) #remaining
p=[]
for i in points:
 p.append(QgsPointXY(i[0],i[1]))
print('Coords of the polygon: ',p)
#------------------------------------------------------------------------------------------
#CREATE ROTATED BOUNDING BOX FROM THESE POINTS
layer = QgsVectorLayer(str('Polygon?crs='+crs), 'polygon' , 'memory')
prov = layer.dataProvider()
feat = QgsFeature()
feat.setGeometry(QgsGeometry.fromPolygonXY([p]))
prov.addFeatures([feat])
layer.updateExtents()
QgsProject.instance().addMapLayers([layer])