Dividir vectoriales tipo línea con puntos mediante algoritmo basado en azimuths con inserción directa del punto de corte en la posición correspondiente

En un post anterior desarrollé un procedimiento para dividir vectoriales tipo línea con puntos mediante un algoritmo basado en azimuths e índices de todos los puntos sin inserción directa de los puntos de corte en la lista final. Aunque funciona de la manera esperada tratando de adaptarlo al lenguaje R me di cuenta que el algoritmo tiene un bucle innecesario y es demasiado complejo para una adaptación expedita.

Durante la simplificación del algoritmo descubrí la manera de insertar los puntos de corte en la lista final en el orden correcto lo que permite un 'slicing' rápido y directo sin necesidad de un uso excesivo de índices.

El código siguiente contiene al inicio las coordenadas de los puntos de la línea y los puntos de corte para probar el algoritmo. Los puntos de corte dividen la línea en 5 partes iguales.

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

from math import fabs points = [QgsPoint(399932.715115, 4447160.8106), QgsPoint(400128.727295, 4447371.90064), QgsPoint(400561.677277, 4447417.13422), QgsPoint(401210.02526, 4447010.032), QgsPoint(401625.743401, 4447548.527), QgsPoint(402291.323224, 4447690.68968), QgsPoint(402409.792125, 4448069.79017)] cut_points = [QgsPoint(400488.108679, 4447409.44795), QgsPoint(401049.003665, 4447111.13858), QgsPoint(401490.677202, 4447373.57079), QgsPoint(402044.667428, 4447638.00592)] tot_points = len(points) + len(cut_points) k = 0 new_points = [] for i in range(len(points) - 1): try: coor_az = points[i].azimuth(points[i+1]) point_az = points[i].azimuth(cut_points[k]) diff_az = fabs(coor_az - point_az) if diff_az > 1e-6: new_points.append(points[i]) new_points.append(points[i+1]) else: new_points.append(cut_points[k]) new_points.append(points[i+1]) k += 1 except IndexError: pass tot_new_points = len(new_points) diff_points = tot_points - tot_new_points n = len(points) add_points = points[(n-diff_points):(n)] for point in add_points: new_points.append(point) k = 0 idx = 0 idxs = [0] for point in new_points: if point == cut_points[k]: idxs.append(idx) if k != len(cut_points) -1: k += 1 else: break idx += 1 idxs.append(len(new_points) - 1) slices = [] for i in range(len(idxs) - 1): slice = new_points[idxs[i]: idxs[i+1]+1] slices.append(slice) epsg = 32612 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(slices)) ] for i, feat in enumerate(feats): feat.setAttributes([i]) feat.setGeometry(QgsGeometry.fromPolyline(slices[i])) prov.addFeatures(feats) QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

La ejecución del código anterior permite dividir, tal como se esperaba, la línea en 5 partes iguales. En la imagen siguiente se observa el resultado obtenido donde se ha seleccionado el cuarto feature objeto de corte.

Si se modifica el código anterior para tener sólo el primer punto de corte, la ejecución del mismo también permite obtener el resultado esperado; tal como se observa a continuación (donde se seleccionó el segundo feature objeto de corte):