I am trying to create a program in Python that will let me check the Geodetic Transformations that the geodesy app in my survey software spits out. I have a test point and know that my program currently returns total rubbish.
The formula that I am using is the 7 Parameter Halmert Transformation:
Vb = T + (1 + s * 10^−6)•R•Va
My code is:
import numpy as np
from math import pow
# Input Parameters
Tx = 89.5 # in metres
Ty = 93.8
Tz = 123.1
Rx = 0 # in arcsecs
Ry = 0
Rz = 0.156
s = -1.2 # ppm
Vax = 58.000086817 # Lat
Vay = 0.723469322 # Long
Vaz = 0 # Height
# Arrays
T = np.array([[Tx], [Ty], [Tz]])
R = np.array([[1, -Rz, Ry],
[Rz, 1, -Rx],
[-Ry, Rx, 1]])
Va = np.array([[Vax], [Vay], [Vaz]])
# Calculations
RVa = R * Va
RVa = np.sum(RVa, axis=1)
RVa = np.array([[RVa[0]], [RVa[1]], [RVa[2]]])
S = 1 + s * pow(10, -6)
SRVa = S * RVa
Vb = T + SRVa
The correct output is:
Lat = 58.0007206472222 Long = 0.72507962777778
I am a long way off that.
I know that there is something funny happening as I tried to transpose RVa after it was summed and nothing happened (hence the manual transposition. I have clearly messed something up badly and appreciate any help.
Cheers,
The Helmert transformation is for mapping a point in cartesian (ECEF, aka earth centred earth fixed) coordinates relative to one coordinate system, to the cartesian coordinates relative to another coordinate system.
If you want to transform lat,long,height coordinates relative to the first coordinate system to lat,long,height coordinates relative to a second, the drill is:
Convert the lat,long,height coordinates to cartesians
Apply Helmert to those cartesians, getting cartesians relative to the second coordinate system
Convert those cartesians to lat,long,height in second coordinate system
A couple of points:
Helmert transformation parameters are usually quoted for ellipsoidal models of the earth. Ignoring this and using a spherical earth will introduce errors, on the scael of tens or even hundreds of metres.
When converting from and to cartesian coordinates the heights are elliposidal heights, that is the distance along the normal to the ellipsoid from the ellipsoid to your point. You could well have a 0 height in one system transform to a non zero point in another. However, if you are close to the earth's surface the error induced by ignoring the heights will be small (typically sub metre).