[gdal-dev] Geopackage precision issues

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

[gdal-dev] Geopackage precision issues

Markus Metz-3
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

                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

_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev
Reply | Threaded
Open this post in threaded view
|

Re: Geopackage precision issues

Even Rouault-2
Makus,

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

Ah indeed. Fixed per
https://github.com/OSGeo/gdal/commit/0f28425fe8a2c70061b8896748f40182aff3f2ee

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

Yes, that's expected: your miny_filter > y_geometry. You probably meant
miny=317579.582420

> Not the reason for this, but another precision issue is at
> https://github.com/OSGeo/gdal/blob/master/gdal/ogr/ogrsf_frmts/gpkg/ogrgeopa
> ckagetablelayer.cpp#L3544
>
>                 SQLEscapeName(pszC).c_str(), sEnvelope.MinX - 1e-11,

I'm not sure if those bbox extension are needed at all. They don't seem
necessary if trying

printf 'id,WKT\n1,"POINT(176418.99585285 317579.62751027)"\n' > point.csv
ogr2ogr point.gpkg point.csv
ogrinfo point.gpkg -al \
    -spat 176418.99585285 317579.62751027 176418.99585285 317579.62751027

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: Geopackage precision issues

Markus Metz-3


On Tue, May 21, 2019 at 10:16 PM Even Rouault <[hidden email]> wrote:

>
> Makus,
>
> > 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.
>
> Ah indeed. Fixed per
> https://github.com/OSGeo/gdal/commit/0f28425fe8a2c70061b8896748f40182aff3f2ee

Thanks for the quick fix! BTW, there are many more instances of %g in other OGR drivers, not sure if some of them need adjustment as well.
>
> > 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.
>
> Yes, that's expected: your miny_filter > y_geometry. You probably meant
> miny=317579.582420

Oops, yes of course, that works, that means it was only the extents suffering from loss of precision.

Markus M


_______________________________________________
gdal-dev mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/gdal-dev