如何批量下载卫星地图

前言

作为城乡规划的学生，有时需要某一个城市片区的卫星地图作为设计作业的底图，地图从何来？百度？高德？谷歌？但是它们都只提供了浏览的服务啊。市面上倒是有不少自动的地图下载器。什么Bigemap啊、水经注啊、91卫图，啥的，但是它们都不是完全免费的（免费版有水印），想要授权就要帮它们推广或者花钱买。蛋疼。。。

然而作为会编程的人，怎么能让这种事情难住？

网上有不少文章分析了地理坐标系、投影坐标系的原理，还有各大地图网站瓦片的分割方式，也有提供下载的网址，只要把瓦片坐标系和缩放级别填到指定的位置，就可以下载到地图的瓦片。

参考文章如下：
客户端地图拼图算法解析
WGS84经纬度坐标与web墨卡托之间的转换【转】
国内主要地图瓦片坐标系定义及计算原理
Google 地图切片URL地址解析
腾讯与百度地图瓦片规则分析

看了这些文章以后，就会理解网络地图的分片原理，并且知道在浏览器中，是可以看到这些瓦片的地址的，只要我们知道从经纬度到瓦片坐标的转换方法，就可以在已知经纬度范围的情况下，下载这一范围内的地图。

谷歌地图

对于坐标偏移问题，国内的经纬度坐标在WGS-84的基础上经过了一些偏移（不是单纯的线性，但是在一定范围内长度面积什么的还是不会变形的，不然国内地图怎么导航？），也就是在高德等地图上采集到的坐标并不是WGS-84的坐标，而是GCJ-02坐标。在参考文章中用的谷歌的源就是国内的源（.cn结尾嘛）。不过对于当底图来看就够用了，实测应该使用http://www.google.cn/maps/采集到的经纬度来进行下载前的坐标获取。

难道瓦片坐标也进行了相应的偏移？
再次详细试验，发现gl参数决定了卫星图是否会偏移。

假设不加这个参数，会按WGS-84来，会和路线图有偏移（因为路线图只有GCJ-02的版本）
加了这个参数（gl=CN）就会使卫星图也变成GCJ-02。
也就是说，如果你使用www.google.cn/maps或者ditu.google.cn来获取坐标，那么最后你还按照这个坐标进行瓦片的获取，同时加上gl参数即可。如果你想按照WGS-84来获取坐标，那么就访问国际版的谷歌地图
，获取瓦片时不加gl参数就好了。

瓦片参数解析：
http://mt{index}.google.com/vt/lyrs={style}&hl=zh-CN&gl=CN&src=app&x={x}&y={y}&z={z}"
{index}为0到3，分别为谷歌的4个服务器，随便取一个就行~

{style}为地图的类型，比如卫星图、路线图什么的.

m：路线图
t：地形图
p：带标签的地形图
s：卫星图
y：带标签的卫星图
h：标签层（路名、地名等）

其中m路线图和s卫星图是我们想要的。
{x}{y}{z} 则分别是瓦片坐标的xy和缩放级别z，z取0-22。但是我测试发现卫星图最大只能取到20。不过即使是路线图，到20级也就足够用了。

对于python3，可以用urllib.request库进行图片的下载，然后用PIL.Image库（pillow）进行图片的合并。
如果想下载快一些，还可以使用多线程。

高德地图

高德地图的坐标是GCJ-02坐标（国内的所有地图都是，有的还进行了二次加偏）
至于怎么在高德地图上看坐标：高德地图怎么看经纬度
还是谷歌方便啊，直接地址栏就有，汗。。。。
高德的瓦片分割方式和谷歌的一样。
这里是瓦片的下载网址：
http://wprd03.is.autonavi.com/appmaptile?style=7&x=427289&y=227618&z=19
经过我的测试
style=7是路线图，6的话是卫星图。其他就不知道了。
实际测试中style=8好像是单独的路线层，背景为透明，结合ltype参数还有不同效果，这个以后可以慢慢摸索。我只下载最基本的图，这里就不深入讨论了。不过也基本够用了。一般不也就要这两个这吗？
x和y自然就是瓦片的坐标，z自然是缩放级别了。高德的z取值范围是[1,19]。不过卫星图实测只能取到18。
wprd03想必是和谷歌一样，有多个服务器提供服务。测试下来可以取到01 到 04。
也就是说，把同样的xyz从谷歌地图（中国版）填到高德里面，可以得到同一块地方的瓦片！

百度地图等

从资料看，百度的坐标在GCJ-02的基础上二次加偏形成了BD-09。百度的瓦片分割方式也和高德谷歌不一样，处理起来略麻烦，就不管了~~
腾讯地图则是瓦片坐标原点在左下角，而谷歌在左上角。xy的输入方式是进行一定的编码。

好麻烦。。。感觉有俩就够用了，其他的就不管了。

python程序

这里贴出我写的程序：

#Longittude 经度
#Latitude   纬度
#Mecator x = y = [-20037508.3427892,20037508.3427892]
#Mecator Latitue = [-85.05112877980659，85.05112877980659]

