Change in coordinate precision with ogr2ogr moving from GDAL 2 to 3

Question

I'm migrating some old code from GDAL 2 to GDAL 3. I use ogr2ogr to import a shapefile into PostGIS like this:

ogr2ogr -f PostgreSQL -nln mwe PG:"host=localhost user=myuser" mwe.shp -overwrite

I've saved a copy of an MWE single feature shapefile here if you'd like to follow along. Plain text would make things easier, but the error only seems to occur when I write data from shapefiles to Postgres.

Expected results - GDAL2

When using gdal2, I can query this data in postgres and see the coordinate precision as expected:

select ST_Y(geom)
from (
  select (ST_DumpPoints(wkb_geometry)).geom
  from mwe
) pts;

yields:

Full info for 'working' stack:

ubuntu 18.04
postgres 11
postgis 3
GDAL 2.2.3, released 2017/11/20

POSTGIS="3.1.1 aaf4c79" [EXTENSION] PGSQL="110" GEOS="3.7.1-CAPI-1.11.1 27a5e771" PROJ="Rel. 4.9.3, 15 August 2016" LIBXML="2.9.4" LIBJSON="0.12
.1" LIBPROTOBUF="1.2.1" WAGYU="0.5.0 (Internal)"

Actually results - GDAL3

When I run this using GDAL 3, some vertices have inflated coordinate precision:

Full info for not working GDAL stack:

ubuntu 20.02
postgres 13
postgis 3
GDAL 3.0.4, released 2020/01/28

 POSTGIS="3.1.1 aaf4c79" [EXTENSION] PGSQL="130" GEOS="3.8.0-CAPI-1.13.1 " PROJ="6.3.1" LIBXML="2.9.10" LIBJSON="0.13.1" LIBPROTOBUF="1.3.3" WAGYU="0.5.0 (Internal)"

ogrinfo coordinate precision

When running ogrinfo on my MWE shapefile using GDAL 3.0.4, I also see the expected coordinate precision:

ogrinfo -al mwe.shp

INFO: Open of `mwe.shp'
      using driver `ESRI Shapefile' successful.
. . .
 POLYGON ((21.04237 39.579149,21.038807 39.577435,21.036219 39.576187,21.033495 39.574875,21.03046 39.573414,21.024775 39.570475,21.023056 39.569588,21.017654 39.563602,21.009989 39.563313,21.008549 39.563259,21.001377 39.565672,21.001286 39.565702,21.004169 39.575203,21.006231 39.581999,21.013037 39.583546,21.021423 39.584133,21.039852 39.585428,21.04237 39.579149))

Questions

• I suspect ogr2ogr is the issue because I'm using PostGIS 3 in both cases - is that reasonable?
• I've tried PG_USE_TEXT YES and PRECISION=NO flags in ogr2ogr to no effect
• Is there any chance that the higher precision coordinates are what's stored in the shapefile, even though ogrinfo with GDAL 3 doesn't report them?

I could certainly truncate / manage this change in PostGIS, but am trying to understand why this happened and if there's a way to get the previous behavior.

Answer 1

I am not sure if this is a GDAL issue. Perhaps it is not an issue at all but rounding errors just belong to conversions between decimal numbers and binary values that computers are using.

I stored your data into PostGIS with GDAL 3.3.0 and the geometry appears to be in the database with unaltered precision when queried this way:

select ST_AsText(wkb_geometry)
from mwe;

POLYGON((21.04237 39.579149,21.038807 39.577435,21.036219 39.576187,21.033495 39.574875,21.03046 39.573414,21.024775 39.570475,21.023056 39.569588,21.017654 39.563602,21.009989 39.563313,21.008549 39.563259,21.001377 39.565672,21.001286 39.565702,21.004169 39.575203,21.006231 39.581999,21.013037 39.583546,21.021423 39.584133,21.039852 39.585428,21.04237 39.579149))

Your DumpPoints query returns me values like 39.570474999999995 so the difference seems to be in ST_DumpPoints but the next test does not confirm that.

select ST_Y(geom)
from (
  select (ST_DumpPoints(
      ST_GeomFromText(
      ' POLYGON((21.04237 39.579149,
          21.038807 39.577435,21.036219 39.576187,
          21.033495 39.574875,21.03046 39.573414,
          21.024775 39.570475,21.023056 39.569588,
          21.017654 39.563602,21.009989 39.563313,
          21.008549 39.563259,21.001377 39.565672,
          21.001286 39.565702,21.004169 39.575203,
          21.006231 39.581999,21.013037 39.583546,
          21.021423 39.584133,21.039852 39.585428,
          21.04237 39.579149))'))).geom) pts
 
Result:

39.579149
39.577435
39.576187
39.574875
39.573414
39.570475
39.569588
39.563602
39.563313
39.563259
39.565672
39.565702
39.575203
39.581999
39.583546
39.584133
39.585428
39.579149

The coordinates are stored as doubles both into shapefiles and PostGIS so you are right that the decimal presentation is not what is used internally. You can test how the computers see the decimal numbers for example with the calculator at https://www.exploringbinary.com/floating-point-converter/.

A reliable way to test if you have identical geometries in both databases is to avoid functions and conversions into decimal numbers and compare the outputs from select wkb_geometry from mwe.

Answer 2

The cause: It seems to be due to a PostgreSQL change that occurred in version 12.0. Namely, the extra_float_digits parameter that is now being used for pg_dump and pg_dumpall functions (see here: https://www.postgresql.org/docs/release/12.0/ and here: https://www.postgresql.org/docs/12/runtime-config-client.html#GUC-EXTRA-FLOAT-DIGITS).

The solution: Set the extra_float_digits parameter to 0 in one of the following ways:

1. In postgresql.conf: uncomment and set the extra_float_digits parameter to 0
2. In psql (for all users): ALTER SYSTEM SET extra_float_digits TO 0;
3. In psql (for specific user): ALTER USER username SET extra_float_digits TO 0;
4. In psql (for specific db): ALTER DATABASE dbname SET extra_float_digits TO 0;

Don't forget to reload the service using sudo systemctl reload postgresql or SELECT pg_reload_conf(); depending on your method above.

I've tested and confirmed this on Ubuntu 20.04 using GDAL 3.0.4 and PostgreSQL 12.6.