[Proj] Computing the Lambert conformal conic projection accurately

Charles Karney ckarney at sarnoff.com
Fri Nov 19 17:27:50 EST 2010

I've posted GeographicLib 1.5 on SourceForge.  See


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

  n = sin(phi1)

if abs(phi2 - phi1) is "sufficiently" small.  However this has the
disadvantage that there's a range of values of (phi2 - phi1) where
neither the limiting form nor the "exact" formula gives accurate results

However, there's a systematic method for dealing with numerical problems
like this, namely the calculus of finite differences.  This method
transforms the expression for n to an equivalent one which can be
evaluated accurately for any phi2 != phi1.  The expression is too
lengthy to include here.  However, I've included fairly extensive
comments within the code for the LambertConformalConic class.  The
technique can also be applied to the projection itself and the code then
naturally deals with the limiting cases of the Mercator and polar
stereographic projections.

A useful resource for divided differences is

  W. M. Kahan and R. J. Fateman,
  Symbolic computation of divided differences,
  SIGSAM Bull. 33(3), 7-28 (1999)

(Kahan, incidentally, the "father" of the IEEE floating point standard.)

Charles Karney <ckarney at sarnoff.com>
Sarnoff Corporation, Princeton, NJ 08543-5300

Tel: +1 609 734 2312
Fax: +1 609 734 2662

More information about the Proj mailing list