There is not an affine relation between geographic (or geodetic) and cartesian coordinates. You can use a projection as an intermediate step, but I usually avoid the projection deformation converting geodetic to cartesian geocentric coordinates, and estimating the affine parameters from geocentric to local (both 3D Cartesian) systems.
I wrote a Python module to estimate those parameters: https://github.com/gabriel-de-luca/simil. In order to install it, just download the simil.py file and make sure that its directory is included in the PYTHONPATH environment variable. simil only depends on numpy: https://numpy.org/.
To transform the geodetic to geocentric coordinates I used pyproj: https://pyproj4.github.io/pyproj/stable/.
To get the coordinates from the LOC file I used lmxl: https://lxml.de/.
I just invented three points in the following some_points.txt file:
<record id="Localization Points" >
<record id="Point 1" >
<value name="Lat" value="0"></value>
<value name="Lon" value="0"></value>
<value name="Ellipsoid_Elv" value="20.1497127691777"></value>
<value name="Local_X" value="883498.7601626500300"></value>
<value name="Local_Y" value="93747.3358902813200"></value>
<value name="Local_Z" value="152.2231851663858"></value>
<value name="HRMS" value="0.019685"></value>
<value name="VRMS" value="0.03937"></value>
<value name="Use_Horizontal" value="No"></value>
<value name="Use_Vertical" value="No"></value>
<value name="Description" value="1050"></value>
</record>
<record id="Point 2" >
<value name="Lat" value="1"></value>
<value name="Lon" value="2"></value>
<value name="Ellipsoid_Elv" value="40.1497127691777"></value>
<value name="Local_X" value="773498.7601626500300"></value>
<value name="Local_Y" value="39747.3358902813200"></value>
<value name="Local_Z" value="512.2231851663858"></value>
<value name="HRMS" value="0.019685"></value>
<value name="VRMS" value="0.03937"></value>
<value name="Use_Horizontal" value="No"></value>
<value name="Use_Vertical" value="No"></value>
<value name="Description" value="1050"></value>
</record>
<record id="Point 3" >
<value name="Lat" value="2"></value>
<value name="Lon" value="1"></value>
<value name="Ellipsoid_Elv" value="2.1497127691777"></value>
<value name="Local_X" value="993498.7601626500300"></value>
<value name="Local_Y" value="77747.3358902813200"></value>
<value name="Local_Z" value="215.2231851663858"></value>
<value name="HRMS" value="0.019685"></value>
<value name="VRMS" value="0.03937"></value>
<value name="Use_Horizontal" value="No"></value>
<value name="Use_Vertical" value="No"></value>
<value name="Description" value="1050"></value>
</record>
</record>
The following script was tested with Python 3.8.8, numpy 1.19.2, pyproj 3.1.0 and lxml 4.6.3:
from lxml import etree
from pyproj import CRS, Transformer
import simil
# Read the text file, create the XML object and get the points lists
with open('some_points.txt') as file:
some_points_xml = etree.fromstring(file.read())
n = 3 # Number of points to get their coordinates
geodetic_points = [[float(some_points_xml[points][coordinates].get("value"))
for coordinates in range(3)]
for points in range(n)]
print('Geodetic points = \n', geodetic_points)
local_points = [[float(some_points_xml[points][coordinates].get("value"))
for coordinates in range(3,6)]
for points in range(n)]
print('Local points = \n', local_points)
# Convert geodetic to geocentric coordinates
geodet_crs = CRS.from_epsg(4979) # Geodetic (lat,lon,h) WGS84 system
geocent_crs = CRS.from_epsg(4978) # Geocentric (X,Y,Z) WGS84 system
geodet_to_geocent = Transformer.from_crs(geodet_crs ,geocent_crs)
geocentric_points = [geodet_to_geocent.transform(p[0],p[1],p[2])
for p in geodetic_points]
print('Geocentric points = \n', geocentric_points)
# calculate Cartesian 3D similitude transformation parameters
# from Geocentric to Local points
m_scalar, r_matrix, t_vector = simil.process(geocentric_points,
local_points)
print('M scalar = ', m_scalar)
print('R Matrix = \n', r_matrix)
print('T Vector = \n', t_vector)
# Define PROJ Affine transformation parameters
print('x_off = ', t_vector[0][0])
print('s11 = ', m_scalar*r_matrix[0][0])
print('s12 = ', m_scalar*r_matrix[0][1])
print('s13 = ', m_scalar*r_matrix[0][2])
print('y_off = ', t_vector[1][0])
print('s21 = ', m_scalar*r_matrix[1][0])
print('s22 = ', m_scalar*r_matrix[1][1])
print('s23 = ', m_scalar*r_matrix[1][2])
print('z_off = ', t_vector[2][0])
print('s31 = ', m_scalar*r_matrix[2][0])
print('s32 = ', m_scalar*r_matrix[2][1])
print('s33 = ', m_scalar*r_matrix[2][2])
And it returns:
Geodetic points =
[[0.0, 0.0, 20.1497127691777], [1.0, 2.0, 40.1497127691777], [2.0, 1.0, 2.1497127691777]]
Local points =
[[883498.76016265, 93747.33589028132, 152.2231851663858], [773498.76016265, 39747.33589028132, 512.2231851663857], [993498.76016265, 77747.33589028132, 215.2231851663858]]
Geocentric points =
[(6378157.149712769, 0.0, 0.0), (6373327.3986458145, 222561.49698355066, 110569.47553367223), (6373308.912659603, 111246.52087806119, 221104.62033501064)]
M scalar = 0.46614318279723016
R Matrix =
[[-0.00251484 -0.61078012 0.79179626]
[ 0.02641048 -0.79156314 -0.61051641]
[ 0.99964802 0.01937637 0.01812165]]
T Vector =
[[ 881845.64010512]
[ 64451.65803851]
[-2972230.43549093]]
x_off = 881845.6401051179
s11 = -0.0011722774288244159
s12 = -0.28471098869694794
s13 = 0.36909043004278536
y_off = 64451.65803850801
s21 = 0.012311065853214564
s22 = -0.3689817608571938
s23 = -0.2845880613810327
z_off = -2972230.4354909346
s31 = 0.4659791092865856
s32 = 0.009032161328442258
s33 = 0.008447285868804786
Note that you must transform with PROJ the geodetic coordinates to geocentric before to apply the affine transformation parameters returned. You can use any other intermediate system in the script by changing the CRS definitions.
