profilesMatlabProcessing/profiles_misc/MITprof_resample.m

function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);
%
%[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);
%
%     resamples a set of fields (fldIn) to profile 3D locations (profIn)
%     and output the result to nc file (filOut) and memory (profOut).
%     using a chosen interpolation method (method).
%
%     profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date
%     fldIn (structure) should contain: fil, name, long_name, units, tim, 
%         missing_value, and FillValue (for nc output).
%
%     fldIn.tim can be set to
%         'const' (for time invariant climatology),
%         'monclim' (for monthly climatology)
%         'monser' (for monthly time series)
%         'monloop' (for cyclic monthly time series)
%
%     method can be set to
%         'polygons' (linear in space)
%         'bindata' (nearest neighbor in space)
%
%note to self: add case when fldIn is a gcmfaces field; nctiles input
%
% Example: 
%  grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
%  profIn=idma_float_plot('4900828');
%  %
%  fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');
%  fldIn.name='prof_Towp';
%  fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
%  fldIn.units='degree C';
%  fldIn.tim='monclim';
%  fldIn.missing_value=-9999.;
%  fldIn.FillValue=-9999.;
%  %
%  profOut=MITprof_resample(profIn,fldIn);


gcmfaces_global;

doOut=~isempty(who('filOut')); doOutInit=false;
if doOut; doOut=~isempty(filOut); end;
if doOut; doOutInit=isempty(dir(filOut)); end;

if isempty(who('method')); method='polygons'; end;

%1) deal with time line
if strcmp(fldIn.tim,'monclim');
    tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];
    tmp1=datevec(profIn.prof_date);
    tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
    tim_prof=(profIn.prof_date-tmp2);
    tim_prof(tim_prof>365)=365;
    tim_prof=tim_prof/365*12;%neglecting differences in months length
elseif strcmp(fldIn.tim,'monloop');
    tmp1=dir(fldIn.fil); 
    nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
    tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1;
    tim_fld=datenum(tmp2);
    tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
    rec_fld=[nt 1:nt 1];
    tmp1=datenum([1992 1 1 0 0 0]);
    tmp2=datenum([1992+nt/12 1 1 0 0 0]);;
    tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);
    %round up tim_prof to prevent interpolation in time:
    %  tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld;
    %  tmp4=sum(tmp3>0,2);
    %  tim_prof=tim_fld(tmp4)';
elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
    tim_fld=[1 2]; rec_fld=[1 1];
    tim_prof=1.5*ones(profIn.np,1);
else;
    error('this case remains to be implemented');
end;

lon=profIn.prof_lon; lat=profIn.prof_lat;
depIn=-mygrid.RC; depOut=profIn.prof_depth;
profOut=NaN*ones(profIn.np,profIn.nr);

%2) loop over record pairs
if ~strcmp(method,'bindata'); gcmfaces_bindata; end;
for tt=1:length(rec_fld)-1;
    tt
    %
    fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));
    fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));
    ndim=length(size(fld0{1}));
    fld=cat(ndim+1,fld0,fld1);
    %
    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
    if ~isempty(ii);
      if ~strcmp(method,'bindata');
        arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
        arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
        %now linear in time:
        k0=floor(tim_prof(ii)); k1=k0+1;
        a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr);
        profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
      else;
        [prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));
        FLD=convert2array(fld(:,:,:,1));
        nk=length(mygrid.RC); kk=ones(1,nk);
        np=length(ii); pp=ones(np,1);
        ind2prof=sub2ind(size(FLD),prof_i*kk,prof_j*kk,pp*[1:nk]);
        arr=FLD(ind2prof);
        arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
        profOut(ii,:)=arr2;
      end;
      %    
      if strcmp(fldIn.tim,'std');
        profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
      end;
    end;
end;

%3) deal with file output
if doOutInit;
    %create a header only file to later append resampled fields
    prof=profIn;
    tmp1=fieldnames(prof);
    nt=length(prof.prof_date);
    nr=length(prof.prof_depth);
    for ii=1:length(tmp1);
        tmp2=prod(size(getfield(prof,tmp1{ii})));
        if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end;
    end;
    MITprof_write(filOut,prof);
end;

if doOut;
    %add the array itelf    
    MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);
    
    %add its attributes
    nc=ncopen(filOut,'write');
    ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name);
    ncaddAtt(nc,fldIn.name,'units',fldIn.units);
    ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value);
    ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue);
    ncclose(nc);
end;

%4) deal with argument output
if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end;


