sábado, 12 de agosto de 2017

Dividir vectoriales tipo línea en partes iguales mediante algoritmo basado en azimuths

En el post anterior se consideró la división de vectoriales tipo línea con puntos mediante un algoritmo basado en azimuths. En éste se va tratar la misma división del vectorial pero en partes iguales. Por tanto, los puntos deberían estar igualmente espaciados.

Como se puede preveer, en este caso no se requiere una capa explícita de puntos sino que se pueden generar mediante el método 'interpolate' de shapely.geometry. El código completo, para una división en 15 partes iguales, es el siguiente:

  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
152
153
154
155
import fiona
from shapely.geometry import shape, mapping, Point, LineString
import math
from fiona.crs import from_epsg

llayer = iface.activeLayer()
provider = llayer.dataProvider()

epsg = llayer.crs().postgisSrid()
crs = from_epsg(epsg)

path = provider.dataSourceUri().split('|')[0]

line = fiona.open(path)

geom = [ feat["geometry"] for feat in line ]

points = [ feat["geometry"]["coordinates"] for feat in line ]

line = shape(geom[0])

coords = points[0]

length = line.length

  
number_parts = 15
interval = length/number_parts

points = [ zip(line.interpolate(interval*i).xy[0], line.interpolate(interval*i).xy[1]) for i in range(1, number_parts) ]
            
new_points = [ point[0] for point in points ]

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 = []

for slice in slices:
    line = LineString(slice)
    lines.append(line)

#creation of the resulting shapefile
schema2 = {'geometry': 'LineString','properties': {'id': 'int'},}

with fiona.open('/home/zeito/pyqgis_data/result2.shp', 'w', 'ESRI Shapefile', schema2, crs)  as output:

    for i, line in enumerate(lines):
        output.write({'geometry':mapping(line),'properties': {'id':i}})

result2 = QgsVectorLayer('/home/zeito/pyqgis_data/result2.shp',
                         'splitted_line',
                         'ogr')

QgsMapLayerRegistry.instance().addMapLayer(result2)

Después de ejecutado en la Python Console de QGIS el resultado se presenta en la imagen a continuación:


Aquí se puede observar, con fines de corroboración de la adecuada ejecución del código, la selección de los features 4 y 8 de la línea dividida. Por otra parte, la línea se dividió hasta en 1000 y 3000 partes y los resultados fueron los esperados.





Publicado por José Guerrero en 10:27

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)