I have reproduced your example with shapefiles.
You can use Shapely and Fiona to solve your problem.
1) Your problem (with a shapely Point):
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)])
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)]
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:
7) but, problem, if you use a a longer radius, the result is different:
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:
I hope that is what you want.
I don't understand your question:
line = LineString([(0, 0), (2, 2)])
# create a point which lies along the line
point = line.interpolate(1)
line.contains(point)
True
You want the two lines which lie either side of the point ?
line1 = LineString([line.coords[0],(point.x, point.y)])
line2 = LineString([(point.x, point.y), line.coords[1]])
Upgrade 1: with a line with multiple vertices
You need to iterate through the segments of the LineString to find the one that contains the point
The LineString must be iterate as pair to divide the line in segments.
def pairs(lst):
for i in range(1, len(lst)):
yield lst[i-1], lst[i]
line = LineString([(0,0),(1,2), (2, 2), (2,3), (4,2),(5,5)])
for pair in pairs(list(line.coords)):
if LineString([pair[0],pair[1]]).contains(point):
print LineString([pair[0],pair[1]])
LINESTRING (2.00 2.00, 2.00 3.00)
And you can use the previous answer: a rapid solution, for example ( this can be done better):
line1 = []
line2 = []
cp = False
for pair in pairs(list(line.coords)):
if cp == False:
line1.append(pair[0])
if cp == True:
line2.append(pair[1])
if LineString([pair[0],pair[1]]).contains(point):
line1.append((point.x,point.y))
line2.append((point.x,point.y))
line2.append(pair[1])
cp = True
line1 = LineString(line1)
line2 = LineString(line2)
Result:
Best Answer
You need to iterate through the edges of the LinearRing of the polygon
Now using Determine if Shapely point is within a LineString/MultiLineString (using the answer of Mike T using the distance with an appropriate threshold because there are floating point precision errors when finding a point on a line)
Control