[gdal-dev] Mismatched output from listgeo/gdalwarp/cs2cs

[Pepijn Van Eeckhoudt]

Hi,

My issue is related to gdal, libgeotiff and proj.4, but the problem 
appeared with gdalwarp so I'm posting here. I'm seeing different output 
from the various tools that I can't really make sense of and I was 
hoping someone here could clarify the meaning of the various output values.

I received a geotiff file from one of my clients that essentially has 
all it's referencing metadata set to 'User-Defined'. Listgeo -proj4 
prints out:

+proj=tmerc +lat_0=0.000000000 +lon_0=15.000000000 +k=0.999900 +x_0=500000.000 +y_0=-5000000.000 +a=6377397.155 +b=6356078.963 +units=m

I've double checked the geotiff data manually and this is indeed the 
information that is in the file. Transverse mercator projection with 
respect to the Bessel 1841 ellipsoid in other words. My client 
reprojected the geotiff to WGS84 using the following gdalwarp input

./gdalwarp.exe -t_srs '+proj=tmerc +lat_0=0 +lon_0=15 +datum=WGS84'

He then loaded both geotiffs in our GIS tool (not based on osgeo 
libraries) and noticed a substantial horizontal shift between the two 
variants of the file. I'm now trying to figure out what is going wrong. 
I listed what I've been able to find out so far below.

Listgeo prints out the following information for the upper left corner 
point:

Upper Left    (  439702.628,   71450.647)  ( 14d13'27.94"E, 45d47' 5.72"N)

If I then take the same projected coordinate and input those into cs2cs 
I get different lon lat values.

cs2cs +proj=tmerc +lat_0=0.000000000 +lon_0=15.000000000 +k=0.999900 +x_0=500000.000 +y_0=-5000000.000 +a=6377397.155 +b=6356078.963 +units=m
14d13'59.5"E    45d8'31.057"N 0.000, 71450.647

Shouldn't I be getting the same values here? If not, then what are the 
two tools actually showing me; it's not clear to me with respect to 
which datum the lon/lat coordinates are defined?

I then tried to get cs2cs to produce WGS84 lat/lon coordinates, but 
oddly enough that produced the same output as the previous command.

cs2cs +proj=tmerc +lat_0=0.000000000 +lon_0=15.000000000 +k=0.999900 +x_0=500000.000 +y_0=-5000000.000 +a=6377397.155 +b=6356078.963 +units=m +to +proj=latlon +datum=WGS84 coords.txt
14d13'59.5"E    45d8'31.057"N 0.000, 71450.647

I was able to get different output by adding an explicit +towgs84 
parameter for the input CRS

cs2cs +proj=tmerc +lat_0=0.000000000 +lon_0=15.000000000 +k=0.999900 +x_0=500000.000 +y_0=-5000000.000 +a=6377397.155 +b=6356078.963 +units=m +towgs84=0,0,0 +to +proj=latlon +datum=WGS84 coords.txt
14d13'59.5"E    45d8'33.211"N -706.488, 71450.647

This leads me to conclude that without an +datum reference or a +towgs84 
parameter value proj (and possibly gdalwarp) no longer performs any 
datum transformation whatsoever. In other words it is not converting 
from one ellipsoid to the other via geocentric. Is that correct? If so, 
what is the rationale behind this choice? Our GIS toolkit takes a 
different approach and defaults to +towgs84=0,0,0 instead. My hunch is 
that this difference in default interpretation is the cause of the 
horizontal shift I'm seeing.

Thanks in advance for any help,

Pepijn Van Eeckhoudt
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20120109/d1ceef2e/attachment.html