How to do stereographic projection conversion (w.r.t. changing point of projection)?

Question

Given a stereographic point P1 projected from point P, and given another point Q, I would like to know how to transform P1's coordinates to that of a stereographic point Q1 projected from point Q.

The Q I'm using is pretty close to P, and the points on the sphere/earth to be projected are pretty near the antipode of P/Q, possibly clustered together, and having a high accuracy on each point in order to differentiate them is important. I'm wondering if a simple adding of an offset value to P1 will provide a good approximation? I need to preserve conformality and relative distances of projected points (since distance is not preserved in a stereographic projection, but I think that relative distances of points are preserved?). P's and Q's latitude and longitude will be known.

Do let me know if this question is better posted at math.stackexchange.com

Answer 1

You appear to be describing rotations of the Riemann Sphere: they are conformal and preserve the spherical metric. There are several simple ways to write them down (not involving any trigonometry!). A good one is to consider coordinates (x,y) in the plane as complex numbers z = x + yI, I^2 = -1. Given any four complex numbers a, b, c, and d, a fractional linear transformation of z is the value (az + b)/(cz + d). Thinking of z as the image of a point on the sphere via stereographic projection, it's easy to work out that for a fractional linear transformation to be a true rotation, a and d must be complex conjugates of each other and b and -c must also be mutual complex conjugates. (This defines the matrix group SU(2,C).)

Let Q be the point (in projected coordinates) you would like to base a new projection on. All we need to do is rotate Q to infinity, because infinity corresponds to the projection's base point (the north pole). This implies the denominator cz + d must equal zero, giving a solution

a = Conjugate of Q,
b = 1,
c = -1,
d = Q.

In other words, the change of coordinates you seek has the formula

z --> (1 + Conjugate(Q)*z) / (Q - z).

As an example, suppose you want to make Q = (5, 5) = 5 + 5I the new origin of projection. Using the rules of complex arithmetic we can work out the transformation in terms of the coordinates (x,y):

z --> (1 + (5 - 5I)*z) / (5 + 5I - z)

 = (1 + (5 - 5I)*(x + yI)) / (5 + 5I - (x + yI))

 = (1 + 5x + 5y + (5y - 5x)I) / (5 - x + (5-y)I)

 = [(5 - x + 50y - 5x^2 - 5y^2) + (-5 - 50x + y + 5x^2 + 5y^2)I]
    / ((5-x)^2 + (5-y)^2).

That is, writing the new coordinates of (x,y) as (x',y'), we have

x' = (5 - x + 50y - 5x^2 - 5y^2) / ((5-x)^2 + (5-y)^2)

and

y' = (-5 - 50x + y + 5x^2 + 5y^2) / ((5-x)^2 + (5-y)^2).

You can follow this up by any rotation around the origin to reorient the new map. A good choice is to multiply the result by Q/Conjugate(Q) = Q^2 / |Q|^2. This is because the transformation has a nice interpretation:

(1 + Conjugate(Q)*z) / (Q - z) * Q / Conjugate(Q)

= (1/Conjugate(Q) + z) (1 + z/Q + (z/Q)^2 + ...).

The first factor merely translates all points by the (small) amount 1/Conjugate(Q). In particular, this keeps the map oriented correctly. The second factor can be ignored when z/Q is very small. Typically, for small rotations, Q is already "near infinity": that is, it is large compared with any point on the map, because it is the projection a point close to the original north pole. This justifies approximating the change in projection by means of a translation and, in addition, it tells us how to measure the error: the size of z/Q (the first, and largest, neglected term in the series on the right) is the ratio of the sizes of z and Q (that is, the ratio of their distances from the map origin). In other words, when the original projection of Q is way beyond the extent of your map, you will likely be ok with the approximate formula.