Getting UTM and stereographic projections to overlap

[Trond Michelsen]

Hi.

I'm having a bit of trouble displaying a shapefile in UTM33-format
properly together with my maps that use polar stereographic
projection.

I've made a small perl-program (included below) that creates two
shapefiles, one with the same projection parameters as my map, and one
in UTM33. Both files simply contains one point for every fifth
latitude and longitude.

I also have a map of Norway as a shapefile in UTM33.

When I plot my basemap together with the two point layers, and the
shapefile of Norway, the stereographic point layer matches the
latitudes and logitudes of the basemap. The UTM33 point layer has an
offset that seems to grow the further south you get. The shapefile of
Norway seems to have pretty much the same offset as the UTM33 point layer

Here are some examples:

http://crusaders.no/~trondmm/mapserv/utm_vs_stere.png
http://crusaders.no/~trondmm/mapserv/utm_vs_stere2.png
http://crusaders.no/~trondmm/mapserv/utm_vs_stere3.png

  red dots    = stereographic point layer
  yellow dots = UTM33 point layer


Now, after some experimenting, I've found out that I can get the two
point layers to overlap completely, if I change the stereographic
shapefile to use the WGS84 ellipsoid. Unfortunately, my basemaps
depend on a spherical earth, and it's not possible for me to convert
them.

So - is there anything I can do in the mapfile to get mapserver to get
the layers to overlap more precisely?


Here are my projection parameters from the mapfile:

The Web object:

PROJECTION
  "proj=stere"
  "lat_0=90"
  "lon_0=0"
  "lat_ts=60"
  "a=6371000"
  "b=6371000"
  "units=m"	
END

The stereographic point layer and all basemaps:

PROJECTION
  "proj=stere"
  "lat_0=90"
  "lon_0=0"
  "lat_ts=60"
  "a=6371000"
  "b=6371000"
  "units=m"	
END

The UTM33-point layer and shapefile of Norway:

PROJECTION
  "proj=utm"
  "zone=33"
  "ellps=WGS84"
  "datum=WGS84"
  "units=m"
END

and here's the program I used to generate the point layers:

--8<--
#!/usr/bin/perl

use strict;
use warnings;
use Geo::Proj4;
use Geo::Shapelib ':constants';

my %proj = (stere => Geo::Proj4->new(proj   => "stere",
                                     lat_0  => 90,
                                     lon_0  => 0,
                                     lat_ts => 60,
                                     a      => 6371000,
                                     b      => 6371000,
                                     units  => "m",
                                     ),
	    utm33 => Geo::Proj4->new(proj   => "utm",
                                     zone   => 33,
                                     ellps  => "WGS84",
                                     datum  => "WGS84",
                                     units  => "m",
                                     ),
	    );
	    
while (my ($name, $proj) = each %proj) {
    my $shape = Geo::Shapelib->new();

    $shape->{Shapetype} = POINT;

    $shape->{FieldNames} = [qw/POS       LAT    LON   /];
    $shape->{FieldTypes} = [qw/String:64 Double Double/];
    
    for (my $lat = 0;$lat <= 90; $lat += 5) {
        for (my $lon = -30; $lon <= 30; $lon += 5) {
            my ($x, $y) = $proj->forward($lat, $lon);
	    
            push @{$shape->{Shapes}}, {SHPType  => POINT,
                                       Vertices => [[$x,$y]]
                                      };
            push @{$shape->{ShapeRecords}}, ["$lat $lon",
                                             $lat,
                                             $lon,
                                            ];
        }
    }

    $shape->save("shapetest_$name");
}
--8<--

-- 
Trond Michelsen