matlab_class/sample_processing/example_griddata.m

function []=example_griddata(choiceV3orV4);
%object:    a griddata example within gcmfaces
%inputs:    choiceV3orV4 ('v3' or 'v4') selects the sample GRID

gcmfaces_global;
if myenv.verbose>0;
    gcmfaces_msg('===============================================');
    gcmfaces_msg(['*** entering example_griddata that will define ' ...
        'a 2x2 gridded variable (an SSH field here), and then ' ...
        'use griddata to map it to the chosen grid. '],'');
end;

%%%%%%%%%%%%%%%%%
%load parameters:
%%%%%%%%%%%%%%%%%
dir0=[myenv.gcmfaces_dir '/sample_input/'];
dirGrid=[dir0 '/GRID' choiceV3orV4  '/'];
dirIn=[dir0 '/SAMPLE' choiceV3orV4  '/'];
if strcmp(choiceV3orV4,'v4'); nF=5; fileFormat='compact'; else; nF=1; fileFormat='straight'; end;
grid_load(dirGrid,nF,fileFormat);

%%%%%%%%%%%%%%%%%%%%%%%
%get sample data: V3 SSH
if myenv.verbose>0;
    gcmfaces_msg('* subsample test-case SSH field it to a 2x2 grid');
end;
dirV3=[myenv.gcmfaces_dir '/sample_input/SAMPLEv3/'];
etan=rdmds([dirV3 'DDetan'],0); etan(etan==0)=NaN;
dirV3=[myenv.gcmfaces_dir '/sample_input/GRIDv3/'];
lon=rdmds([dirV3 'XC']); lat=rdmds([dirV3 'YC']);
lon=lon(1:2:end,21:2:140);
lat=lat(1:2:end,21:2:140);
etan=etan(1:2:end,21:2:140);


%%%%%%%%%%%%%%%%%%%%%%%
%make sure to accomodate the two longitude conventions:
if myenv.verbose>0;
    gcmfaces_msg('* accomodate both lon. conventions, then convert to data vector');
end;
x=[lon-360;lon;lon+360]; y=[lat;lat;lat]; z=[etan;etan;etan];
jj=find(~isnan(z)); x=x(jj); y=y(jj); z=z(jj);

%%%%%%%%%%%%%%%%%%%%%%%
%do the interpolation:
if myenv.verbose>0;
    gcmfaces_msg('* map/interpolate data using griddata');
end;
z_interp=gcmfaces(5);
for ii=1:mygrid.nFaces;
xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
zi = griddata(x',y',z',xi,yi);
z_interp{ii}=zi;
end;

%msk=mygrid.hFacC(:,:,1); z_interp(find(msk==0))=NaN;

%%%%%%%%%%%%%%%%%%%%%%%
%illustrate the result:

figure; x=[lon-360;lon]; y=[lat;lat]; z=[etan;etan];
pcolor(x,y,z); axis([-180 180 -90 90]); shading flat; caxis([-2 1]); colorbar;

figure; [X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,z_interp);
pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat; caxis([-2 1]); colorbar;

if myenv.verbose>0;
    gcmfaces_msg('* note of caution : griddata puts values on land. ');
    gcmfaces_msg('*** leaving example_griddata');
    gcmfaces_msg('===============================================');
end;
1	gforget	1.8	function []=example_griddata(choiceV3orV4);
2			%object: a griddata example within gcmfaces
3			%inputs: choiceV3orV4 ('v3' or 'v4') selects the sample GRID
4	gforget	1.1
5	gforget	1.11	gcmfaces_global;
6			if myenv.verbose>0;
7			gcmfaces_msg('===============================================');
8			gcmfaces_msg(['*** entering example_griddata that will define ' ...
9			'a 2x2 gridded variable (an SSH field here), and then ' ...
10			'use griddata to map it to the chosen grid. '],'');
11			end;
12
13	gforget	1.1	%%%%%%%%%%%%%%%%%
14			%load parameters:
15			%%%%%%%%%%%%%%%%%
16	gforget	1.9	dir0=[myenv.gcmfaces_dir '/sample_input/'];
17	gforget	1.8	dirGrid=[dir0 '/GRID' choiceV3orV4 '/'];
18			dirIn=[dir0 '/SAMPLE' choiceV3orV4 '/'];
19			if strcmp(choiceV3orV4,'v4'); nF=5; fileFormat='compact'; else; nF=1; fileFormat='straight'; end;
20	gforget	1.10	grid_load(dirGrid,nF,fileFormat);
21	gforget	1.1
22			%%%%%%%%%%%%%%%%%%%%%%%
23	gforget	1.2	%get sample data: V3 SSH
24	gforget	1.11	if myenv.verbose>0;
25			gcmfaces_msg('* subsample test-case SSH field it to a 2x2 grid');
26			end;
27	gforget	1.9	dirV3=[myenv.gcmfaces_dir '/sample_input/SAMPLEv3/'];
28	gforget	1.2	etan=rdmds([dirV3 'DDetan'],0); etan(etan==0)=NaN;
29	gforget	1.9	dirV3=[myenv.gcmfaces_dir '/sample_input/GRIDv3/'];
30	gforget	1.2	lon=rdmds([dirV3 'XC']); lat=rdmds([dirV3 'YC']);
31	gforget	1.11	lon=lon(1:2:end,21:2:140);
32			lat=lat(1:2:end,21:2:140);
33			etan=etan(1:2:end,21:2:140);
34
35	gforget	1.1
36			%%%%%%%%%%%%%%%%%%%%%%%
37	gforget	1.11	%make sure to accomodate the two longitude conventions:
38			if myenv.verbose>0;
39			gcmfaces_msg('* accomodate both lon. conventions, then convert to data vector');
40			end;
41			x=[lon-360;lon;lon+360]; y=[lat;lat;lat]; z=[etan;etan;etan];
42	gforget	1.1	jj=find(~isnan(z)); x=x(jj); y=y(jj); z=z(jj);
43
44	gforget	1.11	%%%%%%%%%%%%%%%%%%%%%%%
45			%do the interpolation:
46			if myenv.verbose>0;
47			gcmfaces_msg('* map/interpolate data using griddata');
48			end;
49	gforget	1.1	z_interp=gcmfaces(5);
50	gforget	1.8	for ii=1:mygrid.nFaces;
51	gforget	1.1	xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
52			zi = griddata(x',y',z',xi,yi);
53			z_interp{ii}=zi;
54			end;
55
56			%msk=mygrid.hFacC(:,:,1); z_interp(find(msk==0))=NaN;
57
58			%%%%%%%%%%%%%%%%%%%%%%%
59			%illustrate the result:
60
61	gforget	1.11	figure; x=[lon-360;lon]; y=[lat;lat]; z=[etan;etan];
62	gforget	1.1	pcolor(x,y,z); axis([-180 180 -90 90]); shading flat; caxis([-2 1]); colorbar;
63
64	gforget	1.11	figure; [X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,z_interp);
65	gforget	1.1	pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat; caxis([-2 1]); colorbar;
66
67	gforget	1.11	if myenv.verbose>0;
68			gcmfaces_msg('* note of caution : griddata puts values on land. ');
69			gcmfaces_msg('*** leaving example_griddata');
70			gcmfaces_msg('===============================================');
71			end;