Local Coordinate to Geocentric

Question

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 1

The question asks how to convert local coordinates into geocentric coordinates. As an example,

What would be the geocentric coordinates of a point displaced 250 meters west and 250 meters south of -108.619987 degrees longitude, 36.234600 degrees latitude (at sea level)?

(For the reason why only six decimal places are used here, please see https://gis.stackexchange.com/a/8674.)

This encompasses two geometrical questions:

1. What are the geocentric coordinates corresponding to a given location given in terms of latitude and longitude (at a nominal height of zero)?

2. What rotation will bring the local frame (in which x points east, y points north, and z points up) into the geocentric frame (in which x points towards the Prime Meridian at the Equator, z points at the North Pole, and x-y-z form a positively oriented orthogonal coordinate system)?

Let's solve these geometrically. Because the math for the first one has been done elsewhere on this site, the new part is the second.

To begin with, consider the geocentric frame at the North Pole:

In this and subsequent images, the blue arrow represents the local "x" direction, the red arrow the local "y" direction, and the green arrow the local "z" direction. At the outset, the three local directions correspond with the geocentric axes as shown.

(Sharp-eyed readers will notice the relatively extreme eccentricity of this representation of the ellipsoid, especially in the third and fourth images: the Earth looks somewhat "squashed." This is by design: the software that made these images uses the equations worked out below. Exaggerating the eccentricity emphasizes potentially small errors that can creep into ellipsoidal calculations, making them easier to spot and rectify.)

For reasons that will soon be clear, first rotate the local frame 90 degrees counterclockwise around its z-axis:

Next, rotate this frame around the geocentric y-axis to bring it to the desired latitude. (For this illustration we're headed toward Indonesia, in the southern hemisphere.)

(Notice how the red arrow (local y) is pointing due north and the blue arrow (local x) is pointing due east: it was to establish this x=east, y=north correspondence that the first rotation was done.)

Finally, rotate this frame once more around the geocentric z-axis to bring it to the desired longitude.

These three rotations provide the Euler angles in the Z-Y'-Z'' convention. The amounts of rotation were:

• First, 90 degrees around the z axis from x towards y.

• Second, 90 - 'latitude' degrees around the y axis from z towards x.

• Third, 'longitude' degrees around the z axis from x towards y.

(It should now be evident that these transformations were worked out in reverse: starting with a local east-north-up frame at a given point, a rotation around the earth's axis (geocentric z) can take us to the Prime Meridian, and thence a rotation around the geocentric y axis can take us to the North Pole, where we discover a simple 90 degree rotation around the pole will align the local frame with the geocentric one. The key idea is to use rotations around the geocentric axes to line up the frames: although there are many ways to accomplish this, they will all produce the same answer.)

It remains to convert the geometry into arithmetic. The 3 x 3 matrix representing a rotation of t from axis i towards axis j will have values of cos(t) in its (i,i) and (j,j) positions, sin(t) in the (j,i) position, -sin(t) in the (i,j) position, a 1 in the other diagonal position, and zeros elsewhere. Successive rotations are computed by multiplying the matrices (from right to left). Writing f for the longitude and q for the latitude, the resulting rotation matrix is

-sin(f), -cos(f) sin(q), cos(f) cos(q)
 cos(f), -sin(f) sin(q), sin(f) cos(q)
      0,         cos(q),        sin(q)

For the values of f and q in the question, the matrix is

 0.94765709  0.1887300   -0.2575404
-0.3192899   0.56015335  -0.76438353
 0           0.80660351   0.59109287

The first column is a unit-length vector giving the local x (east) direction in terms of the geocentric coordinates (as a check, notice it has no third component: x is entirely east-west). The second column is a unit vector giving the local y (north) direction and the third column is the local outward unit normal, or "up" direction.

As far as the origin of the local system goes, we can work it out from properties of the Earth's ellipsoid. The WGS 84 ellipsoid, for instance, has a major axis of a = 6378137 meters (that's the equatorial radius) and a minor axis of b = 6356752.3142 meters (that's the polar radius). (The ratio of these axes equals 0.996647189335, which is less than 1 by the "flattening," equal to 1 / 298.257223563.) Consequently, the geocentric coordinates of the origin are given by

a cos(q') cos(f) 
a cos(q') sin(f) 
b sin(q')

where q' = Atan(b/a tan(q)) converts the geodetic latitude q to the angle between the origin and the equator. In the example, these coordinates (in meters) are

-1,644,552.1
-4,881,054.5
 3,749,220.2

Finally, to find the geocentric coordinates of a local displacement from this origin, we form the linear combination in terms of the geocentric frame. E.g., a displacement of 250 meters west is -250 in the local x direction, equivalent to a geocentric displacement of

-250 * (0.94765709, -0.3192899, 0) = (-236.9, 79.8, 0)

and a displacement of 250 meters south is -250 in the local y direction, equivalent to a geocentric displacement of

-250 * (0.1887300, 0.56015335, 0.80660351) = (-47.2, -140.0, -201.7).

Applying these two displacements to the geocentric origin (by adding them component-by-component) yields geocentric coordinates of

(-1644836.2, -4881114.7, 3749018.6).

As the question notes, all this can be compactly represented as a four by four transformation matrix: its upper left three by three block is the rotation matrix, the first three entries in its rightmost column are the geocentric coordinates for the given longitude and latitude, and its bottom row is (0,0,0,1). Local coordinates (x',y',z') are converted to geocentric coordinates by left-multiplying the column vector (x',y',z',1)' by this matrix.

Incidentally, the earth's surface is usually not at sea level at land locations like this one. No problem: the surface elevation at the origin represents a displacement in the local z direction. Simply add that multiple of the local z vector to obtain the final value. In the example, the surface elevation at the origin is approximately 2000 meters. The corresponding geocentric displacement therefore equals

2000 * (-0.2575404, -0.76438353, 0.59109287) = (-515.081, -1528.77, 1182.19)

Adding this to the previous result gives

(-1645351.3, -4882643.5, 3750200.7)

Because non-sea-level contributions are considerable, elevation cannot be ignored when computing geocentric coordinates.

Answer 2

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.

Answer 3

What you essentially have is a plane that you want to rotate and translate so it rests on the surface of a spheroid, with its normal congruent with the sphere's normal at that point.

Your first task therefore is to calculate a rotation matrix so that the normals align. Looking at this Wolfram page, you can see that a R3 rotation can be broken down into 3 separate transforms that can be multiplied together. Now, assuming that you'll be avoiding the poles, and that your forward vector will always point north, you can treat one of the matrices as identity. Which one depends on whether you're using a left or right-handed coordinate system, and whether Z or X points through the north pole.

So your two remaining rotation matrices can be created by substituting the angles represented by alpha, beta, and/or gamma with the latitude and longitude of the point on the surface you want the plane to touch, 36.234600064 and -108.619987456 in your case. Then you need to construct a translation matrix that represents the radius of the globe. Just multiply them all together and there you have it!

Now, this assumes that your globe is spherical, which as long as you steer clear of the poles, deal with small areas, and don't expect millimetre accuracy, is OK. If you want more accuracy, you'll have to take the ellipsoidal nature of the Earth, which can be easily simulated with a slight scaling factor in the north pole axis.

You also need to make sure that the axes you use are the same as those in the Wolfram article; if not, you may have to play about with adding or subtracting 90 degrees to your longitude value (not forgetting to modulus with 360), and possibly choosing the gamma instead of the alpha matrix for instance. If I had a bit more time, I'd work through your example to figure it out, sorry!