[GIS] ST_Intersect with ogr2ogr and Spatialite

intersectionogr2ogrpolygonsqlite

I am trying to intersect Spatialite datasets with ogr2ogr, the OGR Virtual Format and SQL Syntax. I can't seem to locate whether the error is in my SQL query or my VRT file.

Data

rice_temperature.sqlite | Spatialite database containing 1 layer "rice_temp"

rice_elevation.sqlite | Spatialite database containing 1 layer "rice_elev"

rice_vectors.vrt

<OGRVRTDataSource>
<OGRVRTLayer name="temperature">
    <SrcDataSource>rice_temperature.sqlite</SrcDataSource>
</OGRVRTLayer>
<OGRVRTLayer name="elevation">
    <SrcDataSource>rice_elevation.sqlite</SrcDataSource>
</OGRVRTLayer>

All in the same directory I am also executing ogr2ogr from.

Query

ogr2ogr -f "ESRI Shapefile" intersection.shp rice_vectors.vrt -dialect sqlite -sql "SELECT ST_Intersection(rice_temperature.geometry,rice_elevation.geometry) AS geometry FROM rice_temp,rice_elev"

Errors

I tried a number of variations to the VRT file as well as to the SQL query itself. All errors are either no such table: rice temp or Cannot open datasource 'rice_elevation'

I have a feeling that I am close but can't wrap my head around what I am doing wrong – is it the SQL query or my use of OGRs Virtual Format?

Best Answer

With your VRT file GDAL tries to find layers "temperature" and "elevation" from the target databases. Either use the original layer names in VRT in OGRVRTLayer name: "rice_temp" and "rice_elev", of rename them with SrcLayer

<OGRVRTDataSource>
<OGRVRTLayer name="temperature">
    <SrcDataSource>rice_temperature.sqlite</SrcDataSource>
    <SrcLayer>rice_temp</SrcLayer>
</OGRVRTLayer>
<OGRVRTLayer name="elevation">
    <SrcDataSource>rice_elevation.sqlite</SrcDataSource>
    <SrcLayer>rice_elev</SrcLayer>
</OGRVRTLayer>
</OGRVRTDataSource>

Test with ogrinfo if your VRT is correct and two layers are found.

When you have created a VRT what is behind is hidden from ogr2ogr. The only layers you can use in SQL are the ones you have defined in the VRT. GDAL will take care of the rest: you access "temperature" and GDAL knows to read "rice_temp" from "rice_temperature.sqlite".

ogr2ogr -f "ESRI Shapefile" intersection.shp rice_vectors.vrt -dialect sqlite -sql "SELECT ST_Intersection(a.geometry, b.geometry) AS geometry FROM temperature a, elevation b"
Related Question