| 1 |
function [profOut]=MITprof_resample(profIn,fldIn,filOut,method,varargin); |
| 2 |
%object: resample a set of fields in file filFldIn with specified time |
| 3 |
% line timeIn to the positions of profIn and add to file filOut |
| 4 |
%inputs: profIn is a gcmfaces field (nan-masked; up to N3,N4 dimensions) |
| 5 |
% fldIn is a description of the fields being resampled including |
| 6 |
% the corresponding file name and additional specs : |
| 7 |
% fldIn.name, fldIn.long_name, fldIn.units, fldIn.fil (file |
| 8 |
% name) and fldIn.tim (time axis specification). Supported |
| 9 |
% fldIn.tim spec: 'const' (for time invariant climatology), |
| 10 |
% 'monclim' (for monthly climatology), 'monser' (for monthly |
| 11 |
% time series) |
| 12 |
% filOut is the output MITprof file name (if un-specified |
| 13 |
% the resul may only be returned as a function argument) |
| 14 |
% method may be 'polygons' (or 'TriScatteredInterp' ... via |
| 15 |
% gcmfaces_interp_2d in a loop ... to be implemented later) |
| 16 |
%outputs: profOut is the MITprof structure where the interpolated values |
| 17 |
% were appended to profIn (if un-specified the result |
| 18 |
% may only be returned to output file) |
| 19 |
% |
| 20 |
%filFldIn is assumed to be 3D and binary at this point |
| 21 |
|
| 22 |
gcmfaces_global; |
| 23 |
|
| 24 |
doOut=~isempty(who('filOut')); doOutInit=false; |
| 25 |
if doOut; doOut=~isempty(filOut); end; |
| 26 |
if doOut; doOutInit=isempty(dir(filOut)); end; |
| 27 |
|
| 28 |
if isempty(who('method')); method='polygons'; end; |
| 29 |
|
| 30 |
%1) deal with time line |
| 31 |
if strcmp(fldIn.tim,'monclim'); |
| 32 |
tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1]; |
| 33 |
tmp1=datevec(profIn.prof_date); |
| 34 |
tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]); |
| 35 |
tim_prof=(profIn.prof_date-tmp2); |
| 36 |
tim_prof(tim_prof>365)=365; |
| 37 |
tim_prof=tim_prof/365*12;%neglecting differences in months length |
| 38 |
elseif strcmp(fldIn.tim,'monloop'); |
| 39 |
tmp1=dir(fldIn.fil); |
| 40 |
nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4; |
| 41 |
tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1; |
| 42 |
tim_fld=datenum(tmp2); |
| 43 |
tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31]; |
| 44 |
rec_fld=[nt 1:nt 1]; |
| 45 |
tmp1=datenum([1992 1 1 0 0 0]); |
| 46 |
tmp2=datenum([1992+nt/12 1 1 0 0 0]);; |
| 47 |
tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2); |
| 48 |
elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std'); |
| 49 |
tim_fld=[1 2]; rec_fld=[1 1]; |
| 50 |
tim_prof=1.5*ones(profIn.np,1); |
| 51 |
else; |
| 52 |
error('this case remains to be implemented'); |
| 53 |
end; |
| 54 |
|
| 55 |
lon=profIn.prof_lon; lat=profIn.prof_lat; |
| 56 |
depIn=-mygrid.RC; depOut=profIn.prof_depth; |
| 57 |
profOut=NaN*ones(profIn.np,profIn.nr); |
| 58 |
|
| 59 |
%2) loop over record pairs |
| 60 |
if ~strcmp(method,'bindata'); gcmfaces_bindata; end; |
| 61 |
for tt=1:length(rec_fld)-1; |
| 62 |
tt |
| 63 |
% |
| 64 |
fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt)); |
| 65 |
fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1)); |
| 66 |
ndim=length(size(fld0{1})); |
| 67 |
fld=cat(ndim+1,fld0,fld1); |
| 68 |
% |
| 69 |
ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1)); |
| 70 |
if ~isempty(ii); |
| 71 |
if ~strcmp(method,'bindata'); |
| 72 |
arr=gcmfaces_interp(fld,lon(ii),lat(ii),method); |
| 73 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 74 |
%now linear in time: |
| 75 |
k0=floor(tim_prof(ii)); k1=k0+1; |
| 76 |
a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr); |
| 77 |
profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2); |
| 78 |
else; |
| 79 |
[prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii)); |
| 80 |
FLD=convert2array(fld(:,:,:,1)); |
| 81 |
nk=length(mygrid.RC); kk=ones(1,nk); |
| 82 |
np=length(ii); pp=ones(np,1); |
| 83 |
ind2prof=sub2ind(size(FLD),prof_i*kk,prof_j*kk,pp*[1:nk]); |
| 84 |
arr=FLD(ind2prof); |
| 85 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 86 |
profOut(ii,:)=arr2; |
| 87 |
end; |
| 88 |
% |
| 89 |
if strcmp(fldIn.tim,'std'); |
| 90 |
profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:))); |
| 91 |
end; |
| 92 |
end; |
| 93 |
end; |
| 94 |
|
| 95 |
%3) deal with file output |
| 96 |
if doOutInit; |
| 97 |
%create a header only file to later append resampled fields |
| 98 |
prof=profIn; |
| 99 |
tmp1=fieldnames(prof); |
| 100 |
nt=length(prof.prof_date); |
| 101 |
nr=length(prof.prof_depth); |
| 102 |
for ii=1:length(tmp1); |
| 103 |
tmp2=prod(size(getfield(prof,tmp1{ii}))); |
| 104 |
if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end; |
| 105 |
end; |
| 106 |
MITprof_write(filOut,prof); |
| 107 |
end; |
| 108 |
|
| 109 |
if doOut; |
| 110 |
%add the array itelf |
| 111 |
MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut); |
| 112 |
|
| 113 |
%add its attributes |
| 114 |
nc=ncopen(filOut,'write'); |
| 115 |
ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name); |
| 116 |
ncaddAtt(nc,fldIn.name,'units',fldIn.units); |
| 117 |
ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value); |
| 118 |
ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue); |
| 119 |
ncclose(nc); |
| 120 |
end; |
| 121 |
|
| 122 |
%4) deal with argument output |
| 123 |
if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end; |
| 124 |
|
| 125 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph