Python trig precision problem

[Cary West]

Hello all, having a small problem with a trig routine using python math
module. This code is designed to convert geodetic coordinates to lambert
conformal conic coordinates. It implements the formulas found at
http://mathworld.wolfram.com/Lambert...rojection.html . The
problem is that, at least on my machine, the precision is off to the tune of
around 1 km, which is unacceptable. (using the calculator found at
http://www.deh.gov.au/erin/tools/geo2lam-gda.html ) Anyone have any ideas?
here is the offending code.

from math import *

def lambert(lat,lon,ref_lat,ref_lon,st_parallel_1,st_p arallel_2):

earth_radius= 6371.9986
lon*=pi/180.0
lat *=pi/180.0

ref_lon*=pi/180.0
ref_lat *=pi/180.0

st_parallel_1 *=pi/180.0
st_parallel_2 *=pi/180.0

def sec(theta):
return 1.0/cos(theta)

def cot(theta):
return 1.0/tan(theta)

n = log( cos( st_parallel_1 ) * sec( st_parallel_2 ) ) / ( log( tan ( pi
/ 4.0 + st_parallel_2 / 2.0 ) * cot( pi / 4.0 + st_parallel_1 / 2.0 ) ) )

F = ( cos( st_parallel_1 ) * tan ( pi / 4.0 + st_parallel_1 / 2.0 ) **
n ) / n

rho = F * ( cot ( pi / 4.0 + lat / 2.0 ) ** n )

ref_rho = F * ( cot ( pi / 4.0 + ref_lat / 2.0 ) ** n )

x = rho * sin( n * ( lon-ref_lon ) )
y = ref_rho - ( rho * cos( n * ( lon-ref_lon ) ) )

return earth_radius*x,earth_radius*y
#######################################

lat,lon=35,-90

print lambert(lat,lon,40,-97,33,45)

[casevh]

Cary West wrote:

Hello all, having a small problem with a trig routine using python math
module. This code is designed to convert geodetic coordinates to lambert
conformal conic coordinates. It implements the formulas found at
http://mathworld.wolfram.com/Lambert...rojection.html . The
problem is that, at least on my machine, the precision is off to the tune of
around 1 km, which is unacceptable. (using the calculator found at
http://www.deh.gov.au/erin/tools/geo2lam-gda.html ) Anyone have any ideas?
here is the offending code.

I'd assume you are exceeding the precision available with standard
64-bit floating point. Here is a library for extended precision
floating point, including trig functions.

http://calcrpnpy.sourceforge.net/clnumManual.html

casevh

[Roger Miller]

If I were you I would see if I could get the Perl script referred to on
the ERIN web page. You might find that the discrepancy is something as
simple as a slightly different value for the Earth's radius. And by the
way, math.radians() might be a bit clearer than the pi/180 business.

[Robert Kern]

Cary West wrote:

Hello all, having a small problem with a trig routine using python math
module. This code is designed to convert geodetic coordinates to lambert
conformal conic coordinates. It implements the formulas found at
http://mathworld.wolfram.com/Lambert...rojection.html . The
problem is that, at least on my machine, the precision is off to the tune of
around 1 km, which is unacceptable. (using the calculator found at
http://www.deh.gov.au/erin/tools/geo2lam-gda.html ) Anyone have any ideas?
here is the offending code.

If you want to do map projections, I highly recommend using PyProj instead of
implementing the calculations yourself. Geodesy is nontrivial, and MathWorld
doesn't tell you half the details you need to take care of.

http://www.cdc.noaa.gov/people/jeffr...on/pyproj.html

--
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
that is made terrible by our own mad attempt to interpret it as though it had
an underlying truth."
-- Umberto Eco