/[MITgcm]/MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m
ViewVC logotype

Annotation of /MITgcm_contrib/gael/profilesMatlabProcessing/ecco_v4/insitu_cost.m

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


Revision 1.11 - (hide annotations) (download)
Sun Jan 17 16:00:24 2016 UTC (9 years, 6 months ago) by gforget
Branch: MAIN
CVS Tags: checkpoint65t
Changes since 1.10: +1 -1 lines
- revise handling of climdir and griddir (incl. new/old locations)
- update call to grid_load (no more arguments)

1 gforget 1.5 function []=insitu_cost(doComp,varargin);
2 gforget 1.4 %object: computes or displays cost function statistics
3 gforget 1.5 %input: doComp states whether to compute & save (1) or load & display (0)
4 gforget 1.1 %optional must take the following form {'name',param1,param2,...} Active ones are
5 gforget 1.5 % {'dirData',dirData} is the data subdirectory name ('profiles/output/' by default)
6 gforget 1.1 % {'listData',data1,data2,...} where data1 e.g. is 'argo_indian.nc'
7     % {'listVar',var1,var2} where var1 and var2 are 'T' and 'S'
8     % {'years',year0,year1} to limit computation to [year0 year1] interval
9 gforget 1.5 % {'dirMat',dirMat} where dirMat is the output subdirectory ('mat/' by default)
10 gforget 1.4 % {'suffMat',suffMat} is the mat file suffix ('all' by default)
11 gforget 1.5 % {'dirTex',dirTex} where dirTex is the tex directory ('tex/' by default)
12     % {'nameTex',nameTex} is the tex file name ('myPlots' by default)
13 gforget 1.1 % {'addToTex',addToTex} states whether (1) or not (0 default) to
14 gforget 1.5 % augment an ongoing tex file (see write2tex.m, dirTex)
15 gforget 1.1 %
16 gforget 1.4 %example: insitu_cost('./',1,{'listData','argo_in*'});
17 gforget 1.1
18     %load paths and grid:
19     %--------------------
20 gforget 1.11 gcmfaces_global; if isempty(mygrid); grid_load; end;
21 gforget 1.1
22     %set more optional paramaters to default values
23 gforget 1.4 global myparms; year0=myparms.yearInAve(1); year1=myparms.yearInAve(2);
24 gforget 1.6 listData={'argo_*_1992_to_2007*','argo_*_2008_to_2010*',...
25     'argo_*_2011_to_2012*','climode_*',...
26     'ctd_*','itp_*','seals_*','xbt_*'};
27 gforget 1.2 listVar={'T','S'};
28 gforget 1.7 listBas={'atlExt','pacExt','indExt','arct'};
29 gforget 1.9 if sum([90 1170]~=mygrid.ioSize)>0; listBas={}; end;
30 gforget 1.7 listBasTxt=''; for bb=1:length(listBas); listBasTxt=[listBasTxt ' ' listBas{bb} '*']; end;
31 gforget 1.2 %shorter test case: listData={'argo_in*'}; listVar={'T'};
32 gforget 1.5 dirData='/profiles/output/';
33     dirMat='mat/'; suffMat='_all';
34     dirTex='tex/'; addToTex=0; nameTex='myPlots';
35 gforget 1.1 %set more optional paramaters to user defined values
36 gforget 1.5 for ii=1:nargin-1;
37 gforget 1.1 if strcmp(varargin{ii}{1},'listData'); listData={varargin{ii}{2:end}};
38     elseif strcmp(varargin{ii}{1},'listVar'); listVar={varargin{ii}{2:end}};
39     elseif strcmp(varargin{ii}{1},'years'); year0=varargin{ii}{2}; year1=varargin{ii}{3};
40 gforget 1.5 elseif strcmp(varargin{ii}{1},'dirMat'); dirMat=varargin{ii}{2};
41 gforget 1.1 elseif strcmp(varargin{ii}{1},'suffMat'); suffMat=['_' varargin{ii}{2}];
42 gforget 1.5 elseif strcmp(varargin{ii}{1},'dirTex'); dirTex=varargin{ii}{2};
43 gforget 1.1 elseif strcmp(varargin{ii}{1},'addToTex'); addToTex=varargin{ii}{2};
44 gforget 1.5 elseif strcmp(varargin{ii}{1},'nameTex'); nameTex=varargin{ii}{2};
45     elseif strcmp(varargin{ii}{1},'dirData'); dirData=varargin{ii}{2};
46 gforget 1.1 else;
47     warning('inputCheck:insitu_cost',...
48     ['unknown option ''' varargin{ii}{1} ''' was ignored']);
49     end;
50     end;
51    
52     if doComp;
53    
54     %time limits:
55     date0=datenum(year0,1,1); date1=datenum(year1,12,31);
56    
57     for vv=1:length(listVar);
58     varCur=listVar{vv};
59    
60     %get the cost square root:
61 gforget 1.5 [MITprof]=MITprof_stats_load([dirData],listData,varCur);
62 gforget 1.1
63     %mask out values that are not in year range:
64     ii=find(MITprof.prof_date<date0|MITprof.prof_date>date1);
65     MITprof.prof(ii,:)=NaN;
66    
67     %mean and median:
68     costCur=[nanmean(MITprof.prof(:).^2) nanmedian(MITprof.prof(:).^2)];
69     fprintf('mean cost for %s: %0.3g (median: %0.3g) \n',varCur,costCur);
70     eval(['costAve' varCur '.mean=costCur(1);']);
71     eval(['costAve' varCur '.median=costCur(2);']);
72    
73     %depth/time mean cost:
74     tmp_date=(MITprof.prof_date-date0)/365;
75     % [x,y,z,n]=MITprof_stats(MITprof.prof_depth,MITprof.prof_depth,...
76     [x,y,z,n]=MITprof_stats(MITprof.prof_depth,[0:200:1000 1500:500:6000],...
77     tmp_date,[0:1/4:1+(year1-year0)],'mean',MITprof.prof.^2);
78     z(n==0)=NaN; n(n==0)=NaN;
79     eval(['depthTimeCost' varCur '.x=x;']); eval(['depthTimeCost' varCur '.y=y;']);
80     eval(['depthTimeCost' varCur '.z=z;']); eval(['depthTimeCost' varCur '.n=n;']);
81    
82     %latitudinal distribution:
83     [x,y,z,n]=MITprof_stats(MITprof.prof_lat,[-90:5:90],MITprof.prof,[-5:0.25:5]);
84     z(n==0)=NaN; n(n==0)=NaN;
85     eval(['misfitDistrib' varCur '.x=x;']); eval(['misfitDistrib' varCur '.y=y;']);
86     eval(['misfitDistrib' varCur '.z=z;']); eval(['misfitDistrib' varCur '.n=n;']);
87    
88 gforget 1.7 %upper ocean temporal distribution
89     for bb=1:length(listBas);
90     for ll=1:3;
91     bbb=listBas{bb};
92     if ll==1; Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
93     elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
94     elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
95     end;
96     msk=v4_basin(bbb);
97     msk(mygrid.YC<Lmin|mygrid.YC>Lmax)=0;
98     %
99     msk=convert2array(msk);
100     MITprof.prof_msk=msk(MITprof.prof_point);
101     MITprofSub=MITprof_subset(MITprof,'msk',1,'depth',[0 700]);
102     %
103     tmp_date=(MITprofSub.prof_date-date0)/365;
104     [x,y,z,n]=MITprof_stats(tmp_date,[0:1/4:1+(year1-year0)],MITprofSub.prof,[-5:0.25:5]);
105     z(n==0)=NaN; n(n==0)=NaN;
106     eval([txt '.x=x;']); eval([txt '.y=y;']);
107     eval([txt '.z=z;']); eval([txt '.n=n;']);
108     end;
109     end;
110    
111 gforget 1.1 end;
112    
113 gforget 1.10 if ~isdir([dirMat 'cost/']); mkdir([dirMat 'cost/']); end;
114     eval(['save ' dirMat 'cost/insitu_cost' suffMat '.mat costAve* misfitDistrib* depthTimeCost* year* ' ...
115 gforget 1.7 'listData listVar listBas ' listBasTxt ';']);
116 gforget 1.1
117     else;%display result
118    
119 gforget 1.9 if isdir([dirMat 'cost/']); dirMat=[dirMat 'cost/']; end;
120    
121 gforget 1.5 eval(['load ' dirMat 'insitu_cost' suffMat '.mat;']);
122 gforget 1.1
123     figureL;
124    
125     for vv=1:length(listVar);
126     varCur=listVar{vv};
127    
128     eval(['costCur=[costAve' varCur '.mean costAve' varCur '.median];']);
129    
130     eval(['x=depthTimeCost' varCur '.x;']); eval(['y=depthTimeCost' varCur '.y;']);
131     eval(['z=depthTimeCost' varCur '.z;']); eval(['n=depthTimeCost' varCur '.n;']);
132 gforget 1.7 z(n<1e2)=NaN; y=y+year0-1;
133 gforget 1.1 subplot(2,length(listVar),vv); depthStretchPlot('pcolor',{y,x,z},[0:200:1000 1500:500:6000],[0 1000 4000]);
134     caxis([0 8]); shading flat; colorbar; ylabel('depth (in m)');
135     title(sprintf('mean cost for %s: %0.3g (median: %0.3g)',varCur,costCur));
136    
137     %latitudinal distribution:
138     eval(['x=misfitDistrib' varCur '.x;']); eval(['y=misfitDistrib' varCur '.y;']);
139     eval(['z=misfitDistrib' varCur '.z;']); eval(['n=misfitDistrib' varCur '.n;']);
140     tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
141     subplot(2,length(listVar),vv+length(listVar)); contourf(x,y,z,[0:0.05:0.5]); grid on;
142     caxis([0 0.5]); colorbar; ylabel('normalized misfit'); xlabel('latitude'); title('pdf');
143    
144     end;
145    
146 gforget 1.3 myCaption={'Cost function (top) for in situ profiles, as a function of depth and time. ',...
147 gforget 1.4 'Distribution of normalized misfits (bottom) as a function of latitude. For T (left) and S (right).'};
148 gforget 1.5 if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
149 gforget 1.1
150 gforget 1.8 if ~isempty(listBas);
151    
152 gforget 1.7 ii=0;
153     ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=3;
154     ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=3;
155     ii=ii+1; listPanels(ii).bas='arct'; listPanels(ii).ll=3;
156     ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=2;
157     ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=2;
158     ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=2;
159     ii=ii+1; listPanels(ii).bas='atlExt'; listPanels(ii).ll=1;
160     ii=ii+1; listPanels(ii).bas='pacExt'; listPanels(ii).ll=1;
161     ii=ii+1; listPanels(ii).bas='indExt'; listPanels(ii).ll=1;
162    
163     for vv=1:length(listVar);
164     varCur=listVar{vv};
165     figureL;
166     for ii=1:9;
167     bbb=listPanels(ii).bas; ll=listPanels(ii).ll;
168     if ll==1; Lmin=-90; Lmax=-25; txt=[bbb '_90S25S_misDis' varCur];
169     elseif ll==2; Lmin=-25; Lmax=25; txt=[bbb '_25S25N_misDis' varCur];
170     elseif ll==3; Lmin=25; Lmax=90; txt=[bbb '_25N90N_misDis' varCur];
171     end;
172    
173     eval(['x=' txt '.x;']); eval(['y=' txt '.y;']);
174     eval(['z=' txt '.z;']); eval(['n=' txt '.n;']);
175     x=x+year0-1;
176     tmp1=(nansum(n,2)>1e4)*ones(1,size(n,2)); z(tmp1==0)=NaN;
177     subplot(3,3,ii); contourf(x,y,z,[0:0.05:0.5]); grid on;
178     caxis([0 0.5]); colorbar; title(txt,'Interpreter','none');
179     %title([txt ', normalized misfit pdf'],'Interpreter','none');
180     if mod(ii,3)==1; ylabel('normalized misfit'); end;
181     %if ii<=6; set(gca,'XTick',[]); end;
182     end;
183    
184     myCaption={'Distribution of normalized misfits per basin (panel) as ',...
185     'a function of latitude, for ',varCur};
186     if addToTex; write2tex([dirTex '/' nameTex '.tex'],2,myCaption,gcf); end;
187    
188     end;
189    
190 gforget 1.1 end;
191    
192 gforget 1.8 end;
193    
194 gforget 1.1

  ViewVC Help
Powered by ViewVC 1.1.22