绘图提高篇 | Plotnine 核密度空间插值可视化绘制

从本期开始，我会陆续推出系列空间插值的推文教程，包括常见的「Kriging（克里金插值法）、Nearest Neighbor（最近邻点插值法）、Polynomial Regression（多元回归法）、Radial Basis Function（径向基函数法）」 等多种空间插值方法。探索空间可视化带给我们的视觉魅力。

1. 由于自己摸索前进(好在已经走通全程)，更新进度可能会有延迟。

2. 还会继续推出R-Python 的基础图表绘制推文系列。

3. 可能会根据粉丝的需求或者感兴趣图表进行专门的推文教程，大家可以给我发私信，我们会针对需求较多的图表绘制要求进行专门的推文教程。

好了，下面我们就开始今天的推文内容，本期推文主要包括：

• geopandas 绘制空间地图及裁剪操作

• colorbar定制化操作参考代码

• scipy.stats.gaussian_kde()函数进行核密度估计计算

• plotnine 绘制插值结果

geopandas 绘制空间地图及裁剪操作

在上期推文中Python-geopandas 中国地图绘制 中，我们使用了geopandas实现了中国地图的绘制，也相应分享了绘图数据(数据分享链接失效，本期会补上链接)。当然也有人私信我说安装geopandas太麻烦了，下面我们说明一下,如下：

• 针对geopandas的安装问题，最好使用 conda install --channel conda-forge geopandas 进行安装。但考虑到科学上网的问题，这一步就难住了很多人。大多人还是采用pip安装geopandas以及其依赖包，这里给出官网链接，大家可根据需求进行下载安装。geopandas官网。读取geojson 地图文件、散点数据及基础绘图代码如下：

散点数据预览如下：

具体绘图代码如下：

import?numpy?as?np
import?pandas?as?pd
import?geopandas?as?gpd
import?matplotlib.pyplot?as?plt

#统一修改绘图字体
plt.rcParams["font.family"]?=?"Roboto?Condensed"

js?=?gpd.read_file(r"江苏省.json")
nj_data?=?pd.read_excel("pmdata.xlsx")

fig,ax?=?plt.subplots(figsize=(6,4),dpi=130)
cm?=?plt.cm.get_cmap('Spectral_r')

vmin?=?nj_data["PM2.5"].min()
vmax?=?nj_data["PM2.5"].max()
js.plot(fc="none",ec="black",ax=ax,lw=.5,zorder=2)
scatter?=?ax.scatter(nj_data['经度'],nj_data["纬度"],c=nj_data["PM2.5"],s=50,ec="k",lw=0.5,cmap=cm,vmin=vmin,vmax=vmax)
#默认的colorbar?无法满足要求，这里进行定制化操作
scatter_bar?=?plt.colorbar(scatter,shrink=0.75,label="PM2.5")
#设置colorbar的outline(edgecolor)颜色
scatter_bar.outline.set_edgecolor('none')
#s设置网格
ax.grid(which="both",axis="both",c="gray",linewidth=.8,alpha=.6)
ax.set_axisbelow(True)
#轴脊设置
for?spine?in?['top','left','right','bottom']:
????ax.spines[spine].set_visible(None)?#隐去轴脊
ax.text(.5,1.1,"Map?Charts?in?Python?Exercise?01:Map?point",transform?=?ax.transAxes,ha='center',?
????????va='center',fontweight="bold",fontsize=13)
ax.text(.5,1.03,?"processed?map?charts?with?geopandas",
????????transform?=?ax.transAxes,ha='center',?va='center',fontsize?=?10,color='black')
ax.text(.83,-.07,'\nVisualization?by?DataCharm',transform?=?ax.transAxes,
????????ha='center',?va='center',fontsize?=?7,color='black')
plt.savefig(r'scatter_example.png',width=7,height=5,dpi=900,bbox_inches='tight')

可视化效果如下：

colorbar定制化操作参考代码

上面绘图代码中这里我们定制化了colorbar,代码如下：

#默认的colorbar?无法满足要求，这里进行定制化操作
scatter_bar?=?plt.colorbar(scatter,shrink=0.75,label="PM2.5")
#设置colorbar的outline(edgecolor)颜色
scatter_bar.outline.set_edgecolor('none')

此外，我还收集了多个关于设置colorbar的代码语句，如下：

#?COLORBAR
#?set?colorbar?label?plus?label?color
cb.set_label('colorbar?label',?color=fg_color)
#?set?colorbar?tick?color
cb.ax.yaxis.set_tick_params(color=fg_color)
#?set?colorbar?edgecolor?
cb.outline.set_edgecolor(fg_color)
#
cbar?=?plt.colorbar()
cbar.ax.tick_params(color="red",?width=5,?length=10)
#####?two?new?lines?####
color_bar.outline.set_color('w')???????????????????#set?colorbar?box?color
color_bar.ax.yaxis.set_tick_params(color='w')??????#set?colorbar?ticks?color?
#####?two?new?lines?####
clb.set_label('label',?labelpad=-40,?y=1.05,?rotation=0)
#cb=colorbar(cs,?drawedges=True)
cb.outline.set_color('white')
cb.outline.set_linewidth(2)
cb.dividers.set_color('white')
cb.dividers.set_linewidth(2)
cb.outline.set_edgecolor('white')
##?以前文章内容
cbar?=?plt.colorbar(ax=ax,ticks=[0,10,20,30,40],drawedges=False)
#cbar.ax.set_ylabel('Frequency',fontdict=colorbarfontdict)
cbar.ax.set_title('Counts',fontdict=colorbarfontdict,pad=8)
cbar.ax.tick_params(labelsize=12,direction='in')
cbar.ax.set_yticklabels(['0','10','20','30','>40'],family='Times?New?Roman')

