| 1 |
function [Tm,xc,yc,zc] = interp3ddata(xl,yl,zl,TL,xc,yc,zc,msk); |
| 2 |
% |
| 3 |
% [T,x,y,z] = interp3ddata(xl,yl,zl,TL,xc,yc,zc,msk); |
| 4 |
% dataname is a string ('temp' or 'salt') |
| 5 |
|
| 6 |
[nxl,nyl,nzl,nt]=size(TL); |
| 7 |
nxc=size(xc,1); |
| 8 |
nyc=size(yc,2); |
| 9 |
nzc=prod(size(zc)); |
| 10 |
|
| 11 |
% Use NaN's to represent no data |
| 12 |
TL=sq(TL); |
| 13 |
|
| 14 |
% $$$ % interpolate vertically first |
| 15 |
% $$$ disp('Vertical interpolation') |
| 16 |
% $$$ TLL=reshape(TL,[nxl*nyl nzl]); |
| 17 |
% $$$ TN = zeros(nxl*nyl,nzc); |
| 18 |
% $$$ TL = TN*0; |
| 19 |
% $$$ for i=1:size(TLL,1) |
| 20 |
% $$$ iz = find(~isnan(TLL(i,:))); |
| 21 |
% $$$ if length(iz)>1 |
| 22 |
% $$$ TL(i,:)=interp1(zl(iz),TLL(i,iz)',zc,'linear'); |
| 23 |
% $$$ TN(i,:)=interp1(zl(iz),TLL(i,iz)',zc,'nearest','extrap'); |
| 24 |
% $$$ end |
| 25 |
% $$$ end |
| 26 |
% $$$ in = find(isnan(TL)); TL(in) = TN(in); |
| 27 |
% $$$ TL=sq(reshape(TL,[nxl nyl nzc])); |
| 28 |
% $$$ clear TLL TN |
| 29 |
|
| 30 |
|
| 31 |
|
| 32 |
% Create halo in X direction |
| 33 |
disp('Adding wrap around data') |
| 34 |
xl=[xl(end-2:end)-360 xl xl(1:3)+360]; |
| 35 |
for k=1:size(TL,3); |
| 36 |
TLL(:,:,k)=[TL(end-2:end,:,k)' TL(:,:,k)' TL(1:3,:,k)']'; |
| 37 |
end |
| 38 |
TL=TLL; clear TLL |
| 39 |
|
| 40 |
|
| 41 |
% Interpolate horiztonally |
| 42 |
%load MASKS |
| 43 |
|
| 44 |
for k=1:size(TL,3); |
| 45 |
disp( sprintf('Interpolating (x-y) level %i',k) ); |
| 46 |
Tl=interp2(yl,xl,TL(:,:,k),yc,xc,'cubic'); |
| 47 |
Te=interp2(yl,xl,TL(:,:,k),yc,xc,'linear'); |
| 48 |
jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj); |
| 49 |
Te=interp2(yl,xl,TL(:,:,k),yc,xc,'nearest'); |
| 50 |
jj=find( isnan(Tl) ) ; Tl(jj)=Te(jj); |
| 51 |
%Tl(find(isnan(Tl)))=0; |
| 52 |
Tl=xyexpand(Tl.*sq(msk(:,:,1)),11); |
| 53 |
Tm(:,:,k)=Tl.*sq(msk(:,:,1)); |
| 54 |
end |
| 55 |
clear Tl Te jj TL |
| 56 |
|
| 57 |
%load VGRID.mat nzc zc |
| 58 |
|
| 59 |
disp('Vertical interpolation') |
| 60 |
Tl=reshape(Tm,[nxc*nyc nzl]); |
| 61 |
Tm=interp1(zl,Tl',zc,'linear','extrap'); |
| 62 |
Tn=NaN*Tm; |
| 63 |
for i=1:size(Tl,1) |
| 64 |
iz = find(~isnan(Tl(i,:))); |
| 65 |
if length(iz)> 1 |
| 66 |
Tn(:,i)=interp1(zl(iz),Tl(i,iz)',zc,'nearest','extrap'); |
| 67 |
end |
| 68 |
end |
| 69 |
in = find(isnan(Tm)); Tm(in) = Tn(in); |
| 70 |
Tm=reshape(Tm',[nxc nyc nzc]); |
| 71 |
Tm=Tm.*msk; |
| 72 |
clear Tl Tn |
| 73 |
|
| 74 |
Tmsk=mask(Tm); |
| 75 |
it=0; |
| 76 |
while( sum(Tmsk(:)-msk(:)) & it<20) |
| 77 |
it=it+1; |
| 78 |
disp( sprintf('Number of points needing vertical extrapolation = %i', ... |
| 79 |
sum(msk(:)-Tmsk(:)) ) ); |
| 80 |
Te=2*Tm(:,:,[1 1:end-1])-Tm(:,:,[1 1 1:end-2]); |
| 81 |
jj=find( Tmsk==0 & msk==1 ); |
| 82 |
Tm(jj)=Te(jj); |
| 83 |
Tmsk=mask(Tm); |
| 84 |
end |
| 85 |
|
| 86 |
|
| 87 |
% Remove NaN's |
| 88 |
Tm( find(isnan(Tm)) )=0; |
Parent Directory
| 
Revision Log
| 
Revision Graph