Converting points to lines with Python GDAL

Question

I have a dictionary with points. The points are derived from a raster (see image below). The points are not sorted and are not in order. The dictionary looks for example like this:

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

I want to create a multiline. The maximum distance from point to point is 14 m. If the points are further apart, they go to a new line.

So far I have the following code. It works, but the problem is, that the points are connected in the wrong order, as you can see in the image below.

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)

How do I do this in Python, ideally with GDAL/OGR (no ArcPy)?

Answer 1

Try working globally:

1. calculate the distances between all points
2. union the resulting lines pointx - pointy with a distance < 14m

I will use Shapely, much easier for resolving these kinds of problems. You must iterate through all pairs of points to calculate the distance once (as distance point1-point2 = distance point2-point1). There are many solutions in Python and I choose the itertools standard module with combinations.

example:

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

With ogr (look at the 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)

With shapely:

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

So, in your case:

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

And if you want to use the end of your script:

multiline = ogr.CreateGeometryFromWkt(line.wkt)

or using Fiona (an easier Python wrapper of the ogr library)

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)

Result:

But, I do not know if this is what you want.

Answer 2

If the points are in order, you should convert the pointDict keys into a list and then sort it, and use that sorted list to go from one point to the next (dictionary keys are not necessarily iterated in a certain order, so to achieve that order, you will have to provide it yourself:

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