Python trig precision problem

C

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

casevh

Cary said:
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.

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
 
R

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.
 
D

Dennis Lee Bieber

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.
Well... ONE) we don't see the code used by that conversion tool.

SECOND) That tool is specified to be using an Australian datum,
maybe you are expecting WGS-84, and NAD-27 -- since your coordinates are
those of Memphis TN?
--
Wulfraed Dennis Lee Bieber KD6MOG
(e-mail address removed) (e-mail address removed)
HTTP://wlfraed.home.netcom.com/
(Bestiaria Support Staff: (e-mail address removed))
HTTP://www.bestiaria.com/
 
R

Robert Kern

Cary said:
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.

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/jeffrey.s.whitaker/python/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
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Similar Threads

Trig value errors 57
Roots Module 10
atan2 weirdness 3
problem select object in pyode 0
help to convert this fortan program to c++ 2
OpenGL 3D gears demo for Ruby 1
numpy help 2
[ANN] rcalc 2.0 (Ruby Calculator) 11

Members online

No members online now.

Forum statistics

Threads
473,744
Messages
2,569,482
Members
44,901
Latest member
Noble71S45

Latest Threads

Top