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	function []=example_faces2latlon2faces();
2	%object: interpolate back and forth between gcmfaces and lat-lon grid
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_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	%triangulation variable generated by gcmfaces_bindata
25	global mytri;
26
27	%%%%%%%%%%%%%%%%%%%
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	mygrid_latlon.dirGrid='none';
54	mygrid_latlon.fileFormat='straight';
55	mygrid_latlon.ioSize=size(lon);
56
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	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
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	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
122	if myenv.verbose>0;
123	gcmfaces_msg('*** leaving example_faces2latlon2faces');
124	gcmfaces_msg('===============================================');
125	end;