大家好,又见面了,我是你们的朋友全栈君。
0, 说明
关于插值,官网有个小总结,可以直接去参考(从1维到多维),下面是我举的例子。
1, 一维插值interp1(x,y,X1,method)
x = linspace(0,10,11)
y = sin(x)
plot(x,y,'-ro')
插值方法有如下:
method=‘nearest’,‘linear’,‘spline’,‘pchip’,‘cubic’
比如使用三次条样插值spline,则
x = linspace(0,10,11)
y = sin(x)
plot(x,y,'-ro')
xnew = linspace(0,10,101)
f = interp1(x,y,xnew,'spline')
plot(xnew,f)
2,高维插值
2.1 二维插值
使用interp2(X,Y,Z,X1,Y1,method)
函数, method包含:nearest
,‘linear’,‘spline’,‘cubic’
x = linspace(-5,5,11);
y = linspace(-5,5,11);
[X,Y] = meshgrid(x,y);
Z = 3.*sin(X.^2+Y.^2);
xi = linspace(-5,5,100);
yi = linspace(-5,5,100);
[XI,YI]=meshgrid(xi,yi);
ZI = interp2(X,Y,Z,XI,YI,'cubic');
surf(XI,YI,ZI),view(-25,25)
2.2 三维插值
参考官网,用法如下
Vq = interp3(X,Y,Z,V,Xq,Yq,Zq)
Vq = interp3(V,Xq,Yq,Zq)
但这个方法有时候不是很实用,主要是X,Y,Z
都要弄成网格的形式,因此我一般大于等于三维的时候就直接采用interpn
了,见下一节。
2.3 多维插值
大于等于三维一般采用interpn
,可以参考官网用法,
Vq = interpn(X1,X2,...,Xn,V,Xq1,Xq2,...,Xqn)
这里x1,x2,...xn,
以及xq1,xq2,...xqn
都是一维的(注意是单调函数),所以比Interp3简单很多。
方法有如下method: ‘linear’, ‘nearest’, ‘pchip’,‘cubic’, ‘makima’, or ‘spline’. The default method is ‘linear’.
举例:
1)插值一个点
现在有一个高维数据(4维),横坐标是经度,纵坐标是维度,高是海拔,V的值是在这三维中的水汽含量。我现在有了V的数据,这个数据是(37,10,10)的大小,表示高有37层,经纬度分别都是10的大小(因此经纬度构成100的数据网格),现在要计算高500m,经纬度分别为(80,32)的点的值(插值)
data_path = 'C:\Users\Lenovo\Desktop\recent work\matlab_interp\'
el = importdata([data_path,'el_sort.txt']); %海拔
lat = importdata([data_path,'lat_sort.txt']); %纬度
lon = importdata([data_path,'lon_sort.txt']); %经度
V = readNPY([data_path,'s_h_reverse.npy']); %数据
Vq = interpn(el,lon,lat,V,500,80,32,'spline') %插值
2)插值两个点
上面插值只在一个点(500,80,30)上进行,但有时我们要插值的是很多个点构成的数组。
%在(500,80.03,32.3)和(900,81.2,31.4)两个点插值
Vq = interpn(el,lon,lat,V,[500 900],[80.03 81.2],[32.3 31.4],'spline')
- 多个点(数组)
Vq = interpn(el,lon,lat,squeeze(s_h_full(i,:,:,:)),coor(:,1),coor(:,2),coor(:,3),'spline')
这里coor表示坐标,有海拔,经,纬度,都是一列。
完整代码
data_path = 'C:\Users\Lenovo\Desktop\recent work\matlab_interp\';
el = importdata([data_path,'el_sort.txt']); %海拔
lat = importdata([data_path,'lat_sort.txt']); %纬度
lon = importdata([data_path,'lon_sort.txt']); %经度
V = readNPY([data_path,'s_h_full_inv.npy']); %数据
coor = readNPY([data_path,'coor_in_new.npy']);
levell = 50;
hours = 720;
s_h_interp = [];
for i = 1:hours
Vq = interpn(el,lon,lat,squeeze(V(i,:,:,:)),coor(:,1),coor(:,2),coor(:,3),'spline');
s_h_interp = horzcat(s_h_interp,Vq);
end
writeNPY(s_h_interp,[data_path,'s_h_interp_matlab_161.npy'])
发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/136155.html原文链接:https://javaforall.cn
【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛
【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...