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:
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)
Any ideas are appreciated. I want to do this in Python ideally with GDAL/OGR (no ArcPy).
Best Answer
Why not work globally ?
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:
With ogr (look at the Python GDAL/OGR Cookbook!):
With shapely:
So, in your case:
And if you want to use the end of your script:
or using Fiona (an easier Python wrapper of the ogr library)
Result:
But, I do not know if this is what you want.