0

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.

2

The easiest solution is to use union, cascaded_unionor 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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.