[Proj] Computing the Lambert conformal conic projection accurately
OvV_HN
ovv at hetnet.nl
Sun Nov 21 06:24:01 EST 2010
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.....
More information about the Proj
mailing list