[Python知识库]【GIS开发】批量下载和拼接地图瓦片（python）

1、简介

1.1 Web墨卡托投影

墨卡托投影，是正轴等角圆柱投影，又称等角圆柱投影，圆柱投影的一种，由荷兰地图学家墨卡托（G. Mercator）于1569年创拟。为地图投影方法中影响最大的。

用一张纸卷成圆柱，围住地球仪；然后在地球仪的球心放一个发光的灯泡，地球仪上的地图会在圆柱纸上形成投影；把这张纸平展开，就是我们常见平面世界地图。上面所讲的投影法叫墨卡托投影（Mercator Projection）。

地图学是一门严谨的科学，需要数学模型和精确的公式计算，再做近似处理，墨卡托投影数学模型如下：

当然，地形面积按墨卡托投影投射到平面后是有一定变形的，纬度越高，也就是越靠近两极变形越大。

Web墨卡托投影（又称球体墨卡托投影）是墨卡托投影的变种，它接收的输入是Datum为WGS84的经纬度，但在投影时不再把地球当做椭球而当做半径为6378137米的标准球体，以简化计算。

Web墨卡托投影有两个相关的投影标准，经常搞混：

• EPSG4326：Web墨卡托投影后的平面地图，但仍然使用WGS84的经度、纬度表示坐标；
• EPSG3857：Web墨卡托投影后的平面地图，坐标单位为米。

X and Y：

• X goes from 0 (left edge is 180 °W) to 2zoom ? 1 (right edge is 180 °E)
• Y goes from 0 (top edge is 85.0511 °N) to 2zoom ? 1 (bottom edge is 85.0511 °S) in a Mercator projection

1.2 经纬度坐标系

国内主要采用如下三种：

• WGS84：为一种大地坐标系，也是目前广泛使用的GPS全球卫星定位系统使用的坐标系。

• GCJ02：又称火星坐标系，是由中国国家测绘局制订的地理信息系统的坐标系统。由WGS84坐标系经加密后的坐标系。

• BD09：为百度坐标系，在GCJ02坐标系基础上再次加密。其中bd09ll表示百度经纬度坐标，bd09mc表示百度墨卡托米制坐标。

使用OpenStreetMap的坐标为WGS84；使用高德地图、腾讯地图的坐标为GCJ02；使用百度地图的坐标为BD09；谷歌地图和Bing地图的中国部分采用了高德地图的数据，所以坐标为GCJ02。

WGS84的坐标转化为GCJ02的坐标是单向的，即WGS84的坐标能够准确地变换为GCJ02坐标；但GCJ02坐标转换为WGS84时会存在精度损失。

1.3 瓦片定义

瓦片：指将一定范围内的地图按照一定的尺寸和格式，按缩放级别或者比例尺，切成若干行和列的正方形栅格图片，对切片后的正方形栅格图片被形象的称为瓦片（Tile）。

对于经过墨卡托投影为平面的世界地图，在不同的地图分辨率(整个世界地图的像素大小)下，通过切割的方式将世界地图划分为像素为 256 × 256的地图单元，划分成的每一块地图单元称为地图瓦片。

经过Web墨卡托投影后，地图就变为平面的一张地图。为此，我们对这张地图进行等级切分。在最高级（zoom=0），需要的信息最少，只需保留最重要的宏观信息，因此用一张256x256像素的图片表示即可；在下一级（zoom=1），信息量变多，用一张512x512像素的图片表示；以此类推，级别越低的像素越高，下一级的像素是当前级的4倍。这样从最高层级往下到最低层级就形成了一个金字塔坐标体系。

1.4 瓦片编号

• 谷歌XYZ：Z表示缩放层级，Z=zoom；XY的原点在左上角，X从左向右，Y从上向下。
• TMS：开源产品的标准，Z的定义与谷歌相同；XY的原点在左下角，X从左向右，Y从下向上。
• QuadTree：微软Bing地图使用的编码规范，Z的定义与谷歌相同，同一层级的瓦片不用XY两个维度表示，而只用一个整数表示，该整数服从四叉树编码规则
• 百度XYZ：Z从1开始，在最高级就把地图分为四块瓦片；XY的原点在经度为0纬度位0的位置，X从左向右，Y从下向上。

