最近做一个地图可视化的项目需要在地图上画出成都已开通的地铁线路图，中间还是踩了几个小坑，记录一下整个过程。
1. 开发环境
python 3.8
plotly（没安装的自行pip安装）
百度地图API（获取所有线路的地铁站点信息）
2. 效果展示
废话不多说，先上效果图。
3. 实现过程
百度地图开放平台获取站点信息（json格式）
百度坐标系转换到WGS-84坐标系
生成excel或csv表格
plotly绘制线路图
3.1 百度地图开放平台获取站点信息，链接直接上代码：
1
city_code=75 #成都市城市代码
station_info = requests.get('http://map.baidu.com/?qt=bsi&c=%s&t=%s' % (
city_code,
int(time.time() * 1000)
))
station_info_json = eval(station_info.content)
1234567
3.2 为什么要进行坐标转换
plotly用的是WGS-84坐标系，而百度地图所使用的坐标体系，是在火星坐标系的基础上又进行了一次加密处理，百度转WGS-84还稍微有点复杂，不过有开源代码，直接拿来用就可以了。
12
简单科普一下国内常用的坐标系（至于为什么有几种不同坐标系请自行google）：
地球坐标系——WGS84：常见于 GPS 设备，Google 地图等国际标准的坐标系。
火星坐标系——GCJ-02：中国国内使用的被强制加密后的坐标体系，高德坐标就属于该种坐标体系。
百度坐标系——BD-09：百度地图所使用的坐标体系，是在火星坐标系的基础上又进行了一次加密处理。
3.3 将json数据存为excel文件，一遍mapbox调用，excel文件的格式如下：
1
PS：如有需要Python学习资料的小伙伴可以加点击下方链接自行获取
获取站点信息的完整代码如下：
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri Nov 6 11:16:45 2020
Chengdu Metro Lines by plotly with python
@author: 进击的SB
"""
import requests
import time
import numpy as np
import math
# import plotly.offline as py
# import plotly.graph_objs as go
import pandas as pd
PI = math.pi
null = None
city_code=75
station_info = requests.get('http://map.baidu.com/?qt=bsi&c=%s&t=%s' % (
city_code,
int(time.time() * 1000)
))
station_info_json = eval(station_info.content)
# print(station_info_json)
# 解析地铁线路站点信息
for line in station_info_json['content']:
# i = 0
plots = []
plots_name = []
for plot in line['stops']:
plots.append([plot['x'], plot['y']])
plots_name.append(plot['name'])
# print(plots)
plot_mercator = np.array(plots)
def _transformlat(coordinates):
lng = coordinates[ : , 0] - 105
lat = coordinates[ : , 1] - 35
ret = -100 + 2 * lng + 3 * lat + 0.2 * lat * lat + \
0.1 * lng * lat + 0.2 * np.sqrt(np.fabs(lng))
ret += (20 * np.sin(6 * lng * PI) + 20 *
np.sin(2 * lng * PI)) * 2 / 3
ret += (20 * np.sin(lat * PI) + 40 *
np.sin(lat / 3 * PI)) * 2 / 3
ret += (160 * np.sin(lat / 12 * PI) + 320 *
np.sin(lat * PI / 30.0)) * 2 / 3
return ret
def _transformlng(coordinates):
lng = coordinates[ : , 0] - 105
lat = coordinates[ : , 1] - 35
ret = 300 + lng + 2 * lat + 0.1 * lng * lng + \
0.1 * lng * lat + 0.1 * np.sqrt(np.fabs(lng))
ret += (20 * np.sin(6 * lng * PI) + 20 *
np.sin(2 * lng * PI)) * 2 / 3
ret += (20 * np.sin(lng * PI) + 40 *
np.sin(lng / 3 * PI)) * 2 / 3
ret += (150 * np.sin(lng / 12 * PI) + 300 *
np.sin(lng / 30 * PI)) * 2 / 3
return ret
def gcj02_to_wgs84(coordinates):
"""
GCJ-02转WGS-84
:param coordinates: GCJ-02坐标系的经度和纬度的numpy数组
:returns: WGS-84坐标系的经度和纬度的numpy数组
"""
ee = 0.006693421622965943 # 偏心率平方
a = 6378245 # 长半轴
lng = coordinates[ : , 0]
lat = coordinates[ : , 1]
is_in_china= (lng > 73.66) & (lng < 135.05) & (lat > 3.86) & (lat < 53.55)
_transform = coordinates[is_in_china] #只对国内的坐标做偏移
dlat = _transformlat(_transform)
dlng = _transformlng(_transform)
radlat = _transform[ : , 1] / 180 * PI
magic = np.sin(radlat)
magic = 1 - ee * magic * magic
sqrtmagic = np.sqrt(magic)
dlat = (dlat * 180.0) / ((a * (1 - ee)) / (magic * sqrtmagic) * PI)
dlng = (dlng * 180.0) / (a / sqrtmagic * np.cos(radlat) * PI)
mglat = _transform[ : , 1] + dlat
mglng = _transform[ : , 0] + dlng
coordinates[is_in_china] = np.array([
_transform[ : , 0] * 2 - mglng, _transform[ : , 1] * 2 - mglat
]).T
return coordinates
def bd09_to_gcj02(coordinates):
"""
BD-09转GCJ-02
:param coordinates: BD-09坐标系的经度和纬度的numpy数组
:returns: GCJ-02坐标系的经度和纬度的numpy数组
"""
x_pi = PI * 3000 / 180
x = coordinates[ : , 0] - 0.0065
y = coordinates[ : , 1] - 0.006
z = np.sqrt(x * x + y * y) - 0.00002 * np.sin(y * x_pi)
theta = np.arctan2(y, x) - 0.000003 * np.cos(x * x_pi)
lng = z * np.cos(theta)
lat = z * np.sin(theta)
coordinates = np.array([lng, lat]).T
return coordinates
def bd09_to_wgs84(coordinates):
"""
BD-09转WGS-84
:param coordinates: BD-09坐标系的经度和纬度的numpy数组
:returns: WGS-84坐标系的经度和纬度的numpy数组
"""
return gcj02_to_wgs84(bd09_to_gcj02(coordinates))
def mercator_to_bd09(mercator):
"""
BD-09MC转BD-09
:param coordinates: GCJ-02坐标系的经度和纬度的numpy数组
:returns: WGS-84坐标系的经度和纬度的numpy数组
"""
MCBAND = [12890594.86, 8362377.87, 5591021, 3481989.83, 1678043.12, 0]
MC2LL = [[1.410526172116255e-08, 8.98305509648872e-06, -1.9939833816331,
200.9824383106796, -187.2403703815547, 91.6087516669843,
-23.38765649603339, 2.57121317296198, -0.03801003308653,
17337981.2],
[-7.435856389565537e-09, 8.983055097726239e-06, -0.78625201886289,
96.32687599759846, -1.85204757529826, -59.36935905485877,
47.40033549296737, -16.50741931063887, 2.28786674699375,
10260144.86],
[-3.030883460898826e-08, 8.98305509983578e-06, 0.30071316287616,
59.74293618442277, 7.357984074871, -25.38371002664745,
13.45380521110908, -3.29883767235584, 0.32710905363475,
6856817.37],
[-1.981981304930552e-08, 8.983055099779535e-06, 0.03278182852591,
40.31678527705744, 0.65659298677277, -4.44255534477492,
0.85341911805263, 0.12923347998204, -0.04625736007561,
4482777.06],
[3.09191371068437e-09, 8.983055096812155e-06, 6.995724062e-05,
23.10934304144901, -0.00023663490511, -0.6321817810242,
-0.00663494467273, 0.03430082397953, -0.00466043876332,
2555164.4],
[2.890871144776878e-09, 8.983055095805407e-06, -3.068298e-08,
7.47137025468032, -3.53937994e-06, -0.02145144861037,
-1.234426596e-05, 0.00010322952773, -3.23890364e-06,
826088.5]]
x = np.abs(mercator[ : , 0])
y = np.abs(mercator[ : , 1])
coef = np.array([
MC2LL[index] for index in
(np.tile(y.reshape((-1, 1)), (1, 6)) < MCBAND).sum(axis=1)
])
return converter(x, y, coef)
def converter(x, y, coef):
x_temp = coef[ : ,0] + coef[ : ,1] * np.abs(x)
x_n = np.abs(y) / coef[ : ,9]
y_temp = coef[ : ,2] + coef[ : ,3] * x_n + coef[ : ,4] * x_n ** 2 + \
coef[ : ,5] * x_n ** 3 + coef[ : ,6] * x_n ** 4 + coef[ : ,7] * x_n ** 5 + \
coef[ : ,8] * x_n ** 6
x[x < 0] = -1
x[x >= 0] = 1
y[y < 0] = -1
y[y >= 0] = 1
x_temp *= x
y_temp *= y
coordinates = np.array([x_temp, y_temp]).T
return coordinates
data = [] #绘制数据
marked = set()
# cnt = 0
for line in station_info_json['content']:
uid = line['line_uid']
if uid in marked: #由于线路包括了来回两个方向，需要排除已绘制线路的反向线路
continue
plots = [] #站台BD-09MC坐标
plots_name = [] #站台名称
for plot in line['stops']:
plots.append([plot['x'], plot['y']])
plots_name.append(plot['name'])
df1 = pd.DataFrame(columns=['站名'], data=plots_name)
plot_mercator = np.array(plots)
# print(plot_mercator)
plot_coordinates = bd09_to_wgs84(mercator_to_bd09(plot_mercator)) #站台经纬度
list_lonlat = plot_coordinates.tolist()
df2 = pd.DataFrame(columns=['经度', '纬度'], data=list_lonlat)
df3 = pd.concat([df1, df2], axis=1)
df3.to_csv('cdmerto_sites.csv', mode='a', index=False, encoding='gb2312')
csv_to_excel = pd.read_csv('cdmerto_sites.csv', encoding='gb2312')
csv_to_excel.to_excel('cdmerto_allsites.xlsx', sheet_name='Sheet1', index=False)
marked.add(uid) #添加已绘制线路的uid
marked.add(line['pair_line_uid']) #添加已绘制线路反向线路的uid
123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
如果是windows系统，需要增加csv或excel文件的路径。
3.4 使用plotly画出地铁线路图，需要先申请plotly的key，不知道怎么申请的可自行google。
1
import pandas as pd
import plotly.graph_objects as go
import plotly.express as px
prop = pd.read_excel(r'/Users/awesomeo/map/lonlat_gps84.xlsx') #此处为房源可视化数据，可忽略，只关注cd_metro
cd_metro = pd.read_excel(r'/Users/awesomeo/map/cdmerto_allsites.xlsx', encoding='gb2312', sheet_name='Sheet1')
# print(prop.head())
token = 'pk.eyJ1IjoiZm94eGpqIiwiYSI6ImNraDJ1OXNhbzBhYzEydXA2Ymt6N2R0NHAifQ.VkK9tHG3fwSnVr2k2Zxleg'
fig = px.scatter_mapbox(prop,
lon='longitude',
lat='latitude',
size='起拍价/万元',
color='面积/m²',
hover_name='楼盘名称',
hover_data=['district'],
size_max=20,
color_continuous_scale=px.colors.carto.Temps
)
# 使用unique去重得到线路条数
lines = cd_metro['线路'].unique()
# print(lines)
colors =['rgb(255,182,193)', 'rgb(138,43,226)', 'rgb(0,0,255)', 'rgb(0,191,255)', 'rgb(32,178,170)', 'rgb(230,230,250)',
'rgb(255,165,0)', 'rgb(255,69,0)', 'rgb(139,0,0)', 'rgb(240,128,128)', 'rgb(107,142,35)']
# print(list(zip(lines, colors)))
fig.add_traces([go.Scattermapbox(mode='markers + lines',
#取对应每一条线路的经纬度信息
lon = cd_metro.loc[lambda x: x['线路'] ==line]['经度'],
lat = cd_metro.loc[lambda x: x['线路'] ==line]['纬度'],
marker = {'color': color, 'size': 6},
hovertext = cd_metro.loc[lambda x: x['线路'] ==line]['站名'],
hoverinfo = 'text',
showlegend = False) for line, color in list(zip(lines,colors))])
fig.update_layout(mapbox={'accesstoken':token, 'center':{'lon':104.072329, 'lat':30.65342}, 'zoom':11.8},
margin={'l':0, 'r':0, 't':0, 'b':0}
)
fig.write_html(r'cd_metro_all.html')
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647
**重点提醒：
- 一定要进行坐标转换，要不然画出来的地铁线路是偏离的。
