Convert points to lines with Python GDAL

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 got 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:
stop = True

if not stop:
lineDict[i] = ogr.Geometry(ogr.wkbLineString)
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)
``````

Any ideas are appreciated. I want to do this in Python ideally with GDAL/OGR (no ArcPy).

• Your script has errors `y = pointDict[item][2]` -> `y = pointDict[item][1]` and `yPrevious = pointDict[item-1][3] -> yPrevious = pointDict[item-1][1]`or it doesn't work. – gene Nov 22 '13 at 21:59
• @gene: Fixed the mixed up numbers. – ustroetz Nov 24 '13 at 23:01
• The final code `raster2dictionary2lines` can be found here. I solved it with the answer posted by gene. – ustroetz Nov 25 '13 at 17:43

Why not work 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)
distance =  point1.Distance(point2)
line = ogr.Geometry(ogr.wkbLineString)
....
``````

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.

• Excellent! Thank you very much. I solved it with your idea to iterate over all the points and calculate all the distancea. I didn't check out your shapely code, since I want to stay for now with OGR. Two side nodes to the OGR part: Distance in `point1.Distance(point2)` needs to be uppercase. The link to the GDAL/OGR Cookbook doesn't work. – ustroetz Nov 24 '13 at 22:57
• Any ideas on how to modify the code, so that every point has a maximum of two lines (e.g. to make sure that triangles like in your result image are not created). – ustroetz Nov 24 '13 at 22:59
• link to the GDAL/OGR Cookbook corrected – gene Nov 25 '13 at 7:25
• Or, if you like on one heavy line: `cascaded_union([LineString([a, b]) for a, b in itertools.combinations(pointDict.values(), 2) if Point(a).distance(Point(b)) < 14]).wkt` (requires `from shapely.ops import cascaded_union`) – Mike T Nov 25 '13 at 11:47
• Will `cascaded_union` fix the issue that one point has more than two lines connected or is that just another way to create the same code @gene wrote? – ustroetz Nov 25 '13 at 17:38

Are the points in order? If so, 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
``````
• No, they are not in order. – ustroetz Nov 22 '13 at 23:14
• @ustroetz But the problem you stated in your question is that your result is in the wrong order. You will need an order to your points if you want an order to your line verticies. Are you intending to order them by distance? – Cyrus Nov 23 '13 at 15:11
• To order/sort them by distance is the only way I could think of on how to create a line. But I am also open for other ideas. All I want in the end is a line. Don't care how it is done. – ustroetz Nov 24 '13 at 23:00