1.5 瓦片和像素

• Zoom levels

1.6 瓦片计算公式

https://wiki.openstreetmap.org/wiki/Slippy_map_tilenames

• Lon./lat. to tile numbers：

import math
def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
  return (xtile, ytile)

• Tile numbers to lon./lat：

import math
def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

我们这里开始测试一下上面的公式和代码，选择北京大学的未名湖为目标，测试如下：

• （1）查询经纬度
  地球在线，https://map.jiqrxx.com/jingweidu/
  地球在线2，https://www.earthol.com/g/
  
• （2）记录目标的经纬度
  
  未名湖的经纬度（纬度，经度）：39.99476256945049, 116.30985796451569

-（3）计算瓦片坐标编号
计算结果为：215766, 99247

• （4）构造瓦片查询网页地址

这里使用高德地图的服务接口，构造如下：

http://wprd02.is.autonavi.com/appmaptile?x=215766&y=99247&z=18&style=6

• （5）浏览器运行
  
  
  

注意：如果使用上面参数构造OpenStreetMap的瓦片网页地址，结果显示发生明显的偏移，如下：

https://a.tile.openstreetmap.org/18/215766/99247.png

1.7 网络地图服务（WMS）

Open Geospatial Consortium (OGC) home page/
https://www.osgeo.cn/ogc-e-learning/wms/slides/presentation.html
http://www.giswiki.org/wiki/Web_Map_Service

网络地图服务（Web Map Service，WMS）利用具有地理空间位置信息的数据制作地图。其中将地图定义为地理数据可视的表现。能够根据用户的请求返回相应的地图（包括PNG，GIF，JPEG等栅格形式或者是SVG和WEB CGM等矢量形式）。WMS支持网络协议HTTP，所支持的操作是由URL定义的。

1.8 切片地图服务（TMS）

https://wiki.osgeo.org/wiki/Tile_Map_Service_Specification

TMS是纯Restful；而WMTS可以有KVP、SOAP和Restful三种。
TMS瓦片是正方形的，而WMTS是矩形的（正方形是特殊的矩形）。
在纵轴方向上方向相反，TMS瓦片以左下角为原点，WMTS瓦片以左上角为原点。

1.9 切片地图web服务（WMTS）

https://www.ogc.org/standards/wmts

2、各种地图瓦片

腾讯、高德是GCJ02坐标系，百度是BD09坐标系，谷歌、必应是WGS84坐标系，天地图是CGCS2000坐标系，瓦片地图都是平面墨卡托投影。WGS84和CGCS2000坐标系，近似认为它们相等就可以。

2.1 Google

谷歌地图瓦片层级是[0~21]。

// 谷歌矢量：
http://mt2.google.cn/vt/lyrs=m&scale=2&hl=zh-CN&gl=cn&x={x}&y={y}&z={z}
// 谷歌路网：
https://mt1.google.com/vt/lyrs=h&x={x}&y={y}&z={z}
// 谷歌影像：
http://www.google.cn/maps/vt?lyrs=s@189&gl=cn&x={x}&y={y}&z={z}
// 谷歌影像带注记：
https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}
// 谷歌地形：
https://mt1.google.com/vt/lyrs=t&x={x}&y={y}&z={z}
// 谷歌地图矢量带地形渲染：
https://mt1.google.com/vt/lyrs=r&x={x}&y={y}&z={z}
// 其他
http://www.google.cn/maps/vt/pb=!1m4!1m3!1i'+str(zoom)+'!2i'+str(x)+'!3i'+str(y)+'!2m3!1e0!2sm!3i345013117!3m8!2szh-CN!3scn!5e1105!12m4!1e68!2m2!1sset!2sRoadmap!4e0

private static String Google_Satellite_Url = "http://mt0.google.cn/vt/lyrs=y&hl=zh-CN&hl=zh-CN&gl=CN&x={x}&y={y}&z={z}&s=Gali";
private static String Google_Image_Url = "http://mt0.google.cn/vt/lyrs=m&hl=zh-CN&hl=zh-CN&gl=CN&x={x}&y={y}&z={z}&s=Gali";
private static String Google_Terrain_Url = "http://mt0.google.cn/vt/lyrs=p&hl=zh-CN&hl=zh-CN&gl=CN&x={x}&y={y}&z={z}&s=Gali";

