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?

  • 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 – dmci Sep 1 '16 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 – geoeye Sep 1 '16 at 14:25
  • There seems to be some issue with your data but without having a sample it is impossible to say anything more. – user30184 Sep 1 '16 at 21:12
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

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.

  • 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 – geoeye Sep 1 '16 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. – geoeye Sep 1 '16 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. – user30184 Sep 1 '16 at 16:53
  • Updated question with sample output of ogrinfo for one feature. – geoeye Sep 1 '16 at 17:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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