[Proj] Optimal Albers Standard Parallels

[strebe]

Colleagues:

I received an inquiry off list about the posting below from February 2010, so I will update the thread with recent developments.

I have since stumbled onto N. Tissot’s work on this problem* and found that his solution yields the same results as mine. He was a lot smarter than I am, though, so his formulation is much more compact. He rolled the constraints into the Albers formula itself, a move that collapses away much of the calculation.

phiN = northern limit of region of interest
phiS = southern limit of region of interest

      Initialization:
	double		delta (0.5 * (phiN - phiS));
	double 		phi0 (0.5 * (phiN + phiS));
	double		cosDelta = cos (delta);

	n = sin (phi0);
	rN = 1.0 / n;
	C = 0.5 * (n/cosDelta + cosDelta*rN);

     Per point:
	double		rho = sqrt (2.0 * rN * (C - sin (phi)));
	double		theta = n * lambda;

	x = rho * sin (theta);
	y = -rho * cos (theta);


Regards,
— daan Strebe

* 1881, Mémoire sur la représentation des surfaces et les projections des cartes géographiques: Paris, Gauthier Villars, 337 p. + 60 p. tables.

______

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

On Feb 16, 2010, at 3:01:04 AM, "Jan Hartmann" <[hidden email]> 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


-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/proj/attachments/20110624/35580af8/attachment.htm