Implementing Helmert transforms between reference frames / coordinate systems?

Question

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]([email protected])

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]([email protected])

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]([email protected])

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]([email protected])

as opposed to the correct

[3370658.18890,711877.42370,5349787.12430]([email protected])

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?

Answer 1

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.