[GIS] GDAL/OGR ogr.Layer.Intersection fails when producing mixed geometry results

Question

I am trying to intersect layers using the ogr.Layer.Intersection method, in Python. This works fine for pairs of polygons, but I have found intersecting two linestring layers containing matching points can return a combination of linestrings and points, producing an invalid layer for many single-geometry type data sources e.g. Shapefile. How can I avoid the non-matching geometries without going down to the Feature level and filtering them myself?

Simplified example below, which demonstrates a POINT geometry is returned from the two linestrings.

from osgeo import ogr, osr

line1 = "LINESTRING ({0} {0}, {1} {0})".format(0, 1)
line2 = "LINESTRING ({1} {0}, {2} {0})".format(0, 1, 2)

sr = osr.SpatialReference()
sr.ImportFromEPSG(27700)

driver = ogr.GetDriverByName('MEMORY')
driver2 = ogr.GetDriverByName('ESRI Shapefile')

ds1 = driver.CreateDataSource('ds1')
ds2 = driver.CreateDataSource('ds2')
ds3 = driver.CreateDataSource('path to shapefile.shp')

layer1 = ds1.CreateLayer('layer1', sr, ogr.wkbLineString)
layer2 = ds2.CreateLayer('layer2', sr, ogr.wkbLineString)
layer3 = ds3.CreateLayer('layer3', sr, ogr.wkbLineString)

geom1 = ogr.CreateGeometryFromWkt(line1)
feat1 = ogr.Feature(layer1.GetLayerDefn())
feat1.SetGeometry(geom1)

geom2 = ogr.CreateGeometryFromWkt(line2)
feat2 = ogr.Feature(layer2.GetLayerDefn())
feat2.SetGeometry(geom2)

layer1.CreateFeature(feat1)
layer2.CreateFeature(feat2)

layer1.Intersection(layer2, layer3) # This fails, as a point is produced.

UPDATE: The exact error message received if ds3 is referring to an ESRI Shapefile on disk is "ERROR 1: Attempt to write non-linestring (POINT) geometry to ARC type shapefile.". So my main question is: how to avoid the POINT feature being produced at all, or to filter them before the layer is written to the DataSource?

UPDATE #2: Just to be clear, in my real scenario, I am not wanting to intersect two adjacent linestrings, I actually want to clip a roads Shapefile using a polygon of my study area. Coincidentally, it seems one of the road vertex points and the edge of my polygon coincide, producing the error. It worked for all the other layers I've clipped. I do not want the non-matching geometries, and believe the Intersection method should have an option to drop them if the layer has been specified as a particular type.

UPDATE #3: As a working example, consider arcpy.Intersect_analysis:

arcpy.Intersect_analysis(['shape1.shp', 'shape2.shp'], 'shape3.shp')

Even with features that intersect at a point, no multi-geometry outputs are being produced here, they must be filtered out or not produced at all – how can I get the same behavior from Python GDAL/OGR (not shapely or a database) at an ogr.Layer level?

Answer 1

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)'