lunes, 12 de marzo de 2018

Dividir vectoriales tipo línea en partes de igual longitud mediante algoritmo basado en azimuths con PyQGIS 3

En un post anterior se consideró la división de vectoriales tipo línea en partes iguales mediante un algoritmo basado en azimuths. Sin embargo, se usaron módulos third party como fiona que son difíciles de instalar en Windows y ahora, para que funcione sin problemas en QGIS 3 para este sistema operativo, he de adaptarlo completamente a PyQGIS.

El referido código, para una divisón en 5 partes iguales, 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
from shapely.geometry import LineString
import math

registry = QgsProject.instance()

layer = registry.mapLayersByName("new_line")

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

line_geom = feats[0].geometry()
line_qgis = line_geom.asMultiPolyline()

coords = line_qgis[0]

length = line_geom.length()
number_parts = 5
interval = length/number_parts

# create points every x meters along the line
points_qgis = [ line_geom.interpolate(interval*i).asPoint() for i in range(1, number_parts) ]

new_points = points_qgis

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

new_coords = []

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

values = []

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

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

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

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

complete_points = []

if values[0] == 0:

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

    complete_points.append(new_coords[-1][1])
 
else:
    k = -1
    for i, element in enumerate(new_coords):
        complete_points.append(element[0])
        try:
            if i > 0:
                tmp = new_idxs[k]
                for j, item in enumerate(tmp):
                    if j != 0:
                        complete_points.append(item)
        except IndexError:
            pass
        k += 1
        
    complete_points.append(new_coords[-1][1])

sum = 0
j = 1

count = []

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

slices = []

h = 0

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

lines = []

for slice in slices:
    lines.append(LineString(slice).wkt)

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)

Después de ejecutado en la Python Console se obtiene el vectorial de línea de la imagen a continuación (sobreyaciendo al original). Los features se han escogido de manera alternada en la tabla de atributos para que se pueda corroborar que el código funciona de la manera esperada.

Publicado por José Guerrero en 18:21
Etiquetas: ,

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)