matlab_class/gcmfaces_calc/gcmfaces_remap_2d.m

function [fld]=gcmfaces_remap_2d(lon,lat,fld,nDblRes,varargin);
%object:    use bin average to remap ONE lat-lon grid FIELD to a gcmfaces grid
%input:     lon,lat,fld are the gridded product arrays
%           nDblRes is the number of times the input field 
%               resolution should be doubled before bin average.
%optional:  mskExtrap is the mask to which to extrapolate after bin average
%               If it is not specified we do no extrapolation
%note:      nDblRes>0 is useful if the starting resolution is coarse 
%               enough that the bin average would leave many empty points.
%assumption:fld should show NaN for missing values

global mygrid;
if nargin>4; mskExtrap=varargin{1}; else; mskExtrap=[]; end;

%apply meshgrid if not already the case

%switch longitude range to 0-360 if not already the case
if ~isempty(find(lon(:,1)<0));
    ii=max(find(lon(:,1)<0));
    lon=circshift(lon,[-ii 0]);
    lon(end-ii+1:end,:)=lon(end-ii+1:end,:)+360;
    fld=circshift(fld,[-ii 0]);    
end;

%refine grid before bin average:
for ii=1:nDblRes; 
    [lon,lat,fld]=dbl_res(lon,lat,fld,ii); 
end;
%put in vector form:
tmp1=find(~isnan(fld));
lon=lon(tmp1); lat=lat(tmp1); fld=fld(tmp1);
%switch longitude range to -180+180
lon(find(lon>180))=lon(find(lon>180))-360;
%compute bin average:
fld=gcmfaces_bindata(lon,lat,fld);
%potential extrapolation:
if ~isempty(mskExtrap); fld=diffsmooth2D_extrap_inv(fld,mskExtrap); end;
%apply surface land mask:
fld=fld.*mygrid.mskC(:,:,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [lonOut,latOut,obsOut]=dbl_res(lon,lat,obs,extendNaN);
%object:    use neighbor average to double the fields resolution
%inputs:    lon,lat,fld are the gridded product arrays, with 
%               NaN for missing values in fld.
%           extendNaN states whether between a real and a NaN
%               we add a NaN (extendNaN==1) or a real (extendNaN~=1)
%outputs:   lonOut,latOut,fldOut are the x2 resolution arrays
%
%assumptions: 0-360 longitudes and a consistent ordering of arrays

%check longitude range:
if ~isempty(find(lon(:)<0))|~isempty(find(lon(:)>360)); fprintf('problem in longitude specs\n'); return; end;
%check simple ordering:
tmp1=diff(lon(:,1)); if ~isempty(find(tmp1<=0)); fprintf('problem in longitude specs\n'); return; end;

%make sure that we are >=0 and <360:
if lon(1,1)>0&lon(end,1)==360; lon=[lon(end,:)-360;lon(1:end-1,:)];
    lat=[lat(end,:);lat(1:end-1,:)]; obs=[obs(end,:);obs(1:end-1,:)]; end;
%make sure that the last point is not duplicated:
if lon(1,1)==0&lon(end,1)==360; lon=lon(1:end-1,:);
    lat=lat(1:end-1,:); obs=obs(1:end-1,:); end;

%what do we do with last longitude points:
if lon(1,1)>0&lon(end,1)<360;
    addOnePt=1; %add point at the beginning
else;
    addOnePt=2; %add point at the end
end;

for ii=1:3;
    %0) get field
    if ii==1;     tmpA=lon;
    elseif ii==2; tmpA=lat;
    elseif ii==3; tmpA=obs;
    end;
    
    %1) zonal direction:
    tmpB=nanmean2flds(tmpA(1:end-1,:),tmpA(2:end,:),extendNaN);
    tmpC=ones(size(tmpA,1)*2-1,size(tmpA,2));
    tmpC(1:2:end,:)=tmpA;
    tmpC(2:2:end-1,:)=tmpB;
    %treat zonal edge of the domain
    if addOnePt==1; %add one point at the beginning
        if ii==1; offset=-360; else; offset=0; end;
        tmpB=nanmean2flds(tmpA(1,:),offset+tmpA(end,:),extendNaN); tmpC=[tmpB;tmpC];
    end;
    if addOnePt==2; %add one point at the end
        if ii==1; offset=360; else; offset=0; end;
        tmpB=nanmean2flds(offset+tmpA(1,:),tmpA(end,:),extendNaN); tmpC=[tmpC;tmpB];
    end;
    %check that we did not go too far (should not happen)
    if ii==1; tmp1=~isempty(find(lon(:)<0))|~isempty(find(lon(:)>360));
        if tmp1; fprintf('problem in longitude interp\n'); return; end;
    end;
    
    %2) meridional direction:
    tmpB=nanmean2flds(tmpC(:,1:end-1),tmpC(:,2:end),extendNaN);
    tmpA=ones(size(tmpC,1),size(tmpC,2)*2-1);
    tmpA(:,1:2:end)=tmpC;
    tmpA(:,2:2:end-1)=tmpB;
    
    %3) store field:
    if ii==1; lonOut=tmpA;
    elseif ii==2; latOut=tmpA;
    elseif ii==3; obsOut=tmpA;
    end;    
