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 (specified in fldIn) to profile locations
%  (specified in profIn) and output the result either to memory
%  (by default) or to a netcdf file (if filOut is specified) based on
%  on pre-defined interpolation method ('polygons' by default)
%
%     profIn (structure) should contain: prof_depth, prof_lon, prof_lat,
%         and prof_date (serial date number from datenum.m)
%
%     fldIn (structure) should contain: fil, name, and tim (see below
%         for detail and examples), and optionally
%         - long_name, units, missing_value, FillValue; if filOut~=''
%           this information will be used in the netcdf fiel output
%         - fld ([] by default); if provided then it is assumed that
%           user has already read fldIn.fil and stored it to fldIn.fld
%
%     fldIn.tim must be set to one of the following values:
%         'const' (for time invariant climatology),
%         'monclim' (for monthly climatology)
%         'monser' (for monthly time series)
%         'monloop' (for cyclic monthly time series)
%
%     method ('polygons' by default) can be specified as
%         'polygons' (linear in space)
%         'bindata' (nearest neighbor in space)
%
%  Example: (should be revisited)
%
%     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.tim='monclim';
%     %fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
%     %fldIn.units='degree C';
%     %fldIn.missing_value=-9999.;
%     %fldIn.FillValue=-9999.;
%     %fldIn.fld=[];
%     %
%     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;

%0) check for input types 
%   test0=1 <-> binary
%   test1=1 <-> nctiles
%   test2=1 <-> readily available fldIn.fld
test0=isfield(fldIn,'fil');
%
test1=0;
if test0;
  test0=~isempty(dir(fldIn.fil));
  [PATH,NAME,EXT]=fileparts(fldIn.fil);
  fil_nc=fullfile(PATH,NAME,[NAME '.0001.nc']);
  fil_nctiles=fullfile(PATH,NAME,NAME);
  test1=~isempty(dir(fil_nc));
end;
%
if ~isfield(fldIn,'fld'); fldIn.fld=[]; end;
test2=~isempty(fldIn.fld);

