[Proj] DHDN / Gauss Krueger to WGS84

[Frank Warmerdam]

Oliver Eichler wrote:
>> I do have a coordinate in easting/northing [m] from a raster map  using
>> Gauss Krüger projection and Potsdam datum. Now I would like to get the
>> lon/lat coordinate on a WGS84 ellipsoid. And vice versa. I would like to
>> use the proj4 C API not the command line.
>>
>>   
> 
> Ok, further investigations condensed into something like:
> 
> pjgaukru = pj_init_plus("+init=epsg:31467");
> pjwgs84 = pj_init_plus("+init=epsg:4326");
> 
> // easting / northing to lon / lat
> pt = pj_inv(pt,pjgaukru)
> 
> //???
> double z = 0;
> 
> pj_datum_transform(pjgaukru, pjwgs84, 1, 0, &pt.u, &pt.v, &z);
> 
> However the result is still several minutes away from the expected
> coordinate. What do I miss?

Oliver,

How about:

pjgaukru = pj_init_plus("+init=epsg:31467");
pjwgs84 = pj_init_plus("+init=epsg:4326");

//???
double z = 0;

pj_datum_transform(pjgaukru, pjwgs84, 1, 0, &pt.u, &pt.v, &z);

Don't forget that the output is pt.u with the longitude in
radians, and pt.v with the latitude in radians.  The call
you make to pj_inv() is duplicating work that is done by
pj_transform() (corrupting really).

I would add that epsg:31467 expands as:

+proj=tmerc +lat_0=0 +lon_0=9 +k=1.000000 +x_0=3500000 +y_0=0 +ellps=bessel 
+datum=potsdam +units=m +no_defs

Note that this definition is using +datum=potsdam which expands to:

   +ellps=bessel +towgs84=606.0,23.0,413.0

This may or may not be the best datum shift parameters for your area of
work.  You may need to research a better local shift.

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