001 # -*- coding: utf-8 -*-
002
003 from math import floor, pi as PI, sin, log, pow, e as E, asin
004 from urllib import urlopen
005 import socket
006 import time
007
008
009 def dms2dd(data):
010     '''字符串的分度秒转为数字度'''
011     d = data.split(':')
012     dd = float(d[0]) + float(d[1])/60 + float(d[2])/3600
013     return dd
014
015 def dd2dms(data):
016     '''数字度转为字符串的分度秒'''
017     d = floor(data)
018     ms = (data - d)*60
019     m = floor(ms)
020     s = (ms - m)*60
021     return "%d:%d:%.3f"%(d,m,s)
022
023 def lngToPixel(lng, zoom):
024     '''经度到像素X值'''
025     return (lng + 180) * (256 << zoom) / 360
026
027 def pixelToLng(pixelX, zoom):
028     '''像素X到经度'''
029     return pixelX * 360 / (256 << zoom) - 180
030
031 def latToPixel(lat, zoom):
032     '''纬度到像素Y'''
033     siny = sin(lat * PI / 180)
034     y = log((1 + siny) / (1 - siny))
035     return (128 << zoom) * (1 - y / (2 * PI))
036
037 def pixelToLat(pixelY, zoom):
038     '''像素Y到纬度'''
039     y = 2 * PI * (1 - pixelY / (128 << zoom))
040     z = pow(E, y)
041     siny = (z - 1) / (z + 1)
042     return asin(siny) * 180 / PI
043
044 def calibration(rawll, bias, zoom):
045     '''将偏移量bias加到原始经纬度上，得到修正后的经纬度
046     参见：http://blog.csdn.net/dongmeng110/archive/2009/12/13/4997634.aspx'''
047
048     x = lngToPixel(rawll[1], zoom)
049     y = latToPixel(rawll[0], zoom)
050     factor = 18 - zoom
051
052     xx = x - floor(bias[0] / pow(2, factor))
053     yy = y - floor(bias[1] / pow(2, factor))
054
055     return (pixelToLat(yy, zoom), pixelToLng(xx, zoom))
056
057 def getOffset(rawll):
058     '''通过Google服务器得到指定经纬度的偏移量值(只给出最高的18级的偏移量，其他可以简单除2获得)
059     参见：http://blog.csdn.net/dongmeng110/archive/2009/10/31/4750380.aspx'''
060     lat, lng = rawll
061     url = "http://ditu.google.com/maps/vp?spn=0.0,0.0&z=18&vp=%.6f,%.6f" % (lat, lng)
062 #    print url
063     webpage = None
064     retry = 0
065     while not webpage and retry < 5:
066         try:
067             webpage = urlopen(url)
068         except Exception, e:
069             # 处理连接重置
070             print e
071             retry += 1
072             time.sleep(5)
073     if retry == 5:
074         return (0, 0)
075     text = webpage.read()
076 #    print text
077     data = text[text.rfind('[')+1:text.rfind(']')]
078     if not data:
079         return None
080     offset_data = data.split(',')[-2:]
081 #    print offset_data
082     offset = [int(x) for x in offset_data]
083     print offset
084     return offset
085
086
087 def handleRawllFiles(filename, zoom, lat_col, lng_col, has_ret=False):
088     '''校准指定文件中的经纬度值，修正后保存起来。
089     filename是包含原始经纬度值的文件的名字，文件中一行数据就是一个记录，
090     每个记录的不同列用空格分开。
091     zoom是谷歌地图的缩放率，最高18。
092     lat_col和lng_col分别是文件中纬度和经度所在的列号，从0开始算起。
093     has_ret为True返回修正后的经纬度值元组，否则不返回，默认是不返回。
094     '''
095    
096     # timeout in seconds
097     timeout = 10
098     socket.setdefaulttimeout(timeout)
099
100     f = open(filename, 'r')
101     lines = f.readlines()
102     f.close()
103
104     result = []
105     fixedll = []
106     for line in lines:
107         data = line.split(' ')
108         rawll = (float(data[lat_col]), float(data[lng_col]))
109         print rawll
110         bias = getOffset(rawll)
111         if not bias:
112             # 没有获得偏倚值，那就用原来的值啦～
113             result.append(line)
114             fixedll.append(rawll)
115         else:
116             fixll = calibration(rawll, bias, zoom)
117             data[lat_col] = str(fixll[0])
118             data[lng_col] = str(fixll[1])
119             result.append(' '.join(data))
120             fixedll.append(fixll)
121         print result[-1]
122
123     f = open(filename+'.fix', 'w')
124     f.writelines(result)
125     f.close()
126    
127     if has_ret:
128         return fixedll
129
130 def handleMaps(filename, zoom):
131     '''校准指定文件中每个小地图的中心点经纬度值，修正后保存起来'''
132     handleRawllFiles(filename, zoom, 1, 2)
133
134
135 def handleWaypoints(filename, zoom):
136     '''校准指定文件中每个路点的经纬度值，修正后保存起来'''
137     handleRawllFiles(filename, zoom, 1, 2)
138
139 def handleRoute(filename, zoom):
140     '''生成gpsDrive中使用的gpx路径文件'''
141     fixedll = handleRawllFiles(filename, zoom, 1, 2, has_ret=True)
142
143     rtepts = []
144     prefix = '''<?xml version="1.0" encoding="UTF-8"?>
145 <gpx
146   version="1.0"
147   creator="GPSBabel - http://www.gpsbabel.org"
148   xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
149   xmlns="http://www.topografix.com/GPX/1/0"
150   xsi:schemaLocation="http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd">
151 <time>2010-01-09T01:59:52Z</time>
152 <bounds minlat="32.078910000" minlon="118.718760000" maxlat="37.521580000" maxlon="121.561750000"/>
153 <rte>
154 '''
155     suffix = '''</rte>
156 </gpx>'''
157     template = '''<rtept lat="%.9f" lon="%.9f">
158   <ele>0.000000</ele>
159 </rtept>
160 '''
161     for ll in fixedll:
162         rtepts.append(template%ll)
163
164     f = open(filename+'.gpx', 'w')
165     f.write(prefix)
166     f.writelines(rtepts)
167     f.write(suffix)
168     f.close()