matlab_class/sample_processing/example_faces2latlon2faces.m

function []=example_faces2latlon2faces();
%object:    interpolate back and forth between gcmfaces and lat-lon grid

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

gcmfaces_global;

if isempty(mygrid);
  grid_load;
end;

if myenv.verbose>0;
    gcmfaces_msg('===============================================');
    gcmfaces_msg(['*** entering example_faces2latlon2faces that will ' ...
        'interpolate a field on the chosen to a lat-lon grid and back' ...
        'to the original grid'],'');
end;

%store reference grid
mygrid_refgrid=mygrid;

%triangulation variable generated by gcmfaces_bindata
global mytri;

%%%%%%%%%%%%%%%%%%%
%field for testing:
%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* define a sinusoidal field for testing');
end;
fld=sin(2*pi/180*mygrid.XC).*sin(2*pi/180*mygrid.YC);

%%%%%%%%%%%%%%%%%%%%%%%%%%
%lat-lon grid for testing:
%%%%%%%%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* define 2x2 lat-lon grid for testing');
end;
%define lat-lon grid
lon=[-179:2:179]; lat=[-89:2:89]; aa=[-180 180 -90 90];
if max(mygrid_refgrid.XC)>180;
    lon=[1:2:359]; aa=[0 360 -90 90];
end;
% lon=[-179:2:179]; lat=[-30:2:30];
[lat,lon] = meshgrid(lat,lon);
%prepare mygrid for lat-lon with no mask
mygrid_latlon.nFaces=1;
mygrid_latlon.XC=gcmfaces({lon}); mygrid_latlon.YC=gcmfaces({lat});
mygrid_latlon.dirGrid='none';
mygrid_latlon.fileFormat='straight';
mygrid_latlon.ioSize=size(lon);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%interpolate to lat-lon grid:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* interpolate test field to lat-lon grid');
end;
mygrid=mygrid_latlon; gcmfaces_bindata; 
veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);

mygrid=mygrid_refgrid; gcmfaces_bindata; 
vecfld=gcmfaces_interp_2d(fld,veclon,veclat);%go from gcmfaces grid to lat-lon grid

mygrid=mygrid_latlon; gcmfaces_bindata; 
fld_latlon=NaN*convert2array(mygrid.XC); 
fld_latlon(mytri.kk)=vecfld;
fld_latlon=convert2array(fld_latlon); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%interpolate back to ref grid:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* interpolate back to original grid');
end;
mygrid=mygrid_refgrid; gcmfaces_bindata; 
veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);

mygrid=mygrid_latlon; gcmfaces_bindata; 
vecfld=gcmfaces_interp_2d(fld_latlon,veclon,veclat);%go from gcmfaces grid to lat-lon grid

mygrid=mygrid_refgrid; gcmfaces_bindata; 
fld_refgrid=NaN*convert2array(mygrid.XC); 
fld_refgrid(mytri.kk)=vecfld;
fld_refgrid=convert2array(fld_refgrid); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%get and plot the lat-lon array:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* plot lat-lon grid interpolated field');
end;
FLD=fld_latlon{1};
figureL; 
subplot(2,1,1); set(gca,'FontSize',14);
pcolor(lon,lat,FLD); axis([-270 270 -100 100]); 
shading flat; colorbar; title('interpolated lat-lon map');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot difference due to interpolation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if myenv.verbose>0;
    gcmfaces_msg('* plot error due to back and forth interpolation');
end;
mygrid=mygrid_refgrid;
subplot(2,1,2); set(gca,'FontSize',14);
[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,fld_refgrid-fld);
pcolor(X,Y,log10(abs(FLD))); axis([-270 270 -100 100]); 
shading flat; caxis([-5 -2]); colorbar; 
title('log10(error) due to back and forth interpolation');

if myenv.verbose>0;
    gcmfaces_msg('*** leaving example_faces2latlon2faces');
    gcmfaces_msg('===============================================');
