Friday 22 September 2017

gps - Local Coordinate to Geocentric


The ultimate question, I need a transformation matrix to take a point in local space representing a roughly 500m x 500m place in New Mexico centered at -108.619987456 long 36.234600064 lat. The final output needs to be in geocentric coordinates and I can during initialization get any number of points on the map such as the corners, center, western most along the center etc. During run time I would like to have a transformation matrix that can be applied to a position in local space expressed in meters and get back an approximate GCC coordinate.


The center of the point in GCC: -1645380 -4885138 3752889
The bottom left corner of the map in GCC: -1644552, -4881054, 3749220



During run time I need to be able to multiply (-250, -250, 0) with the transformation matrix and get something close to (-1644552, -4881054, 3749220).


I've spent more than a week trying to research a solution for this and so far nothing.


Given an identity matrix as our starting position and orientation. Then using geotrans and the known lat, long, and height of the starting position I get an x,y,z geocentric coordinate. A vector from the origin of (0,0,0) gives both the up and translation for the matrix. However, I need the forward and right so that I can pass a distance in meters from the origin into the transformation matrix and get a roughly accurate GCC. Do I have all of the inputs I need to calculate the right and forward?


Inputs
Origin: 0, 0, 0
Global: -1645380, -4885138, 3752889
Up (Normalized): Global - Origin


Desired Outputs
Right: ? ? ?
Forward: ? ? ?



So with the right and forward added to the up and translation I already calculated I would have a transformation matrix. I could then apply the matrix to a vector of say (50, 50, 0) and get something within about 0-3 cm of the output if I feed the lat/long back into geotrans. This matrix would only ever be used small maps of about 500m x 500m so the curvature of the earth is negligible.


In reply to whuber,


I don't know your level of experience so I will start with the very basics.


An identity matrix is a 4x4 matrix that you can multiply by a vector and it returns the vector unchanged. Such as below.


1 0 0 x=0
0 1 0 y=0
0 0 1 z=0
0 0 0 1


The x, y, and z of the matrix are how much to translate the vector by and the first 3 rows and columns are the rotation. Below is a tranformation matrix that does nothing to the orientation of a vector but does translate. This is what you would want for two axis aligned worlds in different locations.


1 0 0 13

0 1 0 50
0 0 1 -7
0 0 0 1


If you were to multiply a vector of (10, 10, 10) with the above transformation matrix you would get an output of (23, 60, 3).


My problem is the axes are not aligned and the "plane" of the world I am trying to get the local coordinate converted to is projected on the ellipsoid of the earth.


Geotrans is library that you can use to pass a coordinate from one system such as geodetic or GDC (lat, long, height) and get back another such as geocentric or GCC (x,y,z).


For example: If my game world was representing a map centered exactly on the prime meridian and equator the transformation matrix would look something like below. I am just guesstimating the Y and Z rotations and might have them flipped and/or the wrong sign.


0 0 1 6378137
0 1 0 0
1 0 0 0

0 0 0 1



Answer



A third option and the one I ended up using is simple and works for anywhere but with the possible exception of the poles with accuracy for the NM example at less than 2 feet (58 cm).


Given a single input of the geodetic coordinate (GDC)


The first two steps have a lot done "in the background" by geotrans which is available for at least C/C++ and Java. It accounts for the ellipsoidal nature of the earth and can do a number of conversions from one coordinate system to another.
1) Convert GDC to a Geocentric coordinate (GCC) = vector3 Origin
2) Get the GCC for a point 100m above the origin = vector3 AboveOrigin


The "up" component of the matrix is really easy and is just a vector from the point above the origin to the origin. While I was implementing this I thought it was backwards and should be AboveOrigin - Origin but for reasons I don't understand it is as described.
3) vector3 Up = Origin - AboveOrigin
4) Normalize Up



You need a vector that is perpendicular to both your true forward or north vector and your true right. The cross product of your true up and straight up positive Z, which is through the axis in GCC for geotrans, will give you a temp right vector.
5) vector3 temp = Up crossed with vector3 (0,0,1)
6) Normalize temp


This is your final forward or north vector
7) vector3 forward = Up cross temp
8) Normalize forward


Now recalculate the right or east vector
8) vector3 right = Up cross temp
9) Normalize right
10) Initialize your 4x4 matrix with forward, right, up, and origin



For example: On our map we have tags that describe the top, bottom, left, and right edges of the map in decimal degrees. From there we get the center and point above the center.


Center - (-1645379.875, -4885137.5, 3752889)
Above - (-1645406, -4885214, 3752948)


With just those two points you can follow the steps and get a matrix the can be applied to a coordinate in local space and get back a reasonably accurate point in GCC.


To test this I found 9 points across the map through geotrans, the center of the map, the four corners, and the four points along the middle of each edge. Ignoring the curvature of the earth and knowing that this is a small map I just took the average distance from the center of the map in GCC to each of the four edges. I then used that X/Y distance to generate 9 local coordinates to pass through the matrix and compare to the point found through geotrans. In testing the distance of the transformed point against the geotrans point the worst result was for the top middle and top right with 57.3cm of inaccuracy with the left middle being 0 or near 0.


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...