# How to translate coordinates on the earth

I am trying to translate coordinates. I have a set of coordinates for city A which I need to map to city B, given a point P1 in city A maps to P2 in city B. How do I do this? This problem would be similar to a translation transformation but in this case, it's on a sphere and not on a plane.

-
What coordinate systems do the two cities use? – R.K. Feb 8 '12 at 7:38
Jeez, how big is the city? – Hairy Feb 8 '12 at 14:28

On the plane, we do this with least-squares fitting. It's easy there, because the x and y coordinates can be fit separately. On the sphere, using long-lat as coordinates, the interactions among longitude and latitude require simultaneous fitting (one degree of longitude changes length as latitude varies). But the procedure is the same: use least squares. For simplicity--while maintaining high accuracy and automatically coping with the discontinuities at the poles and +-180 meridian--I suggest minimizing squared Euclidean (3D) distances, rather than geodesic distances.

To accomplish this, you will need to code (or have access to) a modestly capable nonlinear regression routine. These are pretty common; you can use `R`'s nonlinear solver, for instance, and I'm sure that Python experts could suggest comparable resources.

Here are some details. The input data consist of "links": each is a (source-->target) pair with both source (in city A) and target (in city B) given as (long, lat). Convert these to Cartesian coordinates. The sphere's radius will not matter, so you might as well code the conversion as

``````(x,y,z) = (Cos[long] Cos[lat], Sin[long] Cos[lat], Sin[lat])
``````

Because you will need it later, the inverse of this is

``````(long, lat) =  (ArcTan[x, y], ArcTan[Sqrt[x^2 + y^2], z])
``````

A general (proper) rotation about the origin in three dimensions can be written in terms of three angular variables t, f, y in 3 by 3 matrix form as

``````Cos(t) Cos(f), -Cos(y) Sin(t) - Cos(t) Sin(f) Sin(y), -Cos(t) Cos(y) Sin(f) + Sin(t) Sin(y)
Cos(f) Sin(t),  Cos(t) Cos(y) - Sin(t) Sin(f) Sin(y), -Cos(y) Sin(t) Sin(f) - Cos(t) Sin(y)
Sin(f),         Cos(f) Sin(y),                         Cos(f) Cos(y)
``````

Call this Rot(t,f,y). To measure how closely Rot(t,f,y) carries the source S of a link to its target T (each expressed in 3D Cartesian coordinates), compute the squared distance between Rot(t,f,y)*S and T. (This is merely the inner product of the difference vector with itself). Sum all these values over all links: that is the objective function to minimize with respect to t, f, and y.

As an example, I generated ten points randomly within a small city-sized region on the earth:

``````{-90.8383, 39.9667}, {-90.9375, 39.9561}, {-90.8831, 39.9842}, {-90.9426, 39.9858}, {-90.7721, 39.7536}, {-90.7525, 39.888}, {-90.9364, 39.9648}, {-90.8006, 39.7983}, {-90.8211, 39.9425}, {-90.795, 39.8444}
``````

I rotated these by a small amount (specifically, t = -15', f = 12', y = 6') and then added normally-distributed errors with standard deviation of 0.5 seconds of arc to both the latitude and longitude. The new points were

``````{-91.2489, 39.859}, {-91.3621, 39.8643}, {-91.3092, 39.8738}, {-91.3679, 39.8896}, {-91.1981, 39.6642}, {-91.1653, 39.7812}, {-91.3533, 39.8602}, {-91.2204, 39.6971}, {-91.2299, 39.8325}, {-91.2015, 39.7403}
``````

Here is a picture (in a Plate Carree projection) of the resulting links, shown as arrows from source points (cyan) to target points (red):

You can see the effects of the errors: the target configuration is a slightly distorted version of the source configuration. Some of the link arrows even cross. (In practice, errors with modern data sources will be much smaller than this, but errors of this magnitude can easily occur when georeferencing historical maps.)

The resulting least-squares rotation matrix is

`````` 0.999932,   0.0108615,  0.00422911
-0.0108536,  0.999939,  -0.00188571
-0.00424934, 0.00183969, 0.999989
``````

The mean square error equals 0.57 seconds of arc, very close to the error built into these data.

Now the arrows show the action of this estimated rotation. The errors are evident in the slight displacements of the red (target) points from the images of the source points (in dark blue).

For such a small distance and such small regions, finding a rotation (as opposed to a translation in a projected map) is overkill. (In the preceding figure, notice how all the arrows appear to be parallel and of the same length: at this scale, a small spherical rotation is very close to a plane translation.) However, fitting a bona fide rotation will work just as well in situations where a map projection may introduce significant errors, such as when the links span the entire globe or the regions are large in extent. For small rotations, you can also read off the rotation parameters of the seven-point Helmert transformation from the result (they are found in the upper right triangle of the matrix). You can also consult textbooks on the theory of least squares (and maximum likelihood estimation) for ways to compute a confidence interval for the rotation matrix. (Providing that information here would take us too far afield and would be better suited for the statistics site.)

-
+1 Great explanation! I disagree though about the "too far afield". The question sounds very similar to that faced when creating Paleogeographic reconstructions which leverage Euler's rotation theorem. – Kirk Kuykendall Feb 8 '12 at 16:53
You're right, @Kirk: a discussion of errors is suitable for GIS, but I suspect most readers would find it too technical and distracting, especially because the question does not specifically request it. In many cases we have the data we have, we do the best we can, and let the chips fall, hoping that any errors in our georeferencing will become obvious through the overlays we create and the maps we draw. (It is always nice, though, to have some sense of the magnitude of possible georeferencing errors, if only to know where to look when layers don't line up right.) – whuber Feb 8 '12 at 17:08
Yeah, I was thinking more in terms of a general solution to title "how to translate coordinates on the earth". By generalizing the problem to "polygon" instead of cities, doesn't Euler suggest a solution could be stated with just 3 numbers: longitude of pole, latitude of pole, and angle of rotation? To me this is non-intuitive. Certainly things get messy quickly when dealing with error. But still, I think a lot of PaleoMaps could benefit from some sort of error analysis. – Kirk Kuykendall Feb 8 '12 at 20:01
There's another idea to add to the list of papers to write "real soon now" :-). – whuber Feb 8 '12 at 20:35