反距离加权插值matlab

%反距离加权插值，可以定义搜索点数
function Z=IDW(data,station,lon,lat,e,num)
%注意data为1*n向量，
%station data 对应坐标的经纬度 n4*2 第一列纬度 第二列经度
%lon lat 需要插值格点的经纬度 数据格式 横着放(横着放，横着放~)，1*n
%e参数 可设为2
%按变量选着 num为离插值点近的点 可设为4~12
%e=2和num=12为arcgis的默认值
%例子：反距离插值
% DATA=IDW2015(P_m,latlon,lon',lat',2,12);
%
[n1 n2]=size(lon);
[n3 n4]=size(data);
% if station(1,2)>station(1,1) %%%%%%%%只适用于特定区域；目的：保证精度为x
x=station(:,2)'; %已知坐标经度
y=station(:,1)'; %已知坐标纬度
% else
% y=station(:,2)'; %已知坐标wei度
% x=station(:,1)'; %已知坐标jing度
% end
z=data; %已知坐标函数值
for i1=1:n1
X=lon(i1,:); %要插值的经度
Y=lat(i1,:); %要插值的纬度
t1=length(x);
t2=length(X);
if(t1~=0&t2~=0)
%距离矩阵r生成 并插值
for j=1:t2
r=sqrt((X(j)-x).^2+(Y(j)-y).^2).^e;
[I J]=sort(r);
m=length(I);
if(m>num)
R=I(1:num);
z11=z(:,J);
z1=z11(:,1:num);
else
display('数据量不足!');
stop;
end
a=sum(z1./repmat(R,n3,1),2); %%%gai
b=sum(1./R);
Z(:,j)=a./repmat(b,n3,1);
% if isnan(Z(:,j))
% Z(:,j)=z(J(1));
% end
end
end
end