基于行政区划的百度个性化地图及遥感影像栅格瓦片下载合成器(可叠加行政区划,适合用作科研遥感影像附图)

文章目录

背景

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

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

其实在早期,我看了一篇知乎的文章,《利用在线地图与审查元素快速获取超高清城市肌理》,这篇文章提到了通过百度地图个性编辑器,加上修改地图元素的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

    原文作者:HouGISer
    原文地址: https://blog.csdn.net/lyandgh/article/details/109964772
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