希望可以大家更好的使用matplotlib魔改颜色条(colorbar)设置。

scipy.stats.gaussian_kde()函数进行核密度估计计算

在系列插值之前，我们先绘制核密度估计的插值图，在Python中物品们可以借用scipy.stats.gaussian_kde()实现空间核密度插值计算，大家也可参考scipy官网关于gaussian_kde() 的用法：高斯核密度估计参考官网。接下来我们使用该函数将散点插值到南京地图的范围之内，这里先给出代码再对应给出解释：

获取地图文件范围

这一步是为了获取插值所需要的范围，使用geopandas的total_bounds()方法即可获取：

js_box?=?js.geometry.total_bounds

gaussian_kde()插值处理

这里直接给出代码，如下：

#生成插值网格
#导入核密度估计包
import?scipy.stats?as?st
#根据获取的地图文件范围设置插值网格的大小范围
xmin?=?js_box[0]
xmax?=?js_box[2]
ymin?=?js_box[1]
ymax?=?js_box[3]

x_valuse?=?nj_data["经度"].values
y_valuse?=?nj_data["纬度"].values
#?这里注意下，也没过多解释，直接看代码即可
X,?Y?=?np.mgrid[xmin:xmax:400j,?ymin:ymax:400j]
positions?=?np.vstack([X.flatten(),Y.flatten()])
values?=?np.vstack([x_valuse,?y_valuse])
kernel?=?st.gaussian_kde(values)
Density?=?kernel(positions)
Density

注意：

1. np.mgrid[xmin:xmax:400j, ymin:ymax:400j]则是更具地图文件的范围插入400x400的网格，当然你也可以设置不同值。

2. X.flatten() 是将数据扁平化处理。

3. np.vstack() 按垂直方向（行顺序）堆叠数组构成一个新的数组,堆叠的数组需要具有相同的维度

接下来，我们将结果转换形状即可：

#reshape
Density_re?=?np.reshape(Density.T,?X.shape)

将结果存成pandas.DataFrame类型

这一步，我们将结果存成pandas的df类型，方便后续的制图操作，如下：

#将插值网格数据整理
df_grid?=pd.DataFrame(dict(long=X.flatten(),lat=Y.flatten()))
df_grid["kde"]?=?Density

结果如下(部分)：

plotnine包可视化展示

这里的可视化绘制，我们直接使用语法和ggplot2相似的python包：plotnine,感兴趣的小伙伴可以自行搜索。绘图代码如下：

import?plotnine
from?plotnine?import?*
plotnine.options.figure_size?=?(5,?4.5)
kde_map_grid?=?(ggplot()?+?
???????????geom_tile(df_grid,aes(x='long',y='lat',fill='kde'),size=0.1)?+
???????????geom_map(js,fill='none',color='gray',size=0.3)?+?
???????????scale_fill_cmap(cmap_name='Spectral_r',name='kde_value',
???????????????????????????breaks=[0.05,0.10,0.15,0.2,0.25],
???????????????????????????)+
???????????scale_x_continuous(breaks=[117,118,119,120,121,122])+
???????????labs(title="Map?Charts?in?Python?Exercise?01:?Map?point?kde",
????????????????)+
???????????#添加文本信息
???????????annotate('text',x=116.5,y=35.3,label="processed?map?charts?with?plotnine",ha="left",
???????????????????size=10)+
???????????annotate('text',x=120,y=30.6,label="Visualization?by?DataCharm",ha="left",size=9)+
???????????theme(
???????????????text=element_text(family="Roboto?Condensed"),
???????????????#修改背景
???????????????panel_background=element_blank(),
???????????????axis_ticks_major_x=element_blank(),
???????????????axis_ticks_major_y=element_blank(),
???????????????axis_text=element_text(size=12),
???????????????axis_title?=?element_text(size=14),
???????????????panel_grid_major_x=element_line(color="gray",size=.5),
???????????????panel_grid_major_y=element_line(color="gray",size=.5),
????????????))
kde_map_grid

最终可视化效果如下：

一般的绘图教程到这里也就结束了，但往往忽略了大多人人关注的“裁剪”操作，在经历过不断探索后，我们最终使用geopandas.clip() 方法完美解决此问题。

geopandas.clip()裁剪操作

在将gaussian_kde()转换成pandas df类型的数据转换成geopandas数据类型后，就可使用geopandas.clip() 方法对geodf对象进行裁剪操作，具体代码如下：

df_grid_geo?=?gpd.GeoDataFrame(df_grid,?geometry=gpd.points_from_xy(df_grid["long"],?df_grid["lat"]),crs="EPSG:4326")
#根据之前的js?geojson格式的地图文件进行裁剪
js_kde_clip?=?gpd.clip(df_grid_geo,js)

注意： 在使用gpd.clip()操作之前，请确保geopandas 安装成功，要不然 crs="EPSG:4326" 无法准确设置，进而导致无法裁剪。个人建议： pyproj must version 2.2.0 or later

再使用plotnine 对裁剪之后的js_kde_clip 数据进行绘图即可，代码和上述绘图代码一样，即数据更改而已，这里就直接放出可视化结果：

注意： 该裁剪方法只限于geopandas + plotnine 组合绘制空间可视化作品。

总结

最为第一篇插值文章，介绍的可能有些啰嗦，后续其他插值的方法我们将更为精简，希望大家可以好好看看本篇文章，下期推文使用Basemap(虽然停止维护，但还有好多优秀功能可以使用，也有对应不同 python 版本的whl文件可供下载安装)绘制此图，当然，也还有「更加实用的裁剪操作方法」。