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):
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)
--
http://mail.python.org/mailman/listinfo/python-list
