作者StubbornLin (Victor)
看板Geography
标题[问题] 为什麽lat/lon转TWD97的lon会有十几公尺的误差
时间Thu Aug 13 17:09:06 2009
我自己照着
http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm
写的公式写了一个从lat/lon转TWD97的python程式
from math import tan, sin, cos, radians
class LatLonToTWD97(object):
"""This object provide method for converting lat/lon coordinate to TWD97
coordinate
the formula reference to
http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm (there is lots of
typo)
http://www.offshorediver.com/software/utm/Converting UTM to Latitude and
Longitude.doc
Parameters reference to
http://rskl.geog.ntu.edu.tw/team/gis/doc/ArcGIS/WGS84%20and%20TM2.htm
http://blog.minstrel.idv.tw/2004/06/taiwan-datum-parameter.html
"""
def __init__(self,
):
# Equatorial radius
self.a = 6378137.0
# Polar radius
self.b = 6356752.314245
# central meridian of zone
self.long0 = radians(121)
# scale along long0
self.k0 = 0.9999
# eccentricity of the earth's elliptical cross-section
self.e = (1-self.b**2/self.a**2)**0.5
# delta x in meter
self.dx = 250000
def convert(self, lat, lon):
"""Convert lat lon to twd97
"""
a = self.a
b = self.b
long0 = self.long0
k0 = self.k0
e = self.e
dx = self.dx
e2 = e**2/(1-e**2)
n = (a-b)/(a+b)
nu = a/(1-(e**2)*(sin(lat)**2))**0.5
p = lon-long0
A = a*(1 - n + (5/4.0)*(n**2 - n**3) + (81/64.0)*(n**4 - n**5))
B = (3*a*n/2.0)*(1 - n + (7/8.0)*(n**2 - n**3) + (55/64.0)*(n**4 -
n**5))
C = (15*a*(n**2)/16.0)*(1 - n + (3/4.0)*(n**2 - n**3))
D = (35*a*(n**3)/48.0)*(1 - n + (11/16.0)*(n**2 - n**3))
E = (315*a*(n**4)/51.0)*(1 - n)
S = A*lat - B*sin(2*lat) + C*sin(4*lat) - D*sin(6*lat) + E*sin(8*lat)
K1 = S*k0
K2 = k0**nu*sin(2*lat)/4.0
K3 = (k0*nu*sin(lat)*(cos(lat)**3)/24.0) * \
(5 - tan(lat)**2 + 9*e2*(cos(lat)**2) + 4*(e2**2)*(cos(lat)**4))
y = K1 + K2*(p**2) + K3*(p**4)
K4 = k0*nu*cos(lat)
K5 = (k0*nu*(cos(lat)**3)/6.0) * \
(1 - tan(lat)**2 + e2*(cos(lat)**2))
x = K4*p + K5*(p**3) + self.dx
return x, y
公式本身我检查好几次了,应该是没有打错了,包括网页里打错的部份
像51被误打成512,然後an被打成tan也都有修正
可是不知道为什麽,算出来的结果lat误差都在公尺之内,而lon总是有几十公尺的误差
像是我用这里面的三角点资料来做测试
http://wiki.osgeo.org/index.php?title=Taiwan_datums/Test_points&uselang=zh-tw
挑个七星山当输入的话
25.17069939, 121.5534442
我的程式算出来的结果是
305787.782217, 2784684.75762
而里面写的资料是
305787.783, 2784799.355
可以发现lat几乎都很精准
而lon却差了-114.59738公尺
不只是那三角点的资料,我也有找其它转换程式去算
结果也都是lat误差在公尺以下,而long误差在十几公尺到上百公尺
请问我可能是哪里弄错了? 参数弄错? 还是?
为什麽会有这long误差这麽大的问题?
因为这要拿来当做导航用的,long误差百公尺可以从路面跑进河里或任何地方 囧
以上,谢谢
--
哇咧咧 创意投票系统
http://walele.com
易记学 程式设计教学
http://ez2learn.com/
易记学 程式设计讨论区
http://forum.ez2learn.com
VICTOR's 个人Blog
http://blog.ez2learn.com/
财报分析王
http://victorlin.serveftp.org/stock/
--
※ 发信站: 批踢踢实业坊(ptt.cc)
◆ From: 218.165.225.76
1F:→ nehex:E = (315*a*(n**4)/512.0)*(1 - n) 08/13 18:05
2F:→ nehex: ^^^ 08/13 18:06
3F:→ nehex:另 geog.ntu 那个网页也不建议使用 @@ 08/13 18:13
4F:→ nehex:噢噢 doc 的 反而又错了 :P 08/13 18:17
5F:→ nehex:llh tm 互转误差是可以无视的... @@ 08/13 18:23
6F:→ nehex:= = 不管是 512 还是 51 都不影响 08/13 18:50