美文网首页CV地图瓦片
如何批量下载卫星地图

如何批量下载卫星地图

作者: fa03399d71ce | 来源:发表于2017-03-12 00:14 被阅读1291次

    前言

    作为城乡规划的学生,有时需要某一个城市片区的卫星地图作为设计作业的底图,地图从何来?百度?高德?谷歌?但是它们都只提供了浏览的服务啊。市面上倒是有不少自动的地图下载器。什么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

    相关文章

      网友评论

      • beiou315:Mac下,pil.new('RGBA', (lenx*256, leny*256))会出错,更改为pil.new('RGB', (lenx*256, leny*256))即可
      • 70683124559c:请教一下,下载的谷歌地图某位置的的地理坐标和谷歌电子地图的坐标有偏差是什么原因啊
        70683124559c:问题已经解决了

      本文标题:如何批量下载卫星地图

      本文链接:https://www.haomeiwen.com/subject/ybuzgttx.html