4

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?

Improve this question
8
  • It is not an error because the intersection is a point. You will get the same result with PostGIS postgis.net/docs/ST_Intersection.html. Filter by geometry type if you can't use format that supports mixed geometries like GML, Spatialite, PostGIS (remember to use general GEOMETRY type), OpenJUMP JML, MapInfo, Oracle etc.
    – user30184
    Apr 27, 2016 at 9:46
  • @user30184 How do I filter by geometry type with the Intersection method?
    – ssast
    Apr 27, 2016 at 9:59
  • @user30184 Also, any idea why this doesn't happen when intersecting two adjacent polygons? e.g. even if they share a vertex, no linestrings or points are being produced by Intersection, only with Clip?
    – ssast
    Apr 27, 2016 at 10:25
  • 1
    I would say that write intersections into memory layer and filter what you want before writing the shapefile. And ogr2ogr has an option -skipfailures, it must be usable somehow with Python as well.
    – user30184
    Apr 27, 2016 at 10:42
  • @user30184 Not too bad a workaround (better than intersecting individual features, anyway). My only concern would be working with two particularly large layers at a later stage, but for my immediate problem, it would suffice.
    – ssast
    Apr 27, 2016 at 10:53

2 Answers 2

Reset to default
3

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

enter image description here

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

enter image description here

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

enter image description here

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

enter image description here

Improve this answer
1
  • The problem isn't so much that a point is produced, I understand that. The issue is that the function described (which is not shapely or PostGIS etc) dutifully produces an invalid layer and fails, with no simple option to avoid it.
    – ssast
    Apr 28, 2016 at 9:15
1

Not a proper answer but does not fit into the comment box. However, at least PostGIS does return a point for two polygons which touch at one point.

SELECT ST_AsText(
ST_Intersection(ST_GeomFromText(
'POLYGON (( 140 360, 140 480, 220 480, 220 360, 140 360 ))'),
ST_GeomFromText(
'POLYGON (( 220 260, 220 360, 300 360, 300 260, 220 260 ))')));

"POINT(220 360)"

Improve this answer
2
  • Same with Spatialite (my preferred choice for SQL). However, with the Python bindings of ogr.Layer.Intersection, an intersection of those two polygons above produces None, but ogr.Layer.Clip yields the point feature. Strangely inconsistent that the intersected point features are apparently filtered from polygons, but not linestrings.
    – ssast
    Apr 27, 2016 at 10:50
  • It can also be just a bug. What "intersection" means is pretty accurately defined.
    – user30184
    Apr 27, 2016 at 11:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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