[Proj] datum conversion, estimating parameters

Tue Mar 22 21:22:04 EST 2005

Robert,

To accomplish what you intend you need the ECEF (geocentric Cartesian) 
coordinates for at least one point in both datums.  Differencing those ECEF 
coordinates will give the the three translations (the first 3 of the 7 
parameters, the rest of the parameters being zero), which is all you can 
hope for with one point.  Solving for the other parameters requires many 
more points over (and this is very important) a LARGE area in order to 
achieve any statistical significance.  The problem is ill-conditioned over a 
small area.

This problem is NOT solved in projected space.

Noel Zinn
Houston

----- Original Message ----- 
From: "Robert Jaeschke" <robert-j at gmx.de>
To: <proj at xserve.flids.com>
Sent: Tuesday, March 22, 2005 3:30 PM
Subject: [Proj] datum conversion, estimating parameters

> Hello,
>
> I'm trying to calculate the seven +towgs84 parameters for an unknown
> datum. Therefore I take known coordinate values to solve an eqation
> system. But test-computation with cs2cs doesn't give correct results.
>
>
> I have the following data:
>
> Salzfurtkapelle 12d10'27.403"E  51d41'38.33"N
> Horstwalde      13d24'24.199"E  52d5'4.063"N
> Jeßnitz         12d18'50.019"E  51d41'7.261"N
>
> which is in latlong/wgs84 format. I transformed it to GKK/wgs84:
>
> awk '{print $2" " $3 " " $1}' < NEWDATA_DMS | proj +init=epsg:31464
>
> 4512049.07      5728719.87 Salzfurtkapelle
> 4596414.79      5773081.35 Horstwalde
> 4521705.77      5727792.03 Jeßnitz
>
> Let's call this NEWDATA.
>
> Furthermore I have six corresponding values (I don't know their datum ...
> and their projection, but let's disregard the last, the goal is, to
> discover their datum, assumed their projection is GKK):
>
> Salzfurtkapelle 4507038 5729278
> Horstwalde      4590582 5774734
> Jeßnitz         4516584 5728546
>
> Let's call this OLDDATA
>
> I used NEWDATA and OLDDATA to solve the following equation system:
>
>          4512049.07 = M*(    4507038 - Rz*5729278 + Ry*0.0) + Dx
>          5728719.87 = M*( Rz*4507038 +    5729278 - Rx*0.0) + Dy
>                0.00 = M*(-Ry*4507038 + Rx*5729278 +    0.0) + Dz
>          4596414.79 = M*(    4590582 - Rz*5774734 + Ry*0.0) + Dx
>          5773081.35 = M*( Rz*4590582 +    5774734 - Rx*0.0) + Dy
>                0.00 = M*(-Ry*4590582 + Rx*5774734 +    0.0) + Dz
>                0.00 = M*(-Ry*4516584 + Rx*5728546 +    0.0) + Dz
>
> Which is basically (2 times + one extra equation):
>
> x_out = M_BF*(       x[io] - Rz_BF*y[io] + Ry_BF*z[io]) + Dx_BF;
> y_out = M_BF*( Rz_BF*x[io] +       y[io] - Rx_BF*z[io]) + Dy_BF;
> z_out = M_BF*(-Ry_BF*x[io] + Rx_BF*y[io] +       z[io]) + Dz_BF;
>
> from http://www.remotesensing.org:16080/proj/gen_parms.html
>
> (I took 7 values für x,y,z so that I have an unique solution and assumed
> z to be 0.0 in both NEWDATA and OLDDATA).
>
>
> The solution is
>
> Dx = -85976.51155
> Dy =  51643.25395
> Dz =      0.0
> Rx =      0.0
> Ry =      0.0
> Rz =     -0.01420808651
> M  =      1.002089055
>
> Then I took these parameters, to transform OLDDATA into NEWDATA (?):
>
> awk '{print $2 " " $3 " 0.0 " }' < OLDDATA | cs2cs \
> +proj=tmerc +lat_0=0 +lon_0=12 +k=1.000000 +x_0=4500000 +y_0=0 
> +ellps=bessel \
> +towgs84=-85976.51155,51643.25395,0.0,0.0,0.0,-0.01420808651,1.002089055 
> +to \
> +proj=tmerc +lat_0=0 +lon_0=12 +k=1.000000 +x_0=4500000 +y_0=0 
> +ellps=bessel
>
> and got these results:
>
> 4575963.52      5787260.83 -44755.54
> 4660087.33      5833043.91 -43436.32
> 4585577.06      5786523.79 -44659.19
>
> What I expected to get, was NEWDATA, which is:
>
> 4512049.07      5728719.87 Salzfurtkapelle
> 4596414.79      5773081.35 Horstwalde
> 4521705.77      5727792.03 Jeßnitz
>
> (with z=0.0)
>
> because OLDDATA together with the above solutions of the equation system
> should give the left hand side of the equations.
> But that didn't happen. So what did I wrong? Where is my error in
> reasoning?
>
> I know, it's very confusing what I did and I'm not sure, what I really 
> did.
> But maybe someone can help me to remove all confusion and mistakes. :-)
>
>
> What I still want to do, is to find out projection/datum of the D-Sat
> data (http://www.math.tu-dresden.de/~jaeschke/dsat/) from where I got
> OLDDATA (NEWDATA contains the real-world coordinates).
>
> But it seems to me, that a pure datum shift gives always constant
> differences in easting/northing deviation. But that is not the case with
> the D-Sat data, as the examples on the homepage and especially the new
> difference vector image show.
>
>
>
> Thank you again for your help and your patience.
>
>
> Robert Jaeschke
>
> -- 
> I really didn't foresee the Internet. But then, neither did the
> computer industry. Not that that tells us very much of course -
> the computer industry didn't even foresee that the century was
> going to end.       --  Douglas Adams, introduction to h2g2.com
> _______________________________________________
> Proj mailing list
> Proj at xserve.flids.com
> http://xserve.flids.com/mailman/listinfo/proj
> 