[Gdal-dev] Lambert to Swiss Obl projection error?

[Jed O. Kaplan]

Hello all,

I am using the latest CVS versions of GDAL and PROJ4 on Mac OS X 10.4

I am trying to warp a small 100x100 raster in a Lambert Azimuthal  
projection to the Swiss Obl. Mercator projection. gdal_warp runs  
successfully, but the output raster is wrong - the grid seems shifted  
to the north by ca. 20km (a simple shift of the grid does not solve  
the problem, however).  You can see the problem if you compare the  
calculated center lon and lat of the input and output rasters - these  
center points should be approximately the same, but are shifted ca.  
00:12' to the north in the output grid.

The relevant output from gdalinfo and the gdal_warp command line, as  
well as the WKT file I use to define the Swiss projection are copied  
below.

Could this be a bug in GDAL, or possibly in PROJ4? Thank you for any  
advice you may be able to give me!

Cheers,

Jed Kaplan

--

jkaplan at kootenay:~/Documents/science_projects/vivaldi/BioClim_SR >  
gdalinfo terrain.tif
Driver: GTiff/GeoTIFF
Size is 100, 100
Coordinate System is:
PROJCS["Lambert Azimuthal Equal Area",
     GEOGCS["unnamed",
         DATUM["unknown",
             SPHEROID["unnamed",6370997,0]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Lambert_Azimuthal_Equal_Area"],
     PARAMETER["latitude_of_center",46.8],
     PARAMETER["longitude_of_center",8.2],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]
Origin = (114914.109000,-17980.043000)
Pixel Size = (250.00000000,-250.00000000)
Metadata:
   AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  114914.109,  -17980.043) (  9d42'18.00"E, 46d37'42.29"N)
Lower Left  (  114914.109,  -42980.043) (  9d41'55.60"E, 46d24'13.08"N)
Upper Right (  139914.109,  -17980.043) ( 10d 1'56.36"E, 46d37'25.12"N)
Lower Right (  139914.109,  -42980.043) ( 10d 1'29.10"E, 46d23'55.99"N)
Center      (  127414.109,  -30480.043) (  9d51'54.75"E, 46d30'49.54"N)
Band 1 Block=100x20 Type=Float32, ColorInterp=Gray

jkaplan at kootenay:~/Documents/science_projects/vivaldi/BioClim_SR >  
gdalwarp -t_srs swiss.wkt -rb terrain.tif terrain-sw.tif
Creating output file that is 101P x 101L.
:0...10...20...30...40...50...60...70...80...90...100 - done.

jkaplan at kootenay:~/Documents/science_projects/vivaldi/BioClim_SR >  
gdalinfo terrain-sw.tif
Driver: GTiff/GeoTIFF
Size is 101, 101
Coordinate System is:
PROJCS["Swiss. Obl. Mercator",
     GEOGCS["bessel",
         DATUM["CH1903",
             SPHEROID["Bessel 1841",6377397.155,299.1528128000008,
                 AUTHORITY["EPSG","7004"]],
             AUTHORITY["EPSG","6149"]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Hotine_Oblique_Mercator"],
     PARAMETER["latitude_of_center",46.95240555555556],
     PARAMETER["longitude_of_center",7.439583333333333],
     PARAMETER["azimuth",90],
     PARAMETER["rectified_grid_angle",90],
     PARAMETER["scale_factor",1],
     PARAMETER["false_easting",600000],
     PARAMETER["false_northing",200000],
     UNIT["metre",1,
         AUTHORITY["EPSG","9001"]]]
Origin = (772925.084811,188209.083260)
Pixel Size = (249.84685851,-249.84685851)
Metadata:
   AREA_OR_POINT=Area
Corner Coordinates:
Upper Left  (  772925.085,  188209.083) (  9d42'22.99"E, 46d49'25.83"N)
Lower Left  (  772925.085,  162974.551) (  9d41'48.72"E, 46d35'48.88"N)
Upper Right (  798159.618,  188209.083) ( 10d 2'13.02"E, 46d49'0.48"N)
Lower Right (  798159.618,  162974.551) ( 10d 1'33.75"E, 46d35'23.64"N)
Center      (  785542.351,  175591.817) (  9d51'59.62"E, 46d42'25.14"N)
Band 1 Block=101x20 Type=Float32, ColorInterp=Gray

jkaplan at kootenay:~/Documents/science_projects/vivaldi/BioClim_SR >  
more swiss.wkt
PROJCS["Swiss. Obl. Mercator",
     GEOGCS["bessel",
         DATUM["CH1903",
             SPHEROID["bessel",6377397.155,299.1528128],
             TOWGS84[674.374,15.056,405.346,0,0,0,0]],
         PRIMEM["Greenwich",0],
         UNIT["degree",0.0174532925199433]],
     PROJECTION["Hotine_Oblique_Mercator"],
     PARAMETER["latitude_of_center",46.95240555555556],
     PARAMETER["longitude_of_center",7.439583333333333],
     PARAMETER["azimuth",90],
     PARAMETER["rectified_grid_angle",90],
     PARAMETER["scale_factor",1],
     PARAMETER["false_easting",600000],
     PARAMETER["false_northing",200000],
     UNIT["metre",1]]