In your solution you mix geometries and predicates of ogr and Shapely
With Shapely only (intersection(), intersects())
from shapely.geometry import Point, LineString, mapping
geom1 = LineString([(0, 0), (0, 1)])
geom2 = LineString([(1, 1), (-1,-1)])
intersection = geom1.intersection(geom2)
# geo_interface -> GeoJSON
mapping(intersection)
{'type': 'Point', 'coordinates': (0.0, 0.0)}
# Shapely geometry
pt1 = [Point(pt) for pt in list(geom1.coords)]
point_start_1 = pt1[0]
.....
# Shapely predicate
intersection.intersects(point_start_1)
True
With ogr only (Intersection(), Intersects()):
from osgeo import ogr
geom1 = ogr.CreateGeometryFromWkt('LINESTRING (0 0, 0 1)')
geom2 = ogr.CreateGeometryFromWkt('LINESTRING (1 1, -1 -1)')
# intersection with ogr
intersection_ogr = geom1.Intersection(geom2)
# geo_interface -> GeoJSON
intersection_ogr.ExportToJson()
'{ "type": "Point", "coordinates": [ 0.0, 0.0 ] }'
# ogr geometry
point_start_1 = ogr.Geometry(ogr.wkbPoint)
point_start_1.AddPoint(geom1.GetPoint_2D(0)[0],geom1.GetPoint_2D(0)[1])
....
# ogr predicate
intersection_ogr.Intersects(point_start_1)
True
In your solution, you use Shapely geometries:
point_start_1 = Point(geom_ogr1.GetPoint_2D(0))
....
with a binary predicate of ogr -> error
intersection_ogr.Intersects(point_start_1)
TypeError: in method 'Geometry_Equal', argument 2 of type 'OGRGeometryShadow *'
What is strange if the intersection result is a point ?
The intersect predicate is
Returns True if the boundary and interior of the object intersect in any way with those of the other.
With a common point between the geometries, the intersects predicate returns TRUE because the boundary of the first geometry intersects the boundary of the second geometry at this point. The result is always this point
from shapely.geometry import LineString, Polygon
from shapely.wkt import loads
Line1 = loads("LINESTRING ({0} {0}, {1} {0})".format(0, 1))
Line2 = loads("LINESTRING ({1} {0}, {2} {0})".format(0, 1, 2))
Line1.intersects(Line2)
True
Line1.touches(Line2)
True
Line1.intersection(Line2).wkt
'POINT (1 0)'
# polygons
polygon1 = loads('POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))')
polygon2 = loads('POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))')
polygon1.intersection(polygon2).wkt
'POINT (220 360)'
In your case the resulting layer is an ogr.wkbPoint
If you change the value of the point, the result is
polygon3 = loads('POLYGON (( 220 260, 220 370, 300 360, 300 260, 220 260 ))')
polygon1.intersection(polygon3).wkt
'LINESTRING (220 370, 220 360)'
Or
polygon4 = loads('POLYGON ((220 260, 200 370, 300 360, 300 260, 220 260))')
polygon1.intersection(polygon4).wkt
'POLYGON ((220 368, 220 360, 201.8181818181818 360, 200 370, 220 368))'
UPDATE 2: With Shapely, there is no problem to clip a roads Shapefile using a polygon even if the road vertex points and the edge of the polygon coincide
roads = loads('LineString (220 360, 120 468)')
polygon1.intersection(roads).wkt
'LINESTRING (220 360, 140 446.4)'
Best Answer
This has been fixed per GDAL ticket https://trac.osgeo.org/gdal/ticket/7091