Geography 板


LINE

我自己照着 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







like.gif 您可能会有兴趣的文章
icon.png[问题/行为] 猫晚上进房间会不会有憋尿问题
icon.pngRe: [闲聊] 选了错误的女孩成为魔法少女 XDDDDDDDDDD
icon.png[正妹] 瑞典 一张
icon.png[心得] EMS高领长版毛衣.墨小楼MC1002
icon.png[分享] 丹龙隔热纸GE55+33+22
icon.png[问题] 清洗洗衣机
icon.png[寻物] 窗台下的空间
icon.png[闲聊] 双极の女神1 木魔爵
icon.png[售车] 新竹 1997 march 1297cc 白色 四门
icon.png[讨论] 能从照片感受到摄影者心情吗
icon.png[狂贺] 贺贺贺贺 贺!岛村卯月!总选举NO.1
icon.png[难过] 羡慕白皮肤的女生
icon.png阅读文章
icon.png[黑特]
icon.png[问题] SBK S1安装於安全帽位置
icon.png[分享] 旧woo100绝版开箱!!
icon.pngRe: [无言] 关於小包卫生纸
icon.png[开箱] E5-2683V3 RX480Strix 快睿C1 简单测试
icon.png[心得] 苍の海贼龙 地狱 执行者16PT
icon.png[售车] 1999年Virage iO 1.8EXi
icon.png[心得] 挑战33 LV10 狮子座pt solo
icon.png[闲聊] 手把手教你不被桶之新手主购教学
icon.png[分享] Civic Type R 量产版官方照无预警流出
icon.png[售车] Golf 4 2.0 银色 自排
icon.png[出售] Graco提篮汽座(有底座)2000元诚可议
icon.png[问题] 请问补牙材质掉了还能再补吗?(台中半年内
icon.png[问题] 44th 单曲 生写竟然都给重复的啊啊!
icon.png[心得] 华南红卡/icash 核卡
icon.png[问题] 拔牙矫正这样正常吗
icon.png[赠送] 老莫高业 初业 102年版
icon.png[情报] 三大行动支付 本季掀战火
icon.png[宝宝] 博客来Amos水蜡笔5/1特价五折
icon.pngRe: [心得] 新鲜人一些面试分享
icon.png[心得] 苍の海贼龙 地狱 麒麟25PT
icon.pngRe: [闲聊] (君の名は。雷慎入) 君名二创漫画翻译
icon.pngRe: [闲聊] OGN中场影片:失踪人口局 (英文字幕)
icon.png[问题] 台湾大哥大4G讯号差
icon.png[出售] [全国]全新千寻侘草LED灯, 水草

请输入看板名称,例如:Soft_Job站内搜寻

TOP