end;

%fix potential truncation errors:
lonOut(lonOut<0)=0; lonOut(lonOut>360)=360; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fld]=nanmean2flds(fld1,fld2,extendNaN);
%object:    compute the average of two fields, accounting for NaNs
%inputs:    fld1 and fld2 are the two fields
%           if extendNaN==1 the result is NaN if either fld1 or fld2 is NaN. 
%               Otherwise the result is real if either fld1 or fld2 is real.

if extendNaN==1;
    fld=(fld1+fld2)/2;
else;
    msk1=~isnan(fld1); fld1(isnan(fld1))=0;
    msk2=~isnan(fld2); fld2(isnan(fld2))=0;
    fld=fld1+fld2;
    msk=msk1+msk2;
    msk(msk==0)=NaN;
    fld=fld./msk;
end;

1	gforget	1.1	function [fld]=gcmfaces_remap_2d(lon,lat,fld,nDblRes,varargin);
2			%object: use bin average to remap ONE lat-lon grid FIELD to a gcmfaces grid
3			%input: lon,lat,fld are the gridded product arrays
4			% nDblRes is the number of times the input field
5			% resolution should be doubled before bin average.
6			%optional: mskExtrap is the mask to which to extrapolate after bin average
7			% If it is not specified we do no extrapolation
8			%note: nDblRes>0 is useful if the starting resolution is coarse
9			% enough that the bin average would leave many empty points.
10			%assumption:fld should show NaN for missing values
11
12			global mygrid;
13			if nargin>4; mskExtrap=varargin{1}; else; mskExtrap=[]; end;
14	gforget	1.2
15			%apply meshgrid if not already the case
16
17			%switch longitude range to 0-360 if not already the case
18			if ~isempty(find(lon(:,1)<0));
19			ii=max(find(lon(:,1)<0));
20			lon=circshift(lon,[-ii 0]);
21			lon(end-ii+1:end,:)=lon(end-ii+1:end,:)+360;
22			fld=circshift(fld,[-ii 0]);
23			end;
24
25	gforget	1.1	%refine grid before bin average:
26			for ii=1:nDblRes;
27			[lon,lat,fld]=dbl_res(lon,lat,fld,ii);
28			end;
29			%put in vector form:
30			tmp1=find(~isnan(fld));
31			lon=lon(tmp1); lat=lat(tmp1); fld=fld(tmp1);
32			%switch longitude range to -180+180
33			lon(find(lon>180))=lon(find(lon>180))-360;
34			%compute bin average:
35			fld=gcmfaces_bindata(lon,lat,fld);
36			%potential extrapolation:
37			if ~isempty(mskExtrap); fld=diffsmooth2D_extrap_inv(fld,mskExtrap); end;
38			%apply surface land mask:
39			fld=fld.*mygrid.mskC(:,:,1);
40
41			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
42
43			function [lonOut,latOut,obsOut]=dbl_res(lon,lat,obs,extendNaN);
44			%object: use neighbor average to double the fields resolution
45			%inputs: lon,lat,fld are the gridded product arrays, with
46			% NaN for missing values in fld.
47			% extendNaN states whether between a real and a NaN
48			% we add a NaN (extendNaN==1) or a real (extendNaN~=1)
49			%outputs: lonOut,latOut,fldOut are the x2 resolution arrays
50			%
51			%assumptions: 0-360 longitudes and a consistent ordering of arrays
52
53			%check longitude range:
54			if ~isempty(find(lon(:)<0))\|~isempty(find(lon(:)>360)); fprintf('problem in longitude specs\n'); return; end;
55			%check simple ordering:
56			tmp1=diff(lon(:,1)); if ~isempty(find(tmp1<=0)); fprintf('problem in longitude specs\n'); return; end;
57
58			%make sure that we are >=0 and <360:
59			if lon(1,1)>0&lon(end,1)==360; lon=[lon(end,:)-360;lon(1:end-1,:)];
60			lat=[lat(end,:);lat(1:end-1,:)]; obs=[obs(end,:);obs(1:end-1,:)]; end;
61			%make sure that the last point is not duplicated:
62			if lon(1,1)==0&lon(end,1)==360; lon=lon(1:end-1,:);
63			lat=lat(1:end-1,:); obs=obs(1:end-1,:); end;
64
65			%what do we do with last longitude points:
66			if lon(1,1)>0&lon(end,1)<360;
67			addOnePt=1; %add point at the beginning
68			else;
69			addOnePt=2; %add point at the end
70			end;
71
72			for ii=1:3;
73			%0) get field
74			if ii==1; tmpA=lon;
75			elseif ii==2; tmpA=lat;
76			elseif ii==3; tmpA=obs;
77			end;
78
79			%1) zonal direction:
80			tmpB=nanmean2flds(tmpA(1:end-1,:),tmpA(2:end,:),extendNaN);
81			tmpC=ones(size(tmpA,1)*2-1,size(tmpA,2));
82			tmpC(1:2:end,:)=tmpA;
83			tmpC(2:2:end-1,:)=tmpB;
84			%treat zonal edge of the domain
85			if addOnePt==1; %add one point at the beginning
86			if ii==1; offset=-360; else; offset=0; end;
87			tmpB=nanmean2flds(tmpA(1,:),offset+tmpA(end,:),extendNaN); tmpC=[tmpB;tmpC];
88			end;
89			if addOnePt==2; %add one point at the end
90			if ii==1; offset=360; else; offset=0; end;
91			tmpB=nanmean2flds(offset+tmpA(1,:),tmpA(end,:),extendNaN); tmpC=[tmpC;tmpB];
92			end;
93			%check that we did not go too far (should not happen)
94			if ii==1; tmp1=~isempty(find(lon(:)<0))\|~isempty(find(lon(:)>360));
95			if tmp1; fprintf('problem in longitude interp\n'); return; end;
96			end;
97
98			%2) meridional direction:
99			tmpB=nanmean2flds(tmpC(:,1:end-1),tmpC(:,2:end),extendNaN);
100			tmpA=ones(size(tmpC,1),size(tmpC,2)*2-1);
101			tmpA(:,1:2:end)=tmpC;
102			tmpA(:,2:2:end-1)=tmpB;
103
104			%3) store field:
105			if ii==1; lonOut=tmpA;
106			elseif ii==2; latOut=tmpA;
107			elseif ii==3; obsOut=tmpA;
108			end;
109			end;
110
111	gforget	1.3	%fix potential truncation errors:
112			lonOut(lonOut<0)=0; lonOut(lonOut>360)=360;
113
114	gforget	1.1	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115
116			function [fld]=nanmean2flds(fld1,fld2,extendNaN);
117			%object: compute the average of two fields, accounting for NaNs
118			%inputs: fld1 and fld2 are the two fields
119			% if extendNaN==1 the result is NaN if either fld1 or fld2 is NaN.
120			% Otherwise the result is real if either fld1 or fld2 is real.
121
122			if extendNaN==1;
123			fld=(fld1+fld2)/2;
124			else;
125			msk1=~isnan(fld1); fld1(isnan(fld1))=0;
126			msk2=~isnan(fld2); fld2(isnan(fld2))=0;
127			fld=fld1+fld2;
128			msk=msk1+msk2;
129			msk(msk==0)=NaN;
130			fld=fld./msk;
131			end;
132