Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It only takes a minute to sign up.
Sign up to join this community
I am working with Python, Shapely and Fiona... Is there any way where I can give in an input of a line shapefile and a polygon shapefile, and get a result of the intersection points and their coordinates? The illustration provides an accurate description....
The intersection of a Polygon and a LineString is a LineString and the intersection of two LineStrings is a Point (or MultiPoint), so you need to transform your Polygon into a LineString -> Shapely: LinearRings
from shapely.geometry import shape
import fiona
# polygon layer
poly = fiona.open("polygons.shp")
# line layer
line = fiona.open("lines.shp")
# First Feature of the shapefiles
s_poly = shape(poly.next()['geometry'])
s_line = shape(line.next()['geometry'])
print s_poly.intersection(s_line)
LINESTRING (360.4742985178883 -286.9847328244275, 450.1982781776156 -140.6494330268984)
# transform the polygon into a LineString
ring = LineString(list(s_poly.exterior.coords))
print ring.intersection(line)
MULTIPOINT (360.4742985178883 -286.9847328244275, 450.1982781776156
# or, more formal
from shapely.geometry.polygon import LinearRing
lring = LinearRing(list(s_poly.exterior.coords))
print lring.intersection(s_line)
MULTIPOINT (360.4742985178883 -286.9847328244275, 450.1982781776156 -140.6494330268984)
If you have many polygons and many polylines:
Multi_pol_ext = MultiLineString([list(shape(pol['geometry']).exterior.coords) for pol in fiona.open("polygons.shp")])
Multi_lines = MultiLineString([shape(line['geometry']) for line in fiona.open("lines.shp")])
Multi_pol_ext.intersection(Multi_lines)
<shapely.geometry.multipoint.MultiPoint object at 0x1091a5210>
