TransformTo precision error when transforming a box that crosses 180 degrees on 32-bit system

[cdestigter]

The following python code reproduces the problem:

from osgeo import osr, ogr

# a polygon that crosses 180 degrees
ogr_geom = ogr.CreateGeometryFromWkt('POLYGON ((165.9134918212895400 -47.7347488403320312, 165.9134918212895400 -34.2310752868652344, 184.6205993652343977 -34.2310752868652344, 184.6205993652343977 -47.7347488403320312, 165.9134918212895400 -47.7347488403320312))')

# from wgs84
src = osr.SpatialReference()
src.ImportFromProj4('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs ')
ogr_geom.AssignSpatialReference(src)

# to google mercator
dest = osr.SpatialReference()
dest.ImportFromProj4('+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs')

ogr_geom.TransformTo(dest)

print ogr_geom.ExportToWkt()

The output is:

POLYGON ((18469405.425279911607504 -6062839.394079929217696,18469405.425279911607504 -4059872.155139415990561,-19523145.574291624128819 -4059872.155139415990561,-19523145.574291624128819 -6062839.394079929217696,18469405.425279915332794 -6062839.394079929217696))

Note that the X coordinates of the first and last points are slightly different. Obviously this causes problems since this is no longer a valid polygon.

This only happens on 32-bit systems. Also, I have tried converting the same sequence of points with pyproj and that works okay, so seems like Proj is not to blame (I could be mistaken about this?)

[Even Rouault]

I don't manage to reproduce your issue with any of GDAL 1.7.0dev, 1.6.1RC2, 1.6.0 on Ubuntu 32bit (with proj 4.6.0). I get the following output each time:

POLYGON ((
18469405.425279915332794 -6062839.394079929217696,
18469405.425279915332794 -4059872.155139415990561,
-19523145.574291624128819 -4059872.155139415990561,
-19523145.574291624128819 -6062839.394079929217696,
18469405.425279915332794 -6062839.394079929217696))

What you get is really surprising. What is your OS and CPU architecture ? Where did you get the binaries for GDAL ? What is your proj.4 version ?

[cdestigter]

Ubuntu Intrepid, i386

proj 4.6.1-4 and gdal 1.6.0-3, both from debian and rebuilt under ubuntu intrepid

The rebuilding happened in a pbuilder (clean) environment

[Even Rouault]

I've tested by running against proj 4.6.1 and I still get correct results.

I, myself, have package gdal 1.6.0 for intrepid (from my hardy, in a intrepid pbuilder environmenet, but I couldn't test). However the same package for hardy doesn't show any issue with your test. My packages for intrepid are here : ​http://even.rouault.free.fr/deb-intrepid-i386-gdal-1.6.0/. You might want to try and test if you can reproduce with it.

A possible way to debug would be to run valgrind --trace-children=yes ./your_script.py and see if anything interesting shows (but I must warn you that there are many "normal" warnings coming from python. So you would have to focus on warnings coming from libgdal itself). You would probably want to compile with debug and without optimizations.

I believe it might be due to your packaging, but I'm not sure how you could have make anything wrong that would result in such a subtle problem.

[Robert Coup]

Bit more detail after some debugging:

• Happens on both 32bit & 64bit systems - the new Debian proj packaging had the proj-data package (gridshifts) as being architecture independent when they're not (fixed now) which caused differences.
• It's a Proj problem: ​http://trac.osgeo.org/proj/ticket/45
• It didn't show up in cs2cs tests because Proj uses a different degrees-to-radians constant which avoids it by coincidence
• The @null gridshift doesn't work on longitudes outside -180º/+180º longitude
• The @null transform isn't quite null, and it bails out part way through because of the above, leaving some points gridshifted and others not (by a very small amount).
• This is why the start-point ends up different from the end-point