[Dart-dev] [4830] DART/trunk/matlab: Support for WRF, streamline some for the other models.
nancy at ucar.edu
nancy at ucar.edu
Mon Mar 28 16:31:53 MDT 2011
Revision: 4830
Author: thoar
Date: 2011-03-28 16:31:52 -0600 (Mon, 28 Mar 2011)
Log Message:
-----------
Support for WRF, streamline some for the other models.
GetWRFInfo.m still needs to remove hardwired d01 bits.
Modified Paths:
--------------
DART/trunk/matlab/CheckModel.m
DART/trunk/matlab/GetWRFInfo.m
DART/trunk/matlab/PlotEnsMeanTimeSeries.m
DART/trunk/matlab/PlotEnsTimeSeries.m
DART/trunk/matlab/plot_ens_mean_time_series.m
DART/trunk/matlab/plot_ens_time_series.m
Added Paths:
-----------
DART/trunk/matlab/get_DARTvars.m
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModel.m
===================================================================
--- DART/trunk/matlab/CheckModel.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/CheckModel.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -270,17 +270,14 @@
case 'wrf'
- % 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)
- % get_varsNdims() ALMOST works ... except for WRF.
+ % requires a 'domain' and 'bottom_top_d01' dimension.
+ % without both of these, it will fail in an ugly fashion.
- varnames = {'U','V','W','PH','MU','QVAPOR','QCLOUD'};
+ varnames = get_DARTvars(fname);
num_vars = length(varnames);
- dinfo = nc_getdiminfo(fname,'domain'); % no graceful error
+ dinfo = nc_getdiminfo(fname,'domain');
num_domains = dinfo.Length;
- dinfo = nc_getdiminfo(fname,'bottom_top_d01'); % no graceful error
+ dinfo = nc_getdiminfo(fname,'bottom_top_d01');
num_levels = dinfo.Length;
vars = struct('model',model, ...
Modified: DART/trunk/matlab/GetWRFInfo.m
===================================================================
--- DART/trunk/matlab/GetWRFInfo.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/GetWRFInfo.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -57,7 +57,7 @@
dID = 1;
if (num_domains > 1)
- dID = GetDomain(num_domains); % Determine prognostic variable
+ dID = GetDomain(num_domains);
end
dn = varget(fname, sprintf( 'DN_d%02d',dID));
@@ -76,21 +76,12 @@
phb = varget(fname, sprintf( 'PHB_d%02d',dID));
hgt = varget(fname, sprintf( 'HGT_d%02d',dID));
-prognostic_bases = {'U', 'V', 'W', 'PH', 'T', ...
- 'MU', 'QVAPOR', 'QCLOUD', 'QRAIN', 'QICE'};
-
-for i = 1:length(prognostic_bases)
- prognostic_vars{i} = sprintf('%s_d%02d',prognostic_bases{i},dID);
-end
-
-varexist(fname, prognostic_vars)
-
switch lower(deblank(routine))
case {'plotbins','plotenserrspread','plotensmeantimeseries','plotenstimeseries'}
- pgvar = GetVarString(prognostic_vars); % Determine prognostic variable
- [level, lvlind] = GetLevel(fname,pgvar); % Determine level and index
+ pgvar = GetVarString(pinfo_in.vars); % Determine prognostic variable
+ [level, lvlind] = GetLevel(fname,pgvar); % Determine level and index
[lat, lon, latind, lonind] = GetLatLon(fname, pgvar);
pinfo.model = model;
@@ -106,14 +97,14 @@
case 'plotcorrel'
disp('Getting information for the ''base'' variable.')
- base_var = GetVarString(prognostic_vars);
+ base_var = GetVarString(pinfo_in.vars);
[base_time, base_tmeind] = GetTime(times);
[base_lvl, base_lvlind] = GetLevel(fname,base_var);
[base_lat, base_latind] = GetLatitude( base_var,TmpJ,VelJ);
[base_lon, base_lonind] = GetLongitude(base_var,TmpI,VelI);
disp('Getting information for the ''comparison'' variable.')
- comp_var = GetVarString(prognostic_vars, base_var);
+ comp_var = GetVarString(pinfo_in.vars, base_var);
[comp_lvl, comp_lvlind] = GetLevel(fname,comp_var,base_lvl);
pinfo.model = model;
@@ -134,14 +125,14 @@
case 'plotvarvarcorrel'
disp('Getting information for the ''base'' variable.')
- base_var = GetVarString(prognostic_vars);
+ base_var = GetVarString(pinfo_in.vars);
[base_time, base_tmeind] = GetTime(times);
[base_lvl , base_lvlind] = GetLevel(fname,base_var);
[base_lat , base_latind] = GetLatitude( base_var,TmpJ,VelJ);
[base_lon , base_lonind] = GetLongitude(base_var,TmpI,VelI);
disp('Getting information for the ''comparison'' variable.')
- comp_var = GetVarString(prognostic_vars, base_var);
+ comp_var = GetVarString(pinfo_in.vars, base_var);
[comp_lvl, comp_lvlind] = GetLevel(fname,comp_var,base_lvl);
[comp_lat, comp_latind] = GetLatitude( comp_var,TmpJ,VelJ, base_lat);
[comp_lon, comp_lonind] = GetLongitude(comp_var,TmpI,VelI, base_lon);
@@ -167,7 +158,7 @@
case 'plotsawtooth'
- pgvar = GetVarString(prognostic_vars);
+ 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);
@@ -192,19 +183,19 @@
case 'plotphasespace'
disp('Getting information for the ''X'' variable.')
- var1 = GetVarString(prognostic_vars);
+ 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);
disp('Getting information for the ''Y'' variable.')
- var2 = GetVarString(prognostic_vars, var1 );
+ 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);
disp('Getting information for the ''Z'' variable.')
- var3 = GetVarString(prognostic_vars, var1 );
+ 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);
@@ -376,7 +367,6 @@
end
-
% Determine a sensible default.
if (nargin > 2)
latind = deflat;
@@ -398,10 +388,10 @@
if ~isempty(varstring)
nums = str2num(varstring);
if (length(nums) ~= 2)
- error('Did not get two numbers for the lat lon pair.')
+ error('Did not get two indices for the lat lon pair.')
end
- lat = nums(1);
- lon = nums(2);
+ latind = nums(1);
+ lonind = nums(2);
end
% Find the closest lat/lon location to the requested one.
@@ -418,11 +408,11 @@
latinds = 1:nlat;
loninds = 1:nlon;
-d = abs(lat - latinds); % crude distance
+d = abs(latind - latinds); % crude distance
ind = find(min(d) == d); % multiple minima possible
latind = ind(1); % use the first one
-d = abs(lon - loninds); % crude distance
+d = abs(lonind - loninds); % crude distance
ind = find(min(d) == d); % multiple minima possible
lonind = ind(1); % use the first one
Modified: DART/trunk/matlab/PlotEnsMeanTimeSeries.m
===================================================================
--- DART/trunk/matlab/PlotEnsMeanTimeSeries.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/PlotEnsMeanTimeSeries.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -305,13 +305,12 @@
plot(pinfo.longitude,pinfo.latitude,'pb','MarkerSize',12,'MarkerFaceColor','b');
axlims = axis;
axlims = axlims + [-20 20 -20 20];
+ grid on
+ axis image
axis(axlims)
if (axlims(2) < 0)
worldmap('hollow','dateline');
else
worldmap('hollow','greenwich');
end
- axis image
- grid on
-
Modified: DART/trunk/matlab/PlotEnsTimeSeries.m
===================================================================
--- DART/trunk/matlab/PlotEnsTimeSeries.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/PlotEnsTimeSeries.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -218,7 +218,7 @@
legend boxoff
end
- case {'fms_bgrid','pe2lyr'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
clf;
@@ -327,8 +327,6 @@
% Subfunctions
%======================================================================
-
-
function x = dim_length(fname,dimname)
y = nc_isvar(fname,dimname);
if (y < 1)
@@ -354,9 +352,6 @@
-
-
-
function var = GetEns(fname, pinfo)
% Gets a time-series of all copies of a prognostic variable
% at a particular 3D location (level, lat, lon).
@@ -372,7 +367,7 @@
disp('To be a valid ensemble member, the CopyMetaData for the member')
disp('must start with the character string ''ensemble member''')
disp('None of them in do in your file.')
- fprintf('%s claims to have %d copies\n',fname, num_copies)
+ fprintf('%s claims to have 0 copies\n',fname)
error('netcdf file has no ensemble members.')
end
ens_num = length(copyindices);
@@ -390,10 +385,14 @@
function PlotLocator(pinfo)
plot(pinfo.longitude,pinfo.latitude,'pb','MarkerSize',12,'MarkerFaceColor','b');
- axis([0 360 -90 90])
- worldmap;
- axis image
+ axlims = axis;
+ axlims = axlims + [-20 20 -20 20];
grid on
+ axis image
+ axis(axlims)
+ if (axlims(2) < 0)
+ worldmap('hollow','dateline');
+ else
+ worldmap('hollow','greenwich');
+ end
-
-
Added: DART/trunk/matlab/get_DARTvars.m
===================================================================
--- DART/trunk/matlab/get_DARTvars.m (rev 0)
+++ DART/trunk/matlab/get_DARTvars.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -0,0 +1,56 @@
+function bob = get_DARTvars(fname)
+%% get_DARTvars returns just the variable names in the netCDF file that are likely to
+%be DART variables.
+%
+% the result is a cell array of strings ... must use {} notation to address elements.
+%
+% EXAMPLE:
+% fname = 'obs_seq.final.nc';
+% DARTvars = get_DARTvars(fname);
+% DARTvars{:}
+% nvars = length(DARTvars);
+% disp(sprintf('first atmospheric variable (of %d) is %s',nvars,DARTvars{1}))
+
+%% 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+fileinfo = nc_info(fname);
+nvars = length(fileinfo.Dataset);
+isDARTvar = zeros(nvars,1);
+
+for i = 1:nvars
+
+ % Anything with a 'copy' dimension is probably a DART state vector variable.
+
+ dimnames = fileinfo.Dataset(i).Dimension;
+
+ if (any(strcmp(dimnames,'copy'))), isDARTvar(i) = 1; end
+
+ % Reject the obvious coordinate variables and some metadata ones
+
+ varname = fileinfo.Dataset(i).Name;
+ if ( nc_iscoordvar(fname,varname)), isDARTvar(i) = 0; end
+ if (strcmp( varname , 'CopyMetaData')), isDARTvar(i) = 0; end
+
+end
+
+if (sum(isDARTvar) == 0)
+ error('No DART state variables in %s',fname)
+end
+
+% coerce just the names into a cell array
+
+varind = 0;
+for i = 1:nvars
+ if (isDARTvar(i) > 0)
+ varind = varind + 1;
+ bob{varind} = fileinfo.Dataset(i).Name;
+ end
+end
Property changes on: DART/trunk/matlab/get_DARTvars.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/matlab/plot_ens_mean_time_series.m
===================================================================
--- DART/trunk/matlab/plot_ens_mean_time_series.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/plot_ens_mean_time_series.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -39,26 +39,38 @@
end
end
-vars = CheckModel(diagn_file); % also gets default values for this model.
+pinfo = CheckModel(diagn_file); % also gets default values for this model.
if (exist(truth_file,'file')==2)
- pinfo = CheckModelCompatibility(truth_file, diagn_file);
+
+ MyInfo = CheckModelCompatibility(truth_file, diagn_file);
+
+ % Combine the information from CheckModel and CheckModelCompatibility
+ mynames = fieldnames(MyInfo);
+
+ for ifield = 1:length(mynames)
+ myname = mynames{ifield};
+ if ( isfield(pinfo,myname) ), warning('%s already exists in pinfo\n'); end
+ eval(sprintf('pinfo.%s = MyInfo.%s;',myname,myname));
+ end
else
- pinfo.truth_file = [];
+ truth_file = [];
end
-truth_file = pinfo.truth_file;
+pinfo.diagn_file = diagn_file;
+pinfo.truth_file = truth_file;
+clear MyInfo mynames myname ifield
-switch lower(vars.model)
+%% For each model, do what needs to be done.
+switch lower(pinfo.model)
+
case {'9var','lorenz_63','lorenz_84','lorenz_96','lorenz_96_2scale', ...
- 'lorenz_04','forced_lorenz_96','ikeda','simple_advection'}
+ 'forced_lorenz_96','lorenz_04','ikeda','simple_advection'}
- varid = SetVariableID(vars); % queries for variable IDs if needed.
- pinfo = setfield(pinfo, 'truth_file', truth_file);
- pinfo = setfield(pinfo, 'diagn_file', diagn_file);
- pinfo = setfield(pinfo, 'var' , varid.var);
- pinfo = setfield(pinfo, 'var_inds' , varid.var_inds);
+ varid = SetVariableID(pinfo); % queries for variable IDs if needed.
+ pinfo.var = varid.var;
+ pinfo.var_inds = varid.var_inds;
fprintf('Comparing %s and \n %s\n', pinfo.truth_file, pinfo.diagn_file)
disp(['Using State Variable IDs ', num2str(pinfo.var_inds)])
@@ -67,46 +79,32 @@
case 'fms_bgrid'
pinfo = GetBgridInfo(pinfo, diagn_file, 'PlotEnsMeanTimeSeries');
- pinfo.truth_file = truth_file; % since it has been verified to be compatible.
- pinfo.diagn_file = diagn_file; % since it has been verified to be compatible.
case 'cam'
- vars.truth_file = truth_file;
- vars.diagn_file = diagn_file;
- vars.prior_file = [];
- vars.posterior_file = [];
- pinfo = GetCamInfo(vars, 'PlotEnsMeanTimeSeries');
- % pinfo.copyindices = SetCopyID2(vars.diagn_file);
- % pinfo.copies = length(pinfo.copyindices);
- pinfo.truth_file = truth_file;
- pinfo.diagn_file = diagn_file;
+ pinfo.prior_file = [];
+ pinfo.posterior_file = [];
+ pinfo = GetCamInfo(pinfo, 'PlotEnsMeanTimeSeries');
case 'wrf'
pinfo = GetWRFInfo(pinfo, diagn_file, 'PlotEnsMeanTimeSeries');
- pinfo.truth_file = truth_file;
- pinfo.diagn_file = diagn_file;
case 'pe2lyr'
pinfo = GetPe2lyrInfo(pinfo, diagn_file, 'PlotEnsMeanTimeSeries');
- pinfo.truth_file = truth_file;
- pinfo.diagn_file = diagn_file;
case 'mitgcm_ocean'
pinfo = GetMITgcm_oceanInfo(pinfo, diagn_file, 'PlotEnsMeanTimeSeries');
- pinfo.truth_file = truth_file;
- pinfo.diagn_file = diagn_file;
otherwise
- error('model %s not implemented yet', vars.model)
+ error('model %s not implemented yet', pinfo.model)
end
pinfo
PlotEnsMeanTimeSeries( pinfo )
-clear vars
+
Modified: DART/trunk/matlab/plot_ens_time_series.m
===================================================================
--- DART/trunk/matlab/plot_ens_time_series.m 2011-03-28 22:19:42 UTC (rev 4829)
+++ DART/trunk/matlab/plot_ens_time_series.m 2011-03-28 22:31:52 UTC (rev 4830)
@@ -1,4 +1,9 @@
-%% DART:plot_ens_time_series - time series of ensemble members, mean and truth
+%% DART:plot_ens_time_series - time series of ensemble and truth (if available)
+%
+% plot_ens_time_series interactively queries for the needed information.
+% Since different models potentially need different pieces of
+% information ... the model types are determined and additional
+% user input may be queried.
%
% Example 2
% truth_file = 'True_State.nc';
@@ -15,15 +20,8 @@
% $Revision$
% $Date$
-if (exist('truth_file','var') ~= 1)
- disp('If the True_State.nc exists, it will be plotted. If not, don''t worry.')
- truth_file = input('Input name of True State file; <cr> for True_State.nc\n','s');
- if isempty(truth_file)
- truth_file = 'True_State.nc';
- end
-end
-
if (exist('diagn_file','var') ~=1)
+ disp(' ')
disp('Input name of prior or posterior diagnostics file;')
diagn_file = input('<cr> for Prior_Diag.nc\n','s');
if isempty(diagn_file)
@@ -31,68 +29,79 @@
end
end
-vars = CheckModel(diagn_file); % also gets default values for this model.
+if (exist('truth_file','var') ~= 1)
+ disp(' ')
+ disp('OPTIONAL: if you have the true state and want it superimposed, provide')
+ disp(' : the name of the input file. If not, enter a dummy filename.')
+ truth_file = input('Input name of True State file; <cr> for True_State.nc\n','s');
+ if isempty(truth_file)
+ truth_file = 'True_State.nc';
+ end
+end
+pinfo = CheckModel(diagn_file); % also gets default values for this model.
+
if (exist(truth_file,'file')==2)
- pinfo = CheckModelCompatibility(truth_file, diagn_file);
+
+ MyInfo = CheckModelCompatibility(truth_file, diagn_file);
+
+ % Combine the information from CheckModel and CheckModelCompatibility
+ mynames = fieldnames(MyInfo);
+
+ for ifield = 1:length(mynames)
+ myname = mynames{ifield};
+ if ( isfield(pinfo,myname) ), warning('%s already exists in pinfo\n'); end
+ eval(sprintf('pinfo.%s = MyInfo.%s;',myname,myname));
+ end
else
- pinfo.truth_file = [];
+ truth_file = [];
end
-truth_file = pinfo.truth_file;
-pinfo.model = vars.model;
+pinfo.diagn_file = diagn_file;
+pinfo.truth_file = truth_file;
-switch lower(vars.model)
+clear MyInfo mynames myname ifield
+%% For each model, do what needs to be done.
+
+switch lower(pinfo.model)
+
case {'9var','lorenz_63','lorenz_84','lorenz_96','lorenz_96_2scale', ...
'forced_lorenz_96','lorenz_04','ikeda','simple_advection'}
- varid = SetVariableID(vars); % queries for variable IDs if needed.
- pinfo = setfield(pinfo,'truth_file', truth_file);
- pinfo = setfield(pinfo,'diagn_file', diagn_file);
- pinfo = setfield(pinfo,'var' , varid.var);
- pinfo = setfield(pinfo,'var_inds' , varid.var_inds);
- %pinfo = struct('truth_file', truth_file, ...
- % 'diagn_file', diagn_file, ...
- % 'var' , varid.var, ...
- % 'var_inds' , varid.var_inds);
+ varid = SetVariableID(pinfo); % queries for variable IDs if needed.
+ pinfo.var = varid.var;
+ pinfo.var_inds = varid.var_inds;
- % disp(sprintf('Comparing %s and \n %s', pinfo.truth_file, pinfo.diagn_file))
- % disp(sprintf('Using Variable %s IDs %s', pinfo.var,num2str(pinfo.var_inds)))
+ fprintf('Comparing %s and \n %s\n', pinfo.truth_file, pinfo.diagn_file)
+ disp(['Using State Variable IDs ', num2str(pinfo.var_inds)])
+ clear varid
case 'fms_bgrid'
- varid = SetVariableID(vars); % queries for variable IDs if needed.
pinfo = GetBgridInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
- pinfo.truth_file = truth_file; % known to be compatible.
- pinfo.diagn_file = diagn_file; % known to be compatible.
case 'cam'
- vars.truth_file = truth_file;
- vars.diagn_file = diagn_file;
- vars.prior_file = [];
- vars.posterior_file = [];
- pinfo = GetCamInfo(vars, 'PlotEnsTimeSeries');
- pinfo.truth_file = truth_file;
- pinfo.diagn_file = diagn_file;
+ pinfo.prior_file = [];
+ pinfo.posterior_file = [];
+ pinfo = GetCamInfo(pinfo, 'PlotEnsTimeSeries');
+ case 'wrf'
+
+ pinfo = GetWRFInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
+
case 'pe2lyr'
- varid = SetVariableID(vars); % queries for variable IDs if needed.
pinfo = GetPe2lyrInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
- pinfo.truth_file = truth_file; % known to be compatible.
- pinfo.diagn_file = diagn_file; % known to be compatible.
case 'mitgcm_ocean'
- varid = SetVariableID(vars); % queries for variable IDs if needed.
- pinfo = GetPe2lyrInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
- pinfo.truth_file = truth_file; % known to be compatible.
- pinfo.diagn_file = diagn_file; % known to be compatible.
+ pinfo = GetMITgcm_oceanInfo(pinfo, diagn_file, 'PlotEnsTimeSeries');
otherwise
- error('model %s not implemented yet', vars.model)
+ error('model %s not implemented yet', pinfo.model)
+
end
pinfo % echo for posterity.
More information about the Dart-dev
mailing list