miércoles, 9 de agosto de 2017

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

En el post anterior se consideró un algoritmo para agrupar puntos por segmentos de recta mediante listas Python. En éste lo vamos a complementar con otro basado en azimuths para dividir rectas mediante puntos.

La lectura del archivo vectorial se va a realizar con fiona y las coordenadas y la geometría se conforman con shapely. El código completo se presenta a continuación e incluye los puntos de corte:

  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
156
157
158
import fiona
from shapely.geometry import shape, mapping, Point, LineString
import math
from fiona.crs import from_epsg
crs = from_epsg(32612)

path = '/home/zeito/pyqgis_data/new_line.shp'

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

new_points = [ (400488.72071200114, 4447409.7246439075), (401049.78326624096, 4447112.397201004)]

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
schema1 = {'geometry': 'Point','properties': {'id': 'int'},}
with fiona.open('/home/zeito/pyqgis_data/result.shp', 'w', 'ESRI Shapefile', schema1, crs)  as output:

    for i, point in enumerate(new_points):
        output.write({'geometry':mapping(Point(point[0], point[1])),'properties': {'id':i}})

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

    for i, point in enumerate(complete_points):
        output.write({'geometry':mapping(Point(point[0], point[1])),'properties': {'id':i}})


#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)

La imagen siguiente contiene la línea y los puntos de corte. La tabla atributiva del vectorial tipo línea señala un sólo feature.


Después de la ejecución del código anterior en la Python Console de QGIS se obtiene:


Los cortes se han producido tal como se esperaba por lo que el vectorial de línea resultante presenta tres features. En la imagen anterior se tiene al segundo rasgo como 'selected feature'.





Publicado por José Guerrero en 22:05
Etiquetas: ,

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)