1

I'm trying to intersect a polyline and polygon and write out the resulting polyline shapefile using ogr2ogr

Here is the command I've been trying: 

 ogr2ogr -dialect SQLITE -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM polygon A, polyline B WHERE ST_Intersects(A.geometry, B.geometry)" test.shp input.vrt

Here are the contents of the input.vrt file:

<OGRVRTDataSource> <OGRVRTLayer name="polygon"> <SrcDataSource>polygon.shp</SrcDataSource> </OGRVRTLayer> <OGRVRTLayer name="polyline"> <SrcDataSource>polyline.shp</SrcDataSource> </OGRVRTLayer> </OGRVRTDataSource>

The resulting output is a .dbf file called "test" with all the desired fields, but no geometries. There is no accompanying shapefile with the same name. I've tried modifications of the command, particularly switching the order of the intersection. I know the features do in fact intersect, as I've done the intersection in ArcGIS and it works, producing 60+ features.

caveat - the polyline has ~1500 features whereas the polygon has been dissolved, and only contains one feature. I don't think this should matter.

I know there are other libraries to accomplish this, but I have an existing backend system built on gdal/ogr and I am looking to reduce dependencies to other libraries.

UPDATE: Debugging based on comments from @User30184.Here are is the output for one feature with the following ogrinfo call:

GEOS error: TopologyException: Input geom 0 is invalid: Self-intersection at or near point 97.562017383948856 16.504292693814882 at 97.562017383948856 16.504292693814882 OGRFeature(SELECT):1 geometry (String) = (null) DN (Real) = 1 osm_id (String) = 27167387 code (Integer) = 5112 fclass (String) = trunk name (String) = (null) ref (String) = 408 oneway (String) = F maxspeed (Integer) = 0 layer (Real) = 0 bridge (String) = F tunnel (String) = F

Looks like the geometry field is not being computed. Any ideas?

Improve this question
3
  • 1
    Could you test your command by modifying A.*, B.* as I think this will select the Geometries (geometry) from both feature types. Perhaps you could try selecting one attribute from each input dataset, e.g. A.attribute_1, B.attribute_1 Commented Sep 1, 2016 at 14:19
  • Just tested and result is the same :-( the only file that is generated is a dbf file, except with only the specified fields instead of all fields Commented Sep 1, 2016 at 14:25
  • There seems to be some issue with your data but without having a sample it is impossible to say anything more. Commented Sep 1, 2016 at 21:12

2 Answers 2

Reset to default
1

Your VRT and command worked for me based on the example features below. I would suggest creating a subset of your inputs to try to further identify the issue.

polygon_polyline

polygon_polyline_intersection

Improve this answer
1

Your SQL gives a bit different result than you believe, because of a.* and b.* the result contains three geometries.

I made a test with these data:

POLYGON (( 220 440, 429 546, 540 220, 200 220, 220 440 )) LINESTRING ( 240 600, 420 120 )

I saved them into polygon.shp and polyline.shp and copied your vrt file. Next a test with ogrinfo:

ogrinfo -dialect sqlite -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, A.*, B.* FROM polygon A, polyline B WHERE ST_Intersects(A.geometry, B.geometry)" input.vrt INFO: Open of `input.vrt' using driver `OGR_VRT' successful. Layer name: SELECT Geometry (geometry): Unknown (any) Geometry (GEOMETRY): Polygon Geometry (GEOMETRY): Line String Feature Count: 1 Extent (geometry): (287.216080, 220.000000) - (382.500000, 474.090452) Extent (GEOMETRY): (200.000000, 220.000000) - (540.000000, 546.000000) Extent (GEOMETRY): (240.000000, 120.000000) - (420.000000, 600.000000) SRS WKT (geometry): (unknown) SRS WKT (GEOMETRY): (unknown) SRS WKT (GEOMETRY): (unknown) Geometry Column 1 = geometry Geometry Column 2 = GEOMETRY Geometry Column 3 = GEOMETRY OGRFeature(SELECT):0 geometry = LINESTRING (287.21608040201 474.090452261307,382.5 220.0) GEOMETRY = POLYGON ((220 440,429 546,540 220,200 220,220 440)) GEOMETRY = LINESTRING (240 600,420 120)

Instead of selecting all attributes with A.* and B.* it might be better to list explicitly the attributes. However, with my data this command does create a shapefile with geometry:

ogr2ogr -f "ESRI Shapefile" -dialect sqlite -sql "SELECT ST_Intersection(A.geometry, B.geometry) AS geometry, a.*, b.* FROM polygon A, po lyline B WHERE ST_Intersects(A.geometry, B.geometry)" test.shp input.vrt

The intersection of a polygon and a linestring is linestring

 LINESTRING (287.21608040201 474.090452261307,382.5 220.0)

I am not sure if that is what you would like to get. If you would like to get polygons split use split function.

Improve this answer
4
  • My desired output is a line, and not a polygon! The lines that @dmci has shown in the answer below are what I'm after. Just tried forcing -f "ESRI Shapefile" as you did in your command, but no luck. More testing underway Commented Sep 1, 2016 at 15:24
  • UPDATE: I regenerated the data, and now performing the intersection gives me the correct .dbf file (matches features generated from testing in ArcGIS). However, there is no accompanying shapefile or other supporting files (.shp, .prj, .shx). The "geometry" field in the .dbf file is blank. Commented Sep 1, 2016 at 16:27
  • 1
    Use ogrinfo as in my example for avoiding possible problems with writing the shapefile out. Adding --debug on may give more info. Use LIMIT 100 because you do not want to get a full list. Commented Sep 1, 2016 at 16:53
  • Updated question with sample output of ogrinfo for one feature. Commented Sep 1, 2016 at 17:06

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.