Python-Basemap核密度空间插值可视化绘制

DataCharm

发布于 2021-02-22 15:32:03

2.2K0

发布于 2021-02-22 15:32:03

上一篇的推文我们使用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包实现空间插值操作，敬请期待

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-12-08，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

python

本文分享自 DataCharm 微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

python

登录后参与评论

0 条评论

热度