[Proj] DHDN / Gauss Krueger to WGS84
Roger Oberholtzer
roger at opq.se
Fri Mar 30 04:48:53 EST 2012
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
More information about the Proj
mailing list