# Finding Nearest Line Segments to Point using shapely?

Background

From a known Point, I require to establish the nearest surrounding "visible perimeter" against a table of MultiLineStrings, as shown on the diagram.

I've searched this site with a number of terms (e.g. minimum edge, minimum perimeter, nearest neighbour, clip, containing polygon, visibility, snap, cut nodes, ray-trace, flood fill, inner boundary, routing, concave hull) but cannot find any previous question which seems to match this scenario.

Diagram

• The green circle is the known Point.
• The black lines are the known MultiLineStrings.
• The grey lines are an indication of a radial sweep from the known Point.
• The red points are the nearest intersection of the radial sweep and the MultiLineStrings.

Parameters

• The Point will never intersect the MultiLineStrings.
• The Point will always be nominally centred within the MultiLineStrings.
• The MultiLineStrings will never fully enclose the Point, therefore the perimeter will be a MultiLineString.
• There will be a table containing approximately 1,000 MultiLineStrings (normally containing a single line of about 100 points).

Methodology Considered

• Undertake a radial sweep by constructing a series of lines from the known Point (at, say, 1 degree increments).
• Establish the nearest intersection point of each radial sweep line with the MultiLineStrings.
• When one of the radial sweep lines doesn't intersect with any of the MultiLineStrings, this would indicate a gap in the perimeter which would be accommodated in the perimeter MultiLineString construction.

Summary

Whilst this technique will find the nearest intersections, it won't necessarily find all the closest perimeter node points, dependent on the resolution of the radial sweep. Can anyone recommend an alternative method to establish all the perimeter points or supplement the radial sweep technique with some form of buffering, sectoring or offsetting?

Software

My preference is to use SpatiaLite and/or Shapely for the solution but would welcome any suggestions which could be implemented using open source software.

Edit: Working Solution (based on answer by @gene)

``````from shapely.geometry import Point, LineString, mapping, shape
from shapely import affinity
import fiona

sweep_res = 10  # sweep resolution (degrees)
focal_pt = Point(0, 0)  # radial sweep centre point

# create the radial sweep lines
line = LineString([(focal_pt.x,focal_pt.y), \

sweep_lines = [affinity.rotate(line, i, (focal_pt.x, focal_pt.y)) \
for i in range(0, 360, sweep_res)]

# load the input lines and combine them into one geometry
input_lines = fiona.open("input_lines.shp")
input_shapes = [shape(f['geometry']) for f in input_lines]

perimeter = []
# traverse each radial sweep line and check for intersection with input lines

if inter.type == "MultiPoint":
# radial line intersects at multiple points
inter_dict = {}
for inter_pt in inter:
inter_dict[focal_pt.distance(inter_pt)] = inter_pt
# save the nearest intersected point to the sweep centre point
perimeter.append(inter_dict[min(inter_dict.keys())])

if inter.type == "Point":
# radial line intersects at one point only
perimeter.append(inter)

if inter.type == "GeometryCollection":
# radial line doesn't intersect, so skip
pass

# combine the nearest perimeter points into one geometry

# save the perimeter geometry
schema = {'geometry': 'MultiPoint', 'properties': {'test': 'int'}}
with fiona.open('perimeter.shp', 'w', 'ESRI Shapefile', schema) as e:
e.write({'geometry':mapping(solution), 'properties':{'test':1}})
``````
• Normally a radial sweep has no meaningful "resolution": it scans from one "event" to the next in order, where events consist of the original nodes of the polylines and of their mutual intersections (which can be found dynamically while sweeping around the original nodes). Its output will be perfectly accurate. Dec 30, 2013 at 17:46

I have reproduced your example with shapefiles.

You can use Shapely and Fiona to solve your problem.

### 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)])
``````

### 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)]
``````

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

``````from shapely.ops import cascaded_union
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]
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:

### 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 = []
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))
• Excellent answer - particularly informative with respect to Fiona usage for input/output via shapefiles. I've appended some code to my question which uses your answer and modified it to reduce the amount of `intersection` calculations required - thank you. Dec 30, 2013 at 21:29