% Program to interpolate temperature fields using 3 sinterpolation
% If you need just one time cressman interp, just repeat using this function
% "sub_cressman"
%=======================================================
%
% User options:
% Reference of R: R=1. The unit of R is degree of longitude.
% Input format: [DATA] = cressman_interp(lon,lat,data_lon,data_lat,data,R,l_w)
%
%=======================================================
% lon : new 1D longitude grid.
% lat : new 1D latitude grid.
% l_w : input('Which weighting function? (1=1,2=Square,3=Exponential) ');
% R : Influenced radius. Reference of R: R=1.
% data_lon,data_lat,data : observations
%% Set up weighting function
function Tem = cressman_interp(lon,lat,data_lon,data_lat,data,R,l_w) % in degree of longitude
% lon : new 1D longitude grid.
% lat : new 1D latitude grid.
% l_w=input('Which weighting function? (1=1,2=Square,3=Exponential) ');
% init_b : we need an initial backgroud data in the grid to be interpolated.
init_b=zeros(length(lon),length(lat));
Rsq=R*R;
anal=zeros(length(lon),length(lat));
if l_w == 1
for ii=1:length(lon)
for jj=1:length(lat)
total_weight=0;
for kk=1:length(data)
dx=abs(lon(ii)-data_lon(kk));
dy=abs(lat(jj)-data_lat(kk));
dis=sqrt(dx.^2+dy.^2);
if (dis<=R)
weight =1;
anal(ii,jj) = anal(ii,jj) + (data(kk)-init_b(ii,jj))*weight;
else
weight=0;
end
total_weight=total_weight+weight;
end
if total_weight == 0
anal(ii,jj)=0;
else
anal(ii,jj)=anal(ii,jj)./total_weight+init_b(ii,jj);
end
end
end
Tem=anal;
elseif l_w == 2
%% Square weighting function
for ii=1:length(lon)
for jj=1:length(lat)
total_weight=0;
for kk=1:length(data)
dx=abs(lon(ii)-data_lon(kk));
dy=abs(lat(jj)-data_lat(kk));
dis=sqrt(dx.^2+dy.^2);
if (dis<=R)
weight = (R.^2-dis.^2) ./ (R.^2+dis.^2);
anal(ii,jj) = anal(ii,jj) + (data(kk)-init_b(ii,jj))*weight;
else
weight=0;
end
total_weight=total_weight+weight;
end
if total_weight == 0
anal(ii,jj)=0;
else
anal(ii,jj)=anal(ii,jj)./total_weight+init_b(ii,jj);
end
end
end
Tem=anal;
%% Exponential weighting function
else
for ii=1:length(lon)
for jj=1:length(lat)
total_weight=0;
for kk=1:length(data)
dx=abs(lon(ii)-data_lon(kk));
dy=abs(lat(jj)-data_lat(kk));
dis=sqrt(dx.^2+dy.^2);
if (dis<=R)
weight = exp(-dis.^2./ (2*(R/2).^2));
anal(ii,jj) = anal(ii,jj) + (data(kk)-init_b(ii,jj))*weight;
else
weight=0;
end
total_weight=total_weight+weight;
end
if total_weight == 0
anal(ii,jj)=0;
else
anal(ii,jj)=anal(ii,jj)./total_weight+init_b(ii,jj);
end
end
end
Tem=anal;
end
end