/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/MITprof_resample.m
ViewVC logotype

Diff of /MITgcm_contrib/gael/profilesMatlabProcessing/profiles_misc/MITprof_resample.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.5 by gforget, Sun Oct 11 13:37:09 2015 UTC revision 1.8 by gforget, Fri Jan 29 15:23:31 2016 UTC
# Line 1  Line 1 
1  function [profOut]=MITprof_resample(profIn,fldIn,filOut,method,varargin);  function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);
 %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)  
2  %  %
3  %filFldIn is assumed to be 3D and binary at this point  %[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, FillValue (for nc output), and fld ([] by default).
12    %         If fld is not [] then it is assumed that user already loaded
13    %         the fields from fldIn.fil.
14    %
15    %     fldIn.tim can be set to
16    %         'const' (for time invariant climatology),
17    %         'monclim' (for monthly climatology)
18    %         'monser' (for monthly time series)
19    %         'monloop' (for cyclic monthly time series)
20    %
21    %     method can be set to
22    %         'polygons' (linear in space)
23    %         'bindata' (nearest neighbor in space)
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    %  fldIn.fld=[];
37    %  %
38    %  profOut=MITprof_resample(profIn,fldIn);
39    
40    
41  gcmfaces_global;  gcmfaces_global;
42    
# Line 27  if doOut; doOutInit=isempty(dir(filOut)) Line 46  if doOut; doOutInit=isempty(dir(filOut))
46    
47  if isempty(who('method')); method='polygons'; end;  if isempty(who('method')); method='polygons'; end;
48    
49    %0) check for file types
50    test0=~isempty(dir(fldIn.fil));
51    [PATH,NAME,EXT]=fileparts(fldIn.fil);
52    fil_nc=fullfile(PATH,NAME,[NAME '.0001.nc']);
53    fil_nctiles=fullfile(PATH,NAME,NAME);
54    test1=~isempty(dir(fil_nc));
55    test2=~isempty(fldIn.fld);
56    
57  %1) deal with time line  %1) deal with time line
58  if strcmp(fldIn.tim,'monclim');  if strcmp(fldIn.tim,'monclim');
59      tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];      tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];
# Line 35  if strcmp(fldIn.tim,'monclim'); Line 62  if strcmp(fldIn.tim,'monclim');
62      tim_prof=(profIn.prof_date-tmp2);      tim_prof=(profIn.prof_date-tmp2);
63      tim_prof(tim_prof>365)=365;      tim_prof(tim_prof>365)=365;
64      tim_prof=tim_prof/365*12;%neglecting differences in months length      tim_prof=tim_prof/365*12;%neglecting differences in months length
65  elseif strcmp(fldIn.tim,'monloop');  elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');
66      tmp1=dir(fldIn.fil);      if test1;
67      nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;        eval(['ncload ' fil_nc ' tim']);
68          nt=length(tim);
69        elseif ~test2;
70          %note: 2D case still needs to be treated here ... or via fldIn.is3d ?
71          tmp1=dir(fldIn.fil);
72          nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
73        else;
74          error('gcmfaces input case is missing here');
75        end;
76      tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1;      tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1;
77      tim_fld=datenum(tmp2);      tim_fld=datenum(tmp2);
78      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
79      rec_fld=[nt 1:nt 1];      rec_fld=[nt 1:nt 1];
80      tmp1=datenum([1992 1 1 0 0 0]);      if strcmp(fldIn.tim,'monloop');
81      tmp2=datenum([1992+nt/12 1 1 0 0 0]);;        tmp1=datenum([1992 1 1 0 0 0]);
82      tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2);        tmp2=datenum([1992+nt/12 1 1 0 0 0]);;
83          tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);
84        else;
85          tim_prof=profIn.prof_date;
86        end;
87        %round up tim_prof to prevent interpolation in time:
88        %  tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld;
89        %  tmp4=sum(tmp3>0,2);
90        %  tim_prof=tim_fld(tmp4)';
91  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
92      tim_fld=[1 2]; rec_fld=[1 1];      tim_fld=[1 2]; rec_fld=[1 1];
93      tim_prof=1.5*ones(profIn.np,1);      tim_prof=1.5*ones(profIn.np,1);
# Line 59  profOut=NaN*ones(profIn.np,profIn.nr); Line 102  profOut=NaN*ones(profIn.np,profIn.nr);
102  %2) loop over record pairs  %2) loop over record pairs
103  if ~strcmp(method,'bindata'); gcmfaces_bindata; end;  if ~strcmp(method,'bindata'); gcmfaces_bindata; end;
104  for tt=1:length(rec_fld)-1;  for tt=1:length(rec_fld)-1;
105      tt    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
106      if ~isempty(ii);
107        %tt
108        %
109        if test2;
110          fld0=fldIn.fld(:,:,:,rec_fld(tt));
111          fld1=fldIn.fld(:,:,:,rec_fld(tt+1));
112        elseif test1;
113          fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));
114          fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));
115        elseif test0;
116          fld0=read_bin(fldIn.fil,rec_fld(tt));
117          fld1=read_bin(fldIn.fil,rec_fld(tt+1));
118        else;
119          error(['file not found:' fldIn.fil]);
120        end;
121      %      %
     fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));  
     fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));  
