Cómo dividir un ángulo en dos mediante PyQGIS: 1ra parte

En los SIGs, probablemente hemos sentido alguna vez la necesidad de encontrar la forma de cómo dividir un ángulo, de un objeto cualquiera que los presente, en dos partes iguales. Por definición, un ángulo es la región del espacio delimitada por dos semirectas que contienen un vértice común. Por tanto, una forma sería encontrar la semirecta (llamada también bisectriz) que produzca la división en partes iguales.

Para ello, si observamos el procedimiento práctico delineado con regla y compás en:

https://es.wikipedia.org/wiki/Bisectriz

se puede intuir que en nuestro SIG una posible operacionalización equivaldría a obtener un buffer centrado en el vértice cómun, encontrar los puntos de intersección del buffer con las semirectas, trazar la semirecta entre estos dos puntos y luego encontrar la menor distancia entre este segmento y el vértice común.

Supongamos que se tiene el shapefile de línea (multiparte) de la imagen siguiente:

line1

el segundo vértice (índice 1) corresponde al común a las dos semirectas y sobre él vamos a obtener el buffer correspondiente a la mitad de la longitud de la menor de ellas. Para producir el buffer se requiere una capa provisional en memoria, tipo punto, sobre la cual extraer su rasgo (feature). Con el método ‘geometry’ de QgsFeature, convertimos en un objeto de QgsGeometry para finalmente obtener el buffer. Este buffer también lo expresamos como una memory layer pero sólo para el anillo exterior. Es necesario de esta manera para que la intersección con las semirectas sólo produzca puntos de intersección; no otra multilínea.

El código es el siguiente:

from math import sqrt

layer = iface.activeLayer()

feat_line = layer.getFeatures().next()

points = feat_line.geometry().asPolyline()

n = len(points)

lengths = [ sqrt(points[i].sqrDist(points[i+1])) 
            for i in range(n-1) ]

epsg = layer.crs().postgisSrid()

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

mem_layer = QgsVectorLayer(uri,
                           "Points",
                           "memory")

#Prepare mem_layer for editing
   
mem_layer.startEditing()

feat =QgsFeature()

feat.setGeometry(QgsGeometry.fromPoint(points[1]))  #Set geometry
feat.setAttributes([1])
mem_layer.addFeature(feat, True)

#stop editing and save changes
mem_layer.commitChanges()

new_feat = mem_layer.getFeatures().next()

polygon_ring = new_feat.geometry().buffer(0.5*min(lengths), -1).asPolygon()

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

mem_layer = QgsVectorLayer(uri,
                           "buffer_ring",
                           "memory")

mem_layer.startEditing()

feat =QgsFeature()

feat.setGeometry(QgsGeometry.fromPolyline(polygon_ring[0]))  #Set geometry
feat.setAttributes([1])
mem_layer.addFeature(feat, True)

#stop editing and save changes
mem_layer.commitChanges()

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

new_feat = mem_layer.getFeatures().next()

print new_feat.geometry().intersection(feat_line.geometry()).exportToWkt()

que ejecutado en la Python Console de QGIS permite obtener lo siguiente:

line2

Los puntos han sido producidos con el plugin QuickWKT a partir de lo que imprime la última instrucción en la Python Console.

Ahora sólo falta la determinación de la bisectriz; tema del próximo post.

Esta entrada fue publicada en PyQGIS, QGIS, SIG, Software Libre. Guarda el enlace permanente.

Una respuesta a Cómo dividir un ángulo en dos mediante PyQGIS: 1ra parte

  1. Pingback: Cómo dividir un ángulo en dos mediante PyQGIS: 2da parte | El Blog de José Guerrero

Responder

Por favor, inicia sesión con uno de estos métodos para publicar tu comentario:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s