1

so I need to perform some kind of transformation. I have a bunch of points like so:

<record id="Localization Points" >
<record id="Point 1" >
<value name="Lat" value="redacted"></value>
<value name="Lon" value="-redacted"></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>

and I have the values that proj needs here:

picture

Source: https://proj.org/operations/transformations/affine.html

How do I figure out the parameters needed for the affine transformation to a local grid with the points like I listed above? That is just one point but each file may have multiple, up to a dozen

3
  • Ok, I'm guessing you've already poked around the proj website but I'll add this link for anyone who may click on this q down the road. proj.org/about.html
    – Lisa
    Jul 6 at 17:37
  • To be able to do any transformation, you need to know where you are and to know the destination point. What you provide us is not enough to give any answer if you really want to make an Affine transformation. From my understanding, what you want to do the job is not an affine transformation but a reprojection. You do not need to try to transform your points yourself using something like an affine transformation but using the state plane definition and the local grid definition, you just do a reprojection. What are the name of your "state plane definition" and of your "local grid"?
    – ThomasG77
    Jul 6 at 20:57
  • Hi @ThomasG77 Thanks for your time. Im getting lat/long & converting that to northing/easting per nad83 datum. On top of that the user has a set of points, like the one I have included in my question, that act to localize or calibrate their site. What I'm not sure about is how I work the values into what equation or transformation that can convert my northing/easting or lat/long into the coordinates that theyve arbitrarily set to localize their job site Jul 7 at 3:22
1

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.

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.