文章目录

• • 背景
  • 效果
  • 项目架构
  • 原理
  • • wgs84经纬度转瓦片号原理
    • wgs84经纬度转像素坐标原理
    • 样式压缩原理
    • style json书写格式
  • 展望

背景

我们知道，地图分为栅格和矢量两种。以往基本都是栅格地图，后来随着技术发展和地理数据不断积累，矢量地图慢慢开始流行。现在网络上看到的百度、高德、腾讯地图等，基本都是基于矢量切片来显示的，而遥感卫星影像，自然还是栅格切片。

当我们做地理相关的研究时，数据获取是一个重要的障碍。土地利用数据一般是很难拿到的，因为这类数据是保密数据。那么我们如果想拿到某个地物的空间分布数据，如何实现呢？

其实在早期，我看了一篇知乎的文章，《利用在线地图与审查元素快速获取超高清城市肌理》，这篇文章提到了通过百度地图个性编辑器，加上修改地图元素的body尺寸，再利用chrome的调试器的截图功能，可以获取到一副范围很大的专题图。如图：

当时的我就震惊了，因为这和我本科的毕设论文有点相关（城市边界识别，基于建设用地栅格数据）。当时我就在想，如果可以通过百度地图的个性化配置，获取建筑二值图，那么我可以把研究扩展到全国范围！（当然，这只是一个想法）

百度地图个性在线编辑器旧版

通过配置地图的style json，我们是可以拿到渲染后的栅格切片的，那么只要我们把所有切片按照x和y进行合并，就可以合成一张大范围的高清的专题图！

基于此，我认为这套想法如果落地，可以满足以下几点需求：

1. 获取web墨卡托坐标系下的某地区的遥感影像图
2. 获取web墨卡托坐标系下的某地区的样式自定义的专题图
3. 通过栅格矢量转换技术，可以实现某一要素的提取
4. 为地理相关研究提供网络地图数据源

因此，我着手开发了百度地图个性地图瓦片下载合成器，开源地址为Baidu-Styled-Tile-Downloader。
先给大家上几个效果图：

效果

以北京市朝阳区为例，下面三幅图依次是遥感影像、绿地和水体。

项目架构

• ChinaAD.py （中国行政区划相关）
• CoordTransform.py （坐标系转换相关）
• LngLatToPixel.py（wgs84经纬度转像素坐标）
• TileDownloader.py （核心类，瓦片下载流程）
• TileXY.py （wgs84经纬度转瓦片号）
• __init__.py （运行主函数）
• config.py （配置文件）
• drawBoundary.py（绘制行政边界）
• getBounds.py （获取区划外接矩形范围）

原理

首先，根据用户输入的行政区划名获取该区域的geojson数据。进一步转化成最小外接矩形的西南（左下）和东北（右上）的经纬度。

其次，根据外接矩形范围获取瓦片号，并存入队列Task。多线程执行爬取任务。

再次，当Task不为空时，根据用户配置的style json，根据压缩规则（processLabel函数）进行压缩，构建待爬瓦片的URL。

然后，爬取URL，将response.content以二进制的格式写入图片文件，并以瓦片号x和y为规则命名图片。

此后，待Task为空，遍历所有图片，解析出瓦片号x和y，并根据x和y，在对应的图片绘制，进行图片的拼接操作，最终输出完整的合成图。

最后，通过行政区划的经纬度坐标，转成以图片左上角为原点的像素坐标，利用opencv-python库绘制在图片上，保存图片。

wgs84经纬度转瓦片号原理

这里面其实涉及两步，第一步是wgs84经纬度转web墨卡托投影坐标。第二步是web墨卡托投影坐标转瓦片坐标。

我们先大体了解下百度的栅格瓦片是如何编号的。借用一张图片表示：

可以看出，百度切片的坐标原点是经纬度都为0的位置。x为东西方向，y为南北方向。可以看到，在level为2的时候，中国基本包含在（0，0）瓦片内，部分在（1，0）中。

其次，我们来了解下wgs84经纬度投影到web墨卡托坐标系下的平面坐标原理。

