matlab_class/ecco_v4/gcmfaces_remap.m

function []=gcmfaces_remap(dirIn,fileIn,dirOut,fileOut,varargin);
%object:    use bin average to remap a lat-lon grid product to a gcmfaces grid
%inputs:    dirIn is the input directory
%           fileIn is the input file name without the final four year characters 
%               (e.g. 'SST_monthly_r2_' to process 'SST_monthly_r2_1992' etc.)
%           dirOut and filOut are the corresponding output names
%optional:  gridIn states the orinitaing grid. It can be set to 
%               1 (default) implying that the grid is [0.5:359.5] & [-89.5:89.5]
%               2 implying that the grid is [0.5:359.5] & [-79.5:79.5]
%               {x,y} where x and y are the position arrays
%
%assumption: mygrid has been read using grid_load

global mygrid mytri;
%create triangulation
gcmfaces_bindata;
%list files to be processed
listFiles=dir([dirIn fileIn '*']);

%input grid (lon must be 0-360 for gcmfaces_remap_one)
if nargin>4; gridIn=varargin{1}; else; gridIn=1; end;
%case of user defined grid
if iscell(gridIn); x=gridIn{1}; y=gridIn{2}; [nx,ny]=size(x); mis=0; gridIn=0; end;
%standard cases
if gridIn==1;
    x=[0.5:359.5]'*ones(1,180);
    y=ones(360,1)*[-89.5:89.5];
    [nx,ny]=size(x);
    mis=0;
else gridIn==2;
    x=[0.5:359.5]'*ones(1,160);
    y=ones(360,1)*[-79.5:79.5];
    [nx,ny]=size(x);
    mis=0;
end;

%process one file after the other
for ii=1:length(listFiles);
    yy=listFiles(ii).name(end-3:end); fprintf(['processing ' fileIn '   for year   ' yy '\n']);
    nt=listFiles(ii).bytes/nx/ny/4; if round(nt)~=nt; error('inconsistent sizes'); end;
    %read data
    fld=reshape(read2memory([dirIn listFiles(ii).name],[nx*ny*nt 1]),[nx ny nt]);
    %mask land
    fld(fld==mis)=NaN;
    %map to v4 grid
    FLD=convert2array(zeros(360,360,nt));
    for tt=1:1; FLD(:,:,tt)=gcmfaces_remap_one(x,y,fld(:,:,tt),3); end;
    %set missing value
    FLD(find(isnan(FLD)))=mis;
    %write data
    write2file([dirOut fileOut yy],convert2gcmfaces(FLD));
end;

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

