floating point rounding errors in OGR

[TobWen]

There seems to be a problem with precision in several OGR-core functions. I'll give an example based on a conversion from GML to GeoJSON. I've already discussed this with other users and it was possible to reproduce the problem.

The following data is released under "CC-BY-SA 4.0, © Metropole Ruhr". It's a mapserver, so the namespace starts with ms, followed by a colon, f.e. ms:poi_hoehlen.

curl -o demo.gml "http://geodaten.metropoleruhr.de/poi/poi?SERVICE=WFS&VERSION=2.0.0&REQUEST=GetFeature&TYPENAMES=ms:poi_hoehlen&STARTINDEX=0&COUNT=100"

Let's check a random POI's coordinate:

$ grep -m 1 '"fid="ms_poi_hoehlen.2"' -A2 demo.gml
<ogr:ms_poi_hoehlen fid="ms_poi_hoehlen.2">
<ogr:geometryProperty><gml:Point srsName="EPSG:25832"><gml:coordinates>414869.165397,5692651.751954</gml:coordinates></gml:Point></ogr:geometryProperty>

And now we convert it to GeoJSON:

$ ogr2ogr -F GeoJSON demo.geojson demo.gml

Now check the coordinates again:

$ grep -m 1 'poi_hoehlen.2' -A2 demo.geojson
"geometry": { "type": "Point", "coordinates": [ 414869.16539699998, 5692651.751954

Ooops. Yeah, we all like floating point errors, don't we? Sure, it's easy to fix the precision to 3 (for millimeters at metric coordinates), but when changing the CRS, f.e. to geographic coordinates, you might run into trouble. Think of files with a mixture of cartesian metric and geographic coordinates ... fixing the precision will be a mess.

From my point of view, when there is no reprojection, the precision should NOT be changed. Internally, those coordinates gets stored as floats, but there is a way to fix this problem. I'm using it in PHP and Python when not using a high-preicision math library:

I'm checking the precision of a number, f.e. 414869.165397 => precision = 6 and I'm storing it in an array: $number = array(414869.165397, 6); When the number is about to be written to the file, I'm rounding the output by the stored value: echo round($number[0][0], $number[0][1]); the result is 414869.165397 again.

[Even Rouault]

That's a known issue with no easy solution. You are converting between text formats but internally coordinates are stored as binary IEEE754 double numbers, so you have no idea of the original precision. As a workaround you can use the COORDINATE_PRECISION layer creation option of the GeoJSON driver if you know the input precision.

[TobWen]

That's a known issue with no easy solution. You are converting between text formats but internally coordinates are stored as binary IEEE754 double numbers, so you have no idea of the original precision.

I've presented an easy to implement and working approach to solve this problem. Since there's no processing, the array (or whatever kind of storage) can be flushed directly after the output. It's damn much easier than a high-precision math engine.

Another easy approach would be to introduce an addition option to store coordinates as a string for GeoJSON. Sure, this will break compatibility at first sight, but one could easily remove the quotation-marks (perhaps your output function is already capable to do so.).

As a workaround you can use the COORDINATE_PRECISION layer creation option of the GeoJSON driver if you know the input precision.

Pardon to be immune to your feedback, but that's not a valid workaround in terms of geoprocessing. Let's say, you're doing batch processing with EPSG:4326 (geographic) and and EPSG:25832 (metric). Cutting the metric coordinates to precision 3 will be fine, but will badly fail on EPSG:4326. Also, when there's no change in projection, I want to copy the raw input to the output without touching the data at all.

[Even Rouault]

This ticket has been automatically closed because Trac is no longer used for GDAL bug tracking, since the project has migrated to GitHub. If you believe this ticket is still valid, you may file it to ​https://github.com/OSGeo/gdal/issues if it is not already reported there.