Python trig precision problem

Cary West cdwest1 at tva.gov
Thu May 18 23:43:06 CEST 2006

```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/LambertConformalConicProjection.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_parallel_2):

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 ) ) )

#######################################

lat,lon=35,-90

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

```