2

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 ?
4
  • I think this is answered here, or at least it will point you to the logic and highlight a couple of issues: gis.stackexchange.com/questions/255254/… – Leigh Bettenay Feb 4 at 23:23
  • Welcome to GIS.SE, Asia. I do it in the following way: convert geographic to geocentric coordinates; find the parameters of a 3D similitude transformation from local cartesian to geocentric system; convert all points to geocentric and then all points to geographic coordinates. If you want to add the three trio of coordinates let me know and I will write an answer. Some Python code involved. – Gabriel De Luca Feb 5 at 6:01
  • @LeighBettenay Thanks, I'll give a look. – Asia Feb 5 at 11:10
  • @GabrielDeLuca I've edited the question, adding a sample dataset. No problem with Python code. Thank you so much! – Asia Feb 5 at 11:11
2

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]]
15
  • Thank you so much! First of all, yes, the dataset is not real. Unfortunately, I need to code in C#, cause the project I'm working on is an Hololens 2 project. By the way, for the conversion between geodetic and ECEF coordinates (and back), I already wrote some functions on my own, so the lack of the pyproj library is not a problem. For the parameters estimation, I will take a look at your library to see if I could write some similar procedure in C#. Again, thanks!! – Asia Feb 5 at 22:56
  • Hi, you are welcome. Glad to know that can be useful, make it open :-) You are welcome again! – Gabriel De Luca Feb 6 at 4:48
  • Gabriel, sorry, I'm here again. :( I tried to convert your python script to C#, using Numpy for .NET, but it gives me a lot of errors, due to different functions in the library. Could you please help me in some way? I tried very hard to understand the maths behind this parameter estimation, but it's far too difficult for me to understand deeply how to do that. Thanks in advance! – Asia Apr 4 at 16:04
  • Hi @Asia, For some of the matrix operations I used numpy's einsum function, and I think the problem may be there. I think operations could be written through other less compact functions. I reproduced the whole paper algorithm by hand, in the sense of writing all the matrices for a simple case of three pairs of coordinates, doing all the operations corresponding to each element of each matrix, in order to understand it and translate it into python. I'm not a great developer or a great mathematician, but as far as I can help you, I will gladly do so. – Gabriel De Luca Apr 4 at 16:52
  • Exactly, I'm struggling on the einsum function. I'm trying to find another way to write all the mathematical operations, without using it. For example, could you please help me understand this line of code (contained in _get_abc_matrices): matrix = np.einsum('i,ijk->jk', alpha_0, np.transpose(m1, (0,2,1)) @ m2) From numpy docs, I understood that the third parameter is an output variable... so why are you outputting the result in a matrix multiplication? How is it possible? Maybe I misunderstood something... Thanks! – Asia Apr 4 at 17:03

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.