Projection roundtip losing too much precision for NAD27 at -128.0 longitude

[joaoportela]

I'm trying to test roundtrip (roundtrip is p1->p2->p1) stability in my projections and found a curious case (a bug?) when trying to convert from WGS84 to NAD27 and back.

We can see the values here:

>>> from pyproj import Proj, transform
>>> wgs84_proj = Proj('+proj=longlat +datum=WGS84 +no_defs')
>>> nad27_proj = Proj('+proj=longlat +datum=NAD27 +no_defs')
>>> wgs84_pos = (-128.0, 50.98)
>>> nad27_pos = transform(wgs84_proj, nad27_proj, wgs84_pos[0], wgs84_pos[1])
>>> nad27_pos
(-127.99845565740134, 50.98025249402805)
>>> wgs84_pos_roundtrip = transform(nad27_proj, wgs84_proj, nad27_pos[0], nad27_pos[1])
>>> wgs84_pos_roundtrip
(-128.00005191895045, 50.98006355013494)

And they are almost 8 meters apart:

>>> from pyproj import Geod
>>> wgs84_geod = Geod(ellps='WGS84')
>>> angle1, angle2, distance = wgs84_geod.inv(wgs84_pos[0], wgs84_pos[1], wgs84_pos_roundtrip[0], wgs84_pos_roundtrip[1])
>>> distance
7.954669999055905

I have tested several positions and I only seem to have problems with longitude at -128.0 and latitudes between 50.98 and 76.98.

For example:

>>> wgs84_pos = (-127.9, 50.98)
>>> nad27_pos = transform(wgs84_proj, nad27_proj, wgs84_pos[0], wgs84_pos[1])
>>> wgs84_pos_roundtrip = transform(nad27_proj, wgs84_proj, nad27_pos[0], nad27_pos[1])
>>> nad27_pos
(-127.89840687790532, 50.98018566035921)
>>> wgs84_pos_roundtrip
(-127.9, 50.98000000000001)
>>> angle1, angle2, distance = wgs84_geod.inv(wgs84_pos[0], wgs84_pos[1], wgs84_pos_roundtrip[0], wgs84_pos_roundtrip[1])
>>> distance
7.069288649822944e-10

Which seems to be just fine!

I also tested this with ArcGis API and it seems to work fine.

[micahcochran]

I tried to replicate this on the command line with cs2cs:

$ cs2cs -f "%.14f" +proj=latlong +datum=WGS84 +to +proj=latlong +datum=NAD27
-127.9 50.98
-127.89840687790534 50.98018566035922 0.00000000000000
$ cs2cs -f "%.14f" +proj=latlong +datum=NAD27 +to +proj=latlong +datum=WGS84
-127.89840687790534 50.98018566035922
-127.90000000000006 50.98000000000003 0.00000000000000

The round trip result was -127.90000000000006 50.98000000000003. Those right-most digits are simply the limits of 64-bit double math.

I used pyproj.Geod for the distance calculation and got: 6.0880917533296244e-09 meters. Close enough to your last result.

This bug could be problematic.

[joaoportela]

That is (kind of) good news.

What else can I do to figure out what is causing this?

[micahcochran]

Sorry, I didn't replicated your example fully. Here is the version with the error. Using -128.0 longitude.

C:\...>cs2cs -f "%.14f" +proj=latlong +datum=WGS84 +to +proj=latlong +datum=NAD27
-128.0  50.98
-127.99845565740137     50.98025249402806 0.00000000000000
C:\...>cs2cs -I -f "%.14f" +proj=latlong +datum=WGS84 +to +proj=latlong +datum=NAD27
# Inverse of last
-127.99845565740137     50.98025249402806
-128.00005191895048     50.98006355013496 0.00000000000000

And the distance calculation:

>>> from pyproj import Geod
>>> wgs84_geod = Geod(ellps='WGS84')
>>> angle1, angle2, distance = wgs84_geod.inv(-128.0, 50.98, -128.00005191895048, 50.98006355013496)
>>> distance
7.954670001528683

That is similar to your results, about 8 meters off.

Base on that I don't think pyproj is at fault, it is just using what the underlying PROJ.4 C library is giving it. Please report this issue to proj.4.

[joaoportela]

Sorry.
I didn't notice that you were using -127.9 in your tests.

From my local tests I already noticed that some values around -128.0 don't trigger the issue.
Examples:

• -127.9 doesn't trigger the issue.
• -127.99 doesn't trigger the issue.
• -127.999 doesn't trigger the issue.
• -127.999999 doesn't trigger the issue.
• -127.99999999999 doesn't trigger the issue.
• -127.999999999999 triggers the issue.
• -128.0 triggers the issue.
• -128.001 triggers the issue.
• -128.01 doesn't trigger the issue.
• -128.1 doesn't trigger the issue.

[joaoportela]

This seems to be an issue of OSGeo/proj.4. Reported here.

Closing.