Projection problems - Mixing maps with different geoids

[Trond Michelsen]

Hi.

I'm having a bit of a problem. 

I'm using mapserver to draw shapefiles on top of output from a
WMS-server. The maps from the remote WMS server uses an unstandard
projection with these parameters:

 +proj=stere
 +lat_0=90
 +lon_0=0
 +lat_ts=60
 +a=6371000
 +b=6371000
 +units=m
 +no_defs

If I set up my WMS server to use this projection as well, and define
all my shapefile layers to use their projection (which is UTM33), the
data from the shapefiles is drawn quite a bit too far south (well, too
close to equator, really). This seems to be because of the difference
in geoids in the two projections.

The cs2cs program demonstrates the offset, here for 60N, 10E:

echo "10 60" \
  | proj  +init=epsg:32633 \
  | cs2cs +init=epsg:32633 +to +proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6371000 +b=6371000 +units=m +no_defs \
  | invproj                    +proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6371000 +b=6371000 +units=m +no_defs

10dE    59d49'59.074"N -8867.78

If I use invproj and proj, instead of cs2cs, everything works as expected:

echo "10 60" \
 | proj    +init=epsg:32633 \
 | invproj +init=epsg:32633 \
 | proj    +proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6371000 +b=6371000 +units=m +no_defs \
 | invproj +proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6371000 +b=6371000 +units=m +no_defs

10dE    60dN

I know this last example seems silly, but the point is that doing an
inverse and forward projection on the coordinates gives a correct
result, while doing cs2cs does not.

Interestingly, the program shpproj, which is included with shapelib,
is apparently doing an inverse and forward projection. So, converting
my shapefiles with the command:

shpproj infile.shp outfile.shp \
 -in=epsg:32633 \
 -out="+proj=stere +lat_0=90 +lon_0=0 +lat_ts=60 +a=6371000 +b=6371000 +units=m +no_defs"

gives me a shapefile I can use.

Converting the shapefiles is a workaround I can live with for a while,
but I was wondering if it's possible to set up mapserver so that I can
use the shapefiles directly. Not too important right now, but once I
start fetching data from remote WFS-servers, it'll probably be a
bigger issue.

-- 
Trond Michelsen