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.
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
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
Because you will need it later, the inverse of this is
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
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:
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
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
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.)