/[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.7 by gforget, Thu Jan 28 23:18:26 2016 UTC revision 1.11 by gforget, Fri Mar 17 18:06:22 2017 UTC
# Line 1  Line 1 
1  function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);  function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);
 %  
2  %[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);  %[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);
3  %  %
4  %     resamples a set of fields (fldIn) to profile 3D locations (profIn)  %  resamples a set of fields (specified in fldIn) to profile locations
5  %     and output the result to nc file (filOut) and memory (profOut).  %  (specified in profIn) and output the result either to memory
6  %     using a chosen interpolation method (method).  %  (by default) or to a netcdf file (if filOut is specified) based on
7  %  %  on pre-defined interpolation method ('polygons' by default)
8  %     profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date  %
9  %     fldIn (structure) should contain: fil, name, long_name, units, tim,  %     profIn (structure) should contain: prof_depth, prof_lon, prof_lat,
10  %         missing_value, and FillValue (for nc output).  %         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  %  %
19  %     fldIn.tim can be set to  %     fldIn.tim must be set to one of the following values:
20  %         'const' (for time invariant climatology),  %         'const' (for time invariant climatology),
21  %         'monclim' (for monthly climatology)  %         'monclim' (for monthly climatology)
22  %         'monser' (for monthly time series)  %         'monser' (for monthly time series)
23  %         'monloop' (for cyclic monthly time series)  %         'monloop' (for cyclic monthly time series)
24  %  %
25  %     method can be set to  %     method ('polygons' by default) can be specified as
26  %         'polygons' (linear in space)  %         'polygons' (linear in space)
27  %         'bindata' (nearest neighbor in space)  %         'bindata' (nearest neighbor in space)
28  %  %
29  %note to self: add case when fldIn is a gcmfaces field; nctiles input  %  Example: (should be revisited)
30  %  %
31  % Example:  %     grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
32  %  grid_load; gcmfaces_global; MITprof_global; addpath matlab/;  %     profIn=idma_float_plot('4900828');
33  %  profIn=idma_float_plot('4900828');  %     %
34  %  %  %     fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');
35  %  fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');  %     fldIn.name='prof_Towp';
36  %  fldIn.name='prof_Towp';  %     fldIn.tim='monclim';
37  %  fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';  %     %fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
38  %  fldIn.units='degree C';  %     %fldIn.units='degree C';
39  %  fldIn.tim='monclim';  %     %fldIn.missing_value=-9999.;
40  %  fldIn.missing_value=-9999.;  %     %fldIn.FillValue=-9999.;
41  %  fldIn.FillValue=-9999.;  %     %fldIn.fld=[];
42  %  %  %     %
43  %  profOut=MITprof_resample(profIn,fldIn);  %     profOut=MITprof_resample(profIn,fldIn);
44    
45    
46  gcmfaces_global;  gcmfaces_global;
# Line 45  if doOut; doOutInit=isempty(dir(filOut)) Line 51  if doOut; doOutInit=isempty(dir(filOut))
51    
52  if isempty(who('method')); method='polygons'; end;  if isempty(who('method')); method='polygons'; end;
53    
54    %0) check for input types
55    %   test0=1 <-> binary
56    %   test1=1 <-> nctiles
57    %   test2=1 <-> readily available fldIn.fld
58    test0=isfield(fldIn,'fil');
59    %
60    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    %
69    if ~isfield(fldIn,'fld'); fldIn.fld=[]; end;
70    test2=~isempty(fldIn.fld);
71    
72  %1) deal with time line  %1) deal with time line
73  if strcmp(fldIn.tim,'monclim');  if strcmp(fldIn.tim,'monclim');
74      tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1];      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      tmp1=datevec(profIn.prof_date);      tmp1=datevec(profIn.prof_date);
80      tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);      tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]);
81      tim_prof=(profIn.prof_date-tmp2);      tim_prof=(profIn.prof_date-tmp2);
82      tim_prof(tim_prof>365)=365;      tim_prof(tim_prof>365)=365;
83      tim_prof=tim_prof/365*12;%neglecting differences in months length      %
84  elseif strcmp(fldIn.tim,'monloop');      if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
85      tmp1=dir(fldIn.fil);  elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');
86      nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;      if test1;
87          eval(['ncload ' fil_nc ' tim']);
88          nt=length(tim);
89        elseif ~test2;
90          warning('Here it is assumed that fldIn.fil contains 3D fields');
91          %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          ndim=length(size(fldIn.fld{1}));
96          fldIs3d=(ndim==4);
97          nt=size(fldIn.fld{1},ndim);
98        end;
99      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;
100      tim_fld=datenum(tmp2);      tim_fld=datenum(tmp2);
101      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
102      rec_fld=[nt 1:nt 1];      rec_fld=[nt 1:nt 1];
103      tmp1=datenum([1992 1 1 0 0 0]);      if strcmp(fldIn.tim,'monloop');
104      tmp2=datenum([1992+nt/12 1 1 0 0 0]);;        tmp1=datenum([1992 1 1 0 0 0]);
105      tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);        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      %round up tim_prof to prevent interpolation in time:      %round up tim_prof to prevent interpolation in time:
111      %  tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld;      %  tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld;
112      %  tmp4=sum(tmp3>0,2);      %  tmp4=sum(tmp3>0,2);
# Line 70  elseif strcmp(fldIn.tim,'monloop'); Line 114  elseif strcmp(fldIn.tim,'monloop');
114  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
115      tim_fld=[1 2]; rec_fld=[1 1];      tim_fld=[1 2]; rec_fld=[1 1];
116      tim_prof=1.5*ones(profIn.np,1);      tim_prof=1.5*ones(profIn.np,1);
117        if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end;
118  else;  else;
119      error('this case remains to be implemented');      error('this case remains to be implemented');
120  end;  end;
# Line 79  depIn=-mygrid.RC; depOut=profIn.prof_dep Line 124  depIn=-mygrid.RC; depOut=profIn.prof_dep
124  profOut=NaN*ones(profIn.np,profIn.nr);  profOut=NaN*ones(profIn.np,profIn.nr);
125    
126  %2) loop over record pairs  %2) loop over record pairs
127  if ~strcmp(method,'bindata'); gcmfaces_bindata; end;  if strcmp(method,'bindata'); gcmfaces_bindata; end;
128  for tt=1:length(rec_fld)-1;  for tt=1:length(rec_fld)-1;
129      tt    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
130      if ~isempty(ii);
131        %tt
132        %
133        if test2&fldIs3d;
134          fld0=fldIn.fld(:,:,:,rec_fld(tt));
135          fld1=fldIn.fld(:,:,:,rec_fld(tt+1));
136        elseif test2&~fldIs3d;
137          fld0=fldIn.fld(:,:,rec_fld(tt));
138          fld1=fldIn.fld(:,:,rec_fld(tt+1));
139        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      %      %
     fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));  
     fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));  
