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.
Some solutions according to the position of an element in a list:
polA = Polygon([(0,0), (3,0), (3,3), (0,3)])
polB = Polygon([(2,-1), (5,-1), (5,2), (2,2)])
polC = Polygon([(5,2), (8,2), (8,5), (5,5)])
collection = [polA, polB, polC]
Iterating by index:
for i in range(len(collection)-1):
print collection[i], collection[i+1], collection[i].touches(collection[i+1])
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
Iterating with itertools (standard module), all combinations without repeated elements
import itertools
for pol in itertools.combinations(collection, 2):
print pol[0],pol[1], pol[0].touches(pol[1])
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
or with the suggestion of Mike T
for pol1, pol2 in itertools.combinations(collection, 2):
print pol1, pol2, pol1.touches(pol2)
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
Compare with
for pol in itertools.permutations(collection, 2):
print pol[0],pol[1], pol[0].touches(pol[1])
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) False
POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) True
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) False
POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) True
Look also at Python - Previous and next values inside a loop
from itertools import tee, islice, chain, izip
def previous_and_next(some_iterable):
prevs, items, nexts = tee(some_iterable, 3)
prevs = chain([None], prevs)
nexts = chain(islice(nexts, 1, None), [None])
return izip(prevs, items, nexts)
for previous, item, nxt in previous_and_next(collection):
print "Item is now", item, "next is", nxt, "previous is", previous
Item is now POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0)) next is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) previous is None
Item is now POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1)) next is POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) previous is POLYGON ((0 0, 3 0, 3 3, 0 3, 0 0))
Item is now POLYGON ((5 2, 8 2, 8 5, 5 5, 5 2)) next is None previous is POLYGON ((2 -1, 5 -1, 5 2, 2 2, 2 -1))
Best Answer
Make concentric circles
Instead of making 360 different line-strings, why not make <360 different concentric circles around the point of interest, until you find the one that intersects land within your desired accuracy.
For example, start with a concentric circle around the point at a distance of $n$. If you have a shapefile in unit distance coordinates, like UTM (as opposed to lon/lat coords), you can take a point
(x, y)
and do:First, determine what you want to be the maximum range at which the function will work. Given the limits of UTM projections, this method won't be too accurate at over 1000km. So, start at 1000 km. If the actual distance is 356 km; then a search to 1 km accuracy would look like 1000->500->250->375->313->344->360->352->356->354->355. There will always be 11 circles (log2 1000 ~= 10 + 1 initial point).
This method will likely generate and call the intersect function on fewer shapes, though the shapes will be more complex. I'm pretty sure this will save you computationally.
Haversine everything
What is the accuracy of your shapes representing land? How many points make up those shapes? What if you just measure the distance of every shape to the point of origin using a haversine formula?
There is a lot of computational complexity building shapes and calculating intersects. haversine is just brute force math, and that tends to run really fast in numpy, in my experience. If your shapes are not that accurate (i.e. shapes with 1km accuracy instead of 20m, or something) then his approach may actually be faster.
While the speed is debatable, it is certainly