/[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.3 - (hide annotations) (download)
Wed Dec 5 18:38:05 2012 UTC (13 years, 11 months ago) by gforget
Branch: MAIN
Changes since 1.2: +3 -0 lines
- update according to latest revision of gcmfaces_bindata.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 gforget 1.3 mygrid_latlon.dirGrid='none';
69     mygrid_latlon.fileFormat='straight';
70     mygrid_latlon.ioSize=size(lon);
71 gforget 1.1
72     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
73     %interpolate to lat-lon grid:
74     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
75    
76     if myenv.verbose>0;
77     gcmfaces_msg('* interpolate test field to lat-lon grid');
78     end;
79     mygrid=mygrid_latlon; gcmfaces_bindata;
80     veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
81     veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
82    
83     mygrid=mygrid_refgrid; gcmfaces_bindata;
84     vecfld=gcmfaces_interp_2d(fld,veclon,veclat);%go from gcmfaces grid to lat-lon grid
85    
86     mygrid=mygrid_latlon; gcmfaces_bindata;
87     fld_latlon=NaN*convert2array(mygrid.XC);
88     fld_latlon(mytri.kk)=vecfld;
89     fld_latlon=convert2array(fld_latlon);
90    
91     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
92     %interpolate back to ref grid:
93     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
94    
95     if myenv.verbose>0;
96     gcmfaces_msg('* interpolate back to original grid');
97     end;
98     mygrid=mygrid_refgrid; gcmfaces_bindata;
99     veclon=convert2array(mygrid.XC); veclon=veclon(mytri.kk);
100     veclat=convert2array(mygrid.YC); veclat=veclat(mytri.kk);
101    
102     mygrid=mygrid_latlon; gcmfaces_bindata;
103     vecfld=gcmfaces_interp_2d(fld_latlon,veclon,veclat);%go from gcmfaces grid to lat-lon grid
104    
105     mygrid=mygrid_refgrid; gcmfaces_bindata;
106     fld_refgrid=NaN*convert2array(mygrid.XC);
107     fld_refgrid(mytri.kk)=vecfld;
108     fld_refgrid=convert2array(fld_refgrid);
109    
110     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111     %get and plot the lat-lon array:
112     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113    
114     if myenv.verbose>0;
115     gcmfaces_msg('* plot lat-lon grid interpolated field');
116     end;
117     FLD=fld_latlon{1};
118     figure; pcolor(lon,lat,FLD); colorbar;
119     axis(aa); title('interpoate lat-lon map');
120    
121     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
122     %plot difference due to interpolation
123     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
124    
125     if myenv.verbose>0;
126     gcmfaces_msg('* plot error due to back and forth interpolation');
127     end;
128     mygrid=mygrid_refgrid;
129     figure;
130 gforget 1.2 aa=which('m_proj');
131     if isempty(aa);
132     [X,Y,FLD]=convert2pcol(mygrid.XC,mygrid.YC,fld_refgrid-fld);
133     pcolor(X,Y,FLD); shading flat; colorbar;
134     else;
135     m_map_gcmfaces(fld_refgrid-fld);
136     end;
137 gforget 1.1 title('error due to back and forth interpolation');
138    
139    
140     if myenv.verbose>0;
141     gcmfaces_msg('*** leaving example_faces2latlon2faces');
142     gcmfaces_msg('===============================================');
143     end;

  ViewVC Help
Powered by ViewVC 1.1.22