| 1 |
function [profOut]=MITprof_resample(profIn,fldIn,filOut,method); |
| 2 |
%[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method); |
| 3 |
% |
| 4 |
% 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; |
| 47 |
|
| 48 |
doOut=~isempty(who('filOut')); doOutInit=false; |
| 49 |
if doOut; doOut=~isempty(filOut); end; |
| 50 |
if doOut; doOutInit=isempty(dir(filOut)); end; |
| 51 |
|
| 52 |
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 |
| 73 |
if strcmp(fldIn.tim,'monclim'); |
| 74 |
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); |
| 80 |
tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]); |
| 81 |
tim_prof=(profIn.prof_date-tmp2); |
| 82 |
tim_prof(tim_prof>365)=365; |
| 83 |
% |
| 84 |
if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end; |
| 85 |
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); |
| 117 |
if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end; |
| 118 |
else; |
| 119 |
error('this case remains to be implemented'); |
| 120 |
end; |
| 121 |
|
| 122 |
lon=profIn.prof_lon; lat=profIn.prof_lat; |
| 123 |
depIn=-mygrid.RC; depOut=profIn.prof_depth; |
| 124 |
profOut=NaN*ones(profIn.np,profIn.nr); |
| 125 |
|
| 126 |
%2) loop over record pairs |
| 127 |
if strcmp(method,'bindata'); gcmfaces_bindata; end; |
| 128 |
for tt=1:length(rec_fld)-1; |
| 129 |
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 |
% |
| 149 |
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); |
| 160 |
% |
| 161 |
if ~strcmp(method,'bindata'); |
| 162 |
arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method); |
| 163 |
if fldIs3d; |
| 164 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 165 |
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 |
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; |
| 199 |
|
| 200 |
%3) deal with file output |
| 201 |
if doOutInit; |
| 202 |
%create a header only file to later append resampled fields |
| 203 |
prof=profIn; |
| 204 |
tmp1=fieldnames(prof); |
| 205 |
nt=length(prof.prof_date); |
| 206 |
nr=length(prof.prof_depth); |
| 207 |
for ii=1:length(tmp1); |
| 208 |
tmp2=prod(size(getfield(prof,tmp1{ii}))); |
| 209 |
if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end; |
| 210 |
end; |
| 211 |
MITprof_write(filOut,prof); |
| 212 |
end; |
| 213 |
|
| 214 |
if doOut; |
| 215 |
if ~fldIs3d; |
| 216 |
dims={'iPROF'}; |
| 217 |
else; |
| 218 |
dims={'iDEPTH','iPROF'}; |
| 219 |
end; |
| 220 |
%add the array itelf |
| 221 |
MITprof_addVar(filOut,fldIn.name,'double',dims,profOut); |
| 222 |
|
| 223 |
%add its attributes |
| 224 |
nc=ncopen(filOut,'write'); |
| 225 |
ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name); |
| 226 |
ncaddAtt(nc,fldIn.name,'units',fldIn.units); |
| 227 |
ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value); |
| 228 |
ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue); |
| 229 |
ncclose(nc); |
| 230 |
end; |
| 231 |
|
| 232 |
%4) deal with argument output |
| 233 |
if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end; |
| 234 |
|
| 235 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph