[Proj] CEA from an ellipse and +R_A
David James Brockley
djb at mssl.ucl.ac.uk
Mon Jul 14 05:09:03 EDT 2008
Thanks for the frank reply. I'm new at map projections!
What I want to do is : call pj_fwd() with a pair of geodetic
coordinates on the WGS84 ellipsoid, have the latitude transformed to
authalic latitude, and get XY in cylindrical equal area with a radius
equal to the authalic radius for WGS84. As far as I can see from the
books, this is the normal way to convert from WGS84 lon/lat to CEA x/
y. Please correct me if I'm hopelessly confused already by this stage.
As I understand it, the Authalic sphere is an intermediate step
between a reference ellipsoid (in my case WGS84) and a map projection
(cylindrical equal area). It allows one to project from the ellipsoid
and causes the result to still be equal area, but you have to replace
geodetic latitude with authalic latitude.
From Map Projections: Theory and Application p102
R^2_A = a^2(1-e^2)( 1 + 2e^2/3 + 3e^4/5 + ... )
sin phi_a = sin phi ( ( 1 + 2e^2 sin phi / 3 + ... ) / ( 1 + 2e^2/3
+ ... )
So, given a reference ellipsoid with known a, proj works out R_A -
this works and is as I expected.
What I would now expect is for proj to take the input Geodetic
latitude and convert it internally into an Authalic latitude before
doing the coordinate transform i.e.
y = R_A * sin( phi_a ) / k
but what actually happens (tracing with GDB) is
y = R_A * sin( phi ) / k
This is when I initialise with +R_A
Is the convention that when using +R_A you have to do the conversion
to Authalic latitude outside of proj?
I can see from PJ_cea.c that for the spheroidal case, pj_qsfn() and
pj_authlat() are used to handle the conversion to and from authalic
latitude, but in this case proj is using WGS84 a for the radius, not
R_A. (If I initialise with +proj=cea +long_0=0 +ellps=WGS84)
When I initialise with +R_A, I get the radius I want, but the
projection switches to using the spherical code so the transformation
to authalic latitude doesn't happen.
Hope this helps explain my confusion. Thanks for your help.
On 11 Jul 2008, at 17:56, Gerald I. Evenden wrote:
> On Friday 11 July 2008 5:57:49 am David James Brockley wrote:
>> I'm trying to convert geodetic lon/lat into CEA xy using the proj API
>> initialised with
>> +proj=cea +long_0=0 +ellps=WGS84 +R_A
>> What I expect to happen is
>> y = a * q / 2
>> where a is the authalic radius calculated from the semi-major of
>> q/2 seems to be a common approximation to q/q_p which is sin
>> ( authalic latitude)
>> What actually happens is that proj uses the equation for a sphere
>> rather than an ellipse.
> That is exactly what is supposed to happen as implied with the
> capital R in
> option's acronym, The R in this case is based upon a sphere of
> surface area
> equal to the ellipse selected. *All* of rhe R_* option series
> create radii
> to be used with the spherical earth projection formulas.
> Quite frankly, I do not understand what you are trying to do and your
> operations do not make sense to me. Authalic projections are based
> upon the
> math of the transformation and not on simple manipulation of the
>> y = a sin( phi )
>> I can get proj to do what I expect by initialising with
>> +proj=cea +long_0=0 +a=6371007 +e=0.0818191
>> but that seems somewhat back to front.
>> Can anyone tell me the correct way to configure proj for what I am
>> trying to do?
>> I'm fairly sure that my expected function for y is correct, since
>> I've both copied it from the literature and derived it by hand from
>> the spherical case by substituting authalic radius for radius and
>> authalic latitude for latitude - I've also realised that the /2 is
>> due to the maximum possible value of q_p being 2 when e=0, and that
>> it doesn't affect the properties of the transform (still equal area).
>> Regards, Brock
>> Proj mailing list
>> Proj at lists.maptools.org
> The whole religious complexion of the modern world is due
> to the absence from Jerusalem of a lunatic asylum.
> -- Havelock Ellis (1859-1939) British psychologist
> Proj mailing list
> Proj at lists.maptools.org
More information about the Proj