气象编程 | Python反距离权重(IDW)插值计算及可视化绘制

气象学家

发布于 2020-12-18 16:21:33

12.4K1

发布于 2020-12-18 16:21:33

文章被收录于专栏：气象学家

前面几篇推文我们分辨介绍了使用Python和R绘制了二维核密度空间插值方法，并使用了Python可视化库plotnine、Basemap以及R的ggplot2完成了相关可视化教程的绘制推文，详细内容如下：

接下来，我们将继续介绍空间插值的其他方法，本期推文，我们将介绍IDW(反距离加权法（Inverse Distance Weighted）) 插值的Python计算方法及插值结果的可视化绘制过程。主要涉及的知识点如下：

IDW简介
自定义Python代码计算空间IDW
分别使用plotnine、Basemap进行IDW插值结果可视化绘制

IDW简介

反距离权重 (IDW) 插值假设：彼此距离较近的事物要比彼此距离较远的事物更相似。当为任何未测量的位置预测值时，反距离权重法会采用预测位置周围的测量值与距离预测位置较远的测量值相比，距离预测位置最近的测量值对预测值的影响更大。反距离权重法假定每个测量点都有一种局部影响，而这种影响会随着距离的增大而减小。由于这种方法为距离预测位置最近的点分配的权重较大，而权重却作为距离的函数而减小，因此称之为反距离权重法。(解释来源于网络)，繁琐的公式也没放，这里我们给出几张示意图即可，原理不解的小伙伴可自行百度。

（基于采样点距离的IDW插值（左）从高程矢量点插值的IDW曲面（右））

自定义Python代码计算空间IDW

我们免去了了繁琐的IDW插值原理部分，这节我们直接根据原理自定义IDW函数，根据已有样例站点位置及对应值，计算IDW结果。在这之前，我们给出所需样例的预览及地图文件的范围(构建插值网格所需)，结果如下：

样例点：

地图文件范围信息：

js_box = js.geometry.total_bounds
js_box
#array([116.36196 ,  30.757975, 121.975185,  35.122924])

小伙伴们对上述计算结果有疑惑的地方可以详细阅读之前的插值文章(文前链接)，或者等我将这系列做完会推出详细的源码及解释文档(目前在整理中)

定义IDW计算函数

这里主要涉及两个计算函数，计算经纬度点转实际距离(km)的haversine方法和计算IDW的函数，定义函数如下：

haversine方法:

import math
import numpy as np
#更换求距离的函数
from math import radians, cos, sin, asin, sqrt

def haversine(lon1, lat1, lon2, lat2):
    R =  6372.8
    dLon = radians(lon2 - lon1)
    dLat = radians(lat2 - lat1)
    lat1 = radians(lat1)
    lat2 = radians(lat2)
    a = sin(dLat/2)**2 + cos(lat1)*cos(lat2)*sin(dLon/2)**2
    c = 2*asin(sqrt(a))
    d = R * c
    return d

def IDW(x, y, z, xi, yi):
    lstxyzi = []
    for p in range(len(xi)):
        lstdist = []
        for s in range(len(x)):
            d = (haversine(x[s], y[s], xi[p], yi[p]))
            lstdist.append(d)
        sumsup = list((1 / np.power(lstdist, 2)))
        suminf = np.sum(sumsup)
        sumsup = np.sum(np.array(sumsup) * np.array(z))
        u = sumsup / suminf
        xyzi = [xi[p], yi[p], u]
        lstxyzi.append(xyzi)
    return(lstxyzi)

计算所需插值的网格

这里直接给出代码，阶段的结果需要更具上面的函数计算对应网格点出的IDW结果，这样就可以实现插值操作，代码如下：

js_box = js.geometry.total_bounds
#还是插入400*400的网格点
grid_lon = np.linspace(js_box[0],js_box[2],400)
grid_lat = np.linspace(js_box[1],js_box[3],400)
xgrid, ygrid = np.meshgrid(grid_lon, grid_lat)

较为简单，这里就不作多解释。

计算IDW结果

结合上面两个部分，我们进行了IDW插值结果，具体计算结果如下：

#将插值网格数据整理
df_grid =pd.DataFrame(dict(long=xgrid.flatten(),lat=ygrid.flatten()))
#这里将数组转成列表
grid_lon_list = df_grid["long"].tolist()
grid_lat_list = df_grid["lat"].tolist()

pm_idw = IDW(know_lon,know_lat,know_z,grid_lon_list,grid_lat_list)
IDW_grid_df = pd.DataFrame(pm_idw,columns=["lon","lat","idw_value"])
IDW_grid_df.head()

