[Proj] Computing the Lambert conformal conic projection accurately

[OvV_HN]

There is probably a more serious problem in one of the other projections, 
namely the Bonne projection.
If lat_0 is close to 0, the Sinuoidal projection is taken.
If one compares in a limiting case the northing, produced by the Sinusoidal, 
and for the same set of parameters the northing, produced by the Bonne 
projection, a difference of perhaps 300 meters will be obtained. This 
observation is valid for both the ellipsoidal and the spherical projections.
This difference can probably not be attributed to instabilities in the 
calculation or a difference between a limiting lat_0 (of 1e-10 deg) and 
zero.
Can somebody enlighten me? Does the Sinusoidal projection differ a bit too 
much from the Bonne (but no one cares)?
Is there something wrong in the proj/libproj code?
Or did I do something wrong?

My test case:

Bessel 1841 ellipsoid
lat=40d 32m; lon=-7d 16m;
lat0=1e-10; lon0=-8.131906111111112;
x0=20000; y0=10000; [meters]
bonne ellipsoidal projection
Result: x = 93299.284; y = 4497999.829; [meters]

for lat0=0 the Sinusoidal projection is taken.
Result: x = 93299.284; y = 4498299.575;
Difference in y: 299.75 m

Oscar van Vlijmen


IN REPLY TO:

From: Charles Karney <ckarney <at> sarnoff.com>
Subject: Computing the Lambert conformal conic projection accurately
Newsgroups: gmane.comp.gis.proj-4.devel
Date: 2010-11-19 22:27:50 GMT (1 day, 12 hours and 30 minutes ago)

I've posted GeographicLib 1.5 on SourceForge.  See

  http://sf.net/projects/geographiclib/files/distrib/

This includes a couple of features of note:

(a) A Matlab interface to the UTMUPS, MGRS, Geoid, Geodesic classes has
    been included.

(b) I've reformulated the Lambert conical projection to improve its
    numerical accuracy.

To follow up on point (b)...  In the 2 parallel case, the Lambert
projection depends on a parameter n, given by (see Snyder, 1987, p. 108)

  n = (log(m1) - log(m2)) / (log(t1) - log(t2))
    = (log(sec(beta2)) - log(sec(beta1))) /
      (asinh(tan(chi2)) - asinh(tan(chi1)))

where {phi,beta,chi}{1,2} are the {geographic,reduced,conformal}
latitudes of the standard parallels.  Obviously, this expression is
inaccurate if the standard parallels are close to one another.  The
normal way of dealing with this (proj4, geotrans) is to substitute the
limiting expression

etc.....