[Proj] DHDN / Gauss Krueger to WGS84

[Roger Oberholtzer]

A while back, the following discussion took place here:


> On Thu, 2007-09-20 at 12:06 -0400, Frank Warmerdam wrote:
> 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.


I find myself needing the exact same thing. So, I thought I would
implement this in the C API, where I usually do these sorts of things. I
usually use pj_transform() instead of pj_datum_transform(). But as the
recipe uses pj_datum_transform(), I thought I would too.

The problem is that the values I get back from pj_datum_transform() are
infinite values (as expressed by printf). If I use pj_transform() the
locations are incorrect - they are located in central Asia.

I am testing with a sample set of coordinates provided later in the
quoted thread:

> it seems to work. A reference point of:
> 
> x = 4459750 m  y = 5448182 m
> 
> is converted to
> 
> E11.446588 N49.169494


If I use pj_transform(), I get:
 
  21.982946, 48.433278

So, perhaps pj_datum_transform() provides the correct calculation. But
as I get infinite values, I can't tell.

I am running proj 4.7.0 on Linux. It works for all other uses I have.

In addition to following the recipe, above I have also tried specifying
the values more directly as:

toProj = pj_init_plus(
                "+proj=tmerc "                  // Projection
                "+ellps=bessel "                // Spheroid
                "+k=1.0 "                       // Scale factor at central meridian
                "+x_0=3500000 "                 // False easting
                "+y_0=0 "                       // False northing
                "+lon_0=9.0 "                   // Longitude of central meridian
                "+lat_0=0.0 "                   // Latitude of origin
                "+towgs84=606.0,23.0,413.0 "    // Datum
                "+units=m "
                "+no_defs");

fromProj = pj_init_plus(
                "+proj=latlong "
                "+datum=WGS84");

double  le = 4459750, // easting,
        ln = 5448182, // northing,
        la = 0.0;     // altitude;

pj_transform(toProj,
             fromProj,
             1, 0,
             &le, &ln, &la)

LONG  = le / DEGREE_TO_RADIAN;
LAT   = ln / DEGREE_TO_RADIAN;
ALT   = la;


In the real code, all return values are checked, and there are no
errors.

My goal is to convert these values to/from WGS84 latitude and
longitudes.

Any useful suggestions are greatly appreciated.

TIA



Yours sincerely,

Roger Oberholtzer

OPQ Systems / Ramböll RST

Office: Int +46 10-615 60 20
Mobile: Int +46 70-815 1696
roger.oberholtzer at ramboll.se
________________________________________

Ramböll Sverige AB
Krukmakargatan 21
P.O. Box 17009
SE-104 62 Stockholm, Sweden
www.rambollrst.se