function [fld]=gcmfaces_remap_one(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
%notes:     nDblRes>0 is useful if the starting resolution is coarse 
%               enough that the bin average would leave many empry points.
%           fld should show NaN for missing values

global mygrid;
if nargin>4; mskExtrap=varargin{1}; else; mskExtrap=[]; 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;

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

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.2	function []=gcmfaces_remap(dirIn,fileIn,dirOut,fileOut,varargin);
2	gforget	1.1	%object: use bin average to remap a lat-lon grid product to a gcmfaces grid
3			%inputs: dirIn is the input directory
4			% fileIn is the input file name without the final four year characters
5			% (e.g. 'SST_monthly_r2_' to process 'SST_monthly_r2_1992' etc.)
6			% dirOut and filOut are the corresponding output names
7	gforget	1.2	%optional: gridIn states the orinitaing grid. It can be set to
8			% 1 (default) implying that the grid is [0.5:359.5] & [-89.5:89.5]
9			% 2 implying that the grid is [0.5:359.5] & [-79.5:79.5]
10			% {x,y} where x and y are the position arrays
11	gforget	1.1	%
12			%assumption: mygrid has been read using grid_load
13
14			global mygrid mytri;
15			%create triangulation
16			gcmfaces_bindata;
17			%list files to be processed
18			listFiles=dir([dirIn fileIn '*']);
19	gforget	1.2
20			%input grid (lon must be 0-360 for gcmfaces_remap_one)
21			if nargin>4; gridIn=varargin{1}; else; gridIn=1; end;
22			%case of user defined grid
23			if iscell(gridIn); x=gridIn{1}; y=gridIn{2}; [nx,ny]=size(x); mis=0; gridIn=0; end;
24			%standard cases
25			if gridIn==1;
26	gforget	1.1	x=[0.5:359.5]'*ones(1,180);
27			y=ones(360,1)*[-89.5:89.5];
28			[nx,ny]=size(x);
29			mis=0;
30	gforget	1.2	else gridIn==2;
31	gforget	1.1	x=[0.5:359.5]'*ones(1,160);
32			y=ones(360,1)*[-79.5:79.5];
33			[nx,ny]=size(x);
34			mis=0;
35			end;
36
37	gforget	1.2	%process one file after the other
38	gforget	1.1	for ii=1:length(listFiles);
39			yy=listFiles(ii).name(end-3:end); fprintf(['processing ' fileIn ' for year ' yy '\n']);
40			nt=listFiles(ii).bytes/nx/ny/4; if round(nt)~=nt; error('inconsistent sizes'); end;
41			%read data
42			fld=reshape(read2memory([dirIn listFiles(ii).name],[nxnynt 1]),[nx ny nt]);
43			%mask land
44			fld(fld==mis)=NaN;
45			%map to v4 grid
46			FLD=convert2array(zeros(360,360,nt));
47			for tt=1:1; FLD(:,:,tt)=gcmfaces_remap_one(x,y,fld(:,:,tt),3); end;
48			%set missing value
49			FLD(find(isnan(FLD)))=mis;
50			%write data
51			write2file([dirOut fileOut yy],convert2gcmfaces(FLD));
52			end;
53
54			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
55
56			function [fld]=gcmfaces_remap_one(lon,lat,fld,nDblRes,varargin);
57			%object: use bin average to remap ONE lat-lon grid FIELD to a gcmfaces grid
58			%input: lon,lat,fld are the gridded product arrays
59			% nDblRes is the number of times the input field
60			% resolution should be doubled before bin average.
61			%optional: mskExtrap is the mask to which to extrapolate after bin average
62			% If it is not specified we do no extrapolation
63			%notes: nDblRes>0 is useful if the starting resolution is coarse
64			% enough that the bin average would leave many empry points.
65			% fld should show NaN for missing values
66
67			global mygrid;
68			if nargin>4; mskExtrap=varargin{1}; else; mskExtrap=[]; end;
69
70			%refine grid before bin average:
71			for ii=1:nDblRes;
72			[lon,lat,fld]=dbl_res(lon,lat,fld,ii);
73			end;
74			%put in vector form:
75			tmp1=find(~isnan(fld));
76			lon=lon(tmp1); lat=lat(tmp1); fld=fld(tmp1);
77			%switch longitude range to -180+180
78			lon(find(lon>180))=lon(find(lon>180))-360;
79			%compute bin average:
80			fld=gcmfaces_bindata(lon,lat,fld);
81			%potential extrapolation:
82			if ~isempty(mskExtrap); fld=diffsmooth2D_extrap_inv(fld,mskExtrap); end;
83			%apply surface land mask:
84			fld=fld.*mygrid.mskC(:,:,1);
85
86			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
87
88			function [lonOut,latOut,obsOut]=dbl_res(lon,lat,obs,extendNaN);
89			%object: use neighbor average to double the fields resolution
90			%inputs: lon,lat,fld are the gridded product arrays, with
91			% NaN for missing values in fld.
92			% extendNaN states whether between a real and a NaN
93			% we add a NaN (extendNaN==1) or a real (extendNaN~=1)
94			%outputs: lonOut,latOut,fldOut are the x2 resolution arrays
95			%
96			%assumptions: 0-360 longitudes and a consistent ordering of arrays
97
98			%check longitude range:
99			if ~isempty(find(lon(:)<0))\|~isempty(find(lon(:)>360)); fprintf('problem in longitude specs\n'); return; end;
100			%check simple ordering:
101			tmp1=diff(lon(:,1)); if ~isempty(find(tmp1<=0)); fprintf('problem in longitude specs\n'); return; end;
102
103			%make sure that we are >=0 and <360:
104			if lon(1,1)>0&lon(end,1)==360; lon=[lon(end,:)-360;lon(1:end-1,:)];
105			lat=[lat(end,:);lat(1:end-1,:)]; obs=[obs(end,:);obs(1:end-1,:)]; end;
106			%make sure that the last point is not duplicated:
107			if lon(1,1)==0&lon(end,1)==360; lon=lon(1:end-1,:);
108			lat=lat(1:end-1,:); obs=obs(1:end-1,:); end;
109
110			%what do we do with last longitude points:
111			if lon(1,1)>0&lon(end,1)<360;
112			addOnePt=1; %add point at the beginning
113			else;
114			addOnePt=2; %add point at the end
115			end;
116
117			for ii=1:3;
118			%0) get field
119			if ii==1; tmpA=lon;
120			elseif ii==2; tmpA=lat;
121			elseif ii==3; tmpA=obs;
122			end;
123
124			%1) zonal direction:
125			tmpB=nanmean2flds(tmpA(1:end-1,:),tmpA(2:end,:),extendNaN);
126			tmpC=ones(size(tmpA,1)*2-1,size(tmpA,2));
127			tmpC(1:2:end,:)=tmpA;
128			tmpC(2:2:end-1,:)=tmpB;
129			%treat zonal edge of the domain
130			if addOnePt==1; %add one point at the beginning
131			if ii==1; offset=-360; else; offset=0; end;
132			tmpB=nanmean2flds(tmpA(1,:),offset+tmpA(end,:),extendNaN); tmpC=[tmpB;tmpC];
133			end;
134			if addOnePt==2; %add one point at the end
135			if ii==1; offset=360; else; offset=0; end;
136			tmpB=nanmean2flds(offset+tmpA(1,:),tmpA(end,:),extendNaN); tmpC=[tmpC;tmpB];
137			end;
138			%check that we did not go too far (should not happen)
139			if ii==1; tmp1=~isempty(find(lon(:)<0))\|~isempty(find(lon(:)>360));
140			if tmp1; fprintf('problem in longitude interp\n'); return; end;
141			end;
142
143			%2) meridional direction:
144			tmpB=nanmean2flds(tmpC(:,1:end-1),tmpC(:,2:end),extendNaN);
145			tmpA=ones(size(tmpC,1),size(tmpC,2)*2-1);
146			tmpA(:,1:2:end)=tmpC;
147			tmpA(:,2:2:end-1)=tmpB;
148
149			%3) store field:
150			if ii==1; lonOut=tmpA;
151			elseif ii==2; latOut=tmpA;
152			elseif ii==3; obsOut=tmpA;
153			end;
154			end;
155
156			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
157
158			function [fld]=nanmean2flds(fld1,fld2,extendNaN);
159			%object: compute the average of two fields, accounting for NaNs
160			%inputs: fld1 and fld2 are the two fields
161			% if extendNaN==1 the result is NaN if either fld1 or fld2 is NaN.
162			% Otherwise the result is real if either fld1 or fld2 is real.
163
164			if extendNaN==1;
165			fld=(fld1+fld2)/2;
166			else;
167			msk1=~isnan(fld1); fld1(isnan(fld1))=0;
168			msk2=~isnan(fld2); fld2(isnan(fld2))=0;
169			fld=fld1+fld2;
170			msk=msk1+msk2;
171			msk(msk==0)=NaN;
172			fld=fld./msk;
173			end;
174