2.2 OpenStreetMap

OpenStreetMap（简称OSM）是一个网上地图协作计划，目标是创造一个内容自由且能让所有人编辑的世界地图。该平台的数据可以自由下载。OpenStreetMap（简称OSM）不仅可以免费在线使用，还可以免费下载原始数据，数据格式有.osm.pbf和shp。

// OSM瓦片网页地址：
https://tile.openstreetmap.org/{z}/{x}/{y}.png
https://c.tile.openstreetmap.org/{z}/{x}/{y}.png
https://tile-b.openstreetmap.fr/hot/{z}/{x}/{y}.png

// OSM经纬度网页地址：
https://www.openstreetmap.org/#map=18/39.99476/116.30986&layers=O

2.3 ArcGIS

// 街道图：
https://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/{z}/{y}/{x}
// 灰色图：
https://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineStreetGray/MapServer/tile/{z}/{y}/{x}
// 深蓝夜色：
https://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineStreetPurplishBlue/MapServer/tile/{z}/{y}/{x}

// 其他
https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}

举例1（北大未名湖）：

https://map.geoq.cn/ArcGIS/rest/services/ChinaOnlineCommunity/MapServer/tile/18/99247/215766
http://wprd02.is.autonavi.com/appmaptile?x=215766&y=99247&z=18&style=7

举例2（北大未名湖）：

https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/18/99247/215766
https://a.tile.openstreetmap.org/18/215766/99247.png

举例3（纽约旁的小岛）：

https://services.arcgisonline.com/arcgis/rest/services/World_Imagery/MapServer/tile/18/98434/77343
https://a.tile.openstreetmap.org/18/77343/98434.png

2.4 高德地图

高德地图瓦片坐标与 Google Map、Open Street Map 相同。高德地图的墨卡托投影截取了纬度（约85.05oS, 约85.05oN）之间部分的地球，使得投影后的平面地图水平方向和垂直方向长度相等。将墨卡托投影地图的左上角作为瓦片坐标系起点，往左方向为X轴，X轴与北纬85.05o重合且方向向左；往下方向为Y轴，Y轴与东经180o（亦为西经180o）重合且方向向下。瓦片坐标最小等级为0级，此时平面地图是一个像素为256*256的瓦片。

高德地图的坐标是GCJ-02坐标（国内的所有地图都是，有的还进行了二次加偏）,高德的瓦片分割方式和谷歌的一样,style=7是路线图，6的话是卫星图。
高德的z取值范围是[1,19]。不过卫星图实测只能取到18。

瓦片请求地址参数说明如下：

// 高德影像
https://webst01.is.autonavi.com/appmaptile?style=6&x={x}&y={y}&z={z}
// 高德路网
https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=2&style=8&ltype=3
// 高德矢量
http://wprd04.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x={x}&y={y}&z={z}

// 高德道路纯图
https://wprd01.is.autonavi.com/appmaptile?lang=zh_cn&size=1&style=7&x=215766&y=99247&z=18&scl=1&ltype=2
// 高德道路简图
http://webrd01.is.autonavi.com/appmaptile?x=215766&y=99247&z=18&size=1&scale=1&style=7
// 高德道路详图
http://webrd01.is.autonavi.com/appmaptile?x=215766&y=99247&z=18&lang=zh_cn&size=1&scale=1&style=8

• 其他

url:"https://webst{s}.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&style=6", // 添加高德影像
subdomains:["01","02","03","04"] //多域名请求

private static String AMap_Satellite_Url = "http://webst02.is.autonavi.com/appmaptile?style=6&x={x}&y={y}&z={z}";
private static String AMap_Cover_Url = "http://webst02.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scale=1&style=8";
private static String AMap_Image_Url = "http://webrd03.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}";

2.5 Bing地图

Bing Maps Tile System，https://docs.microsoft.com/en-us/bingmaps/articles/bing-maps-tile-system?redirectedfrom=MSDN