这里直接贴出leaflet墨卡托投影的源码，已经改为python版，市面上大部分84转web墨卡托的代码我试了下，都和百度的有较大差异，这里我认为百度的投影坐标系可能在web墨卡托基础上做了进一步的修改。其实百度提供了js api中可以将经纬度转平面坐标，但是没有python版的。

  def WGS84_to_WebMercator(self, lng, lat):
        ''' 实现WGS84向web墨卡托的转换 :param lng: WGS84经度 :param lat: WGS84纬度 :return: 转换后的web墨卡托坐标 '''
        # 该代码参考leaflet墨卡托投影源码，与百度地图的最为接近
        # source https://www.jianshu.com/p/778fc3e9f889
        R = 6378137
        R_MINOR = 6356752.314245179
        
        d = self.pi / 180
        r = R
        y = lat * d
        tmp = R_MINOR / r
        e = math.sqrt(1 - tmp * tmp)
        con = e * math.sin(y)
        ts = math.tan(self.pi / 4 - y / 2) / math.pow((1 - con) / (1 + con), e / 2)
        
        y = -r * math.log(max(ts, 1E-10))
        x = lng * d * r
        
        # 纠偏，不加这俩值投影后的坐标和百度的平面坐标是不一样的
        x+=800
        y+=650
        return x, y

获取到平面坐标后，再根据平面坐标转瓦片坐标的原理，得到瓦片坐标。

完整代码如下：

#coding = 'utf-8'
from CoordTransform import LngLatTransfer
import math

# wgs84坐标转百度地图瓦片编号xy
# 公式参考http://cntchen.github.io/2016/05/09/%E5%9B%BD%E5%86%85%E4%B8%BB%E8%A6%81%E5%9C%B0%E5%9B%BE%E7%93%A6%E7%89%87%E5%9D%90%E6%A0%87%E7%B3%BB%E5%AE%9A%E4%B9%89%E5%8F%8A%E8%AE%A1%E7%AE%97%E5%8E%9F%E7%90%86/

def WGS84_to_TILE(lng, lat, level):
	[pointX,pointY] = LngLatTransfer().WGS84_to_WebMercator(lng,lat)
	tileX = math.ceil(pointX*(math.pow(2,level-18))/256)
	tileY = math.ceil(pointY*(math.pow(2,level-18))/256)
	return [tileX,tileY]

wgs84经纬度转像素坐标原理

wgs84转像素坐标代码如下：

def WGS84_to_Pixel(lng, lat, level):
	[pointX,pointY] = LngLatTransfer().WGS84_to_WebMercator(lng,lat)
	pixelX = math.floor(pointX*math.pow(2,level-18)-math.floor(pointX*math.pow(2,level-18)/256)*256)
	pixelY = math.floor(pointY*math.pow(2,level-18)-math.floor(pointY*math.pow(2,level-18)/256)*256)
	return [pixelX,pixelY]

如果直接按照上面的函数操作是不够的，因为得到的是每个瓦片上的像素坐标，我们需要的是在整个拼合图上的像素坐标，因此，还需要利用瓦片号来对像素坐标进一步处理，得到以左上角为原点的像素坐标。函数如下：

# arr是坐标串list，形如[[lng1,lat1],[lng2,lat2],...]
# level是地图的缩放等级
# sw是行政范围西南角瓦片号，[x,y]
# ne是行政范围东北角瓦片号，[x,y]
def lineString(arr,level,sw,ne):
	res = []
	for p in arr:
		[px,py] = WGS84_to_Pixel(p[0], p[1], level)
		[tx,ty] = WGS84_to_TILE(p[0], p[1], level)
		# 实际x坐标=像素x坐标+256*(当前瓦片x-西南瓦片x)
		# 实际y坐标=256*(东北瓦片y-西南瓦片y+1)-(像素y坐标+256*(当前瓦片y-西南瓦片y))
		# y之所以这样算，是因为瓦片坐标原点是左下，但是图像坐标原点是左上，所以要用图片高度-下方的像素数=上方的像素数
		res.append([px+256*(tx-sw[0]),256*(ne[1]-sw[1]+1)-(py+256*(ty-sw[1]))])
	return res

这里还有一个注意点，就是我们获取到的行政区划geojson，featuretype是MultiPolygon的，是四维数组。为什么是四维？其实是因为Polygon是三维的，所以只是多加了一层。
Polygon的geojson形如：

{ 
'coordinates':[
# 第一个polygon
[[114,30],[114.1,30.1],...,[114,30]],
# 第二个polygon
[[114.5,30.5],[114.6,30.6],...,[114.5,30.5]]
]
}

