[gdal-dev] coordinate transformation problem in python

[WolfgangZ]

Hi all,

I was trying to program a short script in python to do coordinate 
transformation from latlong/WGS84 to GaussKrüger Zone 4. An example how 
to do that is on the gdal web page. To keep it simple I used the EPSG 
codes. Here my code (ps I know that it should be "from osgeo import ogr" 
but I run my code on the recent version of FWTools 2.2.2):

import ogr

#original SRS
oSRS=ogr.osr.SpatialReference()
oSRS.ImportFromEPSG(4326)

#target SRS
tSRS=ogr.osr.SpatialReference()
tSRS.ImportFromEPSG(31468)

poCT=ogr.osr.CoordinateTransformation(oSRS,tSRS)

x=12.074760
y=47.893320
res=poCT.TransformPoint(x,y,0.)
print 'x= %10.2f , y= %10.2f' % (res[0],res[1])


The result is
x= 4505693.98 , y= 5306127.81

I checked against another software and got a different result:
4505697, 5306126

Then I tried it with cs2cs.exe ad the result from there is:
cs2cs.exe +proj=latlong +datum=WGS84 +to +proj=tmerc +lat_0=0 +lon_0=12 
k=1.000000 +x_0=4500000 +y_0=0 +ellps=bessel 
+towgs84=591.28,81.35,396.39,1.477,-0.0736,-1.458,9.82 +units=m +no_defs
12.074760 47.893320
4505696.75      5306125.78 -52.43

I suppose the my python code is ignoring the different ellipsoid. Is 
that behaviour desired? What about other gdal software like gdalwarp?

Thanks for any clarification.
Wolfgang