[Dart-dev] [4834] DART/trunk/matlab: Full support for WRF for all routines that do not require a True_State.nc
nancy at ucar.edu
nancy at ucar.edu
Thu Mar 31 17:12:11 MDT 2011
Revision: 4834
Author: thoar
Date: 2011-03-31 17:12:11 -0600 (Thu, 31 Mar 2011)
Log Message:
-----------
Full support for WRF for all routines that do not require a True_State.nc
Those routines are untested until I get someone to perform a perfect model experiment.
The filename is now annotated in the title on a separate line.
The parse.m routine breaks a line into words - and then creates a cell array.
Modified Paths:
--------------
DART/trunk/matlab/CheckModelCompatibility.m
DART/trunk/matlab/GetBgridInfo.m
DART/trunk/matlab/GetWRFInfo.m
DART/trunk/matlab/PlotCorrel.m
DART/trunk/matlab/PlotEnsErrSpread.m
DART/trunk/matlab/PlotEnsMeanTimeSeries.m
DART/trunk/matlab/PlotJeffCorrel.m
DART/trunk/matlab/PlotPhaseSpace.m
DART/trunk/matlab/PlotSawtooth.m
DART/trunk/matlab/PlotTotalErr.m
DART/trunk/matlab/SetCopyID2.m
DART/trunk/matlab/get_copy_index.m
DART/trunk/matlab/plot_correl.m
DART/trunk/matlab/plot_ens_err_spread.m
DART/trunk/matlab/plot_ens_mean_time_series.m
DART/trunk/matlab/plot_ens_time_series.m
DART/trunk/matlab/plot_jeff_correl.m
DART/trunk/matlab/plot_phase_space.m
DART/trunk/matlab/plot_sawtooth.m
DART/trunk/matlab/plot_total_err.m
DART/trunk/matlab/plot_var_var_correl.m
Added Paths:
-----------
DART/trunk/matlab/parse.m
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModelCompatibility.m
===================================================================
--- DART/trunk/matlab/CheckModelCompatibility.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/CheckModelCompatibility.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -22,8 +22,6 @@
if (nargin == 1) % better be a pinfo struct with at least these fields
file1 = arg1.truth_file; % string
file2 = arg1.diagn_file; % string
- time1 = arg1.truth_time; % [a,b] array
- time2 = arg1.diagn_time; % [a,b] array
pinfo_out = arg1;
elseif (nargin == 2) % pair of filenames
file1 = arg1; % truth_file
@@ -38,10 +36,10 @@
if ( exist(file2,'file') ~= 2 ), error('(file2) %s does not exist.',file2); end
% set this up for later
-pinfo_out = setfield(pinfo_out, 'truth_time', [-1,-1]);
-pinfo_out = setfield(pinfo_out, 'diagn_time', [-1,-1]);
+pinfo_out.truth_time = [-1,-1];
+pinfo_out.diagn_time = [-1,-1];
-% Get some information from the file1
+%% Get some information from the file1
tmodel = nc_attget(file1,nc_global,'model');
if (isempty(tmodel))
@@ -50,15 +48,18 @@
tnum_copies = dim_length(file1,'copy');
tnum_times = dim_length(file1,'time');
-ttimes = nc_varget(file1,'time');
+times = nc_varget( file1,'time');
+timeunits = nc_attget( file1,'time','units');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+ttimes = times + timeorigin;
[tnum_vars,tdims] = ModelDimension(file1,tmodel);
if ( tnum_vars <= 0 )
error('Unable to determine resolution of %s.',file1)
end
-
-% Get some information from the file2
+%% Get some information from the file2
dmodel = nc_attget(file1,nc_global,'model');
if (isempty(dmodel))
@@ -67,7 +68,11 @@
dnum_copies = dim_length(file2,'copy');
dnum_times = dim_length(file2,'time');
-dtimes = nc_varget(file2,'time');
+times = nc_varget( file2,'time');
+timeunits = nc_attget( file2,'time','units');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+dtimes = times + timeorigin;
[dnum_vars,ddims] = ModelDimension(file2,dmodel);
if ( dnum_vars <= 0 )
@@ -80,7 +85,7 @@
fprintf('%s has model %s\n',file2,dmodel)
error('no No NO ... models must be the same')
end
-pinfo_out = setfield(pinfo_out, 'model', tmodel);
+pinfo_out.model = tmodel;
if (prod(tnum_vars) ~= prod(dnum_vars))
fprintf('%s has %d state variables\n',file1,prod(tnum_vars))
@@ -223,12 +228,19 @@
function [x,y] = ModelDimension(fname,modelname)
-
+% Check the base geometry of the grid
x = 0;
y = NaN;
switch lower(modelname)
+ case 'wrf'
+ diminfo = nc_getdiminfo(fname, 'west_east_d01'); dnum_lons = diminfo.Length;
+ diminfo = nc_getdiminfo(fname,'south_north_d01'); dnum_lats = diminfo.Length;
+ diminfo = nc_getdiminfo(fname, 'bottom_top_d01'); dnum_lvls = diminfo.Length;
+ x = 3;
+ y = [dnum_lons dnum_lats dnum_lvls];
+
case 'cam'
dnum_lons = dim_length(fname,'lon');
dnum_lats = dim_length(fname,'lat');
Modified: DART/trunk/matlab/GetBgridInfo.m
===================================================================
--- DART/trunk/matlab/GetBgridInfo.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/GetBgridInfo.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -173,9 +173,14 @@
[var3_lat, var3_latind] = GetLatitude( var3, TmpJ, VelJ, var1_lat);
[var3_lon, var3_lonind] = GetLongitude(var3, TmpI, VelI, var1_lon);
- % query for ensemble member
- s1 = input('Input ensemble member metadata STRING. <cr> for ''true state'' ','s');
- if isempty(s1), ens_mem = 'true state'; else ens_mem = s1; end
+ % query for ensemble member string
+ metadata = nc_varget(fname,'CopyMetaData');
+ [N,M] = size(metadata);
+ cell_array = mat2cell(metadata, ones(1,N), M);
+ ens_mem = strtrim(cell_array{1});
+ str1 = sprintf('Input ensemble member metadata STRING. <cr> for ''%s'' ',ens_mem);
+ s1 = input(str1,'s');
+ if ~ isempty(s1), ens_mem = s1; end
% query for line type
s1 = input('Input line type string. <cr> for ''k-'' ','s');
Modified: DART/trunk/matlab/GetWRFInfo.m
===================================================================
--- DART/trunk/matlab/GetWRFInfo.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/GetWRFInfo.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -166,19 +166,17 @@
case 'plotsawtooth'
pgvar = GetVarString(pinfo_in.vars);
- [level, lvlind] = GetLevel(fname,pgvar); % Determine level and index
- [lat , latind] = GetLatitude( pgvar,TmpJ,VelJ);
- [lon , lonind] = GetLongitude(pgvar,TmpI,VelI);
- %[copy , lonind] = GetCopies(pgvar,copy);
- copyindices = SetCopyID(fname);
+ [level, lvlind] = GetLevel(pinfo_in.prior_file, pgvar);
+ [lat, lon, latind, lonind] = GetLatLon(pinfo_in.prior_file, pgvar);
+ [copyindices, copymetadata]= SetCopyID2(pinfo_in.prior_file);
copy = length(copyindices);
pinfo.model = model;
+ pinfo.truth_file = pinfo_in.truth_file;
+ pinfo.prior_file = pinfo_in.prior_file;
+ pinfo.posterior_file = pinfo_in.posterior_file;
pinfo.times = dates;
pinfo.var_names = pgvar;
- %pinfo.truth_file = [];
- %pinfo.prior_file = pinfo.prior_file;
- %pinfo.posterior_file = pinfo.posterior_file;
pinfo.level = level;
pinfo.levelindex = lvlind;
pinfo.latitude = lat;
@@ -188,29 +186,44 @@
pinfo.copies = copy;
pinfo.copyindices = copyindices;
+ % Hui has a habit of cutting out the rest of the ensemble members
+ % from the posterior, but leaving them in the prior. Thanks.
+ % So now I have to figure out if the posterior and prior copy metadata match.
+
+ for i = 1:copy,
+ copyi = get_copy_index(pinfo_in.posterior_file,copymetadata{i});
+ pstruct.postcopyindices = copyi;
+ end
+
+ if ( any(pstruct.postcopyindices < 1) )
+ error('The prior copy does not exist in the posterior file.')
+ end
+
case 'plotphasespace'
disp('Getting information for the ''X'' variable.')
var1 = GetVarString(pinfo_in.vars);
[var1_lvl, var1_lvlind] = GetLevel(fname,var1);
- [var1_lat, var1_latind] = GetLatitude( var1, TmpJ, VelJ);
- [var1_lon, var1_lonind] = GetLongitude(var1, TmpI, VelI);
+ [var1_lat, var1_lon, var1_latind, var1_lonind] = GetLatLon(fname, var1);
disp('Getting information for the ''Y'' variable.')
var2 = GetVarString(pinfo_in.vars, var1 );
[var2_lvl, var2_lvlind] = GetLevel(fname,var2,var1_lvl);
- [var2_lat, var2_latind] = GetLatitude( var2, TmpJ, VelJ, var1_lat);
- [var2_lon, var2_lonind] = GetLongitude(var2, TmpI, VelI, var1_lon);
+ [var2_lat, var2_lon, var2_latind, var2_lonind] = GetLatLon(fname, var2);
disp('Getting information for the ''Z'' variable.')
var3 = GetVarString(pinfo_in.vars, var1 );
[var3_lvl, var3_lvlind] = GetLevel(fname,var3,var1_lvl);
- [var3_lat, var3_latind] = GetLatitude( var3, TmpJ, VelJ, var1_lat);
- [var3_lon, var3_lonind] = GetLongitude(var3, TmpI, VelI, var1_lon);
+ [var3_lat, var3_lon, var3_latind, var3_lonind] = GetLatLon(fname, var3);
- % query for ensemble member
- s1 = input('Input ensemble member metadata STRING. <cr> for ''true state'' ','s');
- if isempty(s1), ens_mem = 'true state'; else ens_mem = s1; end
+ % query for ensemble member string
+ metadata = nc_varget(fname,'CopyMetaData');
+ [N,M] = size(metadata);
+ cell_array = mat2cell(metadata, ones(1,N), M);
+ ens_mem = strtrim(cell_array{1});
+ str1 = sprintf('Input ensemble member metadata STRING. <cr> for ''%s'' ',ens_mem);
+ s1 = input(str1,'s');
+ if ~ isempty(s1), ens_mem = s1; end
% query for line type
s1 = input('Input line type string. <cr> for ''k-'' ','s');
Modified: DART/trunk/matlab/PlotCorrel.m
===================================================================
--- DART/trunk/matlab/PlotCorrel.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotCorrel.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -86,11 +86,11 @@
clf;
contour(correl,-1:0.2:1);
- s1 = sprintf('%s Correlation of variable %s index %d, T = %d of %s', ...
- model, pinfo.base_var, base_var_index, base_time, pinfo.fname);
+ s1 = sprintf('%s Correlation of variable %s index %d, T = %d', ...
+ model, pinfo.base_var, base_var_index, base_time);
s2 = sprintf('against all variables, all times, %d ensemble members', ...
size(state_var,2));
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel('time (timestep #)')
ylabel('state variable (index)')
set(gca,'YTick',1:num_vars)
@@ -143,13 +143,13 @@
contour(lons,lats,correl,-1:0.2:1); hold on;
plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
'MarkerSize',12,'MarkerFaceColor','k');
- s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+ s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
model, pinfo.base_var, pinfo.base_lvl, ...
- pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+ pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
pinfo.comp_var, pinfo.comp_lvl, nmembers);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
worldmap;
@@ -290,13 +290,13 @@
contour(lons,lats,correl,-1:0.2:1); hold on;
plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
'MarkerSize',12,'MarkerFaceColor','k');
- s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+ s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
model, pinfo.base_var, pinfo.base_lvl, ...
- pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+ pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
pinfo.comp_var, pinfo.comp_lvl, nmembers);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
worldmap;
@@ -337,13 +337,13 @@
contour(lons,lats,correl,-1:0.2:1); hold on;
plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
'MarkerSize',12,'MarkerFaceColor','k');
- s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+ s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
model, pinfo.base_var, pinfo.base_lvl, ...
- pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+ pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
pinfo.comp_var, pinfo.comp_lvl, nmembers);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
worldmap;
Modified: DART/trunk/matlab/PlotEnsErrSpread.m
===================================================================
--- DART/trunk/matlab/PlotEnsErrSpread.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotEnsErrSpread.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -107,9 +107,8 @@
subplot(3, 1, j);
plot(times,err, 'b', ...
times,ens_spread(:, ivar), 'r');
- s1 = sprintf('%s model Var %d Ensemble Error Spread for %s', ...
- tmodel, ivar, pinfo.diagn_file);
- title(s1,'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s model Var %d Ensemble Error Spread', tmodel, ivar);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
legend(string1,string2,0)
legend boxoff
xlabel(sprintf('model time (%d timesteps)',tnum_times))
@@ -139,9 +138,8 @@
subplot(length(pinfo.var_inds), 1, iplot);
plot(times,err, 'b', ...
times,ens_spread(:, ivar), 'r');
- s1 = sprintf('%s model Var %d Ensemble Error Spread for %s', ...
- tmodel, ivar, pinfo.diagn_file);
- title(s1,'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s model Var %d Ensemble Error Spread', tmodel, ivar);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
legend(string1,string2,0)
legend boxoff
xlabel(sprintf('model time (%d timesteps)',tnum_times))
@@ -170,15 +168,11 @@
string2 = ['time-mean Ensemble Spread = ' num2str(spreadTotal)];
plot(times,err, 'b', times,ens_spread, 'r');
- s1 = sprintf('%s model Var %s Ensemble Error Spread for %s', ...
- tmodel, pinfo.var, pinfo.diagn_file);
- title(s1,'interpreter','none','fontweight','bold');
- title({ ...
- sprintf('Ensemble Mean Error, Ensemble Spread %s ''%s'' for %s', ...
- tmodel, pinfo.var, pinfo.diagn_file), ...
- sprintf('level %d lat %.2f lon %.2f',pinfo.level, pinfo.latitude, ...
- pinfo.longitude)}, 'interpreter','none','fontweight','bold');
+ s1 = sprintf('Ensemble Mean Error, Ensemble Spread %s ''%s''',tmodel,pinfo.var);
+ s2 = sprintf('level %d lat %.2f lon %.2f', ...
+ pinfo.level, pinfo.latitude, pinfo.longitude);
+ title({s1, s2, pinfo.diagn_file},'interpreter','none','fontweight','bold');
legend(string1,string2,0);
legend boxoff
Modified: DART/trunk/matlab/PlotEnsMeanTimeSeries.m
===================================================================
--- DART/trunk/matlab/PlotEnsMeanTimeSeries.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotEnsMeanTimeSeries.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -131,8 +131,8 @@
legend('Ensemble Mean','True State',0);
end
- title(sprintf('%s Variable %d of %s',pinfo.model,ivar,pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Variable %d',pinfo.model,ivar);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
legend boxoff
end
@@ -178,8 +178,8 @@
hold on; plot(times,truth,'b'); hold off;
legend('Ensemble Mean','True State',0)
end
- title(sprintf('%s Variable %d of %s',pinfo.model,ivar,pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Variable %d',pinfo.model,ivar);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
legend boxoff
end
@@ -205,7 +205,7 @@
pinfo.model, pinfo.var, pinfo.diagn_file);
s2 = sprintf('level %d lat %.2f lon %.2f', ...
pinfo.level, pinfo.latitude, pinfo.longitude);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.diagn_file},'interpreter','none','fontweight','bold')
if (have_truth)
truth = GetCopy(pinfo.truth_file, truth_index, pinfo , ...
@@ -244,8 +244,7 @@
subplot(2,1,2)
plot(times,ens_mean,'r','LineWidth',2);
legend('Ensemble Mean', 0)
- s1 = sprintf('%s model ''%s'' %s Ensemble Mean ', ...
- pinfo.model, pinfo.var, pinfo.diagn_file);
+ s1 = sprintf('%s model ''%s'' Ensemble Mean', pinfo.model, pinfo.var);
s2 = sprintf('level index %d lat %.2f lon %.2f', ...
pinfo.levelindex, pinfo.latitude, pinfo.longitude);
@@ -254,13 +253,13 @@
pinfo.truth_time(1), pinfo.truth_time(2)) ;
hold on; plot(times,truth,'b','LineWidth',2); hold off;
legend('Ensemble Mean','True State',0);
- s1 = sprintf('%s model ''%s'' %s Truth and Ensemble Mean ', ...
- pinfo.model, pinfo.var, pinfo.diagn_file);
+ s1 = sprintf('%s model ''%s'' Truth and Ensemble Mean', ...
+ pinfo.model, pinfo.var);
end
%plot(times,ens_mean,'r','LineWidth',2); % again - on top
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('time (%s) %d timesteps',timeunits, num_times))
ylabel(varunits)
legend boxoff
Modified: DART/trunk/matlab/PlotJeffCorrel.m
===================================================================
--- DART/trunk/matlab/PlotJeffCorrel.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotJeffCorrel.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -46,7 +46,7 @@
switch lower(model)
- case {'fms_bgrid','pe2lyr','mitgcm_ocean'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
clf;
@@ -69,11 +69,11 @@
model, pinfo.base_var, pinfo.base_time, pinfo.base_lvl, ...
pinfo.base_lat, pinfo.base_lon);
- s2 = sprintf('with ''%s'', lvl = %d, lat = %.2f, lon= %.2f, %d ensemble members -- %s', ...
+ s2 = sprintf('with ''%s'', lvl = %d, lat = %.2f, lon= %.2f, %d ensemble members', ...
pinfo.comp_var, pinfo.comp_lvl, pinfo.comp_lat, pinfo.comp_lon, ...
- nmembers, pinfo.fname);
+ nmembers);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel(sprintf('time (%s) %d timesteps',timeunits, num_times))
ylabel('correlation')
@@ -83,7 +83,6 @@
plot([pinfo.base_time pinfo.base_time],[ -1 1 ],'k:', ...
[ax(1) ax(2)],[ 0 0 ],'k:')
-
otherwise
num_vars = dim_length(pinfo.fname,'StateVariable');
@@ -123,8 +122,8 @@
s1 = sprintf('%s Correlation of variable %s %d, T = %d, with variable %s %d', ...
model, pinfo.base_var, pinfo.base_var_index, pinfo.base_time, ...
pinfo.state_var, pinfo.state_var_index);
- s2 = sprintf('%d ensemble members -- %s', nmembers, pinfo.fname);
- title({s1,s2},'interpreter','none','fontweight','bold')
+ s2 = sprintf('%d ensemble members', nmembers);
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
xlabel('time (timestep #)')
ylabel('correlation')
Modified: DART/trunk/matlab/PlotPhaseSpace.m
===================================================================
--- DART/trunk/matlab/PlotPhaseSpace.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotPhaseSpace.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -163,11 +163,8 @@
end
legend boxoff
- case {'fms_bgrid','pe2lyr','mitgcm_ocean'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
- disp('PlotPhaseSpace')
- pinfo
-
ens_mem_id = get_copy_index(pinfo.fname, pinfo.ens_mem); % errors out if no ens_mem
x = Get1Copy(pinfo.fname, pinfo.var1name, ens_mem_id, ...
@@ -190,19 +187,22 @@
% title(sprintf('%s ensemble member %d of %s',model,ens_mem_id,pinfo.fname), ...
% 'interpreter','none','fontweight','bold')
title('The Attractor','fontweight','bold')
- xlabel(sprintf('%s %d %.2f %.2f', ...
- pinfo.var1name, pinfo.var1_lvl, pinfo.var1_lat, pinfo.var1_lon))
- ylabel(sprintf('%s %d %.2f %.2f', ...
- pinfo.var2name, pinfo.var2_lvl, pinfo.var2_lat, pinfo.var2_lon))
- zlabel(sprintf('%s %d %.2f %.2f', ...
- pinfo.var3name, pinfo.var3_lvl, pinfo.var3_lat, pinfo.var3_lon))
+ hx = xlabel(sprintf('%s %d %.2f %.2f', ...
+ pinfo.var1name, pinfo.var1_lvl, pinfo.var1_lat, pinfo.var1_lon));
+ hy = ylabel(sprintf('%s %d %.2f %.2f', ...
+ pinfo.var2name, pinfo.var2_lvl, pinfo.var2_lat, pinfo.var2_lon));
+ hz = zlabel(sprintf('%s %d %.2f %.2f', ...
+ pinfo.var3name, pinfo.var3_lvl, pinfo.var3_lat, pinfo.var3_lon));
+ set(hx,'interpreter','none')
+ set(hy,'interpreter','none')
+ set(hz,'interpreter','none')
s = sprintf('%s %s %s %s %s %s', pinfo.var1name, pinfo.var2name, pinfo.var3name, ...
model, pinfo.fname, pinfo.ens_mem);
h = legend(s,0);
[legh, objh, outh, outm] = legend;
set(objh(1),'interpreter','none')
-
+
else
% Must add salient information to the legend.
% legh handle to the legend axes
Modified: DART/trunk/matlab/PlotSawtooth.m
===================================================================
--- DART/trunk/matlab/PlotSawtooth.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotSawtooth.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -40,7 +40,7 @@
% pinfo.longitude = 45.67;
% PlotSawtooth( pinfo )
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+%% DART software - Copyright � 2004 - 2010 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
% http://www.image.ucar.edu/DAReS/DART/DART_download
%
@@ -54,6 +54,7 @@
pstruct = CheckModelCompatibility(pinfo.prior_file, pinfo.posterior_file);
pinfo.prior_times = pstruct.truth_time;
pinfo.posterior_times = pstruct.diagn_time;
+pinfo.model = pstruct.model;
% Get some information from the truth_file, if it exists.
if ( exist(pinfo.truth_file,'file') == 2 )
@@ -68,20 +69,17 @@
truth = [];
end
-% Get some information from the prior_file
+%% Get some information from the prior_file
+% CheckModelCompatibility ensures that getting this from one is sufficient.
prior.model = nc_attget(pinfo.prior_file, nc_global, 'model');
prior.num_times = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
-prior.times = nc_varget(pinfo.prior_file, 'time', ...
+times = nc_varget(pinfo.prior_file, 'time', ...
pinfo.prior_times(1)-1, prior.num_times);
+timeunits = nc_attget(pinfo.prior_file,'time','units');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+prior.times = times + timeorigin;
-% Get some information from the posterior_file
-post.model = nc_attget(pinfo.posterior_file, nc_global, 'model');
-post.num_times = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
-post.times = nc_varget(pinfo.posterior_file, 'time', ...
- pinfo.posterior_times(1)-1, post.num_times);
-post.timeunits = nc_attget(pinfo.posterior_file, 'time', 'units');
-post.calendar = nc_attget(pinfo.posterior_file, 'time', 'calendar');
-
% Get the indices for the true state and ensemble mean
% The metadata is queried to determine which "copy" is appropriate.
prior.ens_mean_index = get_copy_index(pinfo.prior_file, 'ensemble mean');
@@ -96,217 +94,215 @@
pinfo.xax = x(:);
metadata = nc_varget(pinfo.prior_file,'CopyMetaData'); % get all the metadata
-pinfo
-
if isfield(pinfo,'var_inds')
PlotGivenIndices( pinfo, prior, post, truth, metadata);
elseif isfield(pinfo,'var_names')
- PlotGivenVariable(pinfo, prior, post, truth, metadata);
+ PlotGivenVariable(pinfo, truth, metadata);
end
-%------------------------------------------------------------------------------
-% Plot given an index into the state-space vector
-%------------------------------------------------------------------------------
+
function PlotGivenIndices(pinfo, prior, post, truth, metadata)
- nfigs = length(pinfo.var_inds); % each variable gets its own figure
- iplot = 0;
+%% Plot given an index into the state-space vector
- for ivar = pinfo.var_inds,
+nfigs = length(pinfo.var_inds); % each variable gets its own figure
+iplot = 0;
- % Get the data from the netcdf files
+for ivar = pinfo.var_inds,
- po_ens_mean = get_var_series(pinfo.posterior_file, pinfo.var, ...
- post.ens_mean_index, ivar, ...
- pinfo.posterior_times(1), pinfo.posterior_times(2));
- pr_ens_mean = get_var_series(pinfo.prior_file, pinfo.var, ...
- prior.ens_mean_index, ivar, ...
- pinfo.prior_times(1), pinfo.prior_times(2));
+ % Get the data from the netcdf files
- % Now we paste them together in a clever way to show
- % the effect of the assimilation
- ens_mean( 1,:) = pr_ens_mean;
- ens_mean( 2,:) = po_ens_mean;
- a = ens_mean(:);
+ po_ens_mean = get_var_series(pinfo.posterior_file, pinfo.var, ...
+ post.ens_mean_index, ivar, ...
+ pinfo.posterior_times(1), pinfo.posterior_times(2));
+ pr_ens_mean = get_var_series(pinfo.prior_file, pinfo.var, ...
+ prior.ens_mean_index, ivar, ...
+ pinfo.prior_times(1), pinfo.prior_times(2));
- % Plot the true trajectory if it exists; the ens mean; annotate
+ % Now we paste them together in a clever way to show
+ % the effect of the assimilation
+ ens_mean( 1,:) = pr_ens_mean;
+ ens_mean( 2,:) = po_ens_mean;
+ a = ens_mean(:);
- iplot = iplot + 1;
- figure(iplot); clf;
+ % Plot the true trajectory if it exists; the ens mean; annotate
- if ( exist(pinfo.truth_file,'file') == 2 )
- true_trajectory = get_var_series(pinfo.truth_file, pinfo.var, truth.truth_index, ivar);
+ iplot = iplot + 1;
+ figure(iplot); clf;
+
+ if ( exist(pinfo.truth_file,'file') == 2 )
+ true_trajectory = get_var_series(pinfo.truth_file, pinfo.var, truth.truth_index, ivar);
h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
h = plot(pinfo.xax, a, 'k-','linewidth',2.0);
legend('truth','ensemble mean')
- else
+ else
h = plot(pinfo.xax, a, 'k-','linewidth',2.0); hold on;
legend('ensemble mean')
- end
+ end
- ylabel(sprintf('''%s'' index %d', pinfo.var, ivar ))
- xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
- title(sprintf('%s Trajectories',pinfo.model), ...
- 'interpreter','none','fontweight','bold')
+ ylabel(sprintf('''%s'' index %d', pinfo.var, ivar ))
+ xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
+ title(sprintf('%s Trajectories',pinfo.model), ...
+ 'interpreter','none','fontweight','bold')
- % Now check to see if we are overlaying any individual ensemble members.
- % if pinfo.copyindices = [], nothing happens.
+ % Now check to see if we are overlaying any individual ensemble members.
+ % if pinfo.copyindices = [], nothing happens.
- ens_colors = get(gca,'ColorOrder'); % trying to cycle through colors
- ncolors = size(ens_colors,1) - 1; % last one is black, already used.
+ ens_colors = get(gca,'ColorOrder'); % trying to cycle through colors
+ ncolors = size(ens_colors,1) - 1; % last one is black, already used.
- nmem = 0;
- for imem = pinfo.copyindices
- nmem = nmem + 1;
+ nmem = 0;
+ for imem = pinfo.copyindices
+ nmem = nmem + 1;
%str1 = sprintf('ensemble member %d',imem);
str1 = deblank(metadata(imem,:));
- copy_index = get_copy_index(pinfo.prior_file, str1);
- po_series = get_var_series(pinfo.posterior_file, pinfo.var, ...
- copy_index, ivar, pinfo.posterior_times(1), pinfo.posterior_times(2));
- pr_series = get_var_series(pinfo.prior_file, pinfo.var, ...
- copy_index, ivar, pinfo.prior_times(1), pinfo.prior_times(2));
+ copy_index = get_copy_index(pinfo.prior_file, str1);
+ po_series = get_var_series(pinfo.posterior_file, pinfo.var, ...
+ copy_index, ivar, pinfo.posterior_times(1), pinfo.posterior_times(2));
+ pr_series = get_var_series(pinfo.prior_file, pinfo.var, ...
+ copy_index, ivar, pinfo.prior_times(1), pinfo.prior_times(2));
ens_member(1,:) = pr_series;
ens_member(2,:) = po_series;
b = ens_member(:);
- hold on;
+ hold on;
memcolor = 1 + mod(nmem-1,ncolors); % cycles through colors [1,6]
h = plot(pinfo.xax, b,'linewidth',0.5,'Color',ens_colors(memcolor,:));
[legh, objh, outh, outm] = legend;
nlines = length(outm);
- outm{nlines+1} = str1;
+ outm{nlines+1} = str1;
[legh, objh, outh, outm] = legend([outh; h],outm,0);
- end
- legend boxoff
-
end
+ legend boxoff
-%------------------------------------------------------------------------------
-% Plot given an index into the state-space vector
-%------------------------------------------------------------------------------
+end
-function PlotGivenVariable(pinfo, prior, post, truth, metadata)
- var_names = strread(pinfo.var_names,'%s','delimiter',' ');
- nfigs = length(var_names); % each variable gets its own figure
- iplot = 0;
+function PlotGivenVariable(pinfo, truth, metadata)
+%% Plot given an index into the state-space vector
- for ivar = 1:nfigs
+var_names = parse(pinfo.var_names);
- iplot = iplot + 1;
- figure(iplot); clf;
+nfigs = length(var_names); % each variable gets its own figure
+iplot = 0;
- % Get the data from the netcdf files
+for ivar = 1:nfigs
- vname = var_names{ivar};
-
- [start, count] = GetNCindices(pinfo,'prior',vname);
+ iplot = iplot + 1;
+ figure(iplot); clf;
- ens_colors = get(gca,'ColorOrder'); % trying to cycle through colors
- ncolors = size(ens_colors,1) - 1; % last one is black, already used.
+ % Get the data from the netcdf files
- for i = 1:pinfo.copies
+ vname = var_names{ivar};
+
+ [start, count] = GetNCindices(pinfo,'prior',vname);
- imem = pinfo.copyindices(i);
+ ens_colors = get(gca,'ColorOrder'); % trying to cycle through colors
+ ncolors = size(ens_colors,1) - 1; % last one is black, already used.
- varinfo = nc_getvarinfo(pinfo.prior_file, vname);
+ for i = 1:pinfo.copies
- % Get the Prior
+ imem = pinfo.copyindices(i);
- for j = 1:length(varinfo.Dimension)
- switch( lower(varinfo.Dimension{j}))
- case{'time'}
- start(j) = pinfo.prior_times(1) - 1;
- count(j) = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
- case{'copy'}
- start(j) = imem - 1;
- count(j) = 1;
- otherwise
- end
+ varinfo = nc_getvarinfo(pinfo.prior_file, vname);
+
+ %% Get the Prior
+
+ for j = 1:length(varinfo.Dimension)
+ switch( lower(varinfo.Dimension{j}))
+ case{'time'}
+ start(j) = pinfo.prior_times(1) - 1;
+ count(j) = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
+ case{'copy'}
+ start(j) = imem - 1;
+ count(j) = 1;
+ otherwise
end
+ end
- pr_series = nc_varget(pinfo.prior_file, vname, start, count);
+ pr_series = nc_varget(pinfo.prior_file, vname, start, count);
- % Get the Posterior
+ %% Get the Posterior
- for j = 1:length(varinfo.Dimension)
- switch( lower(varinfo.Dimension{j}))
- case{'time'}
- start(j) = pinfo.posterior_times(1) - 1;
- count(j) = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
- break
- otherwise
- end
+ for j = 1:length(varinfo.Dimension)
+ switch( lower(varinfo.Dimension{j}))
+ case{'time'}
+ start(j) = pinfo.posterior_times(1) - 1;
+ count(j) = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
+ break
+ otherwise
end
+ end
- po_series = nc_varget(pinfo.posterior_file, vname, start, count);
+ po_series = nc_varget(pinfo.posterior_file, vname, start, count);
- % Paste Prior/Posterior into one series
+ %% Paste Prior/Posterior into one series
- ens_member(1,:) = pr_series;
- ens_member(2,:) = po_series;
- b = ens_member(:);
+ ens_member(1,:) = pr_series;
+ ens_member(2,:) = po_series;
+ b = ens_member(:);
- % plot it
+ %% plot it
- hold on;
- memcolor = 1 + mod(i-1,ncolors); % cycles through colors [1,6]
+ hold on;
+ memcolor = 1 + mod(i-1,ncolors); % cycles through colors [1,6]
- str1 = deblank(metadata(imem,:));
+ str1 = deblank(metadata(imem,:));
- if (strcmp(str1,'ensemble mean') == 1)
+ if (strcmp(str1,'ensemble mean') == 1)
h = plot(pinfo.xax, b,'k-','linewidth',2.0);
- else
+ else
h = plot(pinfo.xax, b,'linewidth',0.5,'Color',ens_colors(memcolor,:));
- end
-
- [legh, objh, outh, outm] = legend;
- nlines = length(outm);
- outm{nlines+1} = str1;
- [legh, objh, outh, outm] = legend([outh; h],outm,0);
end
- % Plot the true trajectory if it exists; maybe the ens mean
+ [legh, objh, outh, outm] = legend;
+ nlines = length(outm);
+ outm{nlines+1} = str1;
+ [legh, objh, outh, outm] = legend([outh; h],outm,0);
+ end
- if ( exist(pinfo.truth_file,'file') == 2 )
+ %% Plot the true trajectory if it exists
- varinfo = nc_getvarinfo(pinfo.truth_file, vname);
+ if ( exist(pinfo.truth_file,'file') == 2 )
- for j = 1:length(varinfo.Dimension)
- switch( lower(varinfo.Dimension{j}))
- case{'time'}
- start(j) = 0;
- count(j) = -1;
- case{'copy'}
- start(j) = truth.truth_index - 1;
- count(j) = 1;
- otherwise
- end
+ varinfo = nc_getvarinfo(pinfo.truth_file, vname);
+
+ for j = 1:length(varinfo.Dimension)
+ switch( lower(varinfo.Dimension{j}))
+ case{'time'}
+ start(j) = 0;
+ count(j) = -1;
+ case{'copy'}
+ start(j) = truth.truth_index - 1;
+ count(j) = 1;
+ otherwise
end
+ end
- true_trajectory = nc_varget(pinfo.truth_file, vname, start, count);
+ true_trajectory = nc_varget(pinfo.truth_file, vname, start, count);
- h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
+ h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
- [legh, objh, outh, outm] = legend;
- nlines = length(outm);
- outm{nlines+1} = 'truth';
- [legh, objh, outh, outm] = legend([outh; h],outm,0);
- end
+ [legh, objh, outh, outm] = legend;
+ nlines = length(outm);
+ outm{nlines+1} = 'truth';
+ [legh, objh, outh, outm] = legend([outh; h],outm,0);
+ end
- % annotate
+ % annotate
- ylabel(sprintf('''%s'' (%.3f,%.3fE) level index %d', vname, ...
- pinfo.latitude, pinfo.longitude, pinfo.levelindex ))
- xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
- title(sprintf('%s Trajectories',pinfo.model), ...
- 'interpreter','none','fontweight','bold')
- legend boxoff
+ h = ylabel(sprintf('''%s'' (%.3f,%.3fE) level index %d', vname, ...
+ pinfo.latitude, pinfo.longitude, pinfo.levelindex ));
+ set(h,'Interpreter','none')
+ datetick('x','yyyymmmdd HH:MM')
+ title(sprintf('%s Trajectories',pinfo.model), ...
+ 'Interpreter','none','fontweight','bold')
+ legend boxoff
- end
+end
Modified: DART/trunk/matlab/PlotTotalErr.m
===================================================================
--- DART/trunk/matlab/PlotTotalErr.m 2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotTotalErr.m 2011-03-31 23:12:11 UTC (rev 4834)
@@ -75,9 +75,8 @@
plot(times,err, 'b', times,err_spread, 'r');
legend(string1,string2,0)
legend boxoff
- title(sprintf('%s Total Error over all %d variables for %s',...
- model, num_vars, pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Total Error over all %d variables', model, num_vars);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
ylabel('Total Error')
@@ -116,9 +115,8 @@
plot(times,err, 'b', times,err_spread, 'r');
legend(string1,string2,0)
legend boxoff
- title(sprintf('%s Total Error over all %d variables for %s',...
- model, num_vars, pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Total Error over all %d variables',model,num_vars);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
ylabel('Total Error')
@@ -167,9 +165,8 @@
plot(times,err, 'b', times,err_spread, 'r');
legend(string1,string2,0)
legend boxoff
- title(sprintf('%s Total Error over statevars %d to %d for %s',...
- model, ind1, indN, pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Total Error over statevars %d to %d', model, ind1, indN);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
ylabel('Total Error')
@@ -198,9 +195,8 @@
plot(times,err, 'b', times,err_spread, 'r');
legend(string1,string2,0)
legend boxoff
- title(sprintf('%s Total Error over statevars %d to %d for %s',...
- model, ind1, indN, pinfo.diagn_file), ...
- 'interpreter','none','fontweight','bold')
+ s1 = sprintf('%s Total Error over statevars %d to %d', model, ind1, indN);
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
xlabel(sprintf('model time (%d timesteps)',num_times))
ylabel('Total Error')
@@ -427,7 +423,7 @@
varunits = nc_attget(pinfo.truth_file, pinfo.vars{ivar}, 'units');
- s1 = sprintf('%s %s Ensemble Mean for %s', model,pinfo.vars{ivar},pinfo.diagn_file);
+ s1 = sprintf('%s %s Ensemble Mean', model,pinfo.vars{ivar});
switch lower(pinfo.vars{ivar})
case {'ps'}
@@ -451,7 +447,7 @@
grid on;
xlabel(sprintf('time (%s) %d timesteps',timeunits,num_times))
ylabel(sprintf('global-area-weighted distance (%s)',varunits))
- title(s1,'interpreter','none','fontweight','bold')
+ title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
end
@@ -584,7 +580,7 @@
figure(ivar); clf;
varunits = nc_attget(pinfo.truth_file,pinfo.vars{ivar},'units');
- s1 = sprintf('%s %s Ensemble Mean for %s', model,pinfo.vars{ivar},pinfo.diagn_file);
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list