domingo, 11 de marzo de 2018

Dividir vectoriales tipo línea con puntos mediante algoritmo basado en azimuths con PyQGIS 3

En un post anterior se implementó un código para dividir vectoriales tipo línea con puntos mediante un algoritmo basado en azimuths. El referido código fue creado en Linux con fiona. Este módulo es difícil de instalar en Windows 7, al igual que GeoPandas, debido a la necesidad de compiladores C++ específicos para la configuración de nuestro sistema (ahora difíciles de encontrar porque ya prácticamente no tienen soporte sino se usa Windows 10). Como estoy probando QGIS 3 en este Windows voy a utilizar mucho más código PyQGIS prescindiendo, en lo posible, de los módulos third party.

El código lo produje primero para la versión anterior de QGIS y luego lo adapté a QGIS 3 porque era más sencillo. El número de clases y métodos nuevos no era tan extenso y procuré considerar las geometrías en formato WKT (shapely ayuda mucho en ese sentido). Se presenta a continuación:

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
from shapely.geometry import LineString
import math

registry = QgsProject.instance()

layer = registry.mapLayersByName('new_line')

feats = [ feat for feat in layer[0].getFeatures() ]

feature = feats[0]

line_geom = feature.geometry()
length = line_geom.length()

line_qgis = line_geom.asMultiPolyline()

coords = line_qgis[0]

new_points = [ (400488.826418, 4447409.52296), (401454.700428, 4447326.96875) ]

azimuths = []

for i in range(len(coords)-1):
    point1 = coords[i]
    point2 = coords[i+1]
    azimuths.append(math.atan2(point2[0] - point1[0], point2[1] - point1[1]))

idxs = []

k = 0

for i in range(len(coords)-1):
    for point in new_points:
        new_az = math.atan2(point[0] - coords[i][0], point[1] - coords[i][1])
        for j, item in enumerate(azimuths):
            if math.fabs(item - new_az) < 1e-6:
                idxs.append([j, new_points[k]])
                k +=1

values = []

for item in idxs:
    if item[0] not in values:
        values.append(item[0])

list = [ 1 for i in range(len(values)) ]

for i, item in enumerate(idxs):
    try:
        if idxs[i][0]== idxs[i+1][0] and values[0] == 0:
            list[idxs[i][0]] += 1
        if idxs[i][0]== idxs[i+1][0] and values[0] != 0:
            list[idxs[i][0]-1] += 1

    except IndexError:
        pass

new_idxs = [ [] for i in range(len(values)) ] 

k = 0

for i, item in enumerate(list):
    new_idxs[i].append(idxs[k][0])
    
    for j in range(item):
        new_idxs[i].append(idxs[k][1])
        k +=1

new_coords = []

for i in range(len(coords)-1):
    new_coords.append([coords[i],coords[i+1]])

complete_points = []

for i, element in enumerate(new_coords):
    complete_points.append(new_coords[i][0])
    for item in new_idxs:
        if item[0] == i:
            tmp = item
            for j, element in enumerate(tmp):
                if j != 0:
                    complete_points.append(element)

complete_points.append(new_coords[-1][1])

sum = 0

j = 1

count = []

sum_interv = 0

intervals = []

for point in new_points:
    if point in complete_points:
        l = LineString(complete_points[:complete_points.index(point)+1]).length
        interval = l - sum_interv
        sum_interv += interval
        intervals.append(interval)

for i in range(len(complete_points)-1):
    sum += LineString([complete_points[i], complete_points[i+1]]).length
    j += 1
    for interval in intervals:
        if math.fabs(sum - interval) < 1e-6:
            count.append(j)
            sum = 0
            j = 1

slices = []

h = 0

for i in range(len(new_points)):
    slice = complete_points[h: h + count[i]]
    slices.append(slice)
    h +=  count[i] - 1

slices.append(complete_points[h:len(complete_points)])

lines = []
lines_qgis = [ [] for i in range(len(slices)) ]

for i, slice in enumerate(slices):
    for point in slice:
        lines_qgis[i].append(QgsPoint(point[0], point[1]))

    lines.append(QgsGeometry.fromPolyline(lines_qgis[i]).asWkt())

epsg = layer[0].crs().postgisSrid()

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

mem_layer = QgsVectorLayer(uri,
                           'line',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(len(lines)) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromWkt(lines[i]))

prov.addFeatures(feats)

QgsProject.instance().addMapLayer(mem_layer)

Se va a ejecutar con el shapefile que se visualiza en la imagen a continuación:


Después de ejecutado el código en la Python Console de QGIS 3, se puede observar en la imagen siguiente que el vectorial de línea ha sido dividido, tal como se esparaba, en tres trozos por los puntos especificados dentro del código. En la imagen se selecciona el segundo feature de la línea dividida como evidencia de la ejecución adecuada del código.


Publicado por José Guerrero en 8:56
Etiquetas: ,

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)