[Proj] datum conversion, estimating parameters
Robert Jaeschke
robert-j at gmx.de
Tue Mar 22 16:30:39 EST 2005
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
More information about the Proj
mailing list