| 1 |
%this function does a running mean along dimension dim_cur, |
function [fldOut]=runmean(fldIn,halfWindow,dim,varargin); |
| 2 |
%averaging over indices [-nb_cur:nb_cur] |
%object: compute running mean window ('rmw') over a dimension |
| 3 |
|
%input: fldIn is the field to which the rmw will be applied |
| 4 |
|
% halfWindow is the half width of the rmw |
| 5 |
|
% dim is the dimension over which the rmw will be applied |
| 6 |
|
%optional: doCycle states whether the boundary condition is cyclic (1) or |
| 7 |
|
% not (0; default). If doCycle==0, the no. of averaged points |
| 8 |
|
% decreases from 1+2*halfWindow to halfWindow at the edges. |
| 9 |
|
%output: fldOut is the resulting field |
| 10 |
% |
% |
| 11 |
function [field_out]=runmean(field_in,nb_cur,dim_cur); |
%notes: - NaNs are discarded in the rmw, implying that an average is |
| 12 |
|
% computed if the rmw contains at least 1 valid point. |
| 13 |
|
% - setting halfWindow to 0 implies fldOut=fldIn. |
| 14 |
|
|
| 15 |
|
%determine and check a couple things |
| 16 |
|
fld_isa_gcmfaces=isa(fldIn,'gcmfaces'); |
| 17 |
|
if fld_isa_gcmfaces; |
| 18 |
|
if dim<3; |
| 19 |
|
error('for gcmfaces objects runmean excludes dim=1 or 2'); |
| 20 |
|
end; |
| 21 |
|
end; |
| 22 |
|
|
| 23 |
|
if nargin==3; doCycle=0; else; doCycle=varargin{1}; end; |
| 24 |
|
|
| 25 |
|
%switch to array format if needed |
| 26 |
|
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end; |
| 27 |
|
|
| 28 |
|
%switch dim to first dimension |
| 29 |
|
sizeIn=size(fldIn); |
| 30 |
|
|
| 31 |
|
perm1to2=[1:length(sizeIn)]; |
| 32 |
|
perm1to2=[dim perm1to2(find(perm1to2~=dim))]; |
| 33 |
|
perm2to1=[[1:dim-1]+1 1 [dim+1:length(sizeIn)]]; |
| 34 |
|
sizeIn2=sizeIn(perm1to2); |
| 35 |
|
|
| 36 |
|
fldIn=permute(fldIn,perm1to2); |
| 37 |
|
|
| 38 |
|
sizeCur=size(fldIn); |
| 39 |
|
if ~doCycle; |
| 40 |
|
%add NaNs at both edges |
| 41 |
|
sizeCur(1)=sizeCur(1)+2*halfWindow; |
| 42 |
|
fldIn2=NaN*ones(sizeCur); |
| 43 |
|
fldIn2(halfWindow+1:end-halfWindow,:,:,:)=fldIn; |
| 44 |
|
fldIn=fldIn2; clear fldIn2; |
| 45 |
|
end; |
| 46 |
|
|
| 47 |
|
%create mask and remove NaNs: |
| 48 |
|
fldMsk=~isnan(fldIn); |
| 49 |
|
fldIn(isnan(fldIn))=0; |
| 50 |
|
fldCnt=0*fldIn; |
| 51 |
|
|
| 52 |
|
%apply the running mean |
| 53 |
|
fldOut=zeros(sizeCur); |
| 54 |
|
for tcur=-halfWindow:halfWindow |
| 55 |
|
%To have halfWindow*2 coeffs rather than halfWindow*2+1, centered to the current |
| 56 |
|
%point, it is convenient to reduce the weight of the farthest points to 1/2. |
| 57 |
|
%This became necessary to get proper annual means, from monthly data, with halfWindow=6. |
| 58 |
|
if abs(tcur)==halfWindow; fac=1/2; else; fac=1; end; |
| 59 |
|
tmp1=circshift(fldIn,[tcur zeros(1,length(sizeCur)-1)]); |
| 60 |
|
fldOut=fldOut+fac*tmp1; |
| 61 |
|
tmp1=circshift(fldMsk,[tcur zeros(1,length(sizeCur)-1)]); |
| 62 |
|
fldCnt=fldCnt+fac*tmp1; |
| 63 |
|
end |
| 64 |
|
|
| 65 |
if nb_cur~=0 |
fldCnt(fldCnt==0)=NaN; |
| 66 |
|
fldOut=fldOut./fldCnt; |
| 67 |
|
|
| 68 |
size_cur=size(field_in); |
if ~doCycle; |
| 69 |
|
fldOut=fldOut(halfWindow+1:end-halfWindow,:,:,:); |
| 70 |
|
%consistent with old version bug (one point offset) |
| 71 |
|
% fldOut=fldOut(halfWindow:end-halfWindow,:,:,:); |
| 72 |
|
end; |
| 73 |
|
|
| 74 |
perm1to2=[1:length(size_cur)]; |
%switch dimensions order back to original order |
| 75 |
perm1to2=[dim_cur perm1to2(find(perm1to2~=dim_cur))]; |
fldOut=permute(fldOut,perm2to1); |
|
perm2to1=[[1:dim_cur-1]+1 1 [dim_cur+1:length(size_cur)]]; |
|
|
size_cur2=size_cur(perm1to2); |
|
|
|
|
|
field_in2=permute(field_in,perm1to2); |
|
|
|
|
|
field_out2=zeros(size_cur2); |
|
|
count_out2=zeros(size_cur2(1),2); |
|
|
for tcur=-nb_cur:nb_cur |
|
|
tmp1=[tcur:tcur-1+size_cur2(1)]; |
|
|
tmp2=find(tmp1>=1&tmp1<=size_cur2(1)); |
|
|
field_out2(tmp2,:)=field_out2(tmp2,:)+field_in2(tmp1(tmp2),:); |
|
|
count_out2(tmp2,:)=count_out2(tmp2,:)+1; |
|
|
end |
|
|
for tcur=1:size_cur2(1); |
|
|
field_out2(tcur,:)=field_out2(tcur,:)/count_out2(tcur); |
|
|
end |
|
| 76 |
|
|
| 77 |
%field_out=field_in-permute(field_out2,perm2to1); |
%switch back to gcmfaces format if needed |
| 78 |
field_out=permute(field_out2,perm2to1); |
if fld_isa_gcmfaces; fldOut=convert2array(fldOut); end; |
| 79 |
|
|
|
else |
|
|
field_out=field_in; |
|
|
end%if nb_cur~=0 |
|
| 80 |
|
|
| 81 |
|
|
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch