地图还可以这么画_3D版

背景

前几天有个朋友在微信群里面问，如何画一个这样的图：

其实这种图还是很简单的，我记得我之前在做美赛的时候，就用python的basemap画过。当时也没公众号，代码也丢了，现在basemap也基本上过时了，都是使用cartopy画图。

这种图其实我也不知道叫什么名字，描述一下就是：map add 3d bar

⚠️：代码在文末

相似的图介绍

把这个描述在谷歌里面搜索，其实可以找到很多类似的图：

1. 比如excel可以画：

2. origin可以画：

3. powerBI可以画：

4. tableau可以画：

5. 高德地图、一些地图提供商也可以画：

6. 很多别的网站也提供服务，可以画出这样的：

上面的若干截图只是为了展示图像样式！！！！

但是

1. 很少有看到使用python画的（其实并不少）
2. 也很少能看到画中国的

那么本文，将介绍，如何使用python来画，基于中国地图做一些创作。

代码部分

导入包

# 导入包
from shapely.geometry import Polygon, MultiPolygon
from matplotlib import cm
import matplotlib
from getchinamap.getchinamap import DownloadChmap
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import cartopy.crs as ccrs 
import platform
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import geopandas as gpd
import cartopy.feature
from cartopy.mpl.patch import geos_to_path

import itertools

from mpl_toolkits.mplot3d import Axes3D
# import matplotlib.pyplot as plt
from matplotlib.collections import LineCollection,PolyCollection

# matplotlib 显示中文的问题
if platform.system() == 'Darwin':
    plt.rcParams["font.family"] = 'Arial Unicode MS'
elif platform.system() == 'Windows':
    plt.rcParams["font.family"] = 'SimHei'
else:
    pass

获得中国地图数据

chinamapdata = gpd.read_file("https://geo.datav.aliyun.com/areas_v3/bound/geojson?code=100000_full")
chinamapdata.head(3)

1. 数据来源是阿里云，这种其实只是为了画个轮廓，足够了。
2. 展示前3条，大概了解一下数据情况

样本数据

sampledata = pd.DataFrame({'lon':chinamapdata.centroid.x, 'lat':chinamapdata.centroid.y, 'name':chinamapdata['name']})
sampledata['value'] = np.random.randint(0, 1000, sampledata.shape[0])
sampledata.head(4)

1. 这里基于地图随机生成一系列样本数据，如果你有特定的数据，那么数据中要包括三列。
2. 三列分别为 lon、lat、value，分别是经纬度和对于的值，一般这个值都是正数。
3. 至于name 其实要不要无所谓，你想展示就加，不想展示就算了。反正我这里有（主要是想分享如何在3d空间 添加文本）

地图的box

bounds_box = chinamapdata.bounds
minx = bounds_box['minx'].min()
miny = bounds_box['miny'].min()
maxx = bounds_box['maxx'].max()
maxy = bounds_box['maxy'].max()

1. 计算地图的经纬度的范围而已，为后面画图的时候，限制x、y轴的范围。

数据转换

这个步骤是核心，因为现成的geopands、cartopy、matplotlib组合，是没办法直接画出来这种特定需求的图的。那么就需要对数据做一系列的转换。
主要步骤有：

1. 提取地图数据的geometry。
2. 将geometry转换成成对的数据（类似于线条的数据，要求二维的）。
3. 把数据转换成LineCollection需要的格式。

geoms = chinamapdata.geometry
# target_projection = ccrs.PlateCarree()
# geoms = [target_projection.project_geometry(geom, target_projection)
#          for geom in geoms]

paths = list(itertools.chain.from_iterable(geos_to_path(geom) for geom in geoms))

# At this point, we start working around mpl3d's slightly broken interfaces.
# So we produce a LineCollection rather than a PathCollection.
segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

lc = LineCollection(segments, color='black')

开始画图

1. 如果你是在jupyter里面写的，那这行代码可以加上，不然就可以不加。

%matplotlib widget

2. 主要画图部分

# part 1
segments = []
for path in paths:
    vertices = [vertex for vertex, _ in path.iter_segments()]
    vertices = np.asarray(vertices)
    segments.append(vertices)

# with plt.style.context('fivethirtyeight'):

fig = plt.figure()
ax = Axes3D(fig, xlim=[minx, maxx], ylim=[miny, maxy])
# ax.set_zlim(bottom=0)

lc = LineCollection(segments, color='black',linewidths=1)

ax.add_collection3d(lc)


ax.bar3d(x=sampledata['lon'], y=sampledata['lat'], z=np.zeros_like(sampledata['value']),
        dx=np.ones_like(sampledata['value']), 
        dy=np.ones_like(sampledata['value']), 
        dz=sampledata['value'],alpha=0.8)

for index, iterrow in sampledata.iterrows():
    ax.text(iterrow['lon'], iterrow['lat'], iterrow['value']+2,iterrow['name'], color='green',size=9)

ax.text(80, 30, 1000, '公众号: world of statistics', size=20, color='gray')

ax.set_xlabel('经度')
ax.set_ylabel('维度')
# ax.set_xlim([minx, maxx])
# ax.set_ylim([miny, maxy])
ax.set_zlabel('value')
    # ax.set_title("公众号: world of statistics")


plt.show()

效果如下：

地图区域填充

如果你希望可以对区域填充，只需要稍微更改几行代码即可.

# geoms
# part 2

fig = plt.figure()
ax = Axes3D(fig, xlim=[minx, maxx], ylim=[miny, maxy])


concat = lambda iterable: list(itertools.chain.from_iterable(iterable))
polys = concat(path.to_polygons() for path in paths)

lc = PolyCollection(polys, edgecolor='black',
                    facecolor='green', closed=False, alpha=0.4)

ax.add_collection3d(lc)

ax.bar3d(x=sampledata['lon'], y=sampledata['lat'], z=np.zeros_like(sampledata['value']),
        dx=np.ones_like(sampledata['value']), 
        dy=np.ones_like(sampledata['value']), 
        dz=sampledata['value'],alpha=0.8)


ax.text(80, 30, 1000, '公众号: world of statistics', size=20, color='gray')

ax.set_xlabel('经度')
ax.set_ylabel('维度')
# ax.set_xlim([minx, maxx])
# ax.set_ylim([miny, maxy])
ax.set_zlabel('value')
    # ax.set_title("公众号: world of statistics")


plt.show()

效果如下：

整体的感觉是：如果地图填充了颜色，其实还好看一点。

完整代码

本文章的代码全部在GitHub上免费分享。

1. 链接为：https://github.com/yuanzhoulvpi2017/tiny_python/tree/main/map
2. 文件为01开头的ipynb文件。