[Proj] CEA from an ellipse and +R_A

[David James Brockley]

Hi Gerald,

Yes, that is what I would have expected, but:

msslm4:/Volumes/Space/Slope2 djb$ /usr/local/bin/proj +proj=cea  
+ellps=WGS84
0 10
0.00    1100285.57

so y=1100285.57

from y = r * q / 2

Calculating q for 10 degrees gives

q = 0.345018

rearranging to give r

r = 2 y / q

solving

r = 6378134m

so the calculation has used the Authalic latitude, but not the  
Authalic radius of 6371007m. If we configure proj with +R_A to get  
the authalic radius, we use the spherical equations rather than the  
ellipsoidal - so no conversion to Authalic latitude.

I guess that since it is a linear scaling, and since the calculation  
is already using q / 2 rather than q / q_p which results in another  
linear scaling, it doesn't really matter but it certainly confused me  
while trying to compare the results of proj against the derivation of  
the equations in the books. I'd coded up the equations from the book  
and was trying to use proj as a method of verifying my implementation  
to be correct.

Thanks for your help. I've posted this to the list so that it gets  
archived and hopefully helps anyone suffering the same confusion as  
me in the future.

Regards, David

On 14 Jul 2008, at 18:25, Gerald I. Evenden wrote:

> On Monday 14 July 2008 5:09:03 am you wrote:
>> Hi Gerald
>>
>> 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 far as I can figure out all you need to do is:
>
> lproj +proj=cea +ellps=WGS84
>
> the x-y coming out are determined with authalic latitude.  cea  
> operates in
> both spherical and elliptical mode.
>
>> 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 + ... )
>>
>> and
>>
>> 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.
>>
>> Brock
>>
>> On 11 Jul 2008, at 17:56, Gerald I. Evenden wrote:
>>> On Friday 11 July 2008 5:57:49 am David James Brockley wrote:
>>>> Hello
>>>>
>>>> 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
>>>> WGS84.
>>>>
>>>> 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
>>> elliptic
>>> elements.
>>>
>>>> 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
>>>> http://lists.maptools.org/mailman/listinfo/proj
>>>
>>> --
>>> 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
>>> http://lists.maptools.org/mailman/listinfo/proj
>
>
>
> -- 
> 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