matlab_class/ecco_v4/gcmfaces_remap.m

function []=gcmfaces_remap(dirIn,fileIn,gridIn,dirOut,fileOut);
%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.)
%           gridIn states the originating grid. It can be set to 
%               1 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]
%               3 implying that the grid is [.125:.25:360-.125] & [-90+.125:.25:90-.125]
%               {x,y} where x and y are the position arrays
%           dirOut and filOut are the corresponding output names
%
%assumption: mygrid has been read using grid_load
%            for the input grid, lon must be 0-360 (see gcmfaces_remap_one)

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

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