I am working on a project that involves finding the intersection of linears in a shapefile, if that vertex(the found intersection) doesn't already exist for each(or one) of the features, then that vertex(the found intersection point) should be added to both features. I'm using the shapely and OGR modules for python. Can somebody just give me a general example of adding a vertex to an existing feature in shapefile using either or and saving it. Thanks.
1 Answer
The easiest solution is to use union
, cascaded_union
or unary_union
. All the lines are split at the points of intersection:
from shapely.geometry import LineString
line1 = LineString([(0, 0), (2, 2),(3,1)])
line2 = LineString([(2, 0), (2, 1),(1,2)])
print line1.intersection(line2)
POINT (1.5 1.5)
for line in line1.union(line2):
print line
LINESTRING (0 0, 1.5 1.5)
LINESTRING (1.5 1.5, 2 2, 3 1)
LINESTRING (2 0, 2 1, 1.5 1.5)
LINESTRING (1.5 1.5, 1 2)
# or
from shapely.ops import unary_union
for line in unary_union([line1,line2]):
print line
LINESTRING (0 0, 1.5 1.5)
LINESTRING (1.5 1.5, 2 2, 3 1)
LINESTRING (2 0, 2 1, 1.5 1.5)
LINESTRING (1.5 1.5, 1 2)
The resulting segments 1 and 2 are part of line1 and 3 and 4 of line2: the new lines are obtained with linemerge(lines)
:
from shapely.ops import linemerge
new_line1 = linemerge([line1.union(line2)[0],line1.union(line2)[1]])
new_line2 = linemerge([line1.union(line2)[2],line1.union(line2)[3]])
print new_line1
LINESTRING (0 0, 1.5 1.5, 2 2, 3 1)
print new_line2
LINESTRING (2 0, 2 1, 1.5 1.5, 1 2)
If you use Shapely, it is easiest to use Fiona to read shapefile instead of ogr and if you have many intersections, it is interesting to use a spatial index as rtree or PyQuadTree, for example.