[gdal-dev] Geopackage precision issues

[Markus Metz]

Extents and spatial filter settings with Geopackage are unprecise such that
they are unusable.

For example I created a Geopackage with a single point: 176418.99585285,
317579.62751027
ogrinfo reports extents of
Extent: (176419.000000, 317580.000000) - (176419.000000, 317580.000000)
I think this can happen if the coordinates are converted to float or
converted to a string with something like %.6g or %.g.
Thus the point in the Geopackage is outside the extents of the Geopackage.
If I adjust the extents of the Geopackage to
minx: 176418.823581
maxx: 176419.176419
miny: 317579.682420
maxy: 317580.317580
and use this as a spatial filter, the point is still excluded. A possible
explanation is that the adjusted extents in the spatial filter are
similarly converted, resulting in
minx: 176419.000000
maxx: 176419.000000
miny: 317580.000000
maxy: 317580.000000
which excludes the point in the Geopackage.

Not the reason for this, but another precision issue is at
https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopackagetablelayer.cpp#L3544

                SQLEscapeName(pszC).c_str(), sEnvelope.MinX - 1e-11,

with sEnvelope.MinX = 176418.823581 or similar, subtracting 1e-11 will not
make a difference, it should rather be
                SQLEscapeName(pszC).c_str(), sEnvelope.MinX -
fabs(sEnvelope.MinX * 1e-11),
same for following lines. If this is again converted to float or converted
to a string with something like %.6g or %.g, this will not work and you
need something like
                SQLEscapeName(pszC).c_str(), sEnvelope.MinX -
fabs(sEnvelope.MinX * 2e-6),
to account for single precision floating point limits.

I would like to provide a fix, but I did not find out where this conversion
takes place.

Markus M
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20190521/350cb369/attachment.html>