Django实战（11）——使用folium进行地理数据可视化（准研一项目2）
序
这次的工作是接着上一次的文章做的，上一次的文章主要是实现了“目标定位标绘”这个功能中的“筛选坐标”功能，最终我们将筛选出的坐标打印出来，具体可以看我上一篇文章的内容：
https://timtian.blog.csdn.net/article/details/119353720上一次实际上只是完成了筛选坐标的功能，这一次我们要实现的才是正真意义上的“地理数据可视化“。最终，我们将筛选出来的点显示在地图上，其次是将路径规划的功能和可视化实现了。经过查阅和讨论，我们最终使用了folium这个python三方库，这里附上folium的官方文档地址：
folium的使用
folium是Leaflet.js的Python的API，即可以使用Python语言调用Leaflet的地图可视化能力。
其中，Leaflet是一个非常轻的前端地图可视化库。可以直接使用pip进行安装：
pip install folium我们换用了高德地图的瓦片底图。其中，换用高德地图的瓦片底图的方法，参考了这篇博客。读者可以按照该博客的讲解，改变request的参数以获得不同的效果：
高德WMTS瓦片地图服务地图图源规律https://www.jianshu.com/p/e34f85029fd7![icon-default.png?t=L892]()
https://www.jianshu.com/p/e34f85029fd7
于是，我照猫画虎写出了绘制地图底图的代码drawMapCenter函数，传入参数是一个列表用于选择地图的中心位置的坐标，该坐标使用的是谷歌地图使用的大地坐标系，首先需要转换为高德地图使用的火星坐标系，直接使用百度地图的api实现，这个api是免费的：
# 绘制地图，确定地图的中心位置，地图大小，地图背景。
def drawMapCenter(start_coords):
response = requests.get(
'https://api.map.baidu.com/geoconv/v1/?coords=' + str(start_coords[0]) + ',' + str(
start_coords[1]) + '&from=1&to=3&ak=你的key&output=json')
# print(response.text)
x = json.loads(response.text)["result"][0]['x']
y = json.loads(response.text)["result"][0]['y']
folium_map = folium.Map(
location=[x, y],
zoom_start=15,
tiles='http://webrd02.is.autonavi.com/appmaptile?lang=zh_cn&size=1&scale=1&style=8&x={x}&y={y}&z={z}',
# tiles='https://mt.google.com/vt/lyrs=s&x={x}&y={y}&z={z}',
attr='default'
)
return folium_map接着就是写在地图上绘制点的函数：
def drawMapPoints(points, folium_map):
cnt = 0
for point in points:
print(point)
# 从列表中提取需要的坐标（浮点数）
for p in point:
pattern = re.compile(r'\d+') # 查找数字
result = pattern.findall(p)
print(result)
if len(result) != 0:
x = int(result[0]) + int(result[1]) / 1000 + 2192935.10
y = int(result[2]) + int(result[3]) / 1000 + 19472810.74
# 坐标转换，将输入的CGCS坐标系转为WSG坐标系
tmp_x, tmp_y = CGCS2WSG(x, y)
X = tmp_x[1] + tmp_x[2] / 60 + tmp_x[3] / 3600
Y = tmp_y[1] + tmp_y[2] / 60 + tmp_y[3] / 3600
# 当满足条件的点较多的时候，只在地图上显示前100个点
if cnt < 100:
response = requests.get(
'https://api.map.baidu.com/geoconv/v1/?coords=' + str(X) + ',' + str(
Y) + '&from=1&to=3&ak=你的key&output=json')
xx = json.loads(response.text)["result"][0]['x']
yy = json.loads(response.text)["result"][0]['y']
# 将点加入地图中
folium.Marker([xx, yy]).add_to(folium_map)
cnt += 1
return folium_map传入的参数有两个：points是传入点的列表，我们用for循环一层层解开，然后用正则表达式提取里面的浮点数：
re.compile(r'\d+')这一行的意思是提取出所有的整数，但实际上p这个列表里存放的是一个字符串，形如'[0000.1111, 2222.3333] '这个样子，然后正则提取出来的是0000、1111、2222、3333这样子的结果，于是我们还要把小数部分和整数部分加起来，才是最终的浮点数。由原来统一保留的是3位小数，所以直接除以1000即可得到小数部分。第二个参数folium_map传入的是之前画到一半的图，这里用来继续画下去。
仔细看的话，会看到一个叫做CGCS2WSG的函数，注释也交代的很清楚，是一个用来转换坐标系的函数，原来传入的坐标是CGCS坐标系（2000坐标系）的，然后我们得先把它转为google用的大地坐标系（即WSG坐标系），再用百度api转化为高德用的火星坐标系，这样显示在高德的瓦片底图上才会现实的准。这个CGCS2WSG函数我也不会写，是一个学长之前写过，不过是写的C++版本的，然后他给了我C++版的代码，我直接一行一行的改成了python版的代码，然后用同个数据分别输入到两个版本的程序里，发现返回结果是一样的，就说明我改的没问题。最终发现python和C++版本输入的浮点数保留的小数位数不同，一开始吓了我一跳还以为自己改错了，然后定睛一看只是小数点后的位数不同而已，前面都是一模一样的，这下松了口气。
这个是CGCS2WSG函数的实现，原理是坐标系映射。对于这个地球CGCS和WSG测量出来的椭圆的长轴短轴不一样（提示：地球的截面是一个椭圆，而不是圆），其次是单位精度不一样，映射时候去处理这两个问题即可。代码如下：
# 坐标转化函数
def CGCS2WSG(X, Y):
PI = math.pi
rad_dgr = 180 / PI
rad_min = rad_dgr * 60
rad_sec = rad_min * 60
Maxcnt = 5
a = 6378137.0
b = 6356752.3142
e2 = (a * a - b * b) / a / a
_e2 = (a * a - b * b) / b / b
c = a * a / b
bata0 = 1 - 3.0 / 4 * _e2 + 45. / 64 * _e2 * _e2 - 175.0 / 256 * pow(_e2, 3) + 11025. / 16384 * pow(_e2, 4)
C0 = bata0 * c
m0 = a * (1 - e2)
m2 = 3.0 / 2 * e2 * m0
m4 = 5.0 / 4 * e2 * m2
m6 = 7.0 / 6 * e2 * m4
m8 = 9.0 / 8 * e2 * m6
a2 = m2 / 2 + m4 / 2 + 15.0 / 32 * m6 + 7.0 / 16 * m8
a4 = m4 / 8 + 3.0 / 16 * m6 + 7.0 / 32 * m8
a6 = m6 / 32 + m8 / 16
a8 = m8 / 128
x = X
Daihao = int(Y / 1000000)
y = Y - Daihao * 1000000 - 500000
Bf = X / C0
cnt = 0
Bf1 = Bf
FBf1 = -a2 / 2 * math.sin(2 * Bf1) + a4 / 4 * math.sin(4 * Bf1) - a6 / 6 * math.sin(6 * Bf1)
Bf = (X - FBf1) / C0
cnt += 1
while abs(Bf - Bf1) >= 0.0000001 and cnt < Maxcnt:
Bf1 = Bf
FBf1 = -a2 / 2 * math.sin(2 * Bf1) + a4 / 4 * math.sin(4 * Bf1) - a6 / 6 * math.sin(6 * Bf1)
Bf = (X - FBf1) / C0
cnt += 1
if cnt == Maxcnt and abs(Bf - Bf1) >= 0.0000001:
print("Error in counting Bf!")
return
else:
Vf = math.sqrt(1 + _e2 * math.cos(Bf) * math.cos(Bf))
Nf = c / Vf
Nf2 = Nf * Nf
Mf = c / (Vf * Vf * Vf)
tf2 = math.tan(Bf) * math.tan(Bf)
etarf2 = _e2 * math.cos(Bf) * math.cos(Bf)
y2 = y * y
tb = math.tan(Bf) / (Mf * Nf) * y2
B = []
B.append(Bf - tb / 2 + tb / (24 * Nf2) * y2 * (5 + 3 * tf2 + etarf2 - 9 * etarf2 * tf2) - tb / (
720 * Nf2 * Nf2) * y2 * y2 * (61 + 90 * tf2 + 45 * tf2 * tf2))
B.append(int(B[0] * rad_dgr))
B.append(int(B[0] * rad_min - B[1] * 60))
B.append(B[0] * rad_sec - B[1] * 3600 - B[2] * 60)
L = []
tl = y / Nf / math.cos(Bf)
L.append(tl - tl / (6 * Nf2) * y2 * (1 + 2 * tf2 + etarf2) + tl / (120 * Nf2 * Nf2) * y2 * y2 * (
5 + 28 * tf2 + 24 * tf2 * tf2 + 6 * etarf2 + 8 * etarf2 * tf2))
t = ((6 * Daihao - 3) / rad_dgr + L[0]) * rad_dgr
L.append(int(t))
t = (t - L[1]) * 60
L.append(int(t))
t = (t - L[2]) * 60
L.append(t)
#
# print(B)
# print(L)
return B, L
注意，L或者是B的第1、2、3位（从0开始算）分别为坐标的度、分、秒。使用的时候需要加一下变成小数形式。
然后呢，我把最后的显示效果图展示一下。第一张图是筛选页面，第二张图是选好条件筛选完的结果页面：
图一
图二
路径规划功能以及可视化的实现
接着，就是实现路径规划功能，并给它可视化出来。本来路径规划功能是另一个学长负责的，奈何负责的博士觉得他写的”并不好“，然后客户急着要demo，于是我们还是打算用高德底图的路径规划api先顶一下，以后实在不行的情况下，再考虑自己实现一个路径规划功能。路径规划的核心代码如下所示：
# findWay函数返回值是列表，包含经过的每条路的坐标，第一维度表示单独的一条路，第二维度表示该条路经过的所有坐标，用于画线
def findWay(origin, destination):
origin_x, origin_y = CGCS2GSJ(origin[0] + 2192935.10, origin[1] + 19472810.74)
destination_x, destination_y = CGCS2GSJ(destination[0] + 2192935.10, destination[1] + 19472810.74)
origin_x = round(origin_x, 6)
origin_y = round(origin_y, 6)
destination_x = round(destination_x, 6)
destination_y = round(destination_y, 6)
response = requests.get(
'https://restapi.amap.com/v3/direction/walking?origin=' + str(origin_y) + ',' + str(origin_x) + '&destination=' + str(destination_y) + ',' + str(destination_x) + '&key=95cf688505e14a1d92a0c76a730c1b9d'
)
List = json.loads(response.text)["route"]["paths"][0]["steps"]
# result是最终返回结果
result = []
for line in List:
# one_of_routes是一个列表，每一项都是一条完整的路
# print(line)
one_of_routes = line["polyline"].split(';')
# print(one_of_routes)
route = []
for point_of_one_route in one_of_routes:
point = point_of_one_route.split(',')[::-1]
# print(point)
route.append(castStringListToFloatList(point))
result.append(route)
# print(route)
# print('')
return result首先仍然是将2000坐标系转为火星坐标系，这样才能用，是用CGCS2GSJ函数来实现的：
# 坐标转化
def CGCS2GSJ(x, y):
tmp_x, tmp_y = CGCS2WSG(x, y)
X = tmp_x[1] + tmp_x[2] / 60 + tmp_x[3] / 3600
Y = tmp_y[1] + tmp_y[2] / 60 + tmp_y[3] / 3600
response = requests.get(
'https://api.map.baidu.com/geoconv/v1/?coords=' + str(X) + ',' + str(
Y) + '&from=1&to=3&ak=你自己的key&output=json')
# print(response.text)
xx = json.loads(response.text)["result"][0]['x']
yy = json.loads(response.text)["result"][0]['y']
return xx, yy可以看出，和之前的一样，换汤不换药。先用上面讲过的那个函数将CGCS转为WSG，然后再用百度api去转为火星坐标系。这两行
X = tmp_x[1] + tmp_x[2] / 60 + tmp_x[3] / 3600
Y = tmp_y[1] + tmp_y[2] / 60 + tmp_y[3] / 3600就是之前交代的”坐标的度、分、秒。使用的时候需要加一下变成小数形式“的代码。然后从api返回的json文件中取出要的值就行了。
眼睛尖的同学还会看到这么一个名为castStringListToFloatList函数：
# 将元素类型为string的列表转化为元素类型为float的列表
def castStringListToFloatList(point):
result = []
for tmp in point:
result.append(float(tmp))
return result顾名思义，原来传入castStringListToFloatList的point是一个元素都是string类型的列表。但是呢，我们操作候坐标的时候必须得是浮点类型的，这样我们先把这个列表里的元素全给转为float才行，这个交代下实现起来很简单。
最后我找了两个点实现了下，效果图如下所示：
END
浙公网安备 33010602011771号