Número de inflexiones (bends) de una polilínea mediante PyQGIS

Leyendo esta pregunta en gis.stackexchange.com, en la que se requería obtener el número de inflexiones (bends) de una polilínea para ArcGis, me pregunté que tan difícil sería hacerlo en QGIS a sabiendas que ya tenía parte del código adelantado producto de un post anterior.

En efecto, en mi post se calculaban las longitudes de los segmentos de recta para todos los rasgos de una polilínea. Para determinar estas distancias se requieren las coordenadas de los extremos del segmento; que son las mismas que se necesitan en el método ‘azimuth’ de la clase QgsPoint. En cartografía, un azimuth es el ángulo de una dirección contado en el sentido de las agujas del reloj a partir del norte geográfico. Por tanto, si la diferencia entre los azimuths de dos segmentos consecutivos es cero puede decirse que existen tres puntos colineales en un mismo segmento.

Es sólo cuestión de agregar unas pocas líneas más de código para realizar esa tarea; tal como se tiene a continuación:

from math import sqrt, fabs
import itertools

layer = iface.activeLayer()

features = layer.getFeatures()

lines = [feature.geometry().asPolyline() for feature in features]

k = 0

for points in lines:
    
    n = len(points)
    list = range(n)
    print "line" + str(k) + ", " + str(n) + " points" 

    length_segments = [sqrt(points[i].sqrDist(points[j])) 
                       for i,j in itertools.combinations(list, 2) 
                       if (j - i) == 1]
    
    sum = 0
    
    for length in length_segments:
        i = length_segments.index(length)
        print "segment = %d, length = %.2f" % (i, length)
        sum += length
    
    print "sum = ", sum
    
    k += 1
    
    coor_segments = [[points[i], points[j]]
                     for i,j in itertools.combinations(list, 2)
                     if (j - i) == 1]
    
    points_segments_flat = [item for element in coor_segments for item in element]
    
    n = len(points_segments_flat)
    
    azimuth_list = [ points_segments_flat[i].azimuth(points_segments_flat[i+1]) for i in range (0,n,2) ]
    
    n = len(azimuth_list)
    
    dif_azimuth = [ (azimuth_list[i+1] - azimuth_list[i])  for i in range (n-1) ]
    
    bends = [item for item in dif_azimuth if fabs(item) > 1e-5 ]
    
    print "bends = ", len(bends)

Para probar el código anterior utilicé el siguiente vectorial que tiene 2, 5 y 3 bends, respectivamente, en cada rasgo. Observen que el vectorial está en modo de edición para poder visualizar los nodos que son colineales dentro de cada segmento.

bends1

Cuando se ejecuta el código en la Python Console se obtiene un resultado satisfactorio:

line0, 5 points
segment = 0, length = 21735.11
segment = 1, length = 20208.88
segment = 2, length = 21260.06
segment = 3, length = 16898.42
sum =  80102.4682351
bends =  2
line1, 14 points
segment = 0, length = 18831.03
segment = 1, length = 14788.45
segment = 2, length = 14918.11
segment = 3, length = 20774.56
segment = 4, length = 13850.87
segment = 5, length = 16268.74
segment = 6, length = 10746.22
segment = 7, length = 15430.28
segment = 8, length = 13019.25
segment = 9, length = 16647.36
segment = 10, length = 17445.18
segment = 11, length = 12683.59
segment = 12, length = 9013.99
sum =  194417.62525
bends =  5
line2, 8 points
segment = 0, length = 11789.15
segment = 1, length = 11477.16
segment = 2, length = 10469.81
segment = 3, length = 10716.35
segment = 4, length = 23595.18
segment = 5, length = 16423.37
segment = 6, length = 9394.02
sum =  93865.0436374
bends =  3

donde se contabilizan de la manera esperada el número de puntos, los bends y las longitudes de cada segmento.

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

Una respuesta a Número de inflexiones (bends) de una polilínea mediante PyQGIS

  1. Pingback: Newsletter Octubre 2015 | 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