/[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.10 by gforget, Fri Feb 5 13:01:47 2016 UTC
# Line 8  function [profOut]=MITprof_resample(prof Line 8  function [profOut]=MITprof_resample(prof
8  %  %
9  %     profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date  %     profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date
10  %     fldIn (structure) should contain: fil, name, long_name, units, tim,  %     fldIn (structure) should contain: fil, name, long_name, units, tim,
11  %         missing_value, and FillValue (for nc output).  %         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  %     fldIn.tim can be set to
16  %         'const' (for time invariant climatology),  %         'const' (for time invariant climatology),
# Line 20  function [profOut]=MITprof_resample(prof Line 22  function [profOut]=MITprof_resample(prof
22  %         'polygons' (linear in space)  %         'polygons' (linear in space)
23  %         'bindata' (nearest neighbor in space)  %         'bindata' (nearest neighbor in space)
24  %  %
 %note to self: add case when fldIn is a gcmfaces field; nctiles input  
 %  
25  % Example:  % Example:
26  %  grid_load; gcmfaces_global; MITprof_global; addpath matlab/;  %  grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
27  %  profIn=idma_float_plot('4900828');  %  profIn=idma_float_plot('4900828');
# Line 33  function [profOut]=MITprof_resample(prof Line 33  function [profOut]=MITprof_resample(prof
33  %  fldIn.tim='monclim';  %  fldIn.tim='monclim';
34  %  fldIn.missing_value=-9999.;  %  fldIn.missing_value=-9999.;
35  %  fldIn.FillValue=-9999.;  %  fldIn.FillValue=-9999.;
36    %  fldIn.fld=[];
37  %  %  %  %
38  %  profOut=MITprof_resample(profIn,fldIn);  %  profOut=MITprof_resample(profIn,fldIn);
39    
# Line 45  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=isfield(fldIn,'fil');
51    test1=0;
52    if test0;
53      test0=~isempty(dir(fldIn.fil));
54      [PATH,NAME,EXT]=fileparts(fldIn.fil);
55      fil_nc=fullfile(PATH,NAME,[NAME '.0001.nc']);
56      fil_nctiles=fullfile(PATH,NAME,NAME);
57      test1=~isempty(dir(fil_nc));
58    end;
59    test2=~isempty(fldIn.fld);
60    
61  %1) deal with time line  %1) deal with time line
62  if strcmp(fldIn.tim,'monclim');  if strcmp(fldIn.tim,'monclim');
63      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 53  if strcmp(fldIn.tim,'monclim'); Line 66  if strcmp(fldIn.tim,'monclim');
66      tim_prof=(profIn.prof_date-tmp2);      tim_prof=(profIn.prof_date-tmp2);
67      tim_prof(tim_prof>365)=365;      tim_prof(tim_prof>365)=365;
68      tim_prof=tim_prof/365*12;%neglecting differences in months length      tim_prof=tim_prof/365*12;%neglecting differences in months length
69  elseif strcmp(fldIn.tim,'monloop');      if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
70      tmp1=dir(fldIn.fil);  elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');
71      nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;      if test1;
72          eval(['ncload ' fil_nc ' tim']);
73          nt=length(tim);
74        elseif ~test2;
75          %note: 2D case still needs to be treated here ... or via fldIn.is3d ?
76          tmp1=dir(fldIn.fil);
77          nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
78        else;
79          ndim=length(size(fldIn.fld{1}));
80          fldIs3d=(ndim==4);
81          nt=size(fldIn.fld{1},ndim);
82        end;
83      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;
84      tim_fld=datenum(tmp2);      tim_fld=datenum(tmp2);
85      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];      tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31];
86      rec_fld=[nt 1:nt 1];      rec_fld=[nt 1:nt 1];
87      tmp1=datenum([1992 1 1 0 0 0]);      if strcmp(fldIn.tim,'monloop');
88      tmp2=datenum([1992+nt/12 1 1 0 0 0]);;        tmp1=datenum([1992 1 1 0 0 0]);
89      tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);        tmp2=datenum([1992+nt/12 1 1 0 0 0]);;
90          tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1);
91        else;
92          tim_prof=profIn.prof_date;
93        end;
94      %round up tim_prof to prevent interpolation in time:      %round up tim_prof to prevent interpolation in time:
95      %  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;
96      %  tmp4=sum(tmp3>0,2);      %  tmp4=sum(tmp3>0,2);
# Line 70  elseif strcmp(fldIn.tim,'monloop'); Line 98  elseif strcmp(fldIn.tim,'monloop');
98  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');  elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
99      tim_fld=[1 2]; rec_fld=[1 1];      tim_fld=[1 2]; rec_fld=[1 1];
100      tim_prof=1.5*ones(profIn.np,1);      tim_prof=1.5*ones(profIn.np,1);
101        if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end;
102  else;  else;
103      error('this case remains to be implemented');      error('this case remains to be implemented');
104  end;  end;
# Line 79  depIn=-mygrid.RC; depOut=profIn.prof_dep Line 108  depIn=-mygrid.RC; depOut=profIn.prof_dep
108  profOut=NaN*ones(profIn.np,profIn.nr);  profOut=NaN*ones(profIn.np,profIn.nr);
109    
110  %2) loop over record pairs  %2) loop over record pairs
111  if ~strcmp(method,'bindata'); gcmfaces_bindata; end;  if strcmp(method,'bindata'); gcmfaces_bindata; end;
112  for tt=1:length(rec_fld)-1;  for tt=1:length(rec_fld)-1;
113      tt    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
114      if ~isempty(ii);
115        %tt
116        %
117        if test2&fldIs3d;
118          fld0=fldIn.fld(:,:,:,rec_fld(tt));
119          fld1=fldIn.fld(:,:,:,rec_fld(tt+1));
120        elseif test2&~fldIs3d;
121          fld0=fldIn.fld(:,:,rec_fld(tt));
122          fld1=fldIn.fld(:,:,rec_fld(tt+1));
123        elseif test1;
124          fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));
125          fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));
126        elseif test0;
127          fld0=read_bin(fldIn.fil,rec_fld(tt));
128          fld1=read_bin(fldIn.fil,rec_fld(tt+1));
129        else;
130          error(['file not found:' fldIn.fil]);
131        end;
132      %      %
     fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));  
     fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));  
