domingo, 13 de agosto de 2017

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

Con el advenimiento de QGIS 3 el próximo 29 de septiembre de 2017, si se cumple exactamente el cronograma de liberación, habrán varios cambios que afectarán la compatibilidad de todas las aplicaciones que hasta ahora se han realizado en QGIS 2.x. La razón principal es que QGIS dejará de usar Python 2.7.x en favor de Python 3.x y el framework Qt4/PyQt4 para dar paso a Qt5/PyQt5. Por tanto, conviene estar preparado.

Como en el post anterior se empleó un algoritmo basado en azimuths para dividir vectoriales tipo línea en partes iguales, el cual utilizaba una mezcla de librerías third party de python (fiona y shapely) conjuntamente con PyQGIS, se pensó adaptar dicho script para ser funcional con Python 3.

Como es obvio pensar, tenemos que excluir toda referencia a PyQGIS (todavía basada en Python 2.7.x) y centrarnos sólo en los módulos Python empleados. Como en Debian ya contaba con Python 3, me tocó instalar python3-pip para luego disponer de fiona y shapely, en Python 3, con los siguientes comandos (como superusuario):

1
2
pip3 install fiona
pip3 install shapely

Excluidos los objetos PyQGIS en el código, el primer intento de ejecutarlo arrojó sólo una incompatibilidad:

1
2
3
4
5
6
Traceback (most recent call last):
  File "test2.py", line 29, in <module>
    new_points = [ point[0] for point in points ]
  File "test2.py", line 29, in <listcomp>
    new_points = [ point[0] for point in points ]
TypeError: 'zip' object is not subscriptable

La razón es que en Python 3 zip devuelve un objeto iterable; no una lista como en Python 2. Esto se solucionó con lo señalado aquí.

El código completo 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
import fiona
from shapely.geometry import shape, mapping, Point, LineString
import math
from fiona.crs import from_epsg

epsg = 32612
crs = from_epsg(epsg)

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

  
number_parts = 15
interval = length/number_parts

points = [ list(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/splitted_line.shp', 'w', 'ESRI Shapefile', schema2, crs)  as output:

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

y se ejecutó, sin errores, en consola de bash con:

1
python3 test2.py

Al cargar en QGIS la capa resultante y seleccionar varios rasgos se observó que el script funcionó de la manera esperada; tal como se presenta en la imagen siguiente:


Publicado por José Guerrero en 0:26
Etiquetas: , ,

No hay comentarios:

Publicar un comentario

Suscribirse a: Enviar comentarios (Atom)