[Gdal-dev] OGRSpatialReference: ERROR 1: geocentric transformation missing z or ellps

[Reiner Beck]

Hello,

I try to transform WGS84 Lat/Long coordinates to the swiss coordinate 
system (Swiss Oblique Cylindrical) (see 
www.remotesensing.org/geotiff/proj_list/swiss_oblique_cylindrical.html.)

Using the OGRCoordinateTransformation Object produces the following 
error: ERROR 1: geocentric transformation missing z or ellps.

- the GDAL_DATA environment is set (see below)
- the PROJ_LIB environment is set (see below) (I tried with and without 
this env variable, is it needed?)

- I use GDAL 1.2.3 downloaded last friday (October, 22 2004)
- I use proj4 4.4.8 downloaded today (October, 24 2004) (not sure if 
GDAL uses an intenal copy)
- I tried to link proj.4 statically into GDAL 
(--with-static-proj4=/opt/m4com/util/lib/) and without this option
- Everything is installed into the directory "/opt/m4com/util"

The funny thing is, I upgraded Proj.4 because the second (LV95) 
transformation didn't work (the LV03 did, before upgrading to the latest 
versions ...).
With the new versions, unsetting the GDAL_DATA environment variable will 
produce the correct? output for the LV03 transformation (with the GEOGCS 
part missing in the WKT), while the LV95 still doesn't work. (Result in 
this case: Swiss Geo: LV03: 7.4189 46.6358 0 -> 598416 164808 0ERROR 1: 
latitude or longitude exceeded limits)

Any ideas?
Thanks!

Reiner

Here's the code:
--- code ---
...
OGRCoordinateTransformation* toSwissLV03 = NULL;
OGRCoordinateTransformation* toSwissLV95 = NULL;
OGRSpatialReference wgs84SR;
char* sw;
char* wg;

      wgs84SR.SetWellKnownGeogCS("WGS84");
      wgs84SR.exportToWkt(&wg);
      cerr << "OGRSpatialReference WGS84: " << wg << endl;

      OGRSpatialReference swissSR;
      swissSR.SetProjCS("Schweizer Landeskoordinaten LV03");
      swissSR.SetWellKnownGeogCS("EPSG:4149");
      swissSR.SetHOM(46.95240555555556, 7.439583333333333, 90.0, 90.0, 
1.0, 600000, 200000);

      swissSR.exportToWkt(&sw);
      cerr << "OGRSpatialReference LV03: " << sw << endl;
      toSwissLV03 = OGRCreateCoordinateTransformation(&wgs84SR, &swissSR);

      swissSR.SetProjCS("Schweizer Landeskoordinaten LV95");
      swissSR.SetWellKnownGeogCS("EPSG:4150");
      swissSR.SetHOM(46.95240555555556, 7.439583333333333, 90.0, 90.0, 
1.0, 2600000, 1200000);

      swissSR.exportToWkt(&sw);
      cerr << "OGRSpatialReference LV95: " << sw << endl;
      toSwissLV95 = OGRCreateCoordinateTransformation(&wgs84SR, &swissSR);
  
   VisCoordGeo3d swissGeo = geoOut;
   cerr << "Swiss Geo: LV03: " << swissGeo;
   if ((toSwissLV03) && (toSwissLV03->Transform(1, &swissGeo.geoLong, 
&swissGeo.geoLat, &swissGeo.height) == TRUE))
      cerr << " -> " << swissGeo;
   else
      cerr << "FAIL";
   if ((toSwissLV95) && (toSwissLV95->Transform(1, &swissGeo.geoLong, 
&swissGeo.geoLat, &swissGeo.height) == TRUE))
      cerr << " LV95 -> " << swissGeo;
   else
      cerr << "FAIL";
   cerr << endl;
...
--- code ---
(Note: the structure VisCoordGeo3d contains 3 doubles geoLong, geoLat, 
and height)

The output for the OGRSpatialReference's is:
OGRSpatialReference WGS84:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],
        TOWGS84[0,0,0,0,0,0,0],AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
        
UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9108"]],AXIS["Lat",NORTH],AXIS["Long",EAST],AUTHORITY["EPSG","4326"]]

OGRSpatialReference LV03:
PROJCS["Schweizer Landeskoordinaten LV03",
    GEOGCS["CH1903",
        DATUM["CH1903",
            SPHEROID["Bessel 
1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],AUTHORITY["EPSG","6149"]],
        PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
        
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4149"]],
    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]]

OGRSpatialReference LV95:
PROJCS["Schweizer Landeskoordinaten LV95",
    GEOGCS["CH1903+",
        DATUM["CH1903",
            SPHEROID["Bessel 
1841",6377397.155,299.1528128,AUTHORITY["EPSG","7004"]],
            
TOWGS84[674.374,15.056,405.346,0,0,0,0],AUTHORITY["EPSG","6150"]],
        PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],
        
UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4150"]],
    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",2600000],
        PARAMETER["false_northing",1200000]]

Using the OGRCoordinateTransformation object produces the following error:
Swiss Geo: LV03: 10.5249 43.9997 0
ERROR 1: geocentric transformation missing z or ellps
FAIL
ERROR 1: geocentric transformation missing z or ellps
FAIL

An "ls $GDAL_DATA" produces the follwoing output:
cubewerx_extra.wkt  epsg.wkt        gdal_datum.csv  prime_meridian.csv  
s57expectedinput.csv  seed_3d.dgn
ecw_cs.dat          esri_extra.wkt  gdalicon.png    projop_wparm.csv    
s57objectclasses.csv  stateplane.csv
ellipsoid.csv       gcs.csv         pcs.csv         s57attributes.csv   
seed_2d.dgn           unit_of_measure.csv

A "ls $PROJ_LIB" produces the follwoing output:
epsg  esri  GL27  nad27  nad83  nad.lst  proj_def.dat  world