gdalwarp fails when converting from Mercator spherical (recent break)

[jluis]

FWtools 2.2.0 and trunk (or 1.5.2) versions of gdalwarp are giving different results when converting from a Mercator spherical geotiff file. The FWtools version result is correct, which indicates that the error results from a recent modification

---

gdalwarp -t_srs "+proj=latlong +datum=WGS84" alg_merc.tiff alg_geo.tiff Creating output file that is 565P x 452L. <=== NOTE the 565P

Upper Left ( -8.0862815, 37.3460021) ( 8d 5'10.61"W, 37d20'45.61"N) Lower Left ( -8.0862815, 37.0644230) ( 8d 5'10.61"W, 37d 3'51.92"N) Upper Right ( -7.7343076, 37.3460021) ( 7d44'3.51"W, 37d20'45.61"N) Lower Right ( -7.7343076, 37.0644230) ( 7d44'3.51"W, 37d 3'51.92"N) Center ( -7.9102945, 37.2052125) ( 7d54'37.06"W, 37d12'18.77"N)

---

gdalwarp -t_srs "+proj=latlong +datum=WGS84" alg_merc.tiff alg_geo.tiff Creating output file that is 566P x 452L. <=== BUT NOW 566P

Upper Left ( -8.0862815, 37.1605907) ( 8d 5'10.61"W, 37d 9'38.13"N) Lower Left ( -8.0862815, 36.8792107) ( 8d 5'10.61"W, 36d52'45.16"N) Upper Right ( -7.7339340, 37.1605907) ( 7d44'2.16"W, 37d 9'38.13"N) Lower Right ( -7.7339340, 36.8792107) ( 7d44'2.16"W, 36d52'45.16"N) Center ( -7.9101077, 37.0199007) ( 7d54'36.39"W, 37d 1'11.64"N)

[rouault]

After quite a bit of investigation (trying on gdal trunk, 1.5 branche, 1.4.4 and 1.4.2), I always got the same result. The 565 x 452 case.

I also tried FWTools-2.0.6 on Linux, and I got the 566 x 452 result...

Hum, so I began to suspect the proj4 version to be the cause, and it is ! If I LD_PRELOAD the proj4 version included in FWTools-2.0.6 and run gdalwarp on any of the non FWtools version, I got the 566 x 452 case!

Here's the difference :

even@pc-amd64:~/gdal.svn.commit/gdal$  CPL_DEBUG= apps/gdaltransform -t_srs "+proj=latlong +datum=WGS84" -s_srs "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs"
OGRCT: Source: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs
OGRCT: Target: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
OGRCT: Source: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
OGRCT: Target: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs
-899154.716 4456528.555
-8.08628149867754 37.3460021172736 699.547525894787

even@pc-amd64:~/gdal.svn.commit/gdal$ LD_PRELOAD=../../FWTools-2.0.6/lib/libproj.so CPL_DEBUG= apps/gdaltransform -t_srs "+proj=latlong +datum=WGS84" -s_srs "+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs"
OGRCT: Source: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs
OGRCT: Target: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
OGRCT: Source: +proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs
OGRCT: Target: +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6371008.7714 +b=6371008.7714 +units=m +no_defs
-899154.716 4456528.555
-8.08628149867754 37.160590688729 0

(I got the proj4 definition for the source SRS by spying with CPL_DEBUG when doing gdalwarp.)

The version of proj4 installed on my system is proj 4.5.0. In the readme of FWTools-2.0.6, the version announced is 4.5.0 (CVS)... So I suspected that this means > 4.5.0. I then downloaded and compiled proj 4.6.0, and I get the same results as with FWTools-2.0.6

So the conclusion is that the version of proj installed on your system <= 4.5.0... and that I should consider upgrading to 4.6.0 I guess !

By looking at the Changelog, I suspect that the following entry is at stake :

2007-11-25  Frank Warmerdam  <warmerdam@pobox.com>

        * pj_transform.c: Do ellipsoid comparisons using the _orig ellipse
        values rather than the adjusted one.  Use these original values for
        any conversion to/from geocentric coordinates.

        Also, only do pj_datum_transform if neither the source nor destination
        is PJD_UNKNOWN.  This means we will no longer attempt via-geocentric
        adjustments for coordinate systems lacking a datum definition (having
        only an ellipsoid.

        * projects.h, pj_init.c: added a_orig and es_orig values in the PJ
        structure so we can distinguish between the originally requested
        ellipsoid, and the ellipsoid after adjustment for spherical projections

        Todays changes courtesy of bug 1602.

I suspect you hit paragraph 2. Your source has only an ellipsoid, not a datum.

Hum, maybe Frank can confirm my guesses ?

[jluis]

Even, Nice work. It was it. Building GDAL with external proj4 version > 4.5.0 (which I already had) gives the correct result. Curiously it's the second time I step in this issue (the other was ticket 2025 that originated ticket 1602 of proj4).

Now the question is: How do we update the GDAL internal version of PROJ4?

[rouault]

GDAL doesn't have an internal version of PROJ4. It just relies on the shared object libproj.so that it loads dynamically at runtime from the one it grabs in the path. So you have to update your default proj4 version installed on your system .

An alternative is to define the environment variable PROJSO and makes it point to the full filename of the libproj.so.

[jluis]

Well, I learned that one now. I got mislead by the fact that when keeping commented the PROJ4 section on nmake.opt (I'm on windows) the Dependency Walker doesn't show any dependency of gdal.dll on proj.dll. So I jumped into the conclusion that the thing was dealt internally.