/[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.2 by gforget, Sat Apr 18 12:54:20 2015 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,varargin);  function [profOut]=MITprof_resample(profIn,fldIn,filOut,method);
2  %object:    resample a set of fields in file filFldIn with specified time  %[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method);
 %           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)  
3  %  %
4  %filFldIn is assumed to be 3D and binary at this point  %  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    %
19    %     fldIn.tim must be set to one of the following values:
20    %         '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    %     method ('polygons' by default) can be specified as
26    %         'polygons' (linear in space)
27    %         'bindata' (nearest neighbor in space)
28    %
29    %  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    
45    
46  gcmfaces_global;  gcmfaces_global;
47    
# Line 27  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,'const');      if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
85      tim_fld=[1 2]; rec_fld=[1 2];  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          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;
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        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        %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;
112        %  tmp4=sum(tmp3>0,2);
113        %  tim_prof=tim_fld(tmp4)';
114    elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std');
115        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 47  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;
128  for tt=1:length(rec_fld)-1;  for tt=1:length(rec_fld)-1;
129      fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt));    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
130      fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1));    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        %
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          arr=gcmfaces_interp(fld,lon(ii),lat(ii),'polygons');        if fldIs3d;
164          arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);          arr2=gcmfaces_interp_1d(2,depIn,arr,depOut);
165      end;                else;
166            arr2=arr;
167          end;
168          %now linear in time:
169          a0=(tim_prof(ii)-tim_fld(tt))/(tim_fld(tt+1)-tim_fld(tt));
170          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;
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;
188      %          %    
189      k0=floor(tim_prof(ii)); k1=k0+1;      if strcmp(fldIn.tim,'std');
190      a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr);        profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:)));
191      profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);      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 79  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.2  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22