[Proj] Projecting raster image produces wrong output coordinates

[Frank Warmerdam]

Paul Hiemstra wrote:
> 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.

Paul,

I believe you are seeing PROJ attempt to correct for the change in
ellipsoid which is a big change between the sphere in your
stereo projection and the normal world ellipsoids.  Try appending
" +nadgrids=@null +wktext " to your radar projection when providing
it to gdalwarp.  The +nadgrids=@null directive tells PROJ.4 to
use a "no-op" shift between the datums instead of the default
change of ellipsoid logic.  The +wktext tells GDAL to preserve
the exist PROJ.4 string, otherwise special options like +nadgrids
will be lost in internal conversion to WKT and back.

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org