墨卡托投影， GPS 坐标转像素， GPS 坐标转距离

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:Copygoogle_map', 'x%iy%iz17.jpg'  %(i+ int(x), j + int(y)) )
            urllib.urlretrieve(url, local, scheduler)
#填写你需要从哪张图片开始爬取
#downloader(index_of_x, index_of_y )

　代码写的比较死， 但代码量比较小， 还是非常容易修改的