7

I'm trying to implement a 14-parameter Helmert transform to convert between terrestrial reference frames, and having difficulties with it.

I have successfully followed the example in ITRF to GDA94 coordinate transformations, Dawson & Woods, 2010 (Appendix A); d28rz98at9flks.cloudfront.net/71433/71433.pdf.

Using the coordinate in the sample:

[-4052052.36775,4212836.04113,-2545105.10906](ITRF2005@2010.4559)

I use the following to normalise transform parameters:

const T = { x: tx + ṫx*δt, y: ty + ṫy*δt, z: tz + ṫz*δt };
const R = { x: rx + ṙx*δt, y: ry + ṙy*δt, z: rz + ṙz*δt };
const S = 1 + s + ṡ*δt;

Reference epoch is 1994.0 giving δt=(2010.4559-1994)=16.4559; resulting transform parameters are

T: { x: -0.04270,    y: -0.01706,    z:  0.02881 }
R: { x: -1.17163e-7, y: -1.01576e-7, z: -1.03731e-7 }
S-1: 1.1474e-8

Which match those given in the sample (rotation is negated as “Australia assumes the rotations to be of the coordinate axes”).

Applying the transform as follows:

const x2 = T.x + x1*S   - y1*R.z + z1*R.y;
const y2 = T.y + x1*R.z + y1*S   - z1*R.x;
const z2 = T.z - x1*R.y + y1*R.x + z1*S;

I get

[-4052051.76142,4212836.19454,-2545106.01462](GDA94@2010.456)

which matches (to within a mm) the sample transformation.

However, applying the same calculation to the Onsala observatory example from Proj.4 (github.com/OSGeo/proj.4/blob/2aaf53/test/gie/more_builtins.gie#L357), I don't get the correct result.

From coordinate:

[3370658.37800,711877.31400,5349787.08600](ITRF2000@2017.0)

with a reference epoch of 1988.0 giving δt=29; resulting transform parameters are

T: { x: -0.00714,    y:  0.00007,    z: -0.00383 }
R: { x: -1.73563e-8, y: -2.28347e-8, z:  4.31484e-9 }
S-1: 2.2400e-9

I get

[3370658.25318,711877.42306,5349787.15877](ITRF93@2017.0)

as opposed to the correct

[3370658.18890,711877.42370,5349787.12430](ITRF93@2017.0)

which is wrong by [0.06428,-0.00064,0.03447].

The NOAA HTDP agrees with the Proj.4 result, to 3dp.

Could anyone help with what I'm doing wrong?

  • Can you post the exact input values (14 parameters) that you used in your second example, as you typed them in your code? As well as any formulae used up to the calculation of const T? It could help identify the problem more precisely. – FSimardGIS Feb 6 at 1:12
2
+100

It appears that you performed an extra division by 10 somewhere in your code or your values, because the following parameters in your second example:

T: { x: -0.00714, y: 0.00007, z: -0.00383 }

are all 10 times too small, they should be equal to

T: { x: -0.0714, y: 0.0007, z: -0.0383 }

So, multiply these by 10 and it should work.

With these parameters, I arrive at:

[3370658.18892,711877.42369,5349787.12430]

Which matches the correct set of coordinates within 0.0001.

So perhaps your input values for tx, ty, tz were mistyped in your code, or you have an extra '0' somewhere, or a wrong unit conversion, but this is hard to tell without seeing the full code and values.

Also, usually we multiply the whole rotation matrix by S, however, for such small values, it shouldn't be an issue.

Moreover, the paper that you referenced incorrectly states that a milliarcsecond (mas) is equal to (1 * 10e-3 * pi)/(360 * 180) but it should be (1 * 10e-3 * pi)/(3600 * 180), so this could cause some issues if you used those formulas.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.