https://r1.tiles.ditu.live.com/tiles/r1321001.png?g=1001&mkt=zh-cn
http://dynamic.t2.tiles.ditu.live.com/comp/ch/{quadkey}?it=G,OS,L&mkt=en-us&cstl=w4c&ur=cn

2.6 百度地图

百度图片瓦片的层级是[3~18]。

http://online1.map.bdimg.com/onlinelabel/?qt=tile&x=49310&y=10242&z=18

2.7 搜狗地图

http://p2.go2map.com/seamless1/0/174/717/3/1/744_212.png

2.8 腾讯地图

https://rt1.map.gtimg.com/tile?z=3&x=7&y=5&styleid=1&version=117

2.9 天地图

http://lbs.tianditu.gov.cn/server/MapService.html

http://t0.tianditu.gov.cn/img_w/wmts?SERVICE=WMTS&REQUEST=GetTile&VERSION=1.0.0&LAYER=img&STYLE=default&TILEMATRIXSET=w&FORMAT=tiles&TILEMATRIX=18&TILEROW=99247&TILECOL=215766&tk=您的密钥
https://t0.tianditu.gov.cn/DataServer?T=cia_w&x=16&y=24&l=6&tk=您的密钥

3、代码实现

#***********************************************************************
#   Purpose:   python下载和拼接地图瓦片（高德和google地图）
#   Author:    爱看书的小沐
#   Date:      2022-03-04 
#   Languages: Python
#   Platform:  Python 3.9.7 win64
# ***********************************************************************
# -*-coding:utf-8 -*-

import os
import math
import urllib.request 
import PIL.Image as Image
from PIL import ImageDraw, ImageFont

file_path='d:/maps'
# (1) 全球地图
zoom = 2
xmin = ymin = 0
xmax = ymax = pow(2, zoom)

# (2) 局部地图
# zoom = 18 #19
# xmin = 215768
# ymin = 99253
# xmax = xmin + 4
# ymax = ymin + 4

#################################################################
# 图片拼接
#################################################################
mylog = open(file_path + '/err.log', mode = 'a',encoding='utf-8')

for x in range(xmin, xmax):
    for y in range(ymin, ymax):

        #Google影像瓦片
        # img_url = 'http://mt2.google.cn/vt/lyrs=y&hl=zh-CN&gl=CN&src=app&x='+str(x)+'&y='+str(y)+'&z='+str(zoom)+'&s=G'

        #高德瓦片，wprd03想必是和谷歌一样，有多个服务器提供服务。测试下来可以取到01 到 04。
        img_url = 'https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&style=6'.format(x=x, y=y, z=zoom)
        #img_url = 'https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7&ltype=0'.format(x=x, y=y, z=zoom)
        print(img_url)
        
        file_name=str(y)+'.png'
        path = file_path + "/" +  str(zoom) + "/" + str(x)
        if not os.path.exists(path):
            os.makedirs(path)

        fname = path + '/'+ file_name
        print(fname)

        try:
            urllib.request.urlretrieve(img_url,filename=fname)
        except IOError as e:
            print("IOError:"+str(x)+','+str(y))
            print(img_url, file=mylog)
        except Exception as e:
            print("Exception:"+str(x)+','+str(y))
            print(img_url, file=mylog)

mylog.close()

#################################################################
# 图片拼接
#################################################################

mw = 256 
toImage = Image.new('RGB', (256*(xmax-xmin), 256*(ymax-ymin)) )

def num2deg(xtile, ytile, zoom):
  n = 2.0 ** zoom
  lon_deg = xtile / n * 360.0 - 180.0
  lat_rad = math.atan(math.sinh(math.pi * (1 - 2 * ytile / n)))
  lat_deg = math.degrees(lat_rad)
  return (lat_deg, lon_deg)

