[Gdal-dev] Different projection results between gdalwarp and gdal python bindings

Tue Aug 21 06:36:07 EDT 2007

Hi, 
I am finding some inconsistencies between the gdalwarp utility and the gdal Python bindings (old generation). The aim is to project a raster from UTM30N/ED50 to WGS84 (the projection WKTs are at the end of this message). Using gdalwarp (GDAL 1.4.2.0, released 2007/06/27) works fine, and I get reasonable results. If I then use gdal.CreateAndReprojectImage(), the image has a different GeoTransform. Here are the extents reported by gdalinfo:

* ORIGINAL IMAGE (UTM 30N/ED50)
Pixel Size = (1000.000000000000000,-1000.000000000000000)
Upper Left  (   77000.000, 4309000.000) (  7d52'18.12"W, 38d49'39.10"N)
Lower Left  (   77000.000, 3937000.000) (  7d39'40.95"W, 35d29'9.44"N)
Upper Right (  644000.000, 4309000.000) (  1d20'20.68"W, 38d55'2.82"N)
Lower Right (  644000.000, 3937000.000) (  1d24'39.85"W, 35d33'56.38"N)
Center      (  360500.000, 4123000.000) (  4d34'21.80"W, 37d14'32.98"N)

* GDALWARP OUTPUT IMAGE (WGS84) (this one if fine)
Pixel Size = (0.010672374453967,-0.010672374453967)
Upper Left  (  -7.8716995,  38.9284612) (  7d52'18.12"W, 38d55'42.46"N)
Lower Left  (  -7.8716995,  35.4812842) (  7d52'18.12"W, 35d28'52.62"N)
Upper Right (  -1.3402064,  38.9284612) (  1d20'24.74"W, 38d55'42.46"N)
Lower Right (  -1.3402064,  35.4812842) (  1d20'24.74"W, 35d28'52.62"N)
Center      (  -4.6059529,  37.2048727) (  4d36'21.43"W, 37d12'17.54"N)

* PYTHON OUTPUT IMAGE (WGS84) (this one is shifted/skewed)
Pixel Size = (0.010672660892937,-0.010672660892937)
Upper Left  (  -7.8730950,  38.9280489) (  7d52'23.14"W, 38d55'40.98"N)
Lower Left  (  -7.8730950,  35.4807795) (  7d52'23.14"W, 35d28'50.81"N)
Upper Right (  -1.3414265,  38.9280489) (  1d20'29.14"W, 38d55'40.98"N)
Lower Right (  -1.3414265,  35.4807795) (  1d20'29.14"W, 35d28'50.81"N)
Center      (  -4.6072608,  37.2044142) (  4d36'26.14"W, 37d12'15.89"N)

As you can see, there's a slight displacement. I am using exactly the same WKT string for the gdalwarp and the Python conversions, with explicity +towgs84 parameters. The gdalwarp command was simply:
gdalwarp -s_srs utm30N_ED50_spain.wkt -t_srs wgs84.wkt input.tif output.tif

The Python code looks like this:
g = gdal.OpenShared ("input.tif")
error = gdal.CreateAndReprojectImage ( g, "output.tif", \
		src_wkt=utm30,	dst_wkt= wgs84 )

Any clues?
Many thanks!
Jose

The offending projections:
'''PROJCS["ED50 / UTM zone 30N",GEOGCS["ED50",DATUM["European_Datum_1950",SPHEROID
["International 1924",6378388,297,AUTHORITY["EPSG","7022"]],AUTHORITY["EPSG","62
30"],TOWGS84[-131,-100.3,-163.4,-1.244,-0.020,-1.144,9.39]],PRIMEM["Greenwich",0
,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9
122"]],AUTHORITY["EPSG","4230"]],PROJECTION["Transverse_Mercator"],PARAMETER["la
titude_of_origin",0],PARAMETER["central_meridian",-3],PARAMETER["scale_factor",0
.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["met
re",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","23030"]]'''

'''GEOGCS["WGS 84",
DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.2572235629972,AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]],    PRIMEM["Greenwich",0],    UNIT["degree",0.0174532925199433],    AUTHORITY["EPSG","4326"]]'''

-- 
GMX FreeMail: 1 GB Postfach, 5 E-Mail-Adressen, 10 Free SMS.
Alle Infos und kostenlose Anmeldung: http://www.gmx.net/de/go/freemail