这样就可以得到IDW插值后的DF类型数据了，结果如下(部分）：

可视化绘制

有了规整完的插值结果，那么接下来绘制可视化结果也就非常简单了，方法和之前的几篇推文类似，具体如下：

plotnine绘制

首先，我们还是给出样例点及对应值的映射散点图，绘图过程如下：

「散点图绘制」

import plotnine
from plotnine import *
plotnine.options.figure_size = (5, 4.5)
idw_scatter = (ggplot() +
           geom_map(js,fill='none',color='gray',size=0.4) +
           geom_point(pm,aes(x='经度',y='纬度',fill='PM2.5'),size=5) +
           scale_fill_cmap(cmap_name='Spectral_r',name='PM2.5',
                           breaks=[30,40,60,80]
                           )+
           scale_x_continuous(breaks=[117,118,119,120,121,122])+
           labs(title="Map Charts in Python Exercise 02: Map IDM point",
                )+
           #添加文本信息
           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,weight="bold"),
               panel_grid_major_x=element_line(color="gray",size=.5),
               panel_grid_major_y=element_line(color="gray",size=.5),
            ))
idw_scatter

可视化结果如下：

「IDW插值结果绘制」

idw_scatter_inter = (ggplot() +
           geom_tile(IDW_grid_df,aes(x='lon',y='lat',fill='idw_value'),size=0.1) +
           geom_map(js,fill='none',color='gray',size=0.4) +
           geom_point(pm,aes(x='经度',y='纬度',fill='PM2.5'),size=4,stroke=.3,show_legend=False) +
           scale_fill_cmap(cmap_name='Spectral_r',name='idw_value',
                           breaks=[30,40,60,80]
                           )+
           scale_x_continuous(breaks=[117,118,119,120,121,122])+
           labs(title="Map Charts in Python Exercise 02: Map IDM point",
                )+
           #添加文本信息
           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),
               plot_title=element_text(size=15,weight="bold"),
               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),
            ))
idw_scatter_inter

可视化结果如下：

这里加上了散点是为了更好的对比插值结果，不加的效果如下：

裁剪操作

对研究区域的结果进行裁剪，在之前的推文中我们介绍了很多次，这里主要使用geopandas的clip() 方法进行操作，具体过程不再赘述(可以看我之前的推文教程)，我们直接给出裁剪结果：

Basemap绘制

上面介绍了plotnine包进行绘制的，这里我们再使用Basemap进行绘制，直接给出绘图代码：

from mpl_toolkits.basemap import Basemap

jiangsu_shp = r"F:\DataCharm\shpfile_data\JS\江苏省_行政边界"
fig,ax = plt.subplots(figsize=(6,4.5),dpi=130,facecolor="white")
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)
# lat = np.arange(30,36,1)
# lon = np.arange(116,122,1)
map_base.drawparallels(np.arange(30,36,1), labels=[1,0,0,0],fontsize=12,zorder=0) #画纬度线
map_base.drawmeridians(np.arange(116,122,1), labels=[0,0,0,1],fontsize=12,zorder=0) #画经度线
map_base.readshapefile(shapefile = jiangsu_shp, name = "Js", default_encoding="ISO-8859-1",
                       drawbounds=True)
cp=map_base.pcolormesh(xgrid, ygrid, data=idw_grid,cmap='Spectral_r')  
#ct=map_base.contour(xgrid, ygrid, data=idw_grid,colors='w',linewidths=.7)
#添加散点
vmin = pm["PM2.5"].min()
vmax = pm["PM2.5"].max()
ax.scatter(pm['经度'],pm["纬度"],c=pm["PM2.5"],s=90,ec="k",lw=0.5,cmap="Spectral_r",
           vmin=vmin,vmax=vmax)


colorbar = map_base.colorbar(cp,size='3%',pad="5%",label="IDW_inter")
#设置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 02:Map IDW Grid Charts",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')

可视化结果如下：

裁剪操作

裁剪的操作也在之前的推文中介绍多次，这里我们直接给出结果哈：

当然，你也可以通过basemap.contour() 添加二维等值线，可视化结果如下：

总结

这是IDW插值的第一篇推文教程，好多原理的部分也没有介绍，这里是自定义函数进行插值计算的，当然也是有优秀的第三方包可以完成。下次的R-ggplot2版本的IDW插值我们将使用现有的优秀三方包进行计算操作。文中有很多重复的知识点没有详细介绍，大家可以查看之前的推文，或者等这个系列完成后的详细源码、数据、解释文档的分享哈！(文中出错的地方小伙伴们可以私聊指出或者加群讨论哈)

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

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

编程算法

python

本文分享自气象学家微信公众号，前往查看

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

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

编程算法

python

登录后参与评论

0 条评论

热度