[Dart-dev] [5066] DART/trunk/matlab: more CAM support. Everythin
nancy at ucar.edu
nancy at ucar.edu
Fri Jul 8 16:39:38 MDT 2011
Revision: 5066
Author: thoar
Date: 2011-07-08 16:39:38 -0600 (Fri, 08 Jul 2011)
Log Message:
-----------
more CAM support. Everything but plot_sawtooth.
Modified Paths:
--------------
DART/trunk/matlab/CheckModel.m
DART/trunk/matlab/GetCamInfo.m
DART/trunk/matlab/PlotCorrel.m
DART/trunk/matlab/PlotJeffCorrel.m
DART/trunk/matlab/PlotPhaseSpace.m
DART/trunk/matlab/PlotTotalErr.m
DART/trunk/matlab/PlotVarVarCorrel.m
DART/trunk/matlab/plot_bins.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_sawtooth.m
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModel.m
===================================================================
--- DART/trunk/matlab/CheckModel.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/CheckModel.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -199,13 +199,7 @@
case 'cam'
- % A more robust way would be to use the netcdf low-level ops:
- % bob = var(f); % bob is a cell array of ncvars
- % name(bob{1}) % is the variable name string
- % bob{1}(:) % is the value of the netcdf variable (no offset/scale)
- % have not yet figured out a way to only use non-coordinate variables.
-
- varnames = {'PS','T','U','V','Q','CLDLIQ','CLDICE'};
+ varnames = get_DARTvars(fname);
num_vars = length(varnames);
nlevels = dim_length(fname,'lev'); % determine # of state variables
Modified: DART/trunk/matlab/GetCamInfo.m
===================================================================
--- DART/trunk/matlab/GetCamInfo.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/GetCamInfo.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -1,4 +1,4 @@
-function pinfo = GetCamInfo(pstruct,routine);
+function pinfo = GetCamInfo(pstruct,fname,routine);
%% GetCamInfo prepares a structure of information needed by the subsequent "routine"
% The information is gathered via rudimentary "input" routines.
%
@@ -17,41 +17,35 @@
% $Revision$
% $Date$
-if (isfield(pstruct,'truth_file') && exist(pstruct.truth_file,'file') )
- fname = pstruct.truth_file;
-elseif (isfield(pstruct,'diagn_file') && exist(pstruct.diagn_file,'file') )
- fname = pstruct.diagn_file;
-elseif (isfield(pstruct,'prior_file') && exist(pstruct.prior_file,'file') )
- fname = pstruct.prior_file;
-elseif (isfield(pstruct,'posterior_file') && exist(pstruct.posterior_file,'file') )
- fname = pstruct.posterior_file;
-elseif (isfield(pstruct,'fname') && exist(pstruct.fname,'file') )
- fname = pstruct.fname;
-end
+if (exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
+pinfo = pstruct;
model = nc_attget(fname,nc_global,'model');
if strcmpi(model,'cam') ~= 1
error('Not so fast, this is not a cam model.')
end
-pinfo = pstruct;
+%% Get the domain-independent information.
+varexist(fname, {'copy','time'})
+
copy = nc_varget(fname,'copy');
times = nc_varget(fname,'time');
+
+% Coordinate between time types and dates
+
+timeunits = nc_attget(fname,'time','units');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+dates = times + timeorigin;
+
ilevel = nc_varget(fname,'ilev'); % interfaces
levels = nc_varget(fname, 'lev'); % midpoints
lon = nc_varget(fname, 'lon');
lat = nc_varget(fname, 'lat');
-% A more robust way would be to use the netcdf low-level ops:
-% bob = var(f); % bob is a cell array of ncvars
-% name(bob{1}) % is the variable name string
-% bob{1}(:) % is the value of the netcdf variable (no
-% offset/scale)
-% have not yet figured out a way to only use non-coordinate variables.
-
-prognostic_vars = {'PS','T','U','V','Q','CLDLIQ','CLDICE'};
+prognostic_vars = get_DARTvars(fname);
num_vars = length(prognostic_vars);
switch lower(deblank(routine))
@@ -64,6 +58,7 @@
[lon , lonind] = GetLongitude(pgvar,lon);
pinfo.model = model;
+ pinfo.times = dates;
pinfo.var = pgvar;
pinfo.level = level;
pinfo.levelindex = lvlind;
@@ -77,7 +72,7 @@
disp('Getting information for the ''base'' variable.')
base_var = GetVar(prognostic_vars);
- [base_time, base_tmeind] = GetTime( base_var,times);
+ [base_time, base_tmeind] = GetTime(dates);
[base_lvl, base_lvlind] = GetLevel( base_var,levels);
[base_lat, base_latind] = GetLatitude( base_var,lat);
[base_lon, base_lonind] = GetLongitude(base_var,lon);
@@ -87,6 +82,7 @@
[comp_lvl, comp_lvlind] = GetLevel( comp_var,levels, base_lvlind);
pinfo.model = model;
+ pinfo.times = dates;
pinfo.base_var = base_var;
pinfo.comp_var = comp_var;
pinfo.base_time = base_time;
@@ -105,7 +101,7 @@
disp('Getting information for the ''base'' variable.')
base_var = GetVar(prognostic_vars);
- [base_time, base_tmeind] = GetTime( base_var,times);
+ [base_time, base_tmeind] = GetTime(dates);
[base_lvl , base_lvlind] = GetLevel( base_var,levels);
[base_lat , base_latind] = GetLatitude( base_var,lat);
[base_lon , base_lonind] = GetLongitude(base_var,lon);
@@ -117,6 +113,7 @@
[comp_lon, comp_lonind] = GetLongitude(comp_var,lon,base_lon);
pinfo.model = model;
+ pinfo.times = dates;
pinfo.base_var = base_var;
pinfo.comp_var = comp_var;
pinfo.base_time = base_time;
@@ -144,6 +141,7 @@
% [ lon, lonind] = GetCopies(pgvar,xxx);
pinfo.model = model;
+ pinfo.times = dates;
pinfo.var_names = pgvar;
pinfo.truth_file = [];
pinfo.prior_file = pstruct.prior_file;
@@ -191,6 +189,7 @@
if isempty(s1), ltype = 'k-'; else ltype = s1; end
pinfo.model = model;
+ pinfo.times = dates;
pinfo.var1name = var1;
pinfo.var2name = var2;
pinfo.var3name = var3;
@@ -237,23 +236,41 @@
pgvar(inds) = '';
-function [time, timeind] = GetTime(pgvar, times, deftime)
+
+function [time, timeind] = GetTime(times, deftime)
%----------------------------------------------------------------------
-if (nargin == 3), time = deftime; else time = mean(times); end
+% Query for the time of interest.
-fprintf('Default time is %f, if this is OK, <cr>;\n',time)
-fprintf('If not, enter a time between %.4f and %.4f, we use the closest.\n', ...
- min(times),max(times))
+ntimes = length(times);
+
+% Determine a sensible default.
+if (nargin == 2),
+ time = deftime;
+ tindex = find(times == deftime);
+else
+ if (ntimes < 2)
+ tindex = round(ntimes/2);
+ else
+ tindex = 1;
+ end
+ time = times(tindex);
+end
+
+fprintf('Default time is %s (index %d), if this is OK, <cr>;\n',datestr(time),tindex)
+fprintf('If not, enter an index between %d and %d \n',1,ntimes)
+fprintf('Pertaining to %s and %s \n',datestr(times(1)),datestr(times(ntimes)))
varstring = input('(no syntax required)\n','s');
-if ~isempty(varstring), time = str2num(varstring); end
+if ~isempty(varstring), tindex = str2num(varstring); end
-d = abs(time - times); % crude distance
-ind = find(min(d) == d); % multiple minima possible
-timeind = ind(1); % use the first one
-time = times(timeind);
+timeinds = 1:ntimes;
+d = abs(tindex - timeinds); % crude distance
+ind = find(min(d) == d); % multiple minima possible
+timeind = ind(1); % use the first one
+time = times(timeind);
+
function [level, lvlind] = GetLevel(pgvar, levels, deflevel)
%----------------------------------------------------------------------
% level and lvlind will not be equal for all models, (and probably
@@ -359,3 +376,24 @@
else
dist = reshape(acos(ang1 + ang2)*r,Dlat2);
end
+
+
+function varexist(filename, varnames)
+%% We already know the file exists by this point.
+% Lets check to make sure that file contains all needed variables.
+% Fatally die if they do not exist.
+
+nvars = length(varnames);
+gotone = ones(1,nvars);
+
+for i = 1:nvars
+ gotone(i) = nc_isvar(filename,varnames{i});
+ if ( ~ gotone(i) )
+ fprintf('\n%s is not a variable in %s\n',varnames{i},filename)
+ end
+end
+
+if ~ all(gotone)
+ error('missing required variable ... exiting')
+end
+
Modified: DART/trunk/matlab/PlotCorrel.m
===================================================================
--- DART/trunk/matlab/PlotCorrel.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/PlotCorrel.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -305,7 +305,7 @@
ax = get(h,'Position');
%set(h,'Position',[ax(1) ax(2) ax(3)/2 ax(4)]);
- case 'pe2lyr'
+ case {'pe2lyr','cam'}
% We are going to correlate one var/time/lvl/lat/lon with
% all other lats/lons for a var/time/lvl
Modified: DART/trunk/matlab/PlotJeffCorrel.m
===================================================================
--- DART/trunk/matlab/PlotJeffCorrel.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/PlotJeffCorrel.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -5,25 +5,11 @@
% time and another variable at all times in an ensemble time sequence.
% The correlation is done across ensemble members.
%
-% PlotVarVarCorrel is intended to be called by 'plot_var_var_correl'
+% PlotJeffCorrel is intended to be called by 'plot_jeff_correl'
%
-% USAGE: PlotVarVarCorrel(fname, base_var_index, base_time, state_var_index)
+% USAGE: PlotJeffCorrel(pinfo)
%
-% pinfo is a structure with the following necessary components:
-% fname
-% base_var_index index of one of the state variables to correlate
-% base_time index of the time of interest (timestep)
-% state_var_index index of the other state variable of interest.
-%
-% Example (lorenz 63 model with 1000 timesteps)
-%%--------------------------------------------------------
-% pinfo.fname = 'Posterior_Diag.nc';
-% pinfo.base_var = 'state';
-% pinfo.base_var_index = 2;
-% pinfo.base_time = 500;
-% pinfo.state_var = 'state';
-% pinfo.state_var_index = 1;
-% PlotVarVarCorrel( pinfo )
+% pinfo is a model-dependent structure.
%% DART software - Copyright 2004 - 2011 UCAR. This open source software is
% provided by UCAR, "as is", without charge, subject to all terms of use at
@@ -35,7 +21,7 @@
% $Revision$
% $Date$
-if (exist(pinfo.fname) ~= 2), error('%s does not exist.',pinfo.fname), end
+if (exist(pinfo.fname,'file') ~= 2), error('%s does not exist.',pinfo.fname), end
% Get some file-specific information.
@@ -46,7 +32,7 @@
switch lower(model)
- case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf','cam'}
clf;
@@ -57,13 +43,12 @@
nmembers = size(comp_mem,2);
correl = jeff_correl(base_mem, pinfo.base_tmeind, comp_mem);
- times = nc_varget(pinfo.fname,'time');
subplot(2,1,1)
PlotLocator(pinfo)
subplot(2,1,2)
- plot(times,correl);
+ plot(pinfo.times,correl);
s1 = sprintf('%s Correlation of ''%s'', T = %d, lvl = %d, lat = %.2f, lon=%.2f', ...
model, pinfo.base_var, pinfo.base_time, pinfo.base_lvl, ...
@@ -74,7 +59,7 @@
nmembers);
title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
- xlabel(sprintf('time (%s) %d timesteps',timeunits, num_times))
+ datetick('x','yyyymmmdd HH:MM')
ylabel('correlation')
% call out the time index in question, and put a corr==0 reference line.
@@ -83,6 +68,7 @@
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');
@@ -149,7 +135,7 @@
-function var = GetEns( fname, var, lvlind, latind, lonind)
+function var = GetEns( fname, varname, lvlind, latind, lonind)
% Gets a time-series of all copies of a prognostic variable
% at a particular 3D location (level, lat, lon).
% Determining just the ensemble members (and not mean, spread ...)
@@ -174,9 +160,9 @@
myinfo.levelindex = lvlind;
myinfo.latindex = latind;
myinfo.lonindex = lonind;
-[start, count] = GetNCindices(myinfo,'diagn',var);
+[start, count] = GetNCindices(myinfo,'diagn',varname);
-bob = nc_varget(fname, var, start, count); % 'bob' is only 2D
+bob = nc_varget(fname, varname, start, count); % 'bob' is only 2D
var = bob(:,copyindices);
Modified: DART/trunk/matlab/PlotPhaseSpace.m
===================================================================
--- DART/trunk/matlab/PlotPhaseSpace.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/PlotPhaseSpace.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -163,7 +163,7 @@
end
legend boxoff
- case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf','cam'}
ens_mem_id = get_copy_index(pinfo.fname, pinfo.ens_mem); % errors out if no ens_mem
Modified: DART/trunk/matlab/PlotTotalErr.m
===================================================================
--- DART/trunk/matlab/PlotTotalErr.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/PlotTotalErr.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -261,7 +261,7 @@
otherwise
- disp(sprintf('unknown model %s -- doing nothing',model))
+ disp(sprintf('unsupported model %s -- doing nothing',model))
end
Modified: DART/trunk/matlab/PlotVarVarCorrel.m
===================================================================
--- DART/trunk/matlab/PlotVarVarCorrel.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/PlotVarVarCorrel.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -36,6 +36,7 @@
if (exist(pinfo.fname,'file') ~= 2), error('%s does not exist.',pinfo.fname), end
% Get some file-specific information.
+
model = nc_attget(pinfo.fname, nc_global, 'model');
timeunits = nc_attget(pinfo.fname, 'time', 'units');
@@ -43,7 +44,7 @@
switch lower(model)
- case {'fms_bgrid','pe2lyr','wrf'}
+ case {'fms_bgrid','pe2lyr','wrf','cam'}
clf;
@@ -172,7 +173,7 @@
% is the hard part.
% find which are actual ensemble members
-metadata = nc_varget(fname,'CopyMetaData'); % get all the metadata
+metadata = nc_varget(fname,'CopyMetaData'); % get all the metadata
copyindices = strmatch('ensemble member',metadata); % find all 'member's
if ( isempty(copyindices) )
@@ -185,6 +186,7 @@
end
ens_num = length(copyindices);
+% Get all ensemble members, just return desired ones.
myinfo.diagn_file = fname;
myinfo.levelindex = lvlind;
myinfo.latindex = latind;
Modified: DART/trunk/matlab/plot_bins.m
===================================================================
--- DART/trunk/matlab/plot_bins.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/plot_bins.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -68,7 +68,7 @@
case 'cam'
pinfo = CombineStructs(pinfo,vars);
- pinfo = GetCamInfo(pinfo, 'PlotBins');
+ pinfo = GetCamInfo(pinfo, diagn_file, 'PlotBins');
case 'pe2lyr'
Modified: DART/trunk/matlab/plot_ens_err_spread.m
===================================================================
--- DART/trunk/matlab/plot_ens_err_spread.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/plot_ens_err_spread.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -69,7 +69,7 @@
case 'cam'
pinfo = CombineStructs(pinfo,vars);
- pinfo = GetCamInfo(pinfo, 'PlotEnsErrSpread');
+ pinfo = GetCamInfo(pinfo, truth_file, 'PlotEnsErrSpread');
case 'pe2lyr'
Modified: DART/trunk/matlab/plot_ens_mean_time_series.m
===================================================================
--- DART/trunk/matlab/plot_ens_mean_time_series.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/plot_ens_mean_time_series.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -82,9 +82,7 @@
case 'cam'
- pinfo.prior_file = [];
- pinfo.posterior_file = [];
- pinfo = GetCamInfo(pinfo, 'PlotEnsMeanTimeSeries');
+ pinfo = GetCamInfo(pinfo, diagn_file, 'PlotEnsMeanTimeSeries');
case 'wrf'
Modified: DART/trunk/matlab/plot_ens_time_series.m
===================================================================
--- DART/trunk/matlab/plot_ens_time_series.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/plot_ens_time_series.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -50,7 +50,7 @@
for ifield = 1:length(mynames)
myname = mynames{ifield};
- if ( isfield(pinfo,myname) ), warning('%s already exists in pinfo\n'); end
+ if ( isfield(pinfo,myname) ), warning('plot_ens_time_series: pinfo.%s already exists\n',myname); end
eval(sprintf('pinfo.%s = MyInfo.%s;',myname,myname));
end
else
@@ -84,7 +84,7 @@
pinfo.prior_file = [];
pinfo.posterior_file = [];
- pinfo = GetCamInfo(pinfo, 'PlotEnsTimeSeries');
+ pinfo = GetCamInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
case 'wrf'
Modified: DART/trunk/matlab/plot_sawtooth.m
===================================================================
--- DART/trunk/matlab/plot_sawtooth.m 2011-07-08 20:41:18 UTC (rev 5065)
+++ DART/trunk/matlab/plot_sawtooth.m 2011-07-08 22:39:38 UTC (rev 5066)
@@ -86,7 +86,7 @@
case 'cam'
- pstruct = GetCamInfo(pstruct,'PlotSawtooth');
+ pstruct = GetCamInfo(pstruct, prior_file, 'PlotSawtooth');
pstruct.copyindices = SetCopyID2(pstruct.prior_file);
pstruct.copies = length(pstruct.copyindices);
More information about the Dart-dev
mailing list