/[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.8 by gforget, Fri Jan 29 15:23:31 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, FillValue (for nc output), and fld ([] by default).  %         and prof_date (serial date number from datenum.m)
11  %         If fld is not [] then it is assumed that user already loaded  %
12  %         the fields from fldIn.fil.  %     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  % Example:  %  Example: (should be revisited)
30  %  grid_load; gcmfaces_global; MITprof_global; addpath matlab/;  %
31  %  profIn=idma_float_plot('4900828');  %     grid_load; gcmfaces_global; MITprof_global; addpath matlab/;
32  %  %  %     profIn=idma_float_plot('4900828');
33  %  fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');  %     %
34  %  fldIn.name='prof_Towp';  %     fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin');
35  %  fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';  %     fldIn.name='prof_Towp';
36  %  fldIn.units='degree C';  %     fldIn.tim='monclim';
37  %  fldIn.tim='monclim';  %     %fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)';
38  %  fldIn.missing_value=-9999.;  %     %fldIn.units='degree C';
39  %  fldIn.FillValue=-9999.;  %     %fldIn.missing_value=-9999.;
40  %  fldIn.fld=[];  %     %fldIn.FillValue=-9999.;
41  %  %  %     %fldIn.fld=[];
42  %  profOut=MITprof_resample(profIn,fldIn);  %     %
43    %     profOut=MITprof_resample(profIn,fldIn);
44    
45    
46  gcmfaces_global;  gcmfaces_global;
# Line 46  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 file types  %0) check for input types
55  test0=~isempty(dir(fldIn.fil));  %   test0=1 <-> binary
56  [PATH,NAME,EXT]=fileparts(fldIn.fil);  %   test1=1 <-> nctiles
57  fil_nc=fullfile(PATH,NAME,[NAME '.0001.nc']);  %   test2=1 <-> readily available fldIn.fld
58  fil_nctiles=fullfile(PATH,NAME,NAME);  test0=isfield(fldIn,'fil');
59  test1=~isempty(dir(fil_nc));  %
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);  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        if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end;
85  elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');  elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser');
86      if test1;      if test1;
87        eval(['ncload ' fil_nc ' tim']);        eval(['ncload ' fil_nc ' tim']);
88        nt=length(tim);        nt=length(tim);
89      elseif ~test2;      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 ?        %note: 2D case still needs to be treated here ... or via fldIn.is3d ?
92        tmp1=dir(fldIn.fil);        tmp1=dir(fldIn.fil);
93        nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;        nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4;
94      else;      else;
95        error('gcmfaces input case is missing here');        ndim=length(size(fldIn.fld{1}));
96          fldIs3d=(ndim==4);
97          nt=size(fldIn.fld{1},ndim);
98      end;      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);
# Line 91  elseif strcmp(fldIn.tim,'monloop')|strcm Line 114  elseif strcmp(fldIn.tim,'monloop')|strcm
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 100  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    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));    ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1));
130    if ~isempty(ii);    if ~isempty(ii);
131      %tt      %tt
132      %      %
133      if test2;      if test2&fldIs3d;
134        fld0=fldIn.fld(:,:,:,rec_fld(tt));        fld0=fldIn.fld(:,:,:,rec_fld(tt));
135        fld1=fldIn.fld(:,:,:,rec_fld(tt+1));        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;      elseif test1;
140        fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));        fld0=read_nctiles(fil_nctiles,NAME,rec_fld(tt));
141        fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));        fld1=read_nctiles(fil_nctiles,NAME,rec_fld(tt+1));
# Line 139  for tt=1:length(rec_fld)-1; Line 166  for tt=1:length(rec_fld)-1;
166          arr2=arr;          arr2=arr;
167        end;        end;
168        %now linear in time:        %now linear in time:
169        k0=floor(tim_prof(ii));        a0=(tim_prof(ii)-tim_fld(tt))/(tim_fld(tt+1)-tim_fld(tt));
       a0=tim_prof(ii)-k0;  
170        if fldIs3d;        if fldIs3d;
171          a0=a0*ones(1,profIn.nr);          a0=a0*ones(1,profIn.nr);
172          profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);          profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2);
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22