[Proj] Optimal Albers Standard Parallels
strebe
strebe at aol.com
Fri Jun 24 15:36:00 EST 2011
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
More information about the Proj
mailing list