[gdal-dev] setting a spatial filter on a postgis layer gives wrong results

[Landry Breuil]

Hi,

initially saw this with gdal 2.2.3 on debian, and confirmed the 
behaviour is the same with gdal 3.0 on openbsd.

I'm doing a spatial join between an union of geometries filtered from a 
geopackage and a tileindex layer stored in postgis, and the spatial 
filter set on the postgis layer behaves weird, the provided result is as 
if the bbox of the spatial filter/geometry was used, instead of the 
filter itself.

applying the exact same filter to a copy of the postgis layer stored in 
memory gives correct results.

below my testing code:

from osgeo import gdal
from osgeo import ogr

dbhandle = gdal.OpenEx("PG: host=localhost dbname=gisdata-local 
user=postgis", gdal.OF_VERBOSE_ERROR, allowed_drivers = ['PostgreSQL'])
gpkg = gdal.OpenEx('/tmp/dbdemogeor.gpkg')
communes = gpkg.GetLayerByName('adminexpress_201806_commune')
assert (dbhandle is not None)
layer = dbhandle.GetLayer('rtge.ortho')
print("{} tiles total, {} communes".format(layer.GetFeatureCount(), 
communes.GetFeatureCount()))

outdriver=ogr.GetDriverByName('MEMORY')
source=outdriver.CreateDataSource('memData')
tmp=outdriver.Open('memData',1)
source.CopyLayer(layer, 'out')
layer2 = source.GetLayer('out')

communes.SetAttributeFilter("insee_dep='63'")

multi  = ogr.Geometry(ogr.wkbMultiPolygon)
# compute the multigeometry for all the comms in the area list
for feature in communes:
     g = feature.GetGeometryRef()
     for i in range(0, g.GetGeometryCount()):
         multi.AddGeometry(g.GetGeometryRef(i))
communes.ResetReading()
union = multi.UnionCascaded()
# filter tileindex on a 100m bufferized union of the comms
layer.SetSpatialFilter(union.Buffer(100))

print("{} tiles are inside the geomunion for {} comms ({} 
surf)".format(layer.GetFeatureCount(), communes.GetFeatureCount(), 
union.GetArea()))

print("{} tiles in memory layer".format(layer2.GetFeatureCount()))
layer2.SetSpatialFilter(union.Buffer(100))
print("{} tiles in memory layer after 
filtering".format(layer2.GetFeatureCount()))


And the resulting output:
405557 tiles total, 4092 communes
131055 tiles are inside the geomunion for 467 comms (8001479208.5 surf)
405557 tiles in memory layer
82868 tiles in memory layer after filtering

is my code wrong, or the handling of spatial filters for postgis layers 
should be done in a different way, ie directly using postgis methods 
instead of trying to use gdal APIs?
-- 
Landry Breuil
Administrateur de données du CRAIG

----------------------------------------------------------------------------
Centre Régional Auvergne-Rhône-Alpes de l'Information Géographique
Bât. du CRRI - Campus des Cézeaux
7 avenue Blaise Pascal - CS 60026
63178 Aubière

https://www.craig.fr - @GipCraig
----------------------------------------------------------------------------
Support utilisateurs (tous les jours ouvrés de 8H30 à 12H30) : 04 73 405 405