matlab_class/sample_processing/example_faces2latlon2faces.m

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

if 0;
    clear all; grid_load('./',6,'cube'); gcmfaces_bindata; gcmfaces_global;
    lon=[-179:2:179]; lat=[-30:2:30];
    [lat,lon] = meshgrid(lat,lon);
    fld=sin(2*pi/180*mygrid.XC).*sin(2*pi/180*mygrid.YC);
    fld2=gcmfaces_interp_2d(fld,lon,lat);%go from gcmfaces grid to lat-lon grid
end;

gcmfaces_global; global mytri;

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;

%%%%%%%%%%%%%%%%%
%load parameters:
%%%%%%%%%%%%%%%%%

gcmfaces_global;
dir0=[myenv.gcmfaces_dir '/sample_input/'];
if strcmp(choiceV3orV4,'v4');
    nF=5;
    fileFormat='compact';
    dirGrid=[dir0 '../../GRID/'];
else;
    error('files are missing');
    nF=1;
    fileFormat='straight';
    dirGrid=[dir0 '/GRID' choiceV3orV4  '/'];
end;

mygrid=[]; grid_load(dirGrid,nF,fileFormat,1);

%store reference grid
mygrid_refgrid=mygrid;

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