Coordinate System – Convert X, Y, Z Local Points to WGS84 with Known Coordinates

cartesiancoordinate systemgeometry-conversionwgs84

I have a set of (x, y, z) points on a local coordinates system cartesian plan.
For only 3 of these N points, I have the corresponding WGS84 coordinates (latitude, longitude and altitude).
How can I convert each of the other points in lat/lng coordinates?

EDIT:
Here an example of dataset:

Point # x, y, z lat, lng, h
P1 0, 0, 0 43.88072176381581, 11.072238294111774, 53
P2 0, 12.24, 0 43.88099406334439, 11.072485923398222, 53
P3 32.42, 6.79, 0 43.88080644977896, 11.072808964924704, 53
P4 -12.4, 1, 0.5 ?
P5 55, 8.02, 0.2 ?
P6 0.4, 3, 0 ?
P7 1.03, -187.4, 5 ?

Best Answer

I do it in the following way:

  • Convert geodetic (3D) to geocentric coordinates;
  • Find the parameters of a similitude transformation from local cartesian to geocentric system;
  • Convert all points to geocentric and then all points to geodetic 3D coordinates.

For first and last steps I will use pyproj (2.6.1).
For the second step, I wasn't able to find an open source tool but a great math paper (Huaien Zeng, Xing Fang, Guobin Chang and Ronghua Yang (2018) A dual quaternion algorithm of the Helmert transformation problem. Earth, Planets and Space (2018) 70:26).
I wrote a Python module to partial reproduce the algorithm presented there and published it in my github account (https://github.com/gabriel-de-luca/simil). You can download it to one directory of the Python interpreter Module Search Path, or in the same directory where you save the code that use it.

from pyproj import CRS, Transformer
import numpy as np
import simil

# local cartesian [x, y, z] Control points
local_ctrl_p = [[0, 0, 0],
                [0, 12.24, 0],
                [32.42, 6.79, 0]]

# geodetic [lat,lon,h] Control points
geodet_ctrl_p = [[43.88072176381581, 11.072238294111774, 53],
                 [43.88099406334439, 11.072485923398222, 53],
                 [43.88080644977896, 11.072808964924704, 53]]

# local cartesian [x, y, z] Other points
local_other_p = [[-12.4, 1, 0.5],
                 [55, 8.02, 0.2],
                 [0.4, 3, 0],
                 [1.03, -187.4, 5]]

# CRSes definitions
geodet_crs = CRS.from_epsg(4979) # Geodetic (lat,lon,h) system
geocent_crs = CRS.from_epsg(4978) # Geocentric (X,Y,Z) system

# pyproj transformer object from geodetic to geocentric
geodet_to_geocent = Transformer.from_crs(geodet_crs ,geocent_crs)

# convert geodetic control points to geocentric
geocent_ctrl_p = [geodet_to_geocent.transform(p[0],p[1],p[2])
                  for p in geodet_ctrl_p]

print(np.array(geocent_ctrl_p))

Returns the geodetic control points coordinates converted to geocentric:

[[4518998.16289139  884318.52160765 4398585.26756981]
 [4518973.75937792  884334.02480144 4398607.07516105]
 [4518982.95389214  884362.27851799 4398592.04980611]]

Now let's find the parameters of the 3D linear transformation that better adjust source control points local cartesian coordinates to geocentric ones:

# calculate Cartesian 3D similitude transformation parameters
# from Local Cartesian to Geocentric contol points
m_scalar, r_matrix, t_vector = simil.process(local_ctrl_p,
                                             geocent_ctrl_p)

print('M scalar = ', m_scalar)
print('R Matrix = \n', r_matrix)
print('T Vector = \n',  t_vector)

Returns:

M scalar =  1.2808166135326584
R Matrix = 
 [[-0.00633138 -0.70681673  0.70736838]
 [ 0.98183953  0.12973377  0.13842067]
 [-0.18960762  0.69539863  0.69315922]]
T Vector = 
 [[4518990.78899169]
 [ 884323.63094196]
 [4398591.77207107]]

Control points coordinates must be multiplied by 1.28 to be adjusted to target coordinates. I usually don't have such scale factor (maybe the sample data provided is not real). I usually force the transformation to be adjusted with a fixed scale factor of 1:

# same but force the scale factor to 1
m_scalar, r_matrix, t_vector = simil.process(local_ctrl_p,
                                             geocent_ctrl_p,
                                             scale=False,
                                             lambda_0=1)
print('M scalar = ', m_scalar)
print('R Matrix = \n', r_matrix)
print('T Vector = \n',  t_vector)

Returns:

M scalar =  1.0
R Matrix = 
 [[-0.00634382 -0.7067637   0.70742126]
 [ 0.98183749  0.1297427   0.13842676]
 [-0.18961775  0.69545087  0.69310403]]
T Vector = 
 [[4518989.51051375]
 [ 884326.84158398]
 [4398592.43517149]]

Let's transform the source control points with these parameters. This is easier to me with coordinates instead of points so I transpose:

# transpose source points to get coords
local_ctrl_coords = np.array(local_ctrl_p).T

# transform the control coordinates
transf_ctrl_coords = m_scalar * r_matrix @ local_ctrl_coords + t_vector

# transpose transformed control coordinates to get points
transf_ctrl_p = transf_ctrl_coords.T

print(transf_ctrl_p)

Returns:

[[4518989.51051375  884326.84158398 4398592.43517149]
 [4518980.85972611  884328.42963463 4398600.94749011]
 [4518984.50592159  884359.55370846 4398591.00987537]]

There are (compensated) differences in meters. I hope the sample data provided is just not real at all. Now, transform all other points to geocentric and then to geodetic:

# transform other points coordinates with same parameters
local_other_coords = np.array(local_other_p).T
transf_other_coords = m_scalar * r_matrix @ local_other_coords + t_vector
transf_other_p = transf_other_coords.T

# convert other points to geodetic
geocent_to_geodet = Transformer.from_crs(geocent_crs ,geodet_crs)
geodet_other_p = [geocent_to_geodet.transform(p[0],p[1],p[2])
                  for p in transf_other_p]

print(np.array(geodet_other_p))

Wich returns what we was looking for:

[[43.88084931 11.07221498 53.49981423]
 [43.88075069 11.07304707 53.19980365]
 [43.88083637 11.07237519 52.9998092 ]
 [43.87918162 11.07175954 57.98771276]]
Related Question