IDW空间插值法matlab,基于IDW对PM2.5进行空间插值及可视化
前面几篇推文我们分辨介绍了使用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
IDW
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插值我们将使用现有的优秀三方包进行计算操作。文中有很多重复的知识点没有详细介绍,大家可以查看之前的推文,或者等这个系列完成后的详细源码、数据、解释文档的分享哈!(文中出错的地方小伙伴们可以私聊指出或者加群讨论哈)
推荐阅读
好文!必须在看
IDW空间插值法matlab,基于IDW对PM2.5进行空间插值及可视化相关推荐
- 空间路面matlab,基于Matlab的三维随机路面联合建模与仿真研究
摘要:为了给虚拟样机仿真分析与研究提供复杂三维路面环境,论文基于单点FFT路面不平度时域模型,建立了三维FFT随机路面数学模型,并实现了其在Matlab中的仿真:依据三角网格法理论,建立了生成三维路面 ...
- 空间平面方程matlab求解,向量代数和空间解析几何MATLAB求解.ppt
向量代数和空间解析几何MATLAB求解 第11章 向量代数与空间解析几何MATLAB求解;Outline;11.1 向量及其线性运算;11.2 数量积.向量积与混合积;11.3 曲面及其方程;3.柱面 ...
- matlab patch 六面体,《有限元基础教程》_【MATLAB算例】4.8.2(1) 基于8节点六面体单元的空间块体分析(Hexahedral3D8Node)...
[MATLAB 算例]4.8.2(1) 基于8节点六面体单元的空间块体分析(Hexahedral3D8Node) 如图4-23所示的一个空间块体,在右端部受两个集中力F 作用,其中的参数为: 1051 ...
- 透过散射介质的成像matlab,利用空间光调制器实现基于相位多样性的散射介质成像的制作方法...
本发明属于散射介质成像领域,具体说是一种利用空间光调制器实现基于相位多样性的散射介质成像方法. 背景技术: 传统的光学成像方法通常无法直接获得隐藏在散射介质后方的物体像,因此如何利用光学技术实现穿透散 ...
- 基于 Python 的自然邻域法空间插值的实现与思考
自然邻域法是基于区域大小按比例对这些样本应用权重来进行插值 (Sibson 1981),该插值也称为 Sibson 或"区域占用 (area-stealing)"插值.其基本 ...
- 空间谱估计matlab实现,空间谱估计理论与算法------程序.rar
[实例简介]包含空间谱估计理论与算法(王永良)课本对应各章的matlab程序 MATLAB程序:第2章_空间谱估计基础: 第3章_线性预测算法:第4章_多重信号分类算法:第5章_最大似然及子空间拟合算 ...
- MATLAB基于BP神经网络的手势识别
MATLAB基于BP神经网络的手势识别 摘 要:给出了采用MATLAB来识别五个手指和手背的同心圆距离,并通过颜色肤色的方法来提取手势特征量,同时利用BP神经网络算法进行误差分析来实现手势识别的设计方 ...
- BP神经网络优化 | MATLAB基于飞蛾扑火算法优化BP神经网络(MFO-BP)的预测模型(完整代码在文末)
飞蛾扑火( Moth-flame optimization algorithm,MFO) 是Seyedali Mirjalili等于2015年提出的一种新型智能优化算法.该算法具有并行优化能力强,全局 ...
- matlab 基于 libsvm工具箱的svm分类遇到的问题与解决
matlab 基于 libsvm工具箱的svm分类遇到的问题与解决 参考文章: (1)matlab 基于 libsvm工具箱的svm分类遇到的问题与解决 (2)https://www.cnblogs. ...
- MATLAB实战系列(三十七)-MATLAB基于PQ解耦风电场并网潮流计算
前言 文中涉及代码请参见 电力系统仿真-MATLAB基于PQ解耦风电场并网潮流计算 IEEE30节点.14节点.4节点标准算例,潮流计算的功能:风力发电机组并网潮流计算,并网对大电网的影响. 以下是我 ...
最新文章
- 排序算法的python实现
- 字节跳动AI-Lab算法实习生-敏感文字方向
- 一口气带你踩完五个 List 的大坑,处处坑!| 原力计划
- Mapping Text to Knowledge Graph Entities using Multi-Sense LSTMs
- NetBIOS、NETBEUI、IPX/SPX
- BT5 WIFI破解
- java 一元二次方程_java求解一元二次方程
- ZA7783是一颗将单路MIPI DSI信号转换成单路LVDS/TTL信号的转接芯片
- 18个最受欢迎的低代码开发平台【开源】
- Win11磁盘清理怎么没有了?Win11磁盘清理在哪打开?
- Python中位置参数、关键字参数、默认参数和不定长参数(非固定参数)的简介
- Linux 命令(204)—— ss 命令
- Zabbix主页应用介绍
- Android手机信号
- ovn-controller转换流表
- 错误: Failed to install 'unknown package' from GitHub: schannel: failed to receive handshake, SSL/TL
- batch文件常用命令
- 跟着Cell学作图 | 12.韦恩图(Vennerable包)
- 通过Git同步Obsidian与IOS
- 破解wifi的渗透工具
热门文章
- Django-创建一个完整的项目-详细教程
- arcgis怎么压缩tif文件_PDF文件怎么压缩才能变小?这样压缩,真的很简单!
- 某著名IT公司招聘Axapta/Navision(MBS)顾问
- steam+linux+64+fedora,在fedora中安装steam游戏平台
- 家居品牌如何在小红书上推广?家居产品推广看这里
- Himawari-8葵花数据的python读取和matlab读取
- mappedBy的基本认识
- IDEA2020版本下载、安装
- IDEA官网以往版本下载
- Unity下载安装和Android打包成APK