#程序使用python3哦
from math import floor,pi,log,tan,atan,exp
import urllib.request as ur
import PIL.Image as pil
import io


headers = {'User-Agent':'Mozilla/5.0 (Macintosh; Intel Mac OS X 10_7_5) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/29.0.1547.76 Safari/537.36'}

google_url="http://mt2.google.cn/vt/lyrs={style}&hl=zh-CN&gl=CN&src=app&x={x}&y={y}&z={z}"
amap_url="http://wprd02.is.autonavi.com/appmaptile?style={style}&x={x}&y={y}&z={z}"


# 暂时无用
#WGS-84经纬度转Web墨卡托
def wgs2macator(x,y):
    y = 85.0511287798 if y > 85.0511287798 else y
    y = -85.0511287798 if y < -85.0511287798 else y

    x2 = x * 20037508.34 / 180
    y2 = log(tan((90+y)*pi/360))/(pi/180)
    y2 = y2*20037508.34/180
    return x2, y2

#暂时无用
#Web墨卡托转经纬度
def mecator2wgs(x,y):
    x2 = x / 20037508.34 * 180
    y2 = y / 20037508.34 * 180
    y2= 180/pi*(2*atan(exp(y2*pi/180))-pi/2)
    return x2,y2




'''
东经为正，西经为负。北纬为正，南纬为负
j经度 w纬度 z缩放比例[0-22] ,对于卫星图并不能取到最大，测试值是20最大，再大会返回404.
'''
# 根据WGS-84 的经纬度获取谷歌地图中的瓦片坐标
def getpos(j,w,z):
    isnum=lambda x: isinstance(x,int) or isinstance(x,float)
    if not(isnum(j) and isnum(w)):
        raise TypeError("j and w must be int or float!")
        return None

    if not isinstance(z,int) or z<0:
        raise TypeError("z must be int and between 0 to 22.")
        return None

    if j<0:
        j=180+j
    else:
        j+=180
    j/=360 # make j to (0,1)

    w=85.0511287798 if w>85.0511287798 else w
    w=-85.0511287798 if w<-85.0511287798 else w
    w=log(tan((90+w)*pi/360))/(pi/180)
    w/=180 # make w to (-1,1)
    w=1-(w+1)/2 # make w to (0,1) and left top is 0-point

    num=2**z
    x=floor(j*num)
    y=floor(w*num)
    return x,y


# 根据瓦片坐标获取图像数据
def getdata(x,y,z,source,style='s'):
    '''
    根据瓦片坐标等获取谷歌地图的瓦片的下载地址。
    x y 瓦片坐标
    z   缩放级别
    style:
        s卫星图
        m路线图
    source:
        amap 高德
        google 谷歌
    '''
    if source=="amap":
        style=style.lower()
        style=6 if style=='s' else 7
        source=amap_url
    elif source == "google":
        source=google_url
    else:
        raise Exception("Unknown Map Source ! ")

    furl=source.format(x=x,y=y,z=z,style=style)
    req = ur.Request(url=furl, headers=headers) 
    try:
        data=ur.urlopen(req,timeout=5).read()
    except:
        print("get picture failed!")
        print("url:",furl)
        exit()
    return data


def getpic(x1,y1,x2,y2,z,source='google',outfile="MAP_OUT.png",style='s'):
    '''
    依次输入左上角的经度、纬度，右下角的经度、纬度，缩放级别，地图源，输出文件，影像类型（默认为卫星图）
    获取区域内的瓦片并自动拼合图像。
    '''
    pos1x, pos1y = getpos(x1, y1, z)
    pos2x, pos2y = getpos(x2, y2, z)
    lenx = pos2x - pos1x + 1
    leny = pos2y - pos1y + 1
    print("总数量：{x} X {y}".format(x=lenx,y=leny))
    outpic=pil.new('RGBA',(lenx*256,leny*256))
    datas=[]
    for i in range(lenx):
        datas.append([])
        for j in range(leny):
            print("正在下载：{0},{1}".format(i,j))
            mapdata=getdata(pos1x+i,pos1y+j,z,source,style)
            datas[i].append(mapdata)

    print("下载完成！")
    picio=None
    print('开始拼合图像......')
    for x,sublist in enumerate(datas):
        for y,data in enumerate(sublist):
            picio=io.BytesIO(data)
            small_pic=pil.open(picio)
            outpic.paste(small_pic,(x*256,y*256))
    print('拼合完成！正在导出...')

    outpic.save(outfile)
    print('导出完成！程序退出...')


def getpic_s(x,y,z,source='google',outfile="out_single.png",style="s"):
    #获得单幅瓦片图像
    getpic(x,y,x,y,z,source,outfile,style)


if __name__ == '__main__':
    #下载西安 青龙寺地块 卫星地图
    getpic(108.9797845,34.2356831,108.9949663,34.2275018,16,'amap',outfile="myout.png")

qinglongsi map

后续会增加一些功能，可以访问我的Github：
https://github.com/yuansu2016/pygetmap