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 http://gis.stackexchange.com/a/8674.)
This encompasses two geometrical questions:
What are the geocentric coordinates corresponding to a given location given in terms of latitude and longitude (at a nominal height of zero)?
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)
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
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.