profilesMatlabProcessing/profiles_misc/MITprof_resample.m

function [profOut]=MITprof_resample(profIn,fldIn,filOut,method,varargin);
%object:    resample a set of fields in file filFldIn with specified time
%           line timeIn to the positions of profIn and add to file filOut
%inputs:    profIn is a gcmfaces field (nan-masked; up to N3,N4 dimensions)
%           fldIn is a description of the fields being resampled including
%               the corresponding file name and additional specs :
%               fldIn.name, fldIn.long_name, fldIn.units, fldIn.fil (file 
%               name) and fldIn.tim (time axis specification). Supported 
%               fldIn.tim spec: 'const' (for time invariant climatology),
%               'monclim' (for monthly climatology), 'monser' (for monthly
%               time series)
%           filOut is the output MITprof file name (if un-specified
%               the resul may only be returned as a function argument)
%           method may be 'polygons' (or 'TriScatteredInterp' ... via 
%              gcmfaces_interp_2d in a loop ... to be implemented later)
%outputs:   profOut is the MITprof structure where the interpolated values
%               were appended to profIn (if un-specified the result
%               may only be returned to output file)
%
%filFldIn is assumed to be 3D and binary at this point

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,'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
for tt=1:length(rec_fld)-1;
    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);
        arr=gcmfaces_interp(fld,lon(ii),lat(ii),method);
        arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
        if strcmp(fldIn.tim,'std');
            arr2=arr2.*randn(size(arr2));
        end;
    end;        
    %    
    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);
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	function [profOut]=MITprof_resample(profIn,fldIn,filOut,method,varargin);
2	%object: resample a set of fields in file filFldIn with specified time
3	% line timeIn to the positions of profIn and add to file filOut
4	%inputs: profIn is a gcmfaces field (nan-masked; up to N3,N4 dimensions)
5	% fldIn is a description of the fields being resampled including
6	% the corresponding file name and additional specs :
7	% fldIn.name, fldIn.long_name, fldIn.units, fldIn.fil (file
8	% name) and fldIn.tim (time axis specification). Supported
9	% fldIn.tim spec: 'const' (for time invariant climatology),
10	% 'monclim' (for monthly climatology), 'monser' (for monthly
11	% time series)
12	% filOut is the output MITprof file name (if un-specified
13	% the resul may only be returned as a function argument)
14	% method may be 'polygons' (or 'TriScatteredInterp' ... via
15	% gcmfaces_interp_2d in a loop ... to be implemented later)
16	%outputs: profOut is the MITprof structure where the interpolated values
17	% were appended to profIn (if un-specified the result
18	% may only be returned to output file)
19	%
20	%filFldIn is assumed to be 3D and binary at this point
21
22	gcmfaces_global;
23
24	doOut=~isempty(who('filOut')); doOutInit=false;
25	if doOut; doOut=~isempty(filOut); end;
26	if doOut; doOutInit=isempty(dir(filOut)); end;
27
28	if isempty(who('method')); method='polygons'; end;
29
30	%1) deal with time line
31	if strcmp(fldIn.tim,'monclim');
32	tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];
33	tmp1=datevec(profIn.prof_date);
34	tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
35	tim_prof=(profIn.prof_date-tmp2);
36	tim_prof(tim_prof>365)=365;
37	tim_prof=tim_prof/365*12;%neglecting differences in months length
38	elseif strcmp(fldIn.tim,'const')\|strcmp(fldIn.tim,'std');
39	tim_fld=[1 2]; rec_fld=[1 1];
40	tim_prof=1.5*ones(profIn.np,1);
41	else;
42	error('this case remains to be implemented');
43	end;
44
45	lon=profIn.prof_lon; lat=profIn.prof_lat;
46	depIn=-mygrid.RC; depOut=profIn.prof_depth;
47	profOut=NaN*ones(profIn.np,profIn.nr);
48
49	%2) loop over record pairs
50	for tt=1:length(rec_fld)-1;
51	fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));
52	fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));
53	ndim=length(size(fld0{1}));
54	fld=cat(ndim+1,fld0,fld1);
55	%
56	ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
57	if ~isempty(ii);
58	arr=gcmfaces_interp(fld,lon(ii),lat(ii),method);
59	arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
60	if strcmp(fldIn.tim,'std');
61	arr2=arr2.*randn(size(arr2));
62	end;
63	end;
64	%
65	k0=floor(tim_prof(ii)); k1=k0+1;
66	a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr);
67	profOut(ii,:)=(1-a0).arr2(:,:,1)+a0.arr2(:,:,2);
68	end;
69
70	%3) deal with file output
71	if doOutInit;
72	%create a header only file to later append resampled fields
73	prof=profIn;
74	tmp1=fieldnames(prof);
75	nt=length(prof.prof_date);
76	nr=length(prof.prof_depth);
77	for ii=1:length(tmp1);
78	tmp2=prod(size(getfield(prof,tmp1{ii})));
79	if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end;
80	end;
81	MITprof_write(filOut,prof);
82	end;
83
84	if doOut;
85	%add the array itelf
86	MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);
87
88	%add its attributes
89	nc=ncopen(filOut,'write');
90	ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name);
91	ncaddAtt(nc,fldIn.name,'units',fldIn.units);
92	ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value);
93	ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue);
94	ncclose(nc);
95	end;
96
97	%4) deal with argument output
98	if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end;
99
100