[GIS] Intersect polyline and polygon with ogr

gdalogrpythonsql

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?

Best Answer

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.

Related Solutions

[GIS] Simplify polygon in OGR

You must supply the tolerance, whether in Arc or OGR. It is the measure of your simplification operation and is just a decimal number in the the same units as the SRS of your data.

In Arc ESRI stores an XY tolerance as a separate value. This is proprietary to ESRI and not inherent in the SRS. Either you specify it or ESRI supplies a default - either way it is the same thing as you explicitly supplying a value to OGR. ESRI calculates the default but there is nothing sacred about the value. It's just a suggestion and may not be relevant to your particular data depending on its provenance, collection methods etc.

[GIS] Merging multiple feature classes in file geodatabase with ogr2ogr

To select a single feature class from a gdb, we query for all entities inside it:

-sql "select * from FEATURE_CLASS_NAME"

To generate a list of feature classes following your sample use nested FOR loops:

@echo off
for %%S in (01 02 03) do (
    for %%F in (01 02 03) do (
        echo ogr2ogr out.gdb in.gdb -sql "select * from Sec%%S_Frm%%F"
        )
     )

Emits:

ogr2ogr out.gdb in.gdb -sql "select * from Sec01_Frm01"
ogr2ogr out.gdb in.gdb -sql "select * from Sec01_Frm02"
ogr2ogr out.gdb in.gdb -sql "select * from Sec01_Frm03"
ogr2ogr out.gdb in.gdb -sql "select * from Sec02_Frm01"
...

From there add -nln OUT_LAYER_NAME to say we want it all going to one place, and -nlt multipolygon for geometry type. (NB: nlt is often automatically derived from input, but for some operations like this one geometry is omitted and you get just a table):

ogr2ogr out.gdb in.gdb -nln Form01 -nlt multipolygon -sql ...

And finally we will be updating existing data so:

ogr2ogr out.gdb in.gdb -update -append -nln ...

Useful options:

-config FGDB_BULK_LOAD YES
-skipfailures

How do I include the source filename in a field when merging hundreds of shapefiles (Windows)?

for %f in (*.shp) do (
    ogr2ogr -sql "select *, '%~nf' as s_file from %~nf" -update -append merged\output.shp %f -nln output
    )

Sources

By the way, if you still have the gdb with the "unwieldy number of feature classes" hanging around from previous attempts you might be able to just ogr2ogr -f FileGDB clean.gdb unwieldy.gdb -nln Clean_Form01 -nlt multipolygon in one go.

Code samples are in Windows batch file syntax. Thank you for this question. It prompted me to extend a related project that I've been needing to come back to.

I don't think ogrinfo can write, so that would be the origin of the lock error in the Dec 11 update.

Best Answer

Related Solutions

[GIS] Simplify polygon in OGR

[GIS] Merging multiple feature classes in file geodatabase with ogr2ogr

Related Question