matlab插值计算

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:hoursVq = 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'])