2

I have two datasets, one of which contains millions of postcode points and the other which includes polygon data - county borders in this instance. I'd like to write an ogr2ogr command to export only the points which are contained by a specific polygon.

My source data is within separate SQLite databases, and my output needs to be a ShapeFile.

I've written a VRT file which produces the correct output:

<OGRVRTDataSource>
    <OGRVRTLayer name="pcodes">
        <SrcDataSource>MyPostcodes.sqlite</SrcDataSource>
        <SrcLayer>PostcodeData</SrcLayer>
    </OGRVRTLayer>
    <OGRVRTLayer name="polys">
        <SrcDataSource>Counties.sqlite</SrcDataSource>
        <SrcSql>
            SELECT GEOMETRY
            FROM CountyTable
            WHERE countyname="Nottinghamshire"
        </SrcSql>
    </OGRVRTLayer>
    <OGRVRTLayer name="res">
        <SrcDataSource>test.vrt</SrcDataSource>
        <SrcSql dialect="SQLITE">
            SELECT pcodes.GEOMETRY FROM pcodes, polys
            WHERE Within(pcodes.GEOMETRY, polys.GEOMETRY) = 1
        </SrcSql>
    </OGRVRTLayer>
</OGRVRTDataSource>

My ogr2ogr command is:

ogr2ogr -f "ESRI Shapefile" pointsinpoly.shp test.vrt -overwrite -sql "SELECT * FROM res"

Although this produces the correct ShapeFile, it takes a few minutes to run. I realise that I need to use the spatial indexes, but I'm stumped as to what the correct syntax is.

I have a SpatialIndex table on each of the referenced databases and I've confirmed they're both valid with select CheckSpatialIndex().

  • I'm not sure that SQLite has a powerful enough SQL engine to quickly handle millions of spatial calculations. Have you considered PostgreSQL/PostGIS or something similar? – haff Dec 2 '16 at 19:20
  • Can you make the data (or at least a representative sample) available? Its relatively easy to make a mistake that is easier to debug if we have something to test on. As a start, are your county polygons really within your postcode geometries? – BradHards Dec 4 '16 at 0:09
  • @haff I've also tried SQL Server, but performance was worse. If SQLite doesn't pan out, I'll definitely try PostGIS. Thanks! – Dave R. Dec 5 '16 at 10:33
  • @BradHards Good catch! I've transposed the Within parameters, which were incorrect above. My actual VRT already had the correct ordering, so I'm afraid this isn't the source of the performance issue. It's unlikely that I'll be able to provide an extract of the data, as it's not for public consumption. Would the structure of the two DBs help? In the end, it's just the spatial index sub-query syntax I need. – Dave R. Dec 5 '16 at 13:17
0

The following VRT worked, and the spatial index is now utilised.

<OGRVRTDataSource>
    <OGRVRTLayer name="pcodes">
        <SrcDataSource>MyPostcodes.sqlite</SrcDataSource>
        <SrcLayer>PostcodeData</SrcLayer>
    </OGRVRTLayer>
    <OGRVRTLayer name="polys">
        <SrcDataSource>Counties.sqlite</SrcDataSource>
        <SrcSql>
            SELECT GEOMETRY
            FROM CountyTable
            WHERE countyname="Nottinghamshire"
        </SrcSql>
    </OGRVRTLayer>
    <OGRVRTLayer name="res">
        <SrcDataSource>test.vrt</SrcDataSource>
        <SrcSql dialect="SQLITE">
            SELECT pcodes.GEOMETRY
            FROM pcodes, polys
            WHERE Within(pcodes.GEOMETRY, polys.GEOMETRY) = 1
            AND pcodes.ROWID IN (
                SELECT ROWID
                FROM SpatialIndex
                WHERE f_table_name = 'pcodes'
                AND search_frame = polys.GEOMETRY
            )
        </SrcSql>
    </OGRVRTLayer>
</OGRVRTDataSource>

For my test dataset, the time taken to do the shape export went from nearly 9 minutes to 49 seconds.

However, it appears that this creates an issue with ogr2ogr trying to join on an FID column on the point DB; this doesn't exist on our spatial DBs as we use the default OGC_FID. Adding a new FID column as a copy of OGC_FID fixed the problem, but it's obviously not ideal! I'll raise this as a separate question.

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.