Before:
1. 研究的需要， 在 google map 上爬取了一些的靜態衛星地圖圖片，每張圖片的像素為 256*256
2. 通過 photshop 將這些地圖碎片手動拼成了地圖， 地圖只是覆蓋了學校而已， 還是比較小的
3. 我通過手機采集到一些 GPS 的 trace， 希望在這個離線地圖上畫上這些 GPS 點
4. 最初， 我以為地圖上的經緯度都是等距分布的， 比如一張圖片的頂部緯度是20， 底部緯度為22， 那么圖片中間的橫線緯度是 21. 實際上， 並非如此， 中間橫線的緯度應當是小於 21 的， 想了解具體的原因， 可搜索墨卡托投影
5. 好消息是， 經度是等距分布的， 比如圖片左邊經度是20， 右邊是 22， 那么中間豎線經度是 21。 因此， 我們僅需搞定緯度
解決辦法：
1. 古德曼函數 的輸入是一個點距離赤道的“距離” y，輸出是該點的緯度 lat。 我們的需求恰好相反， 已知 lat, 希望知道 y。 wiki 同時也提供了古德曼函數的反函數， 因此， 直接使用就好
2. 古德曼反函數有多種形式， wiki 的中文和英文版本至少提供了三種實現形式， 選一個就好， 都是等價的
3. 古德曼函數的一段 python 代碼. 注意函數輸入 longitude 不是 30， 60 這種值， 而是 (30/180 * math.pi) 這種類型的輸入
def reverseGoodman(longitude): return math.log((1 + math.sin(longitude)) / (1 - math.sin(longitude)))/2
輸入 0 時，返回 0； 輸入為 85.05112877980659， 返回的是 3.1415...(math.pi)
我本期待輸入 [0, 85.05112877980659] 的輸出應當是 [0, 1]， 所以這個地方疑惑了挺久, 后來發現只需要在原始公式的基礎上除以 math.pi 即可將值域鎖定到 [0,1]
4. 墨卡托投影的假設是地球經度為 [-180, 180], 緯度為[-85.05112877980659，85.05112877980659], 投影出的地球是一個正方形， 長寬分別為 [-20037508.3427892,20037508.3427892]
5. 聯合 (3, 4)， 就能完成一般的計算了. 假如某點的緯度為 31度, 那么該點距離赤道的距離是
20037508.3427892/math.pi*reverseGoodman(31/180*math.pi)
再高級一點， 假如已知兩個點的緯度， 求其在經度方向上的距離差，為此， 我寫了個函數
def reverseGoodman(longitude): return math.log((1 + math.sin(longitude)) / (1 - math.sin(longitude)))/math.pi/2 def convertGeoToMeter(lat1, lat2): lat1 = lat1 * math.pi / 180 lat2 = lat2 * math.pi / 180 worldMapHeight = 20037508.3427892 OffsetY1 = worldMapHeight * reverseGoodman(lat1) OffsetY2 = worldMapHeight * reverseGoodman(lat2) return OffsetY2-OffsetY1
Advanced:
假設已經獲得了 GPS trace, 上面的函數也已經准備好， 現在缺少的只剩靜態地圖的 GPS 信息了， 有了這個信息， 就能實現離線地圖標注
However, google map 的經緯度並非從 0 開始計算， 不能從 zoom level + 地圖碎片的索引 獲得圖片某一個點(如左上角)的GPS
我依然用了比較粗糙辦法： 先找一個地圖碎片，已知其索引和 zoom level， 找到其一個邊緣的某個點， 這個點必須比較好認， 然后打開 web google map, 同樣找到這個點， 右鍵這里是哪， 得到 GPS 信息
得到 GPS 信息后， 需要計算這種辦法的誤差
1. 我爬取的地圖 zoom level 為 17。 根據 google map 的原理， 在 zoom level == 17 這一層， 共有 (4^17) 個靜態地圖碎片。 其中， 橫縱方向， 各有 (2^17) 個地圖碎片
2. 每張地圖碎片的屬性和大家熟悉的地圖相同， 支持分度值， 即橫縱方向上， 每一厘米表示的實際距離是相同的
3. 手動的獲得一張地圖碎片上下邊緣的緯度， 分別為 31.027048, 31.02940。 每張地圖的長度實際為 realDis, 根據經緯度之差計算的來的與真實值之差小於 0.2 米， 准確度很高， 可以使用
realDis = 20037508.3427892/(2**16) print convertGeoToMeter(31.027048, 31.02940)*100/realDis # 99.9327328479% print math.fabs(convertGeoToMeter(31.027048, 31.02940) - realDis) #0.205668048242
4. 經度的誤差計算更加容易， 誤差依然很小
import math realLngDis = 360*1.0/(2**17) observedLngDis = 121.429158-121.426381 print realLatDis print observedLatDis print math.fabs(realLatDis-observedLatDis)/realLatDis*100 #0.00274658203125 #0.00277699999999 #1.10748444425
Material:
1. 爬取的一張地圖碎片
每張地圖碎片含有 x, y, zoom level 三個屬性
2. 地圖碎片的爬取函數， python
import urllib
import os
def scheduler(a, b, c):
per = 100 * a * b / c
if per > 100:
per = 100
print "%.2f%%" %per
def downloader(x, y):
for i in xrange(1,7):
for j in xrange(1,5):
url = 'https://mts2.google.com/vt/lyrs=s@139&hl=en&gl=CN&src=app&x=%i&s=&y=%i&z=17&s=' %(int(i)+ int(x), int(j) + int(y))
local = os.path.join('E:\Copy\google_map', 'x%iy%iz17.jpg' %(i+ int(x), j + int(y)) )
urllib.urlretrieve(url, local, scheduler)
#填寫你需要從哪張圖片開始爬取
#downloader(index_of_x, index_of_y )
代碼寫的比較死， 但代碼量比較小， 還是非常容易修改的