149      ndim=length(size(fld0{1}));      ndim=length(size(fld0{1}));
150        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      fld=cat(ndim+1,fld0,fld1);      fld=cat(ndim+1,fld0,fld1);
160      %      %
161      ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));      if ~strcmp(method,'bindata');
162      if ~isempty(ii);        arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
163        if ~strcmp(method,'bindata');        if fldIs3d;
         arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);  
164          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);  
165        else;        else;
166          [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;  
167        end;        end;
168        %            %now linear in time:
169        if strcmp(fldIn.tim,'std');        a0=(tim_prof(ii)-tim_fld(tt))/(tim_fld(tt+1)-tim_fld(tt));
170          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));        if fldIs3d;
171            a0=a0*ones(1,profIn.nr);
172            profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
173          else;
174            profOut(ii,1)=(1-a0).*arr2(:,1)+a0.*arr2(:,2);
175        end;        end;
176        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_i*kk,prof_j*kk,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;      end;
188        %    
189        if strcmp(fldIn.tim,'std');
190          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
191        end;
192    
193      end;%if ~isempty(ii);
194    end;
195    
196    if ~fldIs3d;
197      profOut=profOut(:,1);
198  end;  end;
199    
200  %3) deal with file output  %3) deal with file output
# Line 129  if doOutInit; Line 212  if doOutInit;
212  end;  end;
213    
214  if doOut;  if doOut;
215        if ~fldIs3d;
216          dims={'iPROF'};
217        else;
218          dims={'iDEPTH','iPROF'};
219        end;
220      %add the array itelf          %add the array itelf    
221      MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);      MITprof_addVar(filOut,fldIn.name,'double',dims,profOut);
222            
223      %add its attributes      %add its attributes
224      nc=ncopen(filOut,'write');      nc=ncopen(filOut,'write');
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22