[GIS] Selecting features within polygon from another layer in python

ogrogr2ogrqgissqlite

Someone asked a nearly identical question about selecting features in a layer by location here Selecting features within polygon from another layer using R or Python?, but the answers are not for general gdal/ogr python users. I'm wondering if anyone has a pythonic solution to this problem that does not require qgis modules. In particular, I'd like to be able to specify the options qgis offers (include input features that touch and/or overlap/cross and/or are completely within and/or are selected).

I've been trying all day to remove polygons from a shapefile that share a border with (but are contained completely within) a polygon in another shapefile. This task takes <5min in QGIS, but after 5+ hours I still haven't found a way to "select by location" and specify "options" using just ogr (I've tried using layer.SetSpatialFilter(polygon_geometry), but the docs specify that this only works with GetNextFeature(), which I am not using to loop through my features, so I was not terribly surprised when it did not produce the desired results. Also, I suspect that features sharing only a border with the selection layer are included rather than excluded). I've also tried ogr2ogr

ogr2ogr -f 'ESRI Shapefile' -dialect sqlite -sql "SELECT * FROM bergs WHERE ST_Intersects(bergs.geometry, AOI_buffered100.geometry)" -overwrite bergs_testborderrem.shp bergs.shp

but cannot get a working SQL statement (I keep getting no such column errors, but all the ogr-sql examples I've found use .geometry, so I have no idea how to query on the geometry otherwise).

Does anyone have any suggestions?

Best Answer

Maybe something like 

ogr2ogr -f 'ESRI Shapefile' -dialect sqlite -sql "SELECT a.* FROM bergs a, 
AOI_buffered100 b WHERE ST_Intersects(a.geometry, b.geometry)" -overwrite 
bergs_testborderrem.shp sub_directory

where the shapefiles are in sub_directory. Also you may want ST_Touches or something instead of ST_Intersects.

The key seems to be the specification of the directory as the source, though relative directories appear to be possible, at least for the second shapefile, as long as there are shapefiles in sub_directory. Finally, make sure that the name of the geometry columns are correct.

Also, see here.

Related Question