for x in range(xmin, xmax):
    for y in range(ymin, ymax):
    
        fname =file_path + '/' + str(zoom) +'/' + str(x)+'/'+str(y)+'.png'
        fromImage = Image.open(fname)
        draw = ImageDraw.Draw(fromImage)
        
        # 添加每个瓦片的文字信息
        fontSize = 18
        setFont = ImageFont.truetype('C:/windows/fonts/arial.ttf', fontSize)
        fillColor = "#ffff00"
        pos = num2deg(x, y, zoom)
        text = 'Tile: ' + str(x)+ "," + str(y)+ "," + str(zoom)
        text2 = 'Lat: %.6f'%pos[0]
        text3 = 'Lon: %.6f'%pos[1]
        draw.text((1,1), text,font=setFont,fill=fillColor)
        draw.text((1,1+fontSize), text2,font=setFont,fill=fillColor)
        draw.text((1,1+2*fontSize), text3,font=setFont,fill=fillColor)
        
        # 添加每个瓦片的边框线条
        draw.rectangle([0, 0, mw-1, mw-1], fill=None, outline='red', width=1)
		
		# 将每个瓦片的小图绘制到大图里面。
        toImage.paste(fromImage, ((x-xmin) * mw, (y-ymin) * mw))

toImage.save(file_path + '/' + str(zoom) +'/preview.png' )
toImage.close()
print("ok")

4、运行结果

4.1 输出全球图

// 地图瓦片坐标相关参数设置
zoom = 5
xmin = ymin = 0
xmax = ymax = pow(2, zoom)

左上角第一行文字：列序号x，行序号y，级别zoom
左上角第二行文字：纬度
左上角第三行文字：经度

• 第1级别：
  

• 第2级别：
  

• 第3级别：
  

4.2 输出局部图

随机取一个位置的经纬度，计算它的瓦片坐标，同时构造瓦片网页地址。

#***********************************************************************
#   Purpose:   python根据经纬度，计算对应的瓦片坐标
#   Author:    爱看书的小沐
#   Date:      2022-03-05
#   Languages: Python
#   Platform:  Python 3.9.7 win64
# ***********************************************************************
# -*-coding:utf-8 -*-
import math

def deg2num(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
  return (xtile, ytile)

def deg2url(lat_deg, lon_deg, zoom):
  lat_rad = math.radians(lat_deg)
  n = 2.0 ** zoom
  xtile = int((lon_deg + 180.0) / 360.0 * n)
  ytile = int((1.0 - math.asinh(math.tan(lat_rad)) / math.pi) / 2.0 * n)
  return 'https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7&ltype=0'.format(x=xtile, y=ytile, z=zoom)

def deg2numrange(lat_deg1, lon_deg1, lat_deg2, lon_deg2, zoom):
  tile1 = deg2num(lat_deg1, lon_deg1, zoom)
  tile2 = deg2num(lat_deg2, lon_deg2, zoom)
  tile1new = min(tile1[0], tile2[0]), min(tile1[1], tile2[1])
  tile2new = max(tile1[0], tile2[0]), max(tile1[1], tile2[1])
  return [tile1new[0], tile1new[1], tile2new[0] - tile1new[0], tile2new[1] - tile1new[1]]

# 输出瓦片坐标
print(deg2num(39.98836718933446, 116.31269474242018, 18))
# 输出瓦片网址
print(deg2url(39.90854390955025, 116.43579204294966, 18))
# 输出两个经纬度之间的瓦片坐标
print(deg2numrange(39.98860557371063, 116.31227095339582, 39.98675601794477, 116.31610114786909, 18))

设置局部地图（以北京大学的足球体育场为例）的瓦片坐标参数如下，进行批量操作。

// 地图瓦片坐标相关参数设置
zoom = 18
xmin = 215768
ymin = 99253
xmax = xmin + 4
ymax = ymin + 4

// 高德卫星影像图
https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&style=6

// 高德交通路网图
https://wprd01.is.autonavi.com/appmaptile?x={x}&y={y}&z={z}&lang=zh_cn&size=1&scl=1&style=7&ltype=0

后记

如果你觉得该方法或代码有一点点用处，可以给作者点个赞、赏杯咖啡；╮(￣▽￣)╭
如果你感觉方法或代码不咋地//(ㄒoㄒ)//，就在评论处留言，作者继续改进。o_O???
谢谢各位童鞋们啦( ′ ▽′ )ﾉ ( ′ ▽′)っ！！！