polygon出现两个点串，只有一种可能，就是孔洞，第一个polygon要减去第二个polygon，才构成了一个polygon。如图

为了得到polygon list，我们需要对行政区划的coordinates的四维数组，去掉最外层的两层，然后再做wgs84转像素坐标的操作，公式如下：

# multipolygon是res['features'][0]['geometry']['coordinates']，res是行政区划接口返回的json
def RegionToPixels(multiPoly,level,sw,ne):
	polys = []
	for poly in multiPoly:
		for singlePoly in poly:
			polys.append(lineString(singlePoly,level,sw,ne))
	return polys

样式压缩原理

我们先来看下一个百度个性化地图瓦片请求的url长什么样。

style = 
{ x: 1584,# 瓦片x
y: 587,# 瓦片y
z: 13,# level
udt: '20181205',# timespan
scale: 1,# 图片缩放比例，例如2就是256*2
ak: '8d6c8b8f3749aed6b1aff3aad6f40e37',# ak，可以不管
styles: 't:all|e:l|v:off,t:road|e:all|v:off,t:background|e:all|c:#ffffffff,t:manmade|e:all|c:#000000ff'# style json
}

大部分参数应该都可以推出含义，就是styles我们看到，并不是完整的json格式，而是有很多缩写的。字段之间用 | 隔开，且每个json之间是用 , 隔开的。具体我就不解释了，直接上代码：

	# 生成压缩style样式字符串
	def generateStyle(self):
		res = ''
		for feature in self.style:
			f_list = []
			for key in feature.keys():
				if(type(feature[key])==dict):
					for inner_key in feature[key].keys():
						f_list.append(self.processLabel(inner_key) + ':' + self.processLabel(feature[key][inner_key]))
				else:
					f_list.append(self.processLabel(key) + ':' + self.processLabel(feature[key]))
			res = res + '|'.join(f_list) + ','
			f_list.clear()
		self.style_str = res
		
	# 压缩规则
	def processLabel(self,name):
		res = str(name)
		if(res=='featureType'):
			res = 't'
		elif(res=='elementType'):
			res = 'e'
		elif(res=='labels'):
			res = 'l'
		elif(res=='geometry'):
			res = 'g'
		elif(res=='visibility'):
			res = 'v'
		elif(res=='color'):
			res = 'c'
		elif(res=='hue'):
			res = 'h'
		elif(res=='weight'):
			res = 'w'
		elif(res=='lightness'):
			res = 'l'
		elif(res=='saturation'):
			res = 's'
		else:
			pass
		return res

style json书写格式

具体参见README.md

展望

目前尚存一些问题，比如部分图片下载成功（http code 200）但是图片合成是全黑的，但是这种情况较少发生。另外，多线程爬取瓦片，显然会对百度的服务器造成一定的压力，因此如果行政区划较大，且level较高，其实是会有很多瓦片要爬取的（上万），这时候也可能存在一些异常情况，比如图片爬取终止等。 （已解决）

后期可能会考虑关闭level的参数入口，直接根据行政区划待爬的瓦片数量来决定level。

另外，考虑到用户的各种情况，后期可能也会开放根据用户提供的矩形范围直接爬取瓦片的接口。

——————————————分割线————————————————

以上是现阶段的一些成果，和我的目标还有一定的距离，我希望后续可以实现以下功能：

1. 瓦片上叠加行政区划的边界线 （已实现）
2. 集成栅格转矢量的算法，根据目标地物直接输出shapefile文件

要实现以上的功能，最主要的难点是：瓦片的左下角平面坐标（原点），以及行政边界在瓦片坐标系的相对坐标（像素为单位）计算，等有空了再实现~

——————————————2020.12.21——————————————
今天把wgs84转web墨卡托的函数改了下，就是给转换后的xy坐标各加上一个值，来平移整个坐标系。虽然估计不是百度实际的偏移方法，但是观察了下武汉市武昌区的行政区划和北京市朝阳区的行政区划，基本和实际吻合，满足作为样例图的需求。至此，一个比较强的需求算是解决了。

本来想用python调用百度js api里面的转换函数，但是发现会用开发者自己的key来转换坐标，很容易耗尽，所以打消了这个念头。

目前这套方法，无需用户提供key，只需要修改配置文件config.py就可以下载到全国各省市的遥感影像和个性化配置地图。
再来附上开源链接，欢迎star

Baidu-Styled-Tile-Downloader