上一篇的推文我们使用geopandas+plotnine 完美绘制高斯核密度插值的空间可视化结果,并提供了一个简单高效的裁剪方法,具体内容点击链接:Python-plotnine 核密度空间插值可视化绘制Python-plotnine 核密度空间插值可视化绘制。
本期推文,我们就使用功能强大但却被人抛弃的Basemap包进行绘制(虽然停止维护,但其空间绘图功能却依旧不能让人忽视,再者,也有对应不同版本的whl文件下载安装),主要涉及的知识点如下:
-
Basemap的pcolormesh()、contour()函数应用
-
fiona、shapely包实现目标区域裁剪操作
-
江苏省shp文件分享
Basemap的pcolormesh()、contour()函数应用
由于上篇推文中已将数据处理完成,这里我们放出处理完的数据格式预览,大家对数据处理不明白的可以参考上篇推文:数据预览如下:
我们可以看出,Density_re 数据为gaussian_kde()处理后并经过reshape操作的核密度估计插值网格数据,接下来,我们就使用Basemap包对该空间插值数据进行可视化展示,我们直接给出绘制代码:
from mpl_toolkits.basemap import Basemap
jiangsu_shp = r"江苏省_行政边界"
fig,ax = plt.subplots(figsize=(6,4.5),dpi=130)
map_base = Basemap(llcrnrlon=js_box[0], urcrnrlon=js_box[2], llcrnrlat=js_box[1],urcrnrlat=js_box[3],projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax)
map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=1) #画纬度线
map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=1) #画经度线
map_base.readshapefile(shapefile = jiangsu_shp, name = "Js", default_encoding="ISO-8859-1",
drawbounds=True)
cp=map_base.pcolormesh(X,Y, data=Density_re,cmap='Spectral_r')
colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="PM2.5_kde")
#设置colorbar
colorbar.outline.set_edgecolor('none')
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 kde point Grid",transform = ax.transAxes,ha='center',
va='center',fontweight="bold",fontsize=14)
ax.text(.5,1.03, "processed map charts with Basemap",
transform = ax.transAxes,ha='center', va='center',fontsize = 10,color='black')
ax.text(.83,-.06,'\nVisualization by DataCharm',transform = ax.transAxes,
ha='center', va='center',fontsize = 8,color='black')
plt.savefig(r'Map_point_kde_Basemap_grid_line.png',width=7,height=5,dpi=900,bbox_inches='tight')
其中:
map_base = Basemap(llcrnrlon=js_box[0], urcrnrlon=js_box[2], llcrnrlat=js_box[1],urcrnrlat=js_box[3],projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax)
构建了Basemap绘图基础,js_box结果如下:
js_box = js.geometry.total_bounds
#array([116.36196 , 30.757975, 121.975185, 35.122924])
map_base.pcolormesh()函数则实现了插值网格数据在地图上的映射效果:
cp=map_base.pcolormesh(X,Y, data=Density_re,cmap='Spectral_r')
最终的可视化效果如下:
从结果中我们可以看到,结果是规整的网格数据,没有根据目标区域(地图文件) 对结果进行裁剪,接下来我们将使用fiona、shapely包 实现对目标区域的裁剪操作。
fiona、shapely包实现目标区域裁剪操作
这里需要用到shapely.geometry的Polygon、Point方法,具体处理代码如下:
import fiona
from shapely.geometry import Polygon,Point
jiangsu_shp = fiona.open(r"江苏省_行政边界.shp")
pol = jiangsu_shp.next()
#next()之后就可以看到具体的属性值(字典类型)
poly_data = pol["geometry"]["coordinates"][0][0]
#构建Polygon(面)对象
shp_ploygeon = Polygon(poly_data)
解释:
-
使用fiona.open()方法打开目标区域(江苏省)的shp文件
-
使用next()查看shp文件的具体属性,结果如下:
-
使用Polygon()方法将其转换成面数据(较重要的一步)
这里我们查看下之前处理好的df_grid 插值网格面数据,如下(部分):
「接下来就是关键的一步操作」: 我们根据df_grid数据中的经纬度信息判断点是否在构建的面(shp_ploygeon)内,不在的点我们赋值为np.nan,在的点不变,这样即可完成“裁剪”操作,具体操作代码如下:
masked_value = [value if shp_ploygeon.contains(Point(long,lat))==True \
else np.nan for long,lat,value in zip(df_grid["long"],df_grid["lat"],df_grid["kde"])]
用到的Python列表表达式,很方便的操作方法,希望大家可以掌握。最终获取的结果如下(部分):
这里还是不太容易看出是否更改结束,我们利用pandas 选取部分数据进行观察,代码如下:
#验证结果
df_grid.loc[(df_grid["long"]>=120)&(df_grid["long"]<=121)&(df_grid["lat"]>=31)&(df_grid["lat"]>=32),\
["long","lat","kde","mask_value"]]
结果如下:
可以看到,已经根据设置更改了,接下来我们更改形状后再使用Basemap进行可视化展示。
Basemap可视化裁剪结果
在可视化之前,我们需要对数据进行reshape操作,代码如下:
mask_value_grid = df_grid["mask_value"].values.reshape(X.shape)
具体的可视化绘制代码如下,这里更改的就是我们转换之后的数据,其他的和上面代码一样,具体如下:
from mpl_toolkits.basemap import Basemap
fig,ax = plt.subplots(figsize=(6,4.5),dpi=130)
map_base = Basemap(llcrnrlon=js_box[0], urcrnrlon=js_box[2], llcrnrlat=js_box[1],urcrnrlat=js_box[3],
projection="cyl",lon_0 = 119,lat_0 = 33,ax = ax)
map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=1) #画纬度线
map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=1) #画经度线
map_base.readshapefile(shapefile = r"F:\DataCharm\shpfile_data\JS\江苏省_行政边界", name = "Js", default_encoding="ISO-8859-1",
drawbounds=True)
cp=map_base.pcolormesh(X,Y, data=mask_value_grid,cmap='Spectral_r')
colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="PM2.5_kde")
#设置colorbar
colorbar.outline.set_edgecolor('none')
#定制化操作
#轴脊设置
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 kde point",transform = ax.transAxes,ha='center',
va='center',fontweight="bold",fontsize=14)
ax.text(.5,1.03, "processed map charts with Basemap",
transform = ax.transAxes,ha='center', va='center',fontsize = 10,color='black')
ax.text(.83,-.06,'\nVisualization by DataCharm',transform = ax.transAxes,
ha='center', va='center',fontsize = 8,color='black')
plt.savefig(r'Map_point_kde_Basemap.png',
width=7,height=5,dpi=900,bbox_inches='tight')
最终裁剪之后的可视化效果如下:
Basemap.contour()绘制二维等高线图
若想在上述的结果中添加等值线,操作也十分简单,这里给出绘制代码:
map_base.contour(X,Y, data=mask_value_grid,colors='w',linewidths=.7)
即可添加成功,我们展示下裁剪后的等值线添加效果,如下:
总结
本期推文我们使用了Basemap绘制了空间插值的可视化效果,虽然这个包停止了维护,但其较为好用的绘图函数还是可以使用的,也别担心安装问题,还是提供不同版本的whl文件进行安装的。下期我们使用R-ggplot+sf包实现空间插值操作,敬请期待啊!