miércoles, 14 de agosto de 2019

Determinando ángulos de rotación de polígonos con base en segmentos rectos de referencia

En los dos posts anteriores se consideraron la rotación de rasgos de una capa vectorial (a partir de un punto de rotación), en sentido horario, un angulo determinado y la obtención de la bounding box del vectorial rotado. El ángulo no era arbitrario porque se determinó con base en una recta de referencia con sólo dos puntos. El objetivo de ello era encontrar la longitud proyectada del polígono sobre las líneas de referencia.

En este post se va a señalar el procedimiento para la determinación de los ángulos referidos. El proceso se resume en el código siguiente:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
registry = QgsProject.instance()

my_line = registry.mapLayersByName('my_line')

feats_line = [ feat for feat in my_line[0].getFeatures() ]
mlines = [ feat.geometry().asMultiPolyline() for feat in feats_line ]
slines = [ line[0] for line in mlines ]

points = [ point for line in slines for point in line ]

azimuth = points[0].azimuth(points[1])
angle = 90 - azimuth

print(angle)

El código anterior se probó con las rectas de diferente pendiente que se encuentran en la imagen siguiente:


Los resultados fueron, respectivamente para my_line y my_line2, de 26.081262365320846 (cuando la pendiente es positiva) y -28.419909027565367 (cuando la pendiente es negativa) grados.

Ahora supongamos que se quiere obtener lo que se conoce como la oriented minimum bounding box. El procedimiento se basa en determinar primero la convex hull del polígono para simplificar la geometría y determinar luego los ángulos de rotación con base en todos los segmentos de dicho polígono; tal como se presenta en la imagen siguiente.


Modificando el código anterior para obtener todos los ángulos basados en los vértices se obtiene lo siguiente:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
registry = QgsProject.instance()

vertices_convex_hull = registry.mapLayersByName('vertices_convex_hull')

feats_vertex = [ feat for feat in vertices_convex_hull[0].getFeatures() ]

points = [ QgsPoint(feat.geometry().asPoint()) for feat in feats_vertex ]

for i in range(len(points)-2):
    print(i, i+1)
    azimuth = points[i].azimuth(points[i+1])
    angle = 90 - azimuth
    print (angle)

cuya ejecución produce los siguientes resultados:

0 1
148.4957332807871
1 2
56.92932217723549
2 3
25.780226574085162
3 4
16.676180725954737
4 5
4.713428974539951
5 6
-43.91907581333771
6 7
255.37912601137148
7 8
243.43494882291841
8 9
231.54629078329296

Uno de los ángulos obtenidos anteriormente permite determinar la oriented minimum bounding box. Para ello hay que rotar con base en los segmentos y calcular el área de cada una de las bounding box con los QgsPointXY(min(x), max(y)) y QgsPointXY(max(x), min(y)). El índice del área mínima permitirá acceder a los valores para confeccionar la bounding box definitiva como una memory layer. Esto se señala en la imagen siguiente para line4; el segmento que produce la oriented minimum bounding box para un ángulo de 16.676180725954737 grados (índice 3). Obsérvese que linea4 es paralela a los dos lados mayores del rectángulo de área mínima.


No hay comentarios: