matlab_class/sample_processing/prep_remap.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%starting point:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gcmfaces_global;
% mygrid_refgrid=mygrid;
mygrid=mygrid_refgrid;

fld=rdmds2gcmfaces([dir1 'ptrac_3d_set2.0000175296'],'rec',2);
fld=fld.*mygrid.mskC;
fld1=fld;
m=repmat(mygrid.mskC(:,:,28),[1 1 50]);%mask coastal region for demo of extrapolation
fld=fld.*m;
fld2=fld;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%create lat-lon grid:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ncload woa13_all_p00_01.nc lon lat depth;
lon=double(lon); lat=double(lat); depth=double(depth);
[lat,lon] = meshgrid(lat,lon);

%prepare mygrid for lat-lon with no mask
mygrid_latlon.nFaces=1;
mygrid_latlon.XC=gcmfaces({lon}); mygrid_latlon.YC=gcmfaces({lat});
mygrid_latlon.dirGrid='none';
mygrid_latlon.fileFormat='straight';
mygrid_latlon.ioSize=size(lon);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%interpolate to lat-lon grid:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii=1:2;
    if ii==1; fld=fld1;
    else; fld=fld2;
    end;
    
    mygrid=mygrid_latlon; gcmfaces_bindata; global mytri;
    veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
    veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
    
    mygrid=mygrid_refgrid; mytri=[]; gcmfaces_bindata;
    vecfld=NaN*veclon*ones(1,50);
    for kk=1:50;
        %go from gcmfaces grid to lat-lon grid
        vecfld(:,kk)=gcmfaces_interp_2d(fld(:,:,kk),veclon,veclat);
    end;
    
    %reformat into a 360*180*50 array
    mygrid=mygrid_latlon; mytri=[]; gcmfaces_bindata;
    fld_latlon=NaN*ones(360*180,50);
    fld_latlon(mytri.kk,:)=vecfld;
    fld_latlon=reshape(fld_latlon,[360 180 50]);
    
    %interpolate vertically
    coeff=interp1(-mygrid_refgrid.RC,[1:50],depth); coeff(1)=1;
    tmp1=find(isnan(coeff)); coeff(tmp1)=50-0.01; %quick fix
    
    %vertical interpolation to mygrid.RC:
    nrOut=length(depth);
    FLD=NaN*zeros(360,180,nrOut);
    for kk=1:nrOut;
        tmp1=coeff(kk); tmp2=floor(tmp1); tmp1=tmp1-tmp2;
        if tmp2==50;
            tmp3=fld_latlon(:,:,tmp2);
        else;
            tmp3=(1-tmp1)*fld_latlon(:,:,tmp2)+tmp1*fld_latlon(:,:,tmp2+1);
        end;
        FLD(:,:,kk)=tmp3;
    end;
    
    if ii==1; FLD1=FLD;
    else; FLD2=FLD;
    end;
end;%for ii=1:2;

% save testcase_remap.mat lon lat depth FLD1 FLD2 fld1 fld2;

mygrid=mygrid_refgrid; mytri=[];

1
2	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3	%starting point:
4	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
5
6	gcmfaces_global;
7	% mygrid_refgrid=mygrid;
8	mygrid=mygrid_refgrid;
9
10	fld=rdmds2gcmfaces([dir1 'ptrac_3d_set2.0000175296'],'rec',2);
11	fld=fld.*mygrid.mskC;
12	fld1=fld;
13	m=repmat(mygrid.mskC(:,:,28),[1 1 50]);%mask coastal region for demo of extrapolation
14	fld=fld.*m;
15	fld2=fld;
16
17	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18	%create lat-lon grid:
19	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
20
21	ncload woa13_all_p00_01.nc lon lat depth;
22	lon=double(lon); lat=double(lat); depth=double(depth);
23	[lat,lon] = meshgrid(lat,lon);
24
25	%prepare mygrid for lat-lon with no mask
26	mygrid_latlon.nFaces=1;
27	mygrid_latlon.XC=gcmfaces({lon}); mygrid_latlon.YC=gcmfaces({lat});
28	mygrid_latlon.dirGrid='none';
29	mygrid_latlon.fileFormat='straight';
30	mygrid_latlon.ioSize=size(lon);
31
32	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33	%interpolate to lat-lon grid:
34	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35
36	for ii=1:2;
37	if ii==1; fld=fld1;
38	else; fld=fld2;
39	end;
40
41	mygrid=mygrid_latlon; gcmfaces_bindata; global mytri;
42	veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
43	veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
44
45	mygrid=mygrid_refgrid; mytri=[]; gcmfaces_bindata;
46	vecfld=NaNveclonones(1,50);
47	for kk=1:50;
48	%go from gcmfaces grid to lat-lon grid
49	vecfld(:,kk)=gcmfaces_interp_2d(fld(:,:,kk),veclon,veclat);
50	end;
51
52	%reformat into a 36018050 array
53	mygrid=mygrid_latlon; mytri=[]; gcmfaces_bindata;
54	fld_latlon=NaNones(360180,50);
55	fld_latlon(mytri.kk,:)=vecfld;
56	fld_latlon=reshape(fld_latlon,[360 180 50]);
57
58	%interpolate vertically
59	coeff=interp1(-mygrid_refgrid.RC,[1:50],depth); coeff(1)=1;
60	tmp1=find(isnan(coeff)); coeff(tmp1)=50-0.01; %quick fix
61
62	%vertical interpolation to mygrid.RC:
63	nrOut=length(depth);
64	FLD=NaN*zeros(360,180,nrOut);
65	for kk=1:nrOut;
66	tmp1=coeff(kk); tmp2=floor(tmp1); tmp1=tmp1-tmp2;
67	if tmp2==50;
68	tmp3=fld_latlon(:,:,tmp2);
69	else;
70	tmp3=(1-tmp1)fld_latlon(:,:,tmp2)+tmp1fld_latlon(:,:,tmp2+1);
71	end;
72	FLD(:,:,kk)=tmp3;
73	end;
74
75	if ii==1; FLD1=FLD;
76	else; FLD2=FLD;
77	end;
78	end;%for ii=1:2;
79
80	% save testcase_remap.mat lon lat depth FLD1 FLD2 fld1 fld2;
81
82	mygrid=mygrid_refgrid; mytri=[];
83