Skip to content

Latest commit

 

History

History
132 lines (83 loc) · 3.99 KB

Find Nearest Line Segments to Point.md

File metadata and controls

132 lines (83 loc) · 3.99 KB

From Find Nearest Line Segments to Point

I have reproduced your example with shapefiles.

You can use Shapely and Fiona to solve your problem.

####1) Your problem (with a shapely Point):

enter image description here

####2) starting with an arbitrary line (with an adequate length):

    from shapely.geometry import Point, LineString
    line = LineString([(point.x,point.y),(final_pt.x,final_pt.y)])

enter image description here

####3) using shapely.affinity.rotate to create the radii (rotating the line from the point, look also the Mike Toews's answer at Python, shapely library: is it possible to do an affine operation on shape polygon?):

    from shapely import affinity
    # Rotate i degrees CCW from origin at point (step 10°)
    radii= [affinity.rotate(line, i, (point.x,point.y)) for i in range(0,360,10)]
    

enter image description here

####4) now, using shapely:cascaded_union (or shapely:unary_union) to get a MultiLineString:

from shapely.ops import cascaded_union
mergedradii = cascaded_union(radii)
print mergedradii.type
MultiLineString

####5) the same with the original lines (shapefile)

    import fiona
    from shapely.geometry import shape
    orlines = fiona.open("orlines.shp")
    shapes = [shape(f['geometry']) for f in orlines]
    mergedlines = cascaded_union(shapes)
    print mergedlines.type
    MultiLineString

####6) the intersection between the two multigeometries is computed and the result is saved to a shapefile:

     points =  mergedlines.intersection(mergedradii)
     print points.type
     MultiPoint
     from shapely.geometry import mapping
     schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
     with fiona.open('intersect.shp','w','ESRI Shapefile', schema) as e:
          e.write({'geometry':mapping(points), 'properties':{'test':1}})

Result:

enter image description here

####7) but, problem, if you use a a longer radius, the result is different:

enter image description here

####8) And if you want to get your result, you need to select only the point with the shortest distance from the original point on a radius:

    points_ok = []
    for line in mergeradii:
       if line.intersects(mergedlines):
           if line.intersection(mergedlines).type == "MultiPoint":
              # multiple points: select the point with the minimum distance
              a = {}
              for pt in line.intersection(merged):
                  a[point.distance(pt)] = pt
              points_ok.append(a[min(a.keys())])
           if line.intersection(mergedlines).type == "Point":
              # ok, only one intersection
              points_ok.append(line.intersection(mergedlines))
    solution = cascaded_union(points_ok)
    schema = {'geometry': 'MultiPoint','properties': {'test': 'int'}}
    with fiona.open('intersect3.shp','w','ESRI Shapefile', schema) as e:
         e.write({'geometry':mapping(solution), 'properties':{'test':1}})

Final result:

enter image description here

I hope that is what you want.