| 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, and FillValue (for nc output). |
| 12 |
% |
| 13 |
% fldIn.tim can be set to |
| 14 |
% 'const' (for time invariant climatology), |
| 15 |
% 'monclim' (for monthly climatology) |
| 16 |
% 'monser' (for monthly time series) |
| 17 |
% 'monloop' (for cyclic monthly time series) |
| 18 |
% |
| 19 |
% method can be set to |
| 20 |
% 'polygons' (linear in space) |
| 21 |
% 'bindata' (nearest neighbor in space) |
| 22 |
% |
| 23 |
%note to self: add case when fldIn is a gcmfaces field; nctiles input |
| 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 |
% % |
| 37 |
% profOut=MITprof_resample(profIn,fldIn); |
| 38 |
|
| 39 |
|
| 40 |
gcmfaces_global; |
| 41 |
|
| 42 |
doOut=~isempty(who('filOut')); doOutInit=false; |
| 43 |
if doOut; doOut=~isempty(filOut); end; |
| 44 |
if doOut; doOutInit=isempty(dir(filOut)); end; |
| 45 |
|
| 46 |
if isempty(who('method')); method='polygons'; end; |
| 47 |
|
| 48 |
%1) deal with time line |
| 49 |
if strcmp(fldIn.tim,'monclim'); |
| 50 |
tim_fld=[-0.5:12.5]; rec_fld=[12 1:12 1]; |
| 51 |
tmp1=datevec(profIn.prof_date); |
| 52 |
tmp2=datenum([tmp1(:,1) ones(profIn.np,2) zeros(profIn.np,3)]); |
| 53 |
tim_prof=(profIn.prof_date-tmp2); |
| 54 |
tim_prof(tim_prof>365)=365; |
| 55 |
tim_prof=tim_prof/365*12;%neglecting differences in months length |
| 56 |
elseif strcmp(fldIn.tim,'monloop'); |
| 57 |
tmp1=dir(fldIn.fil); |
| 58 |
nt=tmp1.bytes/prod(mygrid.ioSize)/length(mygrid.RC)/4; |
| 59 |
tmp1=[1:nt]'; tmp2=ones(nt,1)*[1992 1 15 0 0 0]; tmp2(:,2)=tmp1; |
| 60 |
tim_fld=datenum(tmp2); |
| 61 |
tim_fld=[tim_fld(1)-31 tim_fld' tim_fld(end)+31]; |
| 62 |
rec_fld=[nt 1:nt 1]; |
| 63 |
tmp1=datenum([1992 1 1 0 0 0]); |
| 64 |
tmp2=datenum([1992+nt/12 1 1 0 0 0]);; |
| 65 |
tim_prof=tmp1+mod(profIn.prof_date-tmp1,tmp2-tmp1); |
| 66 |
%round up tim_prof to prevent interpolation in time: |
| 67 |
% tmp3=tim_prof*ones(1,length(tim_fld))-ones(length(tim_prof),1)*tim_fld; |
| 68 |
% tmp4=sum(tmp3>0,2); |
| 69 |
% tim_prof=tim_fld(tmp4)'; |
| 70 |
elseif strcmp(fldIn.tim,'const')|strcmp(fldIn.tim,'std'); |
| 71 |
tim_fld=[1 2]; rec_fld=[1 1]; |
| 72 |
tim_prof=1.5*ones(profIn.np,1); |
| 73 |
else; |
| 74 |
error('this case remains to be implemented'); |
| 75 |
end; |
| 76 |
|
| 77 |
lon=profIn.prof_lon; lat=profIn.prof_lat; |
| 78 |
depIn=-mygrid.RC; depOut=profIn.prof_depth; |
| 79 |
profOut=NaN*ones(profIn.np,profIn.nr); |
| 80 |
|
| 81 |
%2) loop over record pairs |
| 82 |
if ~strcmp(method,'bindata'); gcmfaces_bindata; end; |
| 83 |
for tt=1:length(rec_fld)-1; |
| 84 |
tt |
| 85 |
% |
| 86 |
fld0=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt)); |
| 87 |
fld1=mygrid.mskC.*read_bin(fldIn.fil,rec_fld(tt+1)); |
| 88 |
ndim=length(size(fld0{1})); |
| 89 |
fld=cat(ndim+1,fld0,fld1); |
| 90 |
% |
| 91 |
ii=find(tim_prof>=tim_fld(tt)&tim_prof<tim_fld(tt+1)); |
| 92 |
if ~isempty(ii); |
| 93 |
if ~strcmp(method,'bindata'); |
| 94 |
arr=gcmfaces_interp_2d(fld,lon(ii),lat(ii),method); |
| 95 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 96 |
%now linear in time: |
| 97 |
k0=floor(tim_prof(ii)); k1=k0+1; |
| 98 |
a0=tim_prof(ii)-k0; a0=a0*ones(1,profIn.nr); |
| 99 |
profOut(ii,:)=(1-a0).*arr2(:,:,1)+a0.*arr2(:,:,2); |
| 100 |
else; |
| 101 |
[prof_i,prof_j]=gcmfaces_bindata(lon(ii),lat(ii)); |
| 102 |
FLD=convert2array(fld(:,:,:,1)); |
| 103 |
nk=length(mygrid.RC); kk=ones(1,nk); |
| 104 |
np=length(ii); pp=ones(np,1); |
| 105 |
ind2prof=sub2ind(size(FLD),prof_i*kk,prof_j*kk,pp*[1:nk]); |
| 106 |
arr=FLD(ind2prof); |
| 107 |
arr2=gcmfaces_interp_1d(2,depIn,arr,depOut); |
| 108 |
profOut(ii,:)=arr2; |
| 109 |
end; |
| 110 |
% |
| 111 |
if strcmp(fldIn.tim,'std'); |
| 112 |
profOut(ii,:)=profOut(ii,:).*randn(size(profOut(ii,:))); |
| 113 |
end; |
| 114 |
end; |
| 115 |
end; |
| 116 |
|
| 117 |
%3) deal with file output |
| 118 |
if doOutInit; |
| 119 |
%create a header only file to later append resampled fields |
| 120 |
prof=profIn; |
| 121 |
tmp1=fieldnames(prof); |
| 122 |
nt=length(prof.prof_date); |
| 123 |
nr=length(prof.prof_depth); |
| 124 |
for ii=1:length(tmp1); |
| 125 |
tmp2=prod(size(getfield(prof,tmp1{ii}))); |
| 126 |
if tmp2==nt*nr; prof=rmfield(prof,tmp1{ii}); end; |
| 127 |
end; |
| 128 |
MITprof_write(filOut,prof); |
| 129 |
end; |
| 130 |
|
| 131 |
if doOut; |
| 132 |
%add the array itelf |
| 133 |
MITprof_addVar(filOut,fldIn.name,'double',{'iDEPTH','iPROF'},profOut); |
| 134 |
|
| 135 |
%add its attributes |
| 136 |
nc=ncopen(filOut,'write'); |
| 137 |
ncaddAtt(nc,fldIn.name,'long_name',fldIn.long_name); |
| 138 |
ncaddAtt(nc,fldIn.name,'units',fldIn.units); |
| 139 |
ncaddAtt(nc,fldIn.name,'missing_value',fldIn.missing_value); |
| 140 |
ncaddAtt(nc,fldIn.name,'_FillValue',fldIn.FillValue); |
| 141 |
ncclose(nc); |
| 142 |
end; |
| 143 |
|
| 144 |
%4) deal with argument output |
| 145 |
if nargout>0; profOut=setfield(profIn,fldIn.name,profOut); end; |
| 146 |
|
| 147 |
|
Parent Directory
| 
Revision Log
| 
Revision Graph