%1) deal with time line
if strcmp(fldIn.tim,'monclim');
    tmp1=[1:13]'; tmp2=ones(13,1)*[1991 1 1 0 0 0]; tmp2(:,2)=tmp1;
    tim_fld=datenum(tmp2)-datenum(1991,1,1); 
    tim_fld=1/2*(tim_fld(1:12)+tim_fld(2:13));
    tim_fld=[tim_fld(12)-365 tim_fld' tim_fld(1)+365]; 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;
    %
    if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');
    if test1;
      eval(['ncload ' fil_nc ' tim']);
      nt=length(tim);
    elseif ~test2;
      warning('Here it is assumed that fldIn.fil contains 3D fields');
      %note: 2D case still needs to be treated here ... or via fldIn.is3d ?
      tmp1=dir(fldIn.fil);
      nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
    else;
      ndim=length(size(fldIn.fld{1}));
      fldIs3d=(ndim==4);
      nt=size(fldIn.fld{1},ndim);
    end;
    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];
    if strcmp(fldIn.tim,'monloop');
      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);
    else;
      tim_prof=profIn.prof_date;
    end;
    %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);
    if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end;
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;
  ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
  if ~isempty(ii);
    %tt
    %
    if test2&fldIs3d;
      fld0=fldIn.fld(:,:,:,rec_fld(tt));
      fld1=fldIn.fld(:,:,:,rec_fld(tt+1));
    elseif test2&~fldIs3d;
      fld0=fldIn.fld(:,:,rec_fld(tt));
      fld1=fldIn.fld(:,:,rec_fld(tt+1));
    elseif test1;
      fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));
      fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));
    elseif test0;
      fld0=read_bin(fldIn.fil,rec_fld(tt));
      fld1=read_bin(fldIn.fil,rec_fld(tt+1));
    else;
      error(['file not found:' fldIn.fil]);
    end;
    %
    ndim=length(size(fld0{1}));
    if ndim==2;
      fld0=fld0.*mygrid.mskC(:,:,1);
      fld1=fld1.*mygrid.mskC(:,:,1);
      fldIs3d=0;
    else;
      fld0=fld0.*mygrid.mskC;
      fld1=fld1.*mygrid.mskC;
      fldIs3d=1;
    end;
    fld=cat(ndim+1,fld0,fld1);
    %
    if ~strcmp(method,'bindata');
      arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
      if fldIs3d;
        arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
      else;
        arr2=arr;
      end;
      %now linear in time:
      a0=(tim_prof(ii)-tim_fld(tt))/(tim_fld(tt+1)-tim_fld(tt));
      if fldIs3d;
        a0=a0*ones(1,profIn.nr);
        profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
      else;
        profOut(ii,1)=(1-a0).*arr2(:,1)+a0.*arr2(:,2);
      end;
    elseif fldIs3d;
      [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;
    else;
      error('2D field case is missing here');
    end;
    %    
    if strcmp(fldIn.tim,'std');
      profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
    end;

  end;%if ~isempty(ii);
end;

if ~fldIs3d;
  profOut=profOut(:,1);
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;
    if ~fldIs3d;
      dims={'iPROF'};
    else;
      dims={'iDEPTH','iPROF'};
    end;
    %add the array itelf    
    MITprof_addVar(filOut,fldIn.name,'double',dims,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			%[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);
3			%
4	gforget	1.11	% resamples a set of fields (specified in fldIn) to profile locations
5			% (specified in profIn) and output the result either to memory
6			% (by default) or to a netcdf file (if filOut is specified) based on
7			% on pre-defined interpolation method ('polygons' by default)
8			%
9			% profIn (structure) should contain: prof_depth, prof_lon, prof_lat,
10			% and prof_date (serial date number from datenum.m)
11			%
12			% fldIn (structure) should contain: fil, name, and tim (see below
13			% for detail and examples), and optionally
14			% - long_name, units, missing_value, FillValue; if filOut~=''
15			% this information will be used in the netcdf fiel output
16			% - fld ([] by default); if provided then it is assumed that
17			% user has already read fldIn.fil and stored it to fldIn.fld
18	gforget	1.7	%
19	gforget	1.11	% fldIn.tim must be set to one of the following values:
20	gforget	1.7	% 'const' (for time invariant climatology),
21			% 'monclim' (for monthly climatology)
22			% 'monser' (for monthly time series)
23			% 'monloop' (for cyclic monthly time series)
24			%
25	gforget	1.11	% method ('polygons' by default) can be specified as
26	gforget	1.7	% 'polygons' (linear in space)
27			% 'bindata' (nearest neighbor in space)
28			%
29	gforget	1.11	% Example: (should be revisited)
30			%
31			% grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
32			% profIn=idma_float_plot('4900828');
33			% %
34			% fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');
35			% fldIn.name='prof_Towp';
36			% fldIn.tim='monclim';
37			% %fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
38			% %fldIn.units='degree C';
39			% %fldIn.missing_value=-9999.;
40			% %fldIn.FillValue=-9999.;
41			% %fldIn.fld=[];
42			% %
43			% profOut=MITprof_resample(profIn,fldIn);
44	gforget	1.7
45	gforget	1.1
46			gcmfaces_global;
47
48			doOut=~isempty(who('filOut')); doOutInit=false;
49			if doOut; doOut=~isempty(filOut); end;
50			if doOut; doOutInit=isempty(dir(filOut)); end;
51
52			if isempty(who('method')); method='polygons'; end;
53
54	gforget	1.11	%0) check for input types
55			% test0=1 <-> binary
56			% test1=1 <-> nctiles
57			% test2=1 <-> readily available fldIn.fld
58	gforget	1.9	test0=isfield(fldIn,'fil');
59	gforget	1.11	%
60	gforget	1.9	test1=0;
61			if test0;
62			test0=~isempty(dir(fldIn.fil));
63			[PATH,NAME,EXT]=fileparts(fldIn.fil);
64			fil_nc=fullfile(PATH,NAME,[NAME '.0001.nc']);
65			fil_nctiles=fullfile(PATH,NAME,NAME);
66			test1=~isempty(dir(fil_nc));
67			end;
68	gforget	1.11	%
69			if ~isfield(fldIn,'fld'); fldIn.fld=[]; end;
70	gforget	1.8	test2=~isempty(fldIn.fld);
71
72	gforget	1.1	%1) deal with time line
73			if strcmp(fldIn.tim,'monclim');
74	gforget	1.11	tmp1=[1:13]'; tmp2=ones(13,1)*[1991 1 1 0 0 0]; tmp2(:,2)=tmp1;
75			tim_fld=datenum(tmp2)-datenum(1991,1,1);
76			tim_fld=1/2*(tim_fld(1:12)+tim_fld(2:13));
77			tim_fld=[tim_fld(12)-365 tim_fld' tim_fld(1)+365]; rec_fld=[12 1:12 1];
78			%
79	gforget	1.1	tmp1=datevec(profIn.prof_date);
80			tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
81			tim_prof=(profIn.prof_date-tmp2);
82			tim_prof(tim_prof>365)=365;
83	gforget	1.11	%
84	gforget	1.9	if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
85	gforget	1.8	elseif strcmp(fldIn.tim,'monloop')\|strcmp(fldIn.tim,'monser');
86			if test1;
87			eval(['ncload ' fil_nc ' tim']);
88			nt=length(tim);
89			elseif ~test2;
90	gforget	1.11	warning('Here it is assumed that fldIn.fil contains 3D fields');
91	gforget	1.8	%note: 2D case still needs to be treated here ... or via fldIn.is3d ?
92			tmp1=dir(fldIn.fil);
93			nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
94			else;
95	gforget	1.9	ndim=length(size(fldIn.fld{1}));
96			fldIs3d=(ndim==4);
97			nt=size(fldIn.fld{1},ndim);
98	gforget	1.8	end;
99	gforget	1.5	tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1;
100			tim_fld=datenum(tmp2);
101			tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
102			rec_fld=[nt 1:nt 1];
103	gforget	1.8	if strcmp(fldIn.tim,'monloop');
104			tmp1=datenum([1992 1 1 0 0 0]);
105			tmp2=datenum([1992+nt/12 1 1 0 0 0]);;
106			tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);
107			else;
108			tim_prof=profIn.prof_date;
109			end;
110	gforget	1.6	%round up tim_prof to prevent interpolation in time:
111			% tmp3=tim_profones(1,length(tim_fld))-ones(length(tim_prof),1)tim_fld;
112			% tmp4=sum(tmp3>0,2);
113			% tim_prof=tim_fld(tmp4)';
114	gforget	1.3	elseif strcmp(fldIn.tim,'const')\|strcmp(fldIn.tim,'std');
115			tim_fld=[1 2]; rec_fld=[1 1];
116	gforget	1.1	tim_prof=1.5*ones(profIn.np,1);
117	gforget	1.9	if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end;
118	gforget	1.1	else;
119			error('this case remains to be implemented');
120			end;
121
122			lon=profIn.prof_lon; lat=profIn.prof_lat;
123			depIn=-mygrid.RC; depOut=profIn.prof_depth;
124			profOut=NaN*ones(profIn.np,profIn.nr);
125
126			%2) loop over record pairs
127	gforget	1.10	if strcmp(method,'bindata'); gcmfaces_bindata; end;
128	gforget	1.1	for tt=1:length(rec_fld)-1;
129	gforget	1.8	ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
130			if ~isempty(ii);
131			%tt
132			%
133	gforget	1.9	if test2&fldIs3d;
134	gforget	1.8	fld0=fldIn.fld(:,:,:,rec_fld(tt));
135			fld1=fldIn.fld(:,:,:,rec_fld(tt+1));
136	gforget	1.9	elseif test2&~fldIs3d;
137			fld0=fldIn.fld(:,:,rec_fld(tt));
138			fld1=fldIn.fld(:,:,rec_fld(tt+1));
139	gforget	1.8	elseif test1;
140			fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));
141			fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));
142			elseif test0;
143			fld0=read_bin(fldIn.fil,rec_fld(tt));
144			fld1=read_bin(fldIn.fil,rec_fld(tt+1));
145			else;
146			error(['file not found:' fldIn.fil]);
147			end;
148	gforget	1.5	%
149	gforget	1.1	ndim=length(size(fld0{1}));
150	gforget	1.8	if ndim==2;
151			fld0=fld0.*mygrid.mskC(:,:,1);
152			fld1=fld1.*mygrid.mskC(:,:,1);
153			fldIs3d=0;
154			else;
155			fld0=fld0.*mygrid.mskC;
156			fld1=fld1.*mygrid.mskC;
157			fldIs3d=1;
158			end;
159	gforget	1.1	fld=cat(ndim+1,fld0,fld1);
160			%
161	gforget	1.8	if ~strcmp(method,'bindata');
162			arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
163			if fldIs3d;
164	gforget	1.1	arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
165	gforget	1.8	else;
166			arr2=arr;
167			end;
168			%now linear in time:
169	gforget	1.11	a0=(tim_prof(ii)-tim_fld(tt))/(tim_fld(tt+1)-tim_fld(tt));
170	gforget	1.8	if fldIs3d;
171			a0=a0*ones(1,profIn.nr);
172	gforget	1.5	profOut(ii,:)=(1-a0).arr2(:,:,1)+a0.arr2(:,:,2);
173			else;
174	gforget	1.8	profOut(ii,1)=(1-a0).arr2(:,1)+a0.arr2(:,2);
175	gforget	1.5	end;
176	gforget	1.8	elseif fldIs3d;
177			[prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));
178			FLD=convert2array(fld(:,:,:,1));
179			nk=length(mygrid.RC); kk=ones(1,nk);
180			np=length(ii); pp=ones(np,1);
181			ind2prof=sub2ind(size(FLD),prof_ikk,prof_jkk,pp*[1:nk]);
182			arr=FLD(ind2prof);
183			arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
184			profOut(ii,:)=arr2;
185			else;
186			error('2D field case is missing here');
187			end;
188			%
189			if strcmp(fldIn.tim,'std');
190			profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
191	gforget	1.4	end;
192	gforget	1.8
193			end;%if ~isempty(ii);
194			end;
195
196			if ~fldIs3d;
197			profOut=profOut(:,1);
198	gforget	1.1	end;
199
200			%3) deal with file output
201			if doOutInit;
202			%create a header only file to later append resampled fields
203			prof=profIn;
204			tmp1=fieldnames(prof);
205			nt=length(prof.prof_date);
206			nr=length(prof.prof_depth);
207			for ii=1:length(tmp1);
208			tmp2=prod(size(getfield(prof,tmp1{ii})));
209			if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end;
210			end;
211			MITprof_write(filOut,prof);
212			end;
213
214			if doOut;
215	gforget	1.8	if ~fldIs3d;
216			dims={'iPROF'};
217			else;
218			dims={'iDEPTH','iPROF'};
219			end;
220	gforget	1.1	%add the array itelf
221	gforget	1.8	MITprof_addVar(filOut,fldIn.name,'double',dims,profOut);
222	gforget	1.1
223			%add its attributes
224			nc=ncopen(filOut,'write');
225			ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name);
226			ncaddAtt(nc,fldIn.name,'units',fldIn.units);
227			ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value);
228			ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue);
229			ncclose(nc);
230			end;
231
232			%4) deal with argument output
233			if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end;
234
235