[Proj] datum conversion, estimating parameters

[Robert Jaeschke]

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