I'm trying to intersect a polyline and polygon and write out the resulting polyline shapefile using
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?