/[MITgcm]/MITgcm_contrib/gael/matlab_class/sample_processing/example_faces2latlon2faces.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/matlab_class/sample_processing/example_faces2latlon2faces.m

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.1 - (hide annotations) (download)
Tue Jul 31 19:11:50 2012 UTC (14 years, 3 months ago) by gforget
Branch: MAIN
- example_faces2latlon2faces.m : interpolate back and forth between
  gcmfaces and lat-lon grid using gcmfaces_interp_2d.m.

1 gforget 1.1 function []=example_faces2latlon2faces(choiceGrid);
2     %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(2*pi/180*mygrid.XC).*sin(2*pi/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     dir0=[myenv.gcmfaces_dir '/sample_input/'];
26     dirGrid=[dir0 '/GRID' choiceGrid '/'];
27     dirIn=[dir0 '/SAMPLE' choiceGrid '/'];
28     if strcmp(choiceGrid,'v4');
29     nF=5; fileFormat='compact';
30     elseif strcmp(choiceGrid,'v3');
31     nF=1; fileFormat='straight';
32     elseif strcmp(choiceGrid,'cube');
33     nF=6; fileFormat='cube';
34     else;
35     error('unknow grid');
36     end;
37     grid_load(dirGrid,nF,fileFormat);
38    
39     %store reference grid
40     mygrid_refgrid=mygrid;
41    
42     %%%%%%%%%%%%%%%%%%%
43     %field for testing:
44     %%%%%%%%%%%%%%%%%%%
45    
46     if myenv.verbose>0;
47     gcmfaces_msg('* define a sinusoidal field for testing');
48     end;
49     fld=sin(2*pi/180*mygrid.XC).*sin(2*pi/180*mygrid.YC);
50    
51     %%%%%%%%%%%%%%%%%%%%%%%%%%
52     %lat-lon grid for testing:
53     %%%%%%%%%%%%%%%%%%%%%%%%%%
54    
55     if myenv.verbose>0;
56     gcmfaces_msg('* define 2x2 lat-lon grid for testing');
57     end;
58     %define lat-lon grid
59     lon=[-179:2:179]; lat=[-89:2:89]; aa=[-180 180 -90 90];
60     if max(mygrid_refgrid.XC)>180;
61     lon=[1:2:359]; aa=[0 360 -90 90];
62     end;
63     % lon=[-179:2:179]; lat=[-30:2:30];
64     [lat,lon] = meshgrid(lat,lon);
65     %prepare mygrid for lat-lon with no mask
66     mygrid_latlon.nFaces=1;
67     mygrid_latlon.XC=gcmfaces({lon}); mygrid_latlon.YC=gcmfaces({lat});
68    
69     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70     %interpolate to lat-lon grid:
71     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
72    
73     if myenv.verbose>0;
74     gcmfaces_msg('* interpolate test field to lat-lon grid');
75     end;
76     mygrid=mygrid_latlon; gcmfaces_bindata;
77     veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
78     veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
79    
80     mygrid=mygrid_refgrid; gcmfaces_bindata;
81     vecfld=gcmfaces_interp_2d(fld,veclon,veclat);%go from gcmfaces grid to lat-lon grid
82    
83     mygrid=mygrid_latlon; gcmfaces_bindata;
84     fld_latlon=NaN*convert2array(mygrid.XC);
85     fld_latlon(mytri.kk)=vecfld;
86     fld_latlon=convert2array(fld_latlon);
87    
88     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
89     %interpolate back to ref grid:
90     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
91    
92     if myenv.verbose>0;
93     gcmfaces_msg('* interpolate back to original grid');
94     end;
95     mygrid=mygrid_refgrid; gcmfaces_bindata;
96     veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
97     veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
98    
99     mygrid=mygrid_latlon; gcmfaces_bindata;
100     vecfld=gcmfaces_interp_2d(fld_latlon,veclon,veclat);%go from gcmfaces grid to lat-lon grid
101    
102     mygrid=mygrid_refgrid; gcmfaces_bindata;
103     fld_refgrid=NaN*convert2array(mygrid.XC);
104     fld_refgrid(mytri.kk)=vecfld;
105     fld_refgrid=convert2array(fld_refgrid);
106    
107     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
108     %get and plot the lat-lon array:
109     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
110    
111     if myenv.verbose>0;
112     gcmfaces_msg('* plot lat-lon grid interpolated field');
113     end;
114     FLD=fld_latlon{1};
115     figure; pcolor(lon,lat,FLD); colorbar;
116     axis(aa); title('interpoate lat-lon map');
117    
118     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
119     %plot difference due to interpolation
120     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121    
122     if myenv.verbose>0;
123     gcmfaces_msg('* plot error due to back and forth interpolation');
124     end;
125     mygrid=mygrid_refgrid;
126     figure;
127     m_map_gcmfaces(fld_refgrid-fld); colorbar;
128     title('error due to back and forth interpolation');
129    
130    
131     if myenv.verbose>0;
132     gcmfaces_msg('*** leaving example_faces2latlon2faces');
133     gcmfaces_msg('===============================================');
134     end;

  ViewVC Help
Powered by ViewVC 1.1.22