matlab_class/sample_processing/example_griddata.m

function []=example_griddata();
%object:    a griddata example within gcmfaces

%%%%%%%%%%%%%%%%%
%load grid:
%%%%%%%%%%%%%%%%%

gcmfaces_global;

if isempty(mygrid);
   grid_load;
end;

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

%%%%%%%%%%%%%%%%%%%%%%%
%get sample data: subsampled mygrid.Depth
if myenv.verbose>0;
    gcmfaces_msg('* define test case : subsampled bathymetry');
end;
original_Depth=mygrid.Depth;
original_Depth(original_Depth==0)=NaN;
lon=mygrid.XC; lat=mygrid.YC;
for ff=1:mygrid.nFaces;
    lon{ff}(1:2:end,:)=NaN;
    lat{ff}(1:2:end,:)=NaN;
    original_Depth{ff}(1:2:end,:)=NaN;
    lon{ff}(:,1:2:end)=NaN;
    lat{ff}(:,1:2:end)=NaN;
    original_Depth{ff}(:,1:2:end)=NaN;
end;
lon=convert2vector(lon);
lat=convert2vector(lat);
vec=convert2vector(original_Depth);

%%%%%%%%%%%%%%%%%%%%%%%
%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=[vec;vec;vec];
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;
griddata_Depth=gcmfaces(5);
for ii=1:mygrid.nFaces;
    xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
    zi = griddata(x',y',z',xi,yi);
    griddata_Depth{ii}=zi;
end;

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

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

figureL;
subplot(2,1,1); set(gca,'FontSize',14);
[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,original_Depth);
pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat;
title('input to griddata');
subplot(2,1,2); set(gca,'FontSize',14);
[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,griddata_Depth);
pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat;
title('output of griddata');

if myenv.verbose>0;
    gcmfaces_msg('* note of caution : griddata puts values on land. ');
    gcmfaces_msg('*** leaving example_griddata');
    gcmfaces_msg('===============================================');
end;
1	function []=example_griddata();
2	%object: a griddata example within gcmfaces
3
4	%%%%%%%%%%%%%%%%%
5	%load grid:
6	%%%%%%%%%%%%%%%%%
7
8	gcmfaces_global;
9
10	if isempty(mygrid);
11	grid_load;
12	end;
13
14	if myenv.verbose>0;
15	gcmfaces_msg('===============================================');
16	gcmfaces_msg(['*** entering example_griddata that will define ' ...
17	'a 2x2 gridded variable (an bathymetry field here), and then ' ...
18	'use griddata to map it to the chosen grid. '],'');
19	end;
20
21	%%%%%%%%%%%%%%%%%%%%%%%
22	%get sample data: subsampled mygrid.Depth
23	if myenv.verbose>0;
24	gcmfaces_msg('* define test case : subsampled bathymetry');
25	end;
26	original_Depth=mygrid.Depth;
27	original_Depth(original_Depth==0)=NaN;
28	lon=mygrid.XC; lat=mygrid.YC;
29	for ff=1:mygrid.nFaces;
30	lon{ff}(1:2:end,:)=NaN;
31	lat{ff}(1:2:end,:)=NaN;
32	original_Depth{ff}(1:2:end,:)=NaN;
33	lon{ff}(:,1:2:end)=NaN;
34	lat{ff}(:,1:2:end)=NaN;
35	original_Depth{ff}(:,1:2:end)=NaN;
36	end;
37	lon=convert2vector(lon);
38	lat=convert2vector(lat);
39	vec=convert2vector(original_Depth);
40
41	%%%%%%%%%%%%%%%%%%%%%%%
42	%make sure to accomodate the two longitude conventions:
43	if myenv.verbose>0;
44	gcmfaces_msg('* accomodate both lon. conventions, then convert to data vector');
45	end;
46	x=[lon-360;lon;lon+360]; y=[lat;lat;lat]; z=[vec;vec;vec];
47	jj=find(~isnan(z)); x=x(jj); y=y(jj); z=z(jj);
48
49	%%%%%%%%%%%%%%%%%%%%%%%
50	%do the interpolation:
51	if myenv.verbose>0;
52	gcmfaces_msg('* map/interpolate data using griddata');
53	end;
54	griddata_Depth=gcmfaces(5);
55	for ii=1:mygrid.nFaces;
56	xi=mygrid.XC{ii}; yi=mygrid.YC{ii};
57	zi = griddata(x',y',z',xi,yi);
58	griddata_Depth{ii}=zi;
59	end;
60
61	%msk=mygrid.hFacC(:,:,1); griddata_Depth(find(msk==0))=NaN;
62
63	%%%%%%%%%%%%%%%%%%%%%%%
64	%illustrate the result:
65
66	figureL;
67	subplot(2,1,1); set(gca,'FontSize',14);
68	[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,original_Depth);
69	pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat;
70	title('input to griddata');
71	subplot(2,1,2); set(gca,'FontSize',14);
72	[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,griddata_Depth);
73	pcolor(X,Y,FLD); axis([-180 180 -90 90]); shading flat;
74	title('output of griddata');
75
76	if myenv.verbose>0;
77	gcmfaces_msg('* note of caution : griddata puts values on land. ');
78	gcmfaces_msg('*** leaving example_griddata');
79	gcmfaces_msg('===============================================');
80	end;