| 1 |
function [profOut]=MITprof_resample(profIn,fldIn,filOut,method); |
| 2 |
% |
| 3 |
%[profOut]=MITPROF_RESAMPLE(profIn,fldIn,filOut,method); |
| 4 |
% |
| 5 |
% resamples a set of fields (fldIn) to profile 3D locations (profIn) |
| 6 |
% and output the result to nc file (filOut) and memory (profOut). |
| 7 |
% using a chosen interpolation method (method). |
| 8 |
% |
| 9 |
% profIn (structure) should contain: prof_depth, prof_lon, prof_lat, prof_date |
| 10 |
% fldIn (structure) should contain: fil, name, long_name, units, tim, |
| 11 |
% 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 |
| 16 |
% 'const' (for time invariant climatology), |
| 17 |
% 'monclim' (for monthly climatology) |
| 18 |
% 'monser' (for monthly time series) |
| 19 |
% 'monloop' (for cyclic monthly time series) |
| 20 |
% |
| 21 |
% method can be set to |
| 22 |
% 'polygons' (linear in space) |
| 23 |
% 'bindata' (nearest neighbor in space) |
| 24 |
% |
| 25 |
% Example: |
| 26 |
% grid_load; gcmfaces_global; MITprof_global; addpath matlab/; |
| 27 |
% profIn=idma_float_plot('4900828'); |
| 28 |
% % |
| 29 |
% fldIn.fil=fullfile(myenv.MITprof_climdir,filesep,'T_OWPv1_M_eccollc_90x50.bin'); |
| 30 |
% fldIn.name='prof_Towp'; |
| 31 |
% fldIn.long_name='pot. temp. estimate (OCCA-WOA-PHC combination)'; |
| 32 |
% fldIn.units='degree C'; |
| 33 |
% fldIn.tim='monclim'; |
| 34 |
% fldIn.missing_value=-9999.; |
| 35 |
% fldIn.FillValue=-9999.; |
| 36 |
% fldIn.fld=[]; |
| 37 |
% % |
| 38 |
% profOut=MITprof_resample(profIn,fldIn); |
| 39 |
|
| 40 |
|
| 41 |
gcmfaces_global; |
| 42 |
|
| 43 |
doOut=~isempty(who('filOut')); doOutInit=false; |
| 44 |
if doOut; doOut=~isempty(filOut); end; |
| 45 |
if doOut; doOutInit=isempty(dir(filOut)); end; |
| 46 |
|
| 47 |
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 |
| 62 |
if strcmp(fldIn.tim,'monclim'); |
| 63 |
tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1]; |
| 64 |
tmp1=datevec(profIn.prof_date); |
| 65 |
tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]); |
| 66 |
tim_prof=(profIn.prof_date-tmp2); |
| 67 |
tim_prof(tim_prof>365)=365; |
| 68 |
tim_prof=tim_prof/365*12;%neglecting differences in months length |
| 69 |
if test2; fldIs3d=(length(size(fldIn.fld{1}))==4); end; |
| 70 |
elseif strcmp(fldIn.tim,'monloop')|strcmp(fldIn.tim,'monser'); |
| 71 |
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; |
| 84 |
tim_fld=datenum(tmp2); |
| 85 |
tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31]; |
| 86 |
rec_fld=[nt 1:nt 1]; |
| 87 |
if strcmp(fldIn.tim,'monloop'); |
| 88 |
tmp1=datenum([1992 1 1 0 0 0]); |
| 89 |
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: |
| 95 |
% tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld; |
| 96 |
% tmp4=sum(tmp3>0,2); |
| 97 |
% tim_prof=tim_fld(tmp4)'; |
| 98 |
elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std'); |
| 99 |
tim_fld=[1 2]; rec_fld=[1 1]; |
| 100 |
tim_prof=1.5*ones(profIn.np,1); |
| 101 |
if test2; fldIs3d=(length(size(fldIn.fld{1}))==3); end; |
| 102 |
else; |
| 103 |
error('this case remains to be implemented'); |
| 104 |
end; |
| 105 |
|
| 106 |
lon=profIn.prof_lon; lat=profIn.prof_lat; |
| 107 |
depIn=-mygrid.RC; depOut=profIn.prof_depth; |
| 108 |
profOut=NaN*ones(profIn.np,profIn.nr); |
| 109 |
|
| 110 |
%2) loop over record pairs |
| 111 |
if strcmp(method,'bindata'); gcmfaces_bindata; end; |
| 112 |
for tt=1:length(rec_fld)-1; |
| 113 |
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 |
% |
| 133 |
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); |
| 144 |
% |
| 145 |
if ~strcmp(method,'bindata'); |
| 146 |
arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method); |
| 147 |
if fldIs3d; |
| 148 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 149 |
else; |
| 150 |
arr2=arr; |
| 151 |
end; |
| 152 |
%now linear in time: |
| 153 |
k0=floor(tim_prof(ii)); |
| 154 |
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; |
| 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; |
| 177 |
|
| 178 |
end;%if ~isempty(ii); |
| 179 |
end; |
| 180 |
|
| 181 |
if ~fldIs3d; |
| 182 |
profOut=profOut(:,1); |
| 183 |
end; |
| 184 |
|
| 185 |
%3) deal with file output |
| 186 |
if doOutInit; |
| 187 |
%create a header only file to later append resampled fields |
| 188 |
prof=profIn; |
| 189 |
tmp1=fieldnames(prof); |
| 190 |
nt=length(prof.prof_date); |
| 191 |
nr=length(prof.prof_depth); |
| 192 |
for ii=1:length(tmp1); |
| 193 |
tmp2=prod(size(getfield(prof,tmp1{ii}))); |
| 194 |
if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end; |
| 195 |
end; |
| 196 |
MITprof_write(filOut,prof); |
| 197 |
end; |
| 198 |
|
| 199 |
if doOut; |
| 200 |
if ~fldIs3d; |
| 201 |
dims={'iPROF'}; |
| 202 |
else; |
| 203 |
dims={'iDEPTH','iPROF'}; |
| 204 |
end; |
| 205 |
%add the array itelf |
| 206 |
MITprof_addVar(filOut,fldIn.name,'double',dims,profOut); |
| 207 |
|
| 208 |
%add its attributes |
| 209 |
nc=ncopen(filOut,'write'); |
| 210 |
ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name); |
| 211 |
ncaddAtt(nc,fldIn.name,'units',fldIn.units); |
| 212 |
ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value); |
| 213 |
ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue); |
| 214 |
ncclose(nc); |
| 215 |
end; |
| 216 |
|
| 217 |
%4) deal with argument output |
| 218 |
if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end; |
| 219 |
|
| 220 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph