[Proj] WGS84 to ED50 - SOLVED

[Roger Oberholtzer]

On Thu, 2007-11-29 at 14:55 +0100, Roger Oberholtzer wrote:
> I have a bit more information on the problem with WGS84 lat/longs ->
> ED50(Spain) northings and eastings.
> 
> The software that is being used does not have a 'Spanish' setting. So
> they have defined an ED50 7 parameter transform, with these
> coefficients:
> 
> 
> > International Ellipsoid 1924 definition
> > Mayor semiaxle: 6378388
> > Minor semiaxle:6356911,9461
> > 1/f:297
> > Excentricity: 0.0819910890032286
> > 
> > About the 7 parameters transformation 
> > Translation X (m): 125.098545
> > Translation Y (m): 76.000054
> > Translation Z (m): 156.198703
> > X Rotation(seconds): 0
> > Y Rotation(seconds): 0
> > Z Rotation(seconds): -1.129
> > Scale factor (ppm): -8.30463103
> > 
> > This correspond with applanix header file saying 
> > Mapping frame datum: ED50_ESPAÑA ; a = 6378388.000000; 1/f = 297.000000;
> >  Coordinate transformation from WGS84 to mapping frame datum 
> >     dX = 125.098545; dY = 76.000054; dZ = 156.198703; f = 0.999991695369; 
> >     R1 = 0.000000000000; R2 = 0.000000000000; R3 = -0.000005473550;
> >  sequence of the rotations: x,y,z,1 
> >  Mapping frame projection : TM;
> >  central meridian = -3.000000 deg;
> >  latitude of the grid origin = 0.000000 deg; grid scale factor = 0.999600: 
> >  false easting = 500000.000000 m; false northing = 0.000000 m; 
> > 
> >  
> > Where Rotations are in radians 
> 
> I have tried setting proj to use this, as reported earlier in this thread, with no luck. I tried using the 7 parameters as 
> +towgs84=125.098545,76.000054,156.198703,0.0,0.0,-1.129,-8.30463103, but that made things worse.
> 
> fromProj = pj_init_plus(
>                         "+proj=latlong "
>                         "+ellps=WGS84 "
>                         "+no_defs");
> 
> toProj = pj_init_plus(
>                         "+proj=utm "
>                         "+ellps=intl "
>                         "+k=0.999600 "
>                         "+x_0=500000.0 "
>                         "+y_0=0.0 "
>                         "+zone=30 "
>                         "+units=m"
> 
> "+towgs84=125.098545,76.000054,156.198703,0.0,0.0,-1.129,-8.30463103"
>                         "+no_defs");
> 
> So I am guessing I cannot simply use the 7 parameters as given. Correct?

The towgs84 parameters I had been tying to use were correct - except in
sign. The best settings for matching Applanix/Trimble and a few other
packages is:

fromProj = pj_init_plus( "+proj=latlong "
                                       "+ellps=WGS84 "
                                       "+towgs84=0,0,0 "
                                       "+no_defs");

toProj   = pj_init_plus( "+proj=utm "
                                     "+ellps=intl "
                                     "+k=0.999600 "
                                     "+x_0=500000.0 "
                                     "+y_0=0.0 "
                                     "+zone=30 "
                                     "+units=m "

"+towgs84=-125.098545,-76.000054,-156.198703,0.0,0.0,-1.129,8.30463103 "
                                     "+no_defs");

Note the difference in signs from the earlier one. The earlier one had
at least a 1 meter error. Often more than 2 meters. The one just above
is at worst within 0.03 meters. Quite often it is much better. Often
0.0000 meters :)

I am preparing to be happy.

-- 
Roger Oberholtzer

OPQ Systems / Ramböll RST

Ramböll Sverige AB
Kapellgränd 7
P.O. Box 4205
SE-102 65 Stockholm, Sweden

Office: Int +46 8-615 60 20
Mobile: Int +46 70-815 1696