[Proj] Optimal Albers Standard parallels

[strebe]

[Resending minus tables, since they exceeded the 40K limit for posting.]

On Feb 16, 2010, at 3:01:04 AM, "Jan Hartmann" <j.l.h.hartmann at uva.nl> wrote:
Do you have a programmatic version of the algorithm? It's something that would be really useful for a global map provider. I am working with the Royal Tropical Institute here in Amsterdam (www.kit.nl), on an index for their 10.000s of historical maps from all over the globe. In the long run we would like to offer  some thematic functionality, with an optimal projection
I have no program code (that I can find or recall), though I found the Maple V worksheet I constructed at the time, and from which I made the corrections noted below. I also found tables I calculated giving optimal standard parallels for pairs of bounding parallels.

On Feb 16, 2010, at 3:57:03 AM, Oscar van Vlijmen <ovv at hetnet.nl> wrote:
Very interesting!
But, could you please balance the parentheses in equation 4 (and 3), because 
equation 4 cannot be calculated.

Eek. Yes. Correcting, and using more uniform notation:

1) sin (phi[0]) = A +/- sqrt (A^2 - 1)
2) A = [sin (phi[N]) * cos^2 (phi[S]) - sin (phi[S]) * cos^2 (phi[N])] /
    [cos^2 (phi[S]) - cos^2 (phi[N])]
3) sin (phi[2]) = [sin (phi[1]) - 2 * sin (phi[0]) + sin (phi[1]) * sin^2 (phi[0])] /
    [2 * sin (phi[0]) * sin (phi[1]) - sin^2 (phi[0]) - 1]
4) cos^2 (phi[1]) * sqrt [1 + sin^2 (phi[0]) - 2 * sin (phi[0]) * sin (phi[S])] -
 cos (phi[S]) * [1 + sin^2 (phi[0]) - 2 * sin (phi[0]) * sin (phi[1])] = 0

Equation (4) is of the form

5) a * cos^2 (phi[1]) + b * sin (phi[1]) + c = 0

which has a closed-form solution, though it is a little messy. There

a = sqrt [1 + sin^2 (phi[0]) - 2 * sin (phi[0]) * sin (phi[S])]
b = 2 * cos (phi[S]) * sin (phi[0])
c = - cos (phi[S]) - cos (phi[S]) * sin^2 (phi[0])

with solution

6) phi[1] = Atan2 (-[a*x^2 + c]/b], x)

where

7) x = +/- sqrt (-4*a*c - 2*b^2 +/- 2*z) / (2*a)

and

8) z = sqrt (4*a*b^2*c + b^4 + 4*a^2*b^2)

The sign for (1) is chosen as the opposite of the average of the signs of the northern and southern bounding parallels. Using all combinations of signs in (7) yields four solutions. The two that are between phiS and phiN are the two of interest. (3) becomes superfluous in this scenario, since it yields the same as one of the combinations. I have supplied a few examples below. rn is the inset of the northern standard parallel as a ratio of the latitudinal span of the region of interest; rs is the southern.

Regards,
— daan Strebe


phi[S] = 20 phi[2] = 46.34188   rn = 0.12192
phi[N] = 50 phi[1] = 25.06720   rs = 0.16890

phi[S] = 25 phi[2] = 77.03289   rn = 0.05393
phi[N] = 80 phi[1] = 37.74930   rs = 0.23181

phi[S] = 30 phi[2] = 38.61324   rn = 0.13860
phi[N] = 40 phi[1] = 31.53993   rs = 0.15390

phi[S] = 50 phi[2] = 58.69310   rn = 0.13068
phi[N] = 60 phi[1] = 51.61978   rs = 0.16200


On Feb 15, 2010, at 7:55:38 PM, strebe <strebe at aol.com> wrote:

Jan & colleagues:

Regarding the question of optimal standard parallels for Albers, I found my notes from 8 years ago. To summarize, there is no good rule, so perhaps the exact calculation I describe below has more importance than just academic.

Regards,
— daan Strebe

________________________


The problem: given two latitudinal extents phi[N] and phi[S], find the optimal standard parallels for the Albers projection. Optimal is defined to mean that the maximum angular deformation between the two standard parallels equals the maximum angular deformation at the two latitudinal extents. In this way one is assured that the maximum deviation from true remains minimal throughout the region of interest.