end;
1	gforget	1.5	function []=example_faces2latlon2faces();
2	gforget	1.1	%object: interpolate back and forth between gcmfaces and lat-lon grid
3
4	gforget	1.6	%%%%%%%%%%%%%%%%%
5			%load grid:
6			%%%%%%%%%%%%%%%%%
7
8			gcmfaces_global;
9
10	gforget	1.7	if isempty(mygrid);
11			grid_load;
12			end;
13	gforget	1.1
14			if myenv.verbose>0;
15			gcmfaces_msg('===============================================');
16			gcmfaces_msg(['*** entering example_faces2latlon2faces that will ' ...
17			'interpolate a field on the chosen to a lat-lon grid and back' ...
18			'to the original grid'],'');
19			end;
20
21			%store reference grid
22			mygrid_refgrid=mygrid;
23
24	gforget	1.6	%triangulation variable generated by gcmfaces_bindata
25			global mytri;
26
27	gforget	1.1	%%%%%%%%%%%%%%%%%%%
28			%field for testing:
29			%%%%%%%%%%%%%%%%%%%
30
31			if myenv.verbose>0;
32			gcmfaces_msg('* define a sinusoidal field for testing');
33			end;
34			fld=sin(2pi/180mygrid.XC).sin(2pi/180*mygrid.YC);
35
36			%%%%%%%%%%%%%%%%%%%%%%%%%%
37			%lat-lon grid for testing:
38			%%%%%%%%%%%%%%%%%%%%%%%%%%
39
40			if myenv.verbose>0;
41			gcmfaces_msg('* define 2x2 lat-lon grid for testing');
42			end;
43			%define lat-lon grid
44			lon=[-179:2:179]; lat=[-89:2:89]; aa=[-180 180 -90 90];
45			if max(mygrid_refgrid.XC)>180;
46			lon=[1:2:359]; aa=[0 360 -90 90];
47			end;
48			% lon=[-179:2:179]; lat=[-30:2:30];
49			[lat,lon] = meshgrid(lat,lon);
50			%prepare mygrid for lat-lon with no mask
51			mygrid_latlon.nFaces=1;
52			mygrid_latlon.XC=gcmfaces({lon}); mygrid_latlon.YC=gcmfaces({lat});
53	gforget	1.3	mygrid_latlon.dirGrid='none';
54			mygrid_latlon.fileFormat='straight';
55			mygrid_latlon.ioSize=size(lon);
56	gforget	1.1
57			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58			%interpolate to lat-lon grid:
59			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60
61			if myenv.verbose>0;
62			gcmfaces_msg('* interpolate test field to lat-lon grid');
63			end;
64			mygrid=mygrid_latlon; gcmfaces_bindata;
65			veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
66			veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
67
68			mygrid=mygrid_refgrid; gcmfaces_bindata;
69			vecfld=gcmfaces_interp_2d(fld,veclon,veclat);%go from gcmfaces grid to lat-lon grid
70
71			mygrid=mygrid_latlon; gcmfaces_bindata;
72			fld_latlon=NaN*convert2array(mygrid.XC);
73			fld_latlon(mytri.kk)=vecfld;
74			fld_latlon=convert2array(fld_latlon);
75
76			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
77			%interpolate back to ref grid:
78			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
79
80			if myenv.verbose>0;
81			gcmfaces_msg('* interpolate back to original grid');
82			end;
83			mygrid=mygrid_refgrid; gcmfaces_bindata;
84			veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
85			veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
86
87			mygrid=mygrid_latlon; gcmfaces_bindata;
88			vecfld=gcmfaces_interp_2d(fld_latlon,veclon,veclat);%go from gcmfaces grid to lat-lon grid
89
90			mygrid=mygrid_refgrid; gcmfaces_bindata;
91			fld_refgrid=NaN*convert2array(mygrid.XC);
92			fld_refgrid(mytri.kk)=vecfld;
93			fld_refgrid=convert2array(fld_refgrid);
94
95			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
96			%get and plot the lat-lon array:
97			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
98
99			if myenv.verbose>0;
100			gcmfaces_msg('* plot lat-lon grid interpolated field');
101			end;
102			FLD=fld_latlon{1};
103	gforget	1.4	figureL;
104			subplot(2,1,1); set(gca,'FontSize',14);
105			pcolor(lon,lat,FLD); axis([-270 270 -100 100]);
106			shading flat; colorbar; title('interpolated lat-lon map');
107	gforget	1.1
108			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
109			%plot difference due to interpolation
110			%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111
112			if myenv.verbose>0;
113			gcmfaces_msg('* plot error due to back and forth interpolation');
114			end;
115			mygrid=mygrid_refgrid;
116	gforget	1.4	subplot(2,1,2); set(gca,'FontSize',14);
117			[X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,fld_refgrid-fld);
118			pcolor(X,Y,log10(abs(FLD))); axis([-270 270 -100 100]);
119			shading flat; caxis([-5 -2]); colorbar;
120			title('log10(error) due to back and forth interpolation');
121	gforget	1.1
122			if myenv.verbose>0;
123			gcmfaces_msg('*** leaving example_faces2latlon2faces');
124			gcmfaces_msg('===============================================');
125			end;