1	gforget	1.7	function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);
2	gforget	1.1	%
3	gforget	1.7	%[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);
4			%
5			% resamples a set of fields (fldIn) to profile 3D locations (profIn)
6			% and output the result to nc file (filOut) and memory (profOut).
7			% using a chosen interpolation method (method).
8			%
9			% profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date
10			% fldIn (structure) should contain: fil, name, long_name, units, tim,
11			% missing_value, and FillValue (for nc output).
12			%
13			% fldIn.tim can be set to
14			% 'const' (for time invariant climatology),
15			% 'monclim' (for monthly climatology)
16			% 'monser' (for monthly time series)
17			% 'monloop' (for cyclic monthly time series)
18			%
19			% method can be set to
20			% 'polygons' (linear in space)
21			% 'bindata' (nearest neighbor in space)
22			%
23			%note to self: add case when fldIn is a gcmfaces field; nctiles input
24			%
25			% Example:
26			% grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
27			% profIn=idma_float_plot('4900828');
28			% %
29			% fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');
30			% fldIn.name='prof_Towp';
31			% fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
32			% fldIn.units='degree C';
33			% fldIn.tim='monclim';
34			% fldIn.missing_value=-9999.;
35			% fldIn.FillValue=-9999.;
36			% %
37			% profOut=MITprof_resample(profIn,fldIn);
38
39	gforget	1.1
40			gcmfaces_global;
41
42			doOut=~isempty(who('filOut')); doOutInit=false;
43			if doOut; doOut=~isempty(filOut); end;
44			if doOut; doOutInit=isempty(dir(filOut)); end;
45
46			if isempty(who('method')); method='polygons'; end;
47
48			%1) deal with time line
49			if strcmp(fldIn.tim,'monclim');
50			tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];
51			tmp1=datevec(profIn.prof_date);
52			tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
53			tim_prof=(profIn.prof_date-tmp2);
54			tim_prof(tim_prof>365)=365;
55			tim_prof=tim_prof/365*12;%neglecting differences in months length
56	gforget	1.5	elseif strcmp(fldIn.tim,'monloop');
57			tmp1=dir(fldIn.fil);
58			nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
59			tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1;
60			tim_fld=datenum(tmp2);
61			tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
62			rec_fld=[nt 1:nt 1];
63			tmp1=datenum([1992 1 1 0 0 0]);
64			tmp2=datenum([1992+nt/12 1 1 0 0 0]);;
65	gforget	1.6	tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);
66			%round up tim_prof to prevent interpolation in time:
67			% tmp3=tim_profones(1,length(tim_fld))-ones(length(tim_prof),1)tim_fld;
68			% tmp4=sum(tmp3>0,2);
69			% tim_prof=tim_fld(tmp4)';
70	gforget	1.3	elseif strcmp(fldIn.tim,'const')\|strcmp(fldIn.tim,'std');
71			tim_fld=[1 2]; rec_fld=[1 1];
72	gforget	1.1	tim_prof=1.5*ones(profIn.np,1);
73			else;
74			error('this case remains to be implemented');
75			end;
76
77			lon=profIn.prof_lon; lat=profIn.prof_lat;
78			depIn=-mygrid.RC; depOut=profIn.prof_depth;
79			profOut=NaN*ones(profIn.np,profIn.nr);
80
81			%2) loop over record pairs
82	gforget	1.5	if ~strcmp(method,'bindata'); gcmfaces_bindata; end;
83	gforget	1.1	for tt=1:length(rec_fld)-1;
84	gforget	1.5	tt
85			%
86	gforget	1.2	fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));
87			fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));
88	gforget	1.1	ndim=length(size(fld0{1}));
89			fld=cat(ndim+1,fld0,fld1);
90			%
91			ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
92			if ~isempty(ii);
93	gforget	1.5	if ~strcmp(method,'bindata');
94	gforget	1.7	arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
95	gforget	1.1	arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
96	gforget	1.5	%now linear in time:
97			k0=floor(tim_prof(ii)); k1=k0+1;
98			a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr);
99			profOut(ii,:)=(1-a0).arr2(:,:,1)+a0.arr2(:,:,2);
100			else;
101			[prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));
102			FLD=convert2array(fld(:,:,:,1));
103			nk=length(mygrid.RC); kk=ones(1,nk);
104			np=length(ii); pp=ones(np,1);
105			ind2prof=sub2ind(size(FLD),prof_ikk,prof_jkk,pp*[1:nk]);
106			arr=FLD(ind2prof);
107			arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
108			profOut(ii,:)=arr2;
109			end;
110			%
111			if strcmp(fldIn.tim,'std');
112	gforget	1.4	profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
113	gforget	1.5	end;
114	gforget	1.4	end;
115	gforget	1.1	end;
116
117			%3) deal with file output
118			if doOutInit;
119			%create a header only file to later append resampled fields
120			prof=profIn;
121			tmp1=fieldnames(prof);
122			nt=length(prof.prof_date);
123			nr=length(prof.prof_depth);
124			for ii=1:length(tmp1);
125			tmp2=prod(size(getfield(prof,tmp1{ii})));
126			if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end;
127			end;
128			MITprof_write(filOut,prof);
129			end;
130
131			if doOut;
132			%add the array itelf
133			MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);
134
135			%add its attributes
136			nc=ncopen(filOut,'write');
137			ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name);
138			ncaddAtt(nc,fldIn.name,'units',fldIn.units);
139			ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value);
140			ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue);
141			ncclose(nc);
142			end;
143
144			%4) deal with argument output
145			if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end;
146
147