Convertir puntos en líneas con Python GDAL

Question

Tengo un diccionario con puntos. Los puntos proceden de una trama (véase la imagen de abajo). Los puntos no están ordenados. El diccionario tiene, por ejemplo, este aspecto:

pointDict = {0: (345645.1276541934, 1267223.104499615), 1: (345626.87681620114, 1267223.2540966477), 2: (345645.2772512261, 1267268.581997563), 3: (345617.751397205, 1267223.4036936804), 4: (345654.1034761568, 1267259.306981534), 5: (345636.15183223, 1267231.781127513), 6: (345636.30142926273, 1267268.2828034975), 7: (345626.87681620114, 1267259.306981534), 8: (345617.90099423775, 1267259.306981534), 9: (345608.7755752416, 1267259.6061755994), 10: (345599.7997532782, 1267250.1815625378), 11: (345590.6743342821, 1267250.4807566034)}

Quiero crear una multilínea. La distancia máxima de punto a punto es de 14 m. Si los puntos están más separados, pasan a una nueva línea.

Hasta ahora tengo el siguiente código. Funciona, pero el problema es que los puntos están conectados en el orden equivocado, como se puede ver en la imagen de abajo.

import ogr, gdal, os
from math import sqrt

pointDict = {0: (345645.1276541934, 1267223.104499615), 1: (345626.87681620114, 1267223.2540966477), 2: (345645.2772512261, 1267268.581997563), 3: (345617.751397205, 1267223.4036936804), 4: (345654.1034761568, 1267259.306981534), 5: (345636.15183223, 1267231.781127513), 6: (345636.30142926273, 1267268.2828034975), 7: (345626.87681620114, 1267259.306981534), 8: (345617.90099423775, 1267259.306981534), 9: (345608.7755752416, 1267259.6061755994), 10: (345599.7997532782, 1267250.1815625378), 11: (345590.6743342821, 1267250.4807566034)}

multiline = ogr.Geometry(ogr.wkbMultiLineString)

i = 0
lineDict = {}
for item in pointDict:
    stop = False

    x = pointDict[item][0]
    y = pointDict[item][1]

    if item != 0:
        xPrevious = pointDict[item-1][0]
        yPrevious = pointDict[item-1][1]
        distance = sqrt((y-yPrevious)**2+(x-xPrevious)**2)

    for line in multiline:
        if line.GetPointCount() > 0:
            j = 0 
            for j in range(0, line.GetPointCount()):
                point = line.GetPoint(j)
                xExisting = point[0]
                yExisting  = point[1]
                distance = sqrt((y-yExisting)**2+(x-xExisting)**2)
                j += 1
                if distance < 14:
                    line.AddPoint(x,y)
                    stop = True

    if not stop:
        lineDict[i] = ogr.Geometry(ogr.wkbLineString) 
        lineDict[i].AddPoint(x,y)  
        multiline.AddGeometry(lineDict[i])
        i += 1

for line in multiline:
    print line

outSHPfn = 'test1.shp'   
shpDriver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outSHPfn):
    shpDriver.DeleteDataSource(outSHPfn)
outDataSource = shpDriver.CreateDataSource(outSHPfn)
outLayer = outDataSource.CreateLayer(outSHPfn, geom_type=ogr.wkbMultiLineString )
featureDefn = outLayer.GetLayerDefn()
outFeature = ogr.Feature(featureDefn)
outFeature.SetGeometry(multiline)
outLayer.CreateFeature(outFeature)

Cualquier idea es bienvenida. Quiero hacer esto en Python idealmente con GDAL/OGR (no ArcPy).

Answer 1

¿Por qué no trabajar globalmente?

1. calcular las distancias entre todos los puntos
2. unión de las líneas resultantes puntox - punta con una distancia < 14m

Utilizaré Shapely mucho más fácil para resolver este tipo de problemas. Debes iterar a través de todos los pares de puntos para calcular la distancia una vez (como distancia punto1-punto2 = distancia punto2-punto1). Hay muchas soluciones en Python y yo elijo el módulo estándar itertools con combinaciones .

ejemplo:

