I'm using GDAL (3.3.1) Python (3.8.11) bindings to process vector layers and am trying to use OGR SQL dialect to filter the layer, like this:
sql = f"SELECT * from {out_lyr.GetName()} WHERE OGR_GEOM_AREA < {area_threshold}"
out_lyr_selection = out_ds.ExecuteSQL(sql)
This works fine on FileGDB and Shapefile, when I try to use it on GeoPackage it gives the following error.
RuntimeError: In ExecuteSQL(): sqlite3_prepare_v2(SELECT * from merged WHERE OGR_GEOM_AREA < 20.0):
no such column: OGR_GEOM_AREA
This post suggests I should use "indirect_sqlite" dialect, so I checked the OGR doc page. It suggested this would work:
sql = f"SELECT * from {out_lyr.GetName()} WHERE ST_Area(GEOMETRY) < {area_threshold}"
out_lyr_selection = out_ds.ExecuteSQL(sql, dialect="indirect_sqlite")
However, this time the error was:
RuntimeError: In ExecuteSQL(): sqlite3_prepare_v2(SELECT * from merged WHERE ST_Area(GEOMETRY) < 20.0):
no such column: GEOMETRY
Am I doing something wrong, or is it not possible to use OGR SQL dialect to filter by area with GeoPackage?
Best Answer
For OGR SQL use
dialect="OGRSQL"
. It is the default SQL dialect but only when the datasource has no support for a native SQL dialect. Therefore if you read PostGIS, SpatiaLite, GeoPackage, Oracle etc. you must explicitly select the ORG SQL dialect.When using indirect_sqlite dialect pay attention to error message
no such column: GEOMETRY
. Check the real name of the geometry column with ogrinfo. It can be whatever but quite often you will seeGeometry Column = geom