[Proj] Projecting raster image produces wrong output coordinates

[Paul Hiemstra]

Dear proj-list,

I am trying to convert a radar rainfall image from one projection to 
another. The projection of the input file is:

"+proj=stere +x_0=0 +y_0=0 +lat_0=90 +lon_0=0 +lat_ts=60 +a=6378.388 
+b=6356.906"

and I am trying to project to "EPSG:28992" (Dutch CRS). But when I use 
gdal to perform the projection the output raster file is no longer on 
the correct location but has shifted thousands of kilometers. I used 
both Linux (Debian Etch, proj 4.4.9, GDAL 1.4.0.0) and Windows (Proj 
4.5.0, GDAL 1.5dev, FWTools 2.0.0) and ran into the same problems. The 
problems also occurs when the target projection is not "EPSG:28992", 
e.g. "EPSG:4326" or "EPSG:3035".

Transforming from "EPSG:28992" to the projection of the input file 
yields good results.

kind regards,
Paul Hiemstra

Gdal output:
 > gdalinfo test.tif
Driver: GTiff/GeoTIFF
Files: test.tif
       test.aux
Size is 256, 256
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378.388,296.9177916395094]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Polar_Stereographic"],
    PARAMETER["latitude_of_origin",60],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (0.000000000000000,-3727.291000000000200)
Pixel Size = (2.500000000000000,-2.500000000000000)
Metadata:
  AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (       0.000,   -3727.291) (  0d 0'0.01"E, 55d28'14.64"N)
Lower Left  (       0.000,   -4367.291) (  0d 0'0.01"E, 49d58'1.45"N)
Upper Right (     640.000,   -3727.291) (  9d44'34.97"E, 54d59'42.06"N)
Lower Right (     640.000,   -4367.291) (  8d20'13.23"E, 49d34'22.17"N)
Center      (     320.000,   -4047.291) (  4d31'14.52"E, 52d35'19.33"N)
Band 1 Block=256x8 Type=Float32, ColorInterp=Gray
 >gdalwarp -t_srs "+init=EPSG:28992" test.tif output.tif
Creating output file that is 304P x 259L.
Processing input file test_1800_stere.tif.
:0...10...20...30...40...50...60...70...80...90...100 - done.
 >gdalinfo output.tif
Driver: GTiff/GeoTIFF
Files: output.tif
Size is 304, 259
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["Bessel 1841",
        DATUM["unknown",
            SPHEROID["unnamed",6377397.155,299.1528128000009]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Oblique_Stereographic"],
    PARAMETER["latitude_of_origin",52.15616055555555],
    PARAMETER["central_meridian",5.38763888888889],
    PARAMETER["scale_factor",0.999908],
    PARAMETER["false_easting",155000],
    PARAMETER["false_northing",463000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (98602.580578793495000,4311602.125703314300000)
Pixel Size = (318.397229988391810,-318.397229988391810)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (   98602.581, 4311602.126) (  0d44'48.56"W, 85d39'38.22"N)
Lower Left  (   98602.581, 4229137.243) (  0d 3'42.14"E, 84d59'6.70"N)
Upper Right (  195395.338, 4311602.126) (  9d47'21.85"E, 85d40'20.58"N)
Lower Right (  195395.338, 4229137.243) (  9d12'27.11"E, 84d59'43.47"N)
Center      (  146998.960, 4270369.684) (  4d34'34.58"E, 85d20'43.24"N)
Band 1 Block=304x6 Type=Float32, ColorInterp=Gray

-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone: 	+31302535773
Fax:	+31302531145
http://intamap.geo.uu.nl/~paul