Cómo dividir un ángulo en dos mediante PyQGIS: 2da parte

En el post anterior se consideró la primera parte de cómo dividir un ángulo en dos mediante PyQGIS. Esto abarcó hasta la obtención de los puntos de intersección de las semirectas del ángulo con el anillo exterior del buffer centrado en su vértice común.

Lo que resta ahora es producir el segmento de recta que une ambos puntos y luego utilizar el método ‘closestSegmentWithContext’ de QgsGeometry para encontrar el punto de dicho segmento que se encuentra más cercano al vértice común; que precisamente corresponde a su punto medio. Con este último y el vértice común se obtiene la bisectriz. El código completo se incluye a continuación:

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_points = mem_layer.getFeatures().next()

polygon_ring = new_feat_points.geometry().buffer(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()

new_feat = mem_layer.getFeatures().next()

segment = new_feat.geometry().intersection(feat_line.geometry()).asMultiPoint()

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

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

mem_layer.startEditing()

feat =QgsFeature()

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

#stop editing and save changes
mem_layer.commitChanges()

feat = mem_layer.getFeatures().next()

closest_point= feat.geometry().closestSegmentWithContext(QgsPoint(points[1]))[1]

bisector = [points[1], closest_point]

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

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

mem_layer.startEditing()

feat =QgsFeature()

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

#stop editing and save changes
mem_layer.commitChanges()

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

Observe que a pesar de que se usan diferentes memory layer (todas con el mismo nombre) sólo la última, que es la bisectriz, es la que se añade al registro para su visualización.

El resultado de ejecución del código anterior para el shapefile ya considerado es el siguiente:

line3

es decir, el esperado.

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

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