[Proj] wrong computation of meridinal distance / pj mlfn.c
Peter Albert
peter.albert at gmx.de
Tue Jul 12 11:24:10 EDT 2005
Dear all,
I am new to proj and also relatively new to projection in general;however,
I think I found a problem in how the meridinal distance is approximated in
proj (using the subroutine in pj_mlfn.c).
Actually, it seems to me as if the meridinal distance on a given ellipsoid
is not calculated from the equator to the given (geodetic) latitude but
between the Pole and (90° - the appropriate geocentric latitude)!
Let me give you an example: For the sake of simplicity, let the major axis
of our ellipsoid be a = 1.0 and the squared eccentricity es = 0.007.
Let the geodetic latitude be 10N, then the appropriate geocentric latitude
calculates to 9.931398465N.
(Following this transformation:
float(arctan((a-es*a)/a*tan(10*PI/180))*180/PI)
)
Now, if I understand the sinusoidal projection correctly, then the
y-coordinate of a given point in this projection should just equal its
meridinal distance from the equator, in our case:
echo "0 10" | proj +proj=sinu +a=1.0 +es=0.007 -f %.8f | cut -f 2
The result is 0.17332956, which generally makes sense given the major axis
being set to 1. (If I choose a realistic ellipsoid, the resulting value is
just 0.17332956 times a.)
Just to be sure, let's cross-check my calculation of geocentric latitude by
adding "+geoc" to the previous call and using the geocentric latitude
instead of 10N:
echo "0 9.931398465" | proj +geoc +proj=sinu +a=1.0 +es=0.007 -f %.8f |
cut -f 2
gives 0.17332956, no problems so far.
However, now I tried to calculate the elliptical integral myself, i.e.
int(sqrt(a-es*(cos(t)^2)) dt) from the equator up to the given geocentric
latitude. I did this using mupad, a computer algebra system:
numeric::int(sqrt(1.0-0.007*(cos(x))^2), x=0..9.931398465*PI/180)
The result now is 0.1727339231, which significantly differs from the proj
output.
However, if you integrate from 90° - the geocentric latitude up to the pole
like this:
numeric::int(sqrt(1.0-0.007*(cos(x))^2),
x=(90-9.931398465)*PI/180..90*PI/180)
then you get 0.1733295629, which perfectly matches proj's results.
I would greatly appreciate if someone could tell me whether I really got
something wrong here, or if there actually is a problem in proj. I looked
into the code, and the meridinal distance is unfortunately approximated
using a weird 8th order polynomial which I can't say anything about ...
Thanks in advance and kind regards,
Peter
P.S.: BTW, I am using proj-4.4.9
--
5 GB Mailbox, 50 FreeSMS http://www.gmx.net/de/go/promail
+++ GMX - die erste Adresse für Mail, Message, More +++
More information about the Proj
mailing list