A confocal projection is a projection which maps confocal ellipses on the sphere onto confocal ellipses on the plane. Thus it takes as parameters the coordinates of the two focus points, or, equivalently, the center of the projection plus the direction and distance of one of the focus points from the center (either way, four scalar degrees of freedom).
Confocal projections can be seen as a generalization of azimuthal projections, which map concentric circles onto concentric circles. Indeed, in the limit case for a confocal projection as the two focus points approach the same point is an azimuthal projection.
The best-known confocal projection is the two-point equidistant projection (generalizing the azimuthal equidistant projection), whose confocal property is obvious from its definition.
Because they have the property of mapping all ellipses with the same focus points as ellipses, confocal projections will in general not have optimal distortion for any specific ellipse. Indeed, some fairly simple reasoning can show that confocal projections in general cannot be optimal. (Note that on the sphere, for any two focus points, one of the "ellipses" around those focus points is in fact a great circle. Since it is a circle, obviously the optimal way to project it is, for nearly any optimality metric you can name, an azimuthal projection onto a planar circle. But confocal projections would project it as an ellipse, not a circle, thereby proving suboptimality.)
Recently, I worked out the formula for the conformal confocal projection. (Yes, that's a bit of a tongue-twister, sorry. You could call it "two-point conformal" by analogy with the two-point equidistant, but then it's not like it's only conformal at two points.)
...Well, probably.
Poking around with a search engine, I found a formula in a few places for conformally mapping a planar circle onto a planar ellipse:
ellipse = sin(pi/2/K(s^2) * asn(circle/sqrt(s)), s^2)
Or in reverse direction:
circle = sqrt(s) * sn(asin(ellipse) / (pi/2/K(s^2)), s^2)
Here sn is a Jacobi elliptic function (defined in terms of m = k^2, not k as used in some other sources), asn is its inverse function, equivalent to one form of the incomplete elliptic integral of the first kind (but it's more commonly implemented in a slightly different form that satisfies F(asin(x),m) = asn(x,m), so using the asn terminology is less ambiguous), and K is the complete elliptic integral of the first kind, K(m) = asn(1,m), and for that matter, pi/2/K(m) = agm(1,sqrt(1-m)), where agm is the arithmetic-geometric mean, the most common way of implementing the elliptic integral in practice (so if you're computing it that way, it's actually easier for this formula to just never bother with the pi/2 constant at all!).
Then, of course, the problem becomes projecting a spherical ellipse onto a planar circle, so that you can then apply the above formula to take it to its final form. That might sound difficult, but conveniently, there's a trick: like I mentioned before, every great circle is also an ellipse! Since a confocal projection, per definition, preserves all ellipses with the same focus points, we might as well pick the one that's easiest to use and set the projection in terms of that. A great circle (or any other circle) can, of course, be projected onto a circle using the stereographic projection.
Next, one has to detemine the value of s. For this, we observe the obvious constraint that the focus points on the sphere must, in fact, be mapped onto the focus points (of the great circle as mapped above) on the plane. This turns out to be conveniently easy, and corresponds to a value of s = tan(focus/2)^2, where "focus" is half the distance between the focus points, or the distance between a focus point and the center of the projection. (On the plane, this is called linear eccentricity, but I'm not sure if the concept of "eccentricity" is really meaningful on the sphere.)
Thus, we obtain the following final formula for the projection:
oblated_stereographic = sin(pi/2/K(tan(focus/2)^4) * asn(stereographic/tan(focus/2), tan(focus/2)^4))
And its inverse projection, which is what I actually use:
stereographic = tan(focus/2) * sn(asin(oblated_stereographic) / (pi/2/K(tan(focus/2)^4)), tan(focus/2)^4)
Here "stereographic" is the complex number denoting the stereographic projection with scale 1/2 at the center, and with the focus points lying on the real number line.
Now, because of the general complication of the formula involved (or even of drawing an ellipse on a sphere), I haven't actually proven that this projection satisfies the confocal property. However, by the uniqueness theorem on conformal projections, it's easy to prove that if any conformal projection satisfies the confocal property, it must be this one.
In the given normalization for the projection, the projected focus points are always at +1 and -1 (regardless of the distance of the actual focus points on the sphere), and its projected linear scale at the center equals pi/4/tan(focus/2)/K(tan(focus/2)^4). Furthermore, the projected ellipses representing a hemisphere (viewed as an ellipse of semiminor axis pi/2=90°) and the full sphere (viewed as an ellipse of semiminor axis pi=180°) have particularly easy formulae in terms of only the complete rather than incomplete elliptic function of the first kind:
projected semimajor axis of hemisphere = cosh(pi/4/K(tan(focus/2)^4)*K(1-tan(focus/2)^4))
projected semiminor axis of hemisphere = sinh(pi/4/K(tan(focus/2)^4)*K(1-tan(focus/2)^4))
projected semimajor axis of sphere = cosh(pi/2/K(tan(focus/2)^4)*K(1-tan(focus/2)^4))
projected semiminor axis of sphere = sinh(pi/2/K(tan(focus/2)^4)*K(1-tan(focus/2)^4))
This projection may, or may not, in fact be the same as the Miller oblated stereographic projection, which also seems to do something with ellipses. It's hard to tell, since it doesn't get discussed much. Only a few places mention it, none that I've found give formulae, and many of the places that do mention it seem to use the name to refer to the one specific aspect favored by Miller (for projecting Europe and Africa without Asia) rather than a general parametrizable projection... but don't mention what the parameters for that aspect are, so I can't even compare the pictures to see if they look the same.
I managed to cobble up a working implementation of the Jacobi elliptic functions, and so here is an example of this projection, with the focus points being 147°W 66°N (in Alaska) and 54°W and 45°S (in the Atlantic Ocean off the coast of South America - I picked this value by eyeballing what seemed to produce a good-looking projection, rather than any theoretical motivation).
- conformal.png (687.6 KiB) Viewed 2651 times
- conformal60.png (409.39 KiB) Viewed 2651 times
Now you might recall from above that I mentioned that an optimal projection would, in at least some circumstances, take the form of a circle. Certainly this would be true for a hemispherical projection, and also for a full-spherical projection, where the latter would correspond to the modified Lagrange projections I described here. Given these examples, it's probably true at least for all elliptic regions larger than a hemisphere. It's probably not true for smaller elliptic regions (which are the ones you're more likely to be interested in), but it's hard to say where exactly the cutoff lies. Whenever the optimal projection is a circle, though, we can in theory compute it by using the above formula... to map an ellipse onto a circle! So we project the hemisphere onto a circle, map the circle onto an ellipse, and then map an ellipse of a different size back onto a circle. (If this is giving you flashbacks to the Lagrange, Hammer, or Aitoff projections' projection-resize-unprojection shenanigans, you're not the only one, although in this case we only do one actual sphere-to-plane projection and do the rest of the work in the plane.)
projection = sqrt(s) * sn(asn(stereographic/tan(focus/2), tan(focus/2)^4) / K(tan(focus/2)^4) * K(s^2), s^2)
stereographic = tan(focus/2) * sn(asn(projection/sqrt(s), s^2) / K(s^2) * K(tan(focus/2)^4), tan(focus/2)^4)
Where s is determined (probably numerically) by:
asn(tan(semimajor/2)/tan(focus/2), tan(focus/2)^4) / K(tan(focus/2)^4) = asn(1/sqrt(s), s^2) / K(s^2) = 1 - K(1-s^2)/K(s^2)*i/2
Unfortunately, I don't currently have an implementation of the asn function for complex numbers (I've worked out how to compute it for arbitrary real or imaginary numbers, but unfortunately, asn doesn't have a convenient addition theorem the way sn does, so this does not suffice to compute it for general complex numbers), so I can't test it. Even if I could, actually proving or disproving optimality would be difficult.
While I was at it, I also included an implementation of the confocal equidistant projection (better known as the two-point equidistant projection), shown here:
- equidistant.png (619.3 KiB) Viewed 2651 times
Because this forum only allows three images per post, you don't get a cropped version of the equidistant projection, but then it has far less area exaggeration than the conformal projection, so it doesn't really need it.
I wanted to also compute a confocal equal-area projection, but although I managed to derive some useful formulae, the actual solution proved beyond me and I eventually gave up.
One other thing of possible mild interest is that along the way, I worked out formulae for determining the center point between the two focus points. Ways of computing this have, of course, been around for a long, long time, but unlike the formulae I've been able to find anywhere else, the ones I came up with satisfy the condition of being completely symmetric in the two focus points. (Wikipedia's approach, for example, involves a roundabout way that first calculates a "node" point on the equator and then works from there.)
center_longitude = atan2(sin((focus2_longitude-focus1_longitude)/2) * (cos(focus2_latitude)-cos(focus1_latitude)), cos((focus2_longitude-focus1_longitude)/2) * (cos(focus2_latitude)+cos(focus1_latitude))) + (focus1_longitude+focus2_longitude)/2
center_latitude = asin((sin(focus1_latitude)+sin(focus2_latitude)) / (2*cos(focus_distance)))
bearing = atan2(cos(focus1_latitude) * cos(focus2_latitude) * sin(focus2_longitude-focus1_longitude), (sin(focus2_latitude)-sin(focus1_latitude)) * cos(focus_distance))
The downside of this approach, of course, being that it only computed the exact halfway point and not arbitrary fractions. If you think you can solve the latter case too, feel free to try, but this was hard enough for me. Another viable formula, but less favorable because it isn't longitude-neutral, is:
center_longitude = atan2(cos(focus1_latitude)*sin(focus1_longitude) + cos(focus2_latitude)*sin(focus2_longitude), cos(focus1_latitude)*cos(focus1_longitude) + cos(focus2_latitude)*cos(cos_longitude2))