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

classic Classic list List threaded Threaded
3 messages Options
Reply | Threaded
Open this post in threaded view
|

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

Landry Breuil-2
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
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: setting a spatial filter on a postgis layer gives wrong results

Even Rouault-2
On vendredi 5 juillet 2019 10:12:41 CEST Landry Breuil wrote:

> 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.

This indeed exactly this. The PG driver implements spatial filtering only
based on the envelope of the filter geometry, whereas the MEM driver uses the
full geometry

Quoting
https://gdal.org/api/
ogrlayer_cpp.html#_CPPv4N8OGRLayer16SetSpatialFilterEP11OGRGeometry
"""
Currently this test is may be inaccurately implemented, but it is guaranteed
that all features whose envelope (as returned by OGRGeometry::getEnvelope())
overlaps the envelope of the spatial filter will be returned. This can result
in more shapes being returned that should strictly be the case.
"""

I guess the PG driver should be fixed. This is probably just an historical
artifact.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: setting a spatial filter on a postgis layer gives wrong results

jratike80
Hi,

It should be already possible to make accurate filtering from PostGIS by
executing SQL with ST_Intersects.

-Jukka Rahkonen-


Even Rouault-2 wrote

> On vendredi 5 juillet 2019 10:12:41 CEST Landry Breuil wrote:
>> 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.
>
> This indeed exactly this. The PG driver implements spatial filtering only
> based on the envelope of the filter geometry, whereas the MEM driver uses
> the
> full geometry
>
> Quoting
> https://gdal.org/api/
> ogrlayer_cpp.html#_CPPv4N8OGRLayer16SetSpatialFilterEP11OGRGeometry
> """
> Currently this test is may be inaccurately implemented, but it is
> guaranteed
> that all features whose envelope (as returned by
> OGRGeometry::getEnvelope())
> overlaps the envelope of the spatial filter will be returned. This can
> result
> in more shapes being returned that should strictly be the case.
> """
>
> I guess the PG driver should be fixed. This is probably just an historical
> artifact.
>
> Even
>
> --
> Spatialys - Geospatial professional services
> http://www.spatialys.com
> _______________________________________________
> gdal-dev mailing list

> gdal-dev@.osgeo

> https://lists.osgeo.org/mailman/listinfo/gdal-dev





--
Sent from: http://osgeo-org.1560.x6.nabble.com/GDAL-Dev-f3742093.html
_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev