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_2d)

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_2d(x,y,fld(:,:,tt),3); end;
    %set missing value
    FLD(find(isnan(FLD)))=mis;
    %write data
    write2file([dirOut fileOut yy],convert2gcmfaces(FLD));
end;

1	function []=gcmfaces_remap(dirIn,fileIn,gridIn,dirOut,fileOut);
2	%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	% 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	% 2 implying that the grid is [0.5:359.5] & [-79.5:79.5]
9	% 3 implying that the grid is [.125:.25:360-.125] & [-90+.125:.25:90-.125]
10	% {x,y} where x and y are the position arrays
11	% dirOut and filOut are the corresponding output names
12	%
13	%assumption: mygrid has been read using grid_load
14	% for the input grid, lon must be 0-360 (see gcmfaces_remap_2d)
15
16	global mygrid mytri;
17	%create triangulation
18	gcmfaces_bindata;
19	%list files to be processed
20	listFiles=dir([dirIn fileIn '*']);
21
22	%case of user defined grid
23	if iscell(gridIn); x=gridIn{1}; y=gridIn{2}; [nx,ny]=size(x); mis=0;
24	else;
25	%standard cases
26	if gridIn==1;%gloabl 1 degree grid
27	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	elseif gridIn==2;%ECCO 1 degree grid
32	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	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	end;
45
46	%process one file after the other
47	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	for tt=1:nt; FLD(:,:,tt)=gcmfaces_remap_2d(x,y,fld(:,:,tt),3); end;
57	%set missing value
58	FLD(find(isnan(FLD)))=mis;
59	%write data
60	write2file([dirOut fileOut yy],convert2gcmfaces(FLD));
61	end;
62