122      ndim=length(size(fld0{1}));      ndim=length(size(fld0{1}));
123        if ndim==2;
124          fld0=fld0.*mygrid.mskC(:,:,1);
125          fld1=fld1.*mygrid.mskC(:,:,1);
126          fldIs3d=0;
127        else;
128          fld0=fld0.*mygrid.mskC;
129          fld1=fld1.*mygrid.mskC;
130          fldIs3d=1;
131        end;
132      fld=cat(ndim+1,fld0,fld1);      fld=cat(ndim+1,fld0,fld1);
133      %      %
134      ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));      if ~strcmp(method,'bindata');
135      if ~isempty(ii);        arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
136        if ~strcmp(method,'bindata');        if fldIs3d;
         arr=gcmfaces_interp(fld,lon(ii),lat(ii),method);  
137          arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);          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);  
138        else;        else;
139          [prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));          arr2=arr;
         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;  
140        end;        end;
141        %            %now linear in time:
142        if strcmp(fldIn.tim,'std');        k0=floor(tim_prof(ii));
143          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));        a0=tim_prof(ii)-k0;
144          if fldIs3d;
145            a0=a0*ones(1,profIn.nr);
146            profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
147          else;
148            profOut(ii,1)=(1-a0).*arr2(:,1)+a0.*arr2(:,2);
149        end;        end;
150        elseif fldIs3d;
151          [prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));
152          FLD=convert2array(fld(:,:,:,1));
153          nk=length(mygrid.RC); kk=ones(1,nk);
154          np=length(ii); pp=ones(np,1);
155          ind2prof=sub2ind(size(FLD),prof_i*kk,prof_j*kk,pp*[1:nk]);
156          arr=FLD(ind2prof);
157          arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
158          profOut(ii,:)=arr2;
159        else;
160          error('2D field case is missing here');
161      end;      end;
162        %    
163        if strcmp(fldIn.tim,'std');
164          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
165        end;
166    
167      end;%if ~isempty(ii);
168    end;
169    
170    if ~fldIs3d;
171      profOut=profOut(:,1);
172  end;  end;
173    
174  %3) deal with file output  %3) deal with file output
# Line 107  if doOutInit; Line 186  if doOutInit;
186  end;  end;
187    
188  if doOut;  if doOut;
189        if ~fldIs3d;
190          dims={'iPROF'};
191        else;
192          dims={'iDEPTH','iPROF'};
193        end;
194      %add the array itelf          %add the array itelf    
195      MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);      MITprof_addVar(filOut,fldIn.name,'double',dims,profOut);
196            
197      %add its attributes      %add its attributes
198      nc=ncopen(filOut,'write');      nc=ncopen(filOut,'write');
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22