In a breathtakingly brief treatment, Melluish touches on the issue by describing how to select standard parallels such that the "error" on the extreme parallels is equal in magnitude and opposite in sign to the error of the central parallel. The distortion measurement is equivalent to mine, but the central parallel (the average of the north and south extreme parallels) generally does not carry the greatest distortion, so his method does not yield an optimal map. Hinks discusses the problem precisely, even noting that the central parallel is not the parallel of greatest distortion... but his treatment applies to the equidistant conic, not the equal-area, contrary to Snyder's comment in USGS P.P 1395.

We define a family of Albers projections to be all Albers projections characterized by a common 'central' latitude at which angular deformation reaches a peak between the two standard parallels, or sinks to zero in the case of one standard parallel. I have placed 'central' between quote marks because it is not any simple average of the standard parallels. We can demonstrate that the ratio of the angular deformation of any two points on the map remains constant for every member of the family. This is the beautiful and useful characteristic of the family. The property arises out of a happy mathematical coincidence.

First I declare the member of the family with only one standard parallel to be the canonical form. Because of the family's property of invariant scale ratios, we can first reduce the problem of choosing an optimal Albers projection to the single standard parallel case, which is much more tractable than trying to deal with two free parameters. This gives us a projection with zero distortion along the canonical parallel and identical distortion at the northern and southern extremities. Once we have found the canonical projection, we can choose the member of its family that produces the same magnitude of distortion at the canonical parallel as at the extreme parallels. We are guaranteed that the projection that fulfills this condition is one from the family of the canonical projection.

Let phi[S] denote the southern parallel of interest; phi[N] the northern parallel of interest; phi[0] the canonical parallel; phi[1] a standard parallel; and phi[2] the other standard parallel.

phi[S] and phi[N] are given. From them we must determine phi[0], phi[1], and phi[2], such that the magnitude of the angular deformation at phi[0] is the same as the magnitude of the angular deformation at phi[N] and phi[S]. Furthermore, there must be no greater angular deformation registered anywhere between the phi[S] and phi[N]. Given those conditions:

1) sin (phi[0]) = A +/- sqrt (A^2 - 1)

where

2) A = (sin (phi[N]) * cos^2 (phi[S]) - sin (phi[S]) * cos^2 (phi[N])) / (cos^2 (phi[S] - cos^2 (phi[N]))

That establishes the canonical parallel of the canonical form, and therefore the family, wherein lies the solution.

If we arbitrarily select some phi[1], one may ask how to arrive at the correct phi[2] within the family, given phi[0].

3) sin (phi[2]) = [sin (phi[1]) - 2 sin (phi[N]) + sin (phi[1]) sin^2 (phi[N])] / [2 sin (phi[N]) * sin (phi[1]) - sin^2 (phi[N]) - 1]

However, we do not yet know the proper phi[1]. With much manipulation, we discover the iterative equation:

4) cos^2 (phi[1]) * sqrt (1 + sin^2 (phi[0]) - 2 sin (phi[0]) * sin (phi[S])) - cos (phi[S]) (1 + sin^2 (phi[0] - 2 sin (phi[0]) sin (phi[1]) = 0

We solve for phi[1] with (4) and then use (3) to compute phi[2]. Now have the optimal Albers for any circumstance.

The literature fixates quite a bit on rules of thumb for choosing the standard parallels, with some texts recommending inward 1/5th the meridional separation of the extreme parallels. Other texts recommend 1/7th, and yet others recommend 1/6th. My own calculations based on the above criteria show that the correct values range from 5% to 25% for reasonable conics, and the symmetry between north and south is not necessarily very good. Hence, no rule of thumb is very good. The higher the latitude of the northern parallel, the greater the asymmetry. The asymmetry remains remarkable even if both latitudes are high. For instance, extreme latitudes at 50° and 70° result in standard parallels placed optimally at roughly 10% below phi[N] and 18% above phi[S].

Deetz's and Adams's recommendation of 29°30' and 45°30' for the US are pretty good, of course. Given maximum extents of 24°34' and 49°25', the optimal standard parallels are 28°42' and 46°19'. If we just can't bring ourselves to treat Lake of the Woods equitably then the upper standard parallel drops to 45°56'. The lower parallel only drops by 5'. The former map yield a maximum east-west scale distortion of 1.19%; the latter yields 1.15% (excluding Lake of the Woods, of course). Deetz's and Adams's recommendation yields 0.98% in the heartland; 1.25% at the northern border; 1.45% at the northern tip of the Lake of the Woods boundary; and 1.37% at Key West.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20100216/2127ead0/attachment.html>