Modifier l'algorithme de boîte de délimitation minimale

J'essaie de créer un algorithme similaire à la zone de délimitation minimale (bien qu'il puisse finir par ne rien ressembler). Dans ce cas, l'angle sera passé en paramètre et étant donné l'angle, j'ai besoin du plus petit rectangle couvrant tous mes points / polygones. Jusqu'à présent, ma ligne de pensée est de trouver le centre de mes points (algorithme centroïde), et de là créer deux lignes parallèles avec le même angle que l'angle paramétrique, et deux autres lignes perpendiculaires à celles-ci. Ensuite, en utilisant l'itération, déplacez ces lignes vers l'extérieur (dans des directions opposées) jusqu'à ce qu'elles contiennent tous les points. Il n'a pas non plus besoin d'être une zone de délimitation minimale exacte, des travaux approximatifs (je suppose que cela dépend de la taille de chaque étape d'itération).

Voici mon code jusqu'à présent. J'ai déjà dissous tous mes polygones en un seul. Je prends ensuite une coque convexe pour réduire les sommets. J'ai ensuite mis tous les sommets dans une liste - je ne sais pas si cela aide encore ...

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)

Remarques: À un moment donné, je dois passer l'angle de la boîte. J'utilise QGIS 3 et je dois le créer en Python. La couche «couche» a une géométrie, le polygone dissous de tous les autres polygones - peut-être qu'une itération n'est pas nécessaire pour y accéder.

Veuillez me faire savoir si je devrais transmettre plus de détails / informations.

Voici le code complet. Il contient trop de lignes (bien plus que nécessaire pour sûr) mais cela fonctionne. Vous pouvez maintenant le nettoyer si vous le souhaitez.

En résumé, l'algorithme calcule la distance maximale entre les lignes parallèles qui ont la pente définie par le paramètre de rotation et passent par les points. Pour chaque point, une ligne «horizontale» et «verticale» sera créée. Ces noms sont juste orientatifs car ils sont définis à la position 0 (rotation = 0). Ainsi, pour chaque point externe, il sera créé ces 2 lignes posibles puis, de manière itérative, le poligon sera créé sur la base des 4 externes, ou dit autrement, où la distance des lignes parallèles est maximale.

Une dernière chose: il est conçu pour être utilisé dans QGIS 3.8 avec de l'herbe.

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])