133      ndim=length(size(fld0{1}));      ndim=length(size(fld0{1}));
134        if ndim==2;
135          fld0=fld0.*mygrid.mskC(:,:,1);
136          fld1=fld1.*mygrid.mskC(:,:,1);
137          fldIs3d=0;
138        else;
139          fld0=fld0.*mygrid.mskC;
140          fld1=fld1.*mygrid.mskC;
141          fldIs3d=1;
142        end;
143      fld=cat(ndim+1,fld0,fld1);      fld=cat(ndim+1,fld0,fld1);
144      %      %
145      ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));      if ~strcmp(method,'bindata');
146      if ~isempty(ii);        arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);
147        if ~strcmp(method,'bindata');        if fldIs3d;
         arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method);  
148          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);  
149        else;        else;
150          [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;  
151        end;        end;
152        %            %now linear in time:
153        if strcmp(fldIn.tim,'std');        k0=floor(tim_prof(ii));
154          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));        a0=tim_prof(ii)-k0;
155          if fldIs3d;
156            a0=a0*ones(1,profIn.nr);
157            profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
158          else;
159            profOut(ii,1)=(1-a0).*arr2(:,1)+a0.*arr2(:,2);
160        end;        end;
161        elseif fldIs3d;
162          [prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii));
163          FLD=convert2array(fld(:,:,:,1));
164          nk=length(mygrid.RC); kk=ones(1,nk);
165          np=length(ii); pp=ones(np,1);
166          ind2prof=sub2ind(size(FLD),prof_i*kk,prof_j*kk,pp*[1:nk]);
167          arr=FLD(ind2prof);
168          arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
169          profOut(ii,:)=arr2;
170        else;
171          error('2D field case is missing here');
172        end;
173        %    
174        if strcmp(fldIn.tim,'std');
175          profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
176      end;      end;
177    
178      end;%if ~isempty(ii);
179    end;
180    
181    if ~fldIs3d;
182      profOut=profOut(:,1);
183  end;  end;
184    
185  %3) deal with file output  %3) deal with file output
# Line 129  if doOutInit; Line 197  if doOutInit;
197  end;  end;
198    
199  if doOut;  if doOut;
200        if ~fldIs3d;
201          dims={'iPROF'};
202        else;
203          dims={'iDEPTH','iPROF'};
204        end;
205      %add the array itelf          %add the array itelf    
206      MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut);      MITprof_addVar(filOut,fldIn.name,'double',dims,profOut);
207            
208      %add its attributes      %add its attributes
209      nc=ncopen(filOut,'write');      nc=ncopen(filOut,'write');
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.10
  ViewVC Help
Powered by ViewVC 1.1.22