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