myPointDict = {0:(1,1), 1:(2,2), 2:(3,3),3:(4,4),4:(5,5)}
import itertools
for i in  itertools.combinations(PointDict.values(), 2):
   print i
((1, 1), (2, 2))
((1, 1), (3, 3))
((1, 1), (4, 4))
((1, 1), (5, 5))
((2, 2), (3, 3))
((2, 2), (4, 4))
((2, 2), (5, 5))
((3, 3), (4, 4))
((3, 3), (5, 5))
((4, 4), (5, 5))

Con ogr (mire ¡el libro de cocina Python GDAL/OGR Cookbook! ):

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(x,y)
distance =  point1.Distance(point2)
line = ogr.Geometry(ogr.wkbLineString) 
    line.AddPoint(x1, y1)
    ....
    line.AddPoint(xn,yn)

Con forma:

point = Point(x,y)
distance = Point(x1,y1).distance(Point(x2,y2)
linestring = LineString([point1,..., pointn] 

Así que, en tu caso:

from shapely.geometry import Point, LineString
# creation of a empty line for unioning the resulting geometries
line = LineString()
for i in  itertools.combinations(pointDict.values(), 2):
     # if distance < 14m union the line ptx-pty to line
     if Point(i[0]).distance(Point(i[1])) < 14:
            line = line.union(LineString([(Point(i[0]).x, Point(i[0]).y), (Point(i[1]).x, Point(i[1]).y)]))
     # result 
     print line.wkt
     'MULTILINESTRING ((345672.493225679441821 1267286.555012494325638,345681.57590266619809 1267286.555012494325638),(345672.493225679441821 1267286.555012494325638,345663.410548692685552 1267277.472335507394746),(345672.493225679441821 1267286.555012494325638,345681.57590266619809 1267277.472335507394746),(345681.57590266619809 1267286.555012494325638,345681.57590266619809 1267277.472335507394746),(345654.327871705929283 1267277.472335507394746,345663.410548692685552 1267277.472335507394746),(345654.327871705929283 1267277.472335507394746,345645.245194719173014 1267268.389658520696685),(345681.57590266619809 1267277.472335507394746,345690.658579652954359 1267268.389658520696685),(345636.162517732358538 1267268.389658520696685,345645.245194719173014 1267268.389658520696685),(345636.162517732358538 1267268.389658520696685,345627.079840745602269 1267259.306981533998623),(345690.658579652954359 1267268.389658520696685,345681.57590266619809 1267259.306981533998623),(345608.914486772089731 1267259.306981533998623,345617.997163758846 1267259.306981533998623),(345608.914486772089731 1267259.306981533998623,345599.831809785333462 1267250.224304547300562),(345617.997163758846 1267259.306981533998623,345627.079840745602269 1267259.306981533998623),(345681.57590266619809 1267259.306981533998623,345672.493225679441821 1267250.224304547300562),(345590.749132798577193 1267250.224304547300562,345599.831809785333462 1267250.224304547300562),(345672.493225679441821 1267250.224304547300562,345663.410548692685552 1267241.14162756036967))'

Y si quieres usar el final de tu guión:

multiline = ogr.CreateGeometryFromWkt(line.wkt)

o utilizando Fiona (una envoltura Python más sencilla de la biblioteca ogr)

import fiona
from shapely.geometry import mapping
# schema of the shapefile
schema = {'geometry': 'MultiLineString','properties': {'test': 'int'}}
with fiona.open('myshp3.shp','w','ESRI Shapefile', schema) as c:
       record = {'geometry':mapping(line), 'properties':{'test':1}}
       c.write(record)

Resultado:

Pero no sé si esto es lo que quieres.

Answer 2

¿Están los puntos en orden? Si es así, deberías convertir las claves del pointDict en una lista y luego ordenarla, y usar esa lista ordenada para ir de un punto al siguiente (las claves del diccionario no se iteran necesariamente en un orden determinado, así que para conseguir ese orden, tendrás que proporcionarlo tú mismo:

listofkeys - pointDict.keys()
listofkeys.sort()
for item in listofkeys:
     #code