[Dart-dev] [6125] DART/branches/development: These routines support the SQG model and use (only) the snctools toolbox.
nancy at ucar.edu
nancy at ucar.edu
Fri May 10 09:30:41 MDT 2013
Revision: 6125
Author: thoar
Date: 2013-05-10 09:30:38 -0600 (Fri, 10 May 2013)
Log Message:
-----------
These routines support the SQG model and use (only) the snctools toolbox.
I have started to replace the calls to deprecated Matlab functions.
The startup.m function is much more graceful and readable.
It will not error out if you are not using DART, for example.
There is a new check in the plot_profile.m function to detect
duplicate 'level' items - which cause a logical conundrum.
The old version of obs_diag.f90 did not check for this - problematic.
Modified Paths:
--------------
DART/branches/development/diagnostics/matlab/get_varnames.m
DART/branches/development/diagnostics/matlab/obs_num_time.m
DART/branches/development/diagnostics/matlab/plot_profile.m
DART/branches/development/matlab/CheckModel.m
DART/branches/development/matlab/CheckModelCompatibility.m
DART/branches/development/matlab/GetCamInfo.m
DART/branches/development/matlab/GetNCindices.m
DART/branches/development/matlab/PlotBins.m
DART/branches/development/matlab/PlotCorrel.m
DART/branches/development/matlab/PlotEnsErrSpread.m
DART/branches/development/matlab/PlotEnsMeanTimeSeries.m
DART/branches/development/matlab/PlotEnsTimeSeries.m
DART/branches/development/matlab/PlotJeffCorrel.m
DART/branches/development/matlab/PlotTotalErr.m
DART/branches/development/matlab/PlotVarVarCorrel.m
DART/branches/development/matlab/plot_bins.m
DART/branches/development/matlab/plot_correl.m
DART/branches/development/matlab/plot_ens_err_spread.m
DART/branches/development/matlab/plot_ens_mean_time_series.m
DART/branches/development/matlab/plot_ens_time_series.m
DART/branches/development/matlab/plot_jeff_correl.m
DART/branches/development/matlab/plot_phase_space.m
DART/branches/development/matlab/plot_sawtooth.m
DART/branches/development/matlab/plot_total_err.m
DART/branches/development/matlab/plot_var_var_correl.m
DART/branches/development/matlab/startup.m
Added Paths:
-----------
DART/branches/development/matlab/GetSqgInfo.m
-------------- next part --------------
Modified: DART/branches/development/diagnostics/matlab/get_varnames.m
===================================================================
--- DART/branches/development/diagnostics/matlab/get_varnames.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/diagnostics/matlab/get_varnames.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -46,7 +46,7 @@
end
if (isempty(inds))
- error('No atmospheric variables in %s',name(f))
+ error('No atmospheric variables in %s',fname)
end
% coerce just the names into a cell array
Modified: DART/branches/development/diagnostics/matlab/obs_num_time.m
===================================================================
--- DART/branches/development/diagnostics/matlab/obs_num_time.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/diagnostics/matlab/obs_num_time.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -45,7 +45,7 @@
% Get attributes from obs_diag run.
%----------------------------------------------------------------------
-if ( exist(datafile) == 2 )
+if ( exist(datafile,'file') == 2 )
eval(datafile)
@@ -63,7 +63,7 @@
end
else
- error(sprintf('%s cannot be found.', datafile))
+ error('%s cannot be found.', datafile)
end
% set up a structure with all static plotting components
@@ -128,7 +128,7 @@
figure(page2); clf;
h = plot( xax, NbyRegion(:,1), ptypes{1}, 'LineWidth', 2.0);
- [legh, objh, outh, outm] = legend(Regions{1});
+ [~, ~, outh, outm] = legend(Regions{1});
hold on;
for iregion = 2:length(Regions),
@@ -140,7 +140,7 @@
% outm cell array for the text in the legend
nlines = length(outm);
outm{nlines + 1} = Regions{iregion};
- [legh, objh, outh, outm] = legend([outh; h],outm,0);
+ [~, ~, outh, outm] = legend([outh; h],outm,0);
end
legend boxoff
@@ -243,7 +243,7 @@
subplot('position',[0.48 0.01 0.04 0.04])
axis off
bob = which(main);
-[pathstr,name,ext,versn] = fileparts(bob);
+[pathstr,~,~] = fileparts(bob);
h = text(0.0,0.5,pathstr);
set(h,'HorizontalAlignment','center', ...
'VerticalAlignment','middle',...
Modified: DART/branches/development/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/branches/development/diagnostics/matlab/plot_profile.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/diagnostics/matlab/plot_profile.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -8,7 +8,7 @@
%
% fname : netcdf file produced by 'obs_diag'
% copystring : 'copy' string == quantity of interest. These
-% can be any of the ones available in the netcdf
+% can be any of the ones available in the netcdf
% file 'CopyMetaData' variable.
% (ncdump -v CopyMetaData obs_diag_output.nc)
%
@@ -90,7 +90,7 @@
plotdat.nvars = length(plotdat.varnames);
-plotdat.copyindex = get_copy_index(fname,copystring);
+plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -99,16 +99,16 @@
%----------------------------------------------------------------------
for ivar = 1:plotdat.nvars
-
+
% create the variable names of interest.
-
+
plotdat.myvarname = plotdat.varnames{ivar};
plotdat.guessvar = sprintf('%s_VPguess',plotdat.varnames{ivar});
plotdat.analyvar = sprintf('%s_VPanaly',plotdat.varnames{ivar});
% remove any existing postscript file - will simply append each
% level as another 'page' in the .ps file.
-
+
psfname = sprintf('%s_%s_profile.ps',plotdat.varnames{ivar},plotdat.copystring);
fprintf('Removing %s from the current directory.\n',psfname)
system(sprintf('rm %s',psfname));
@@ -119,10 +119,11 @@
analydims = nc_var_dims( fname, plotdat.analyvar);
varinfo = nc_getvarinfo(fname, plotdat.analyvar);
- if ( findstr('surface',guessdims{2}) > 0 )
+ if (~ isempty(strfind(guessdims{2},'surface')))
fprintf('%s is a surface field.\n',plotdat.guessvar)
fprintf('Cannot display a surface field this way.\n')
- elseif ( findstr('undef',guessdims{2}) > 0 )
+ end
+ if (~ isempty(strfind(guessdims{2},'undef')))
fprintf('%s has no vertical definition.\n',plotdat.guessvar)
fprintf('Cannot display this field this way.\n')
end
@@ -144,15 +145,18 @@
plotdat.YDir = 'normal';
end
[levels, indices] = sort(plotdat.level_org);
- plotdat.level = levels;
+ plotdat.level = unique(levels);
+ if (length(plotdat.level) ~= length(levels))
+ error('There is a duplicated value in the array specifying the levels - must change your input.nml and rerun obs_diag')
+ end
plotdat.indices = indices;
level_edges = sort(plotdat.level_edges);
plotdat.level_edges = level_edges;
-
- guess = nc_varget(fname, plotdat.guessvar);
- analy = nc_varget(fname, plotdat.analyvar);
+
+ guess = nc_varget(fname, plotdat.guessvar);
+ analy = nc_varget(fname, plotdat.analyvar);
n = size(analy);
-
+
% singleton dimensions are auto-squeezed - which is unfortunate.
% We want these things to be 3D. [copy-level-region]
% Sometimes there is one region, sometimes one level, ...
@@ -175,7 +179,7 @@
guess = bob; clear bob
analy = ted; clear ted
end
-
+
% check to see if there is anything to plot
nposs = sum(guess(plotdat.Npossindex,:,:));
@@ -197,7 +201,7 @@
clf; orient tall
for iregion = 1:plotdat.nregions
- plotdat.region = iregion;
+ plotdat.region = iregion;
plotdat.myregion = deblank(plotdat.region_names(iregion,:));
myplot(plotdat);
@@ -238,8 +242,8 @@
% Determine some quantities for the legend
nobs = sum(nobs_used);
if ( nobs > 1 )
- other_guess = mean(CG(isfinite(CG)));
- other_analy = mean(CA(isfinite(CA)));
+ other_guess = mean(CG(isfinite(CG)));
+ other_analy = mean(CA(isfinite(CA)));
else
other_guess = NaN;
other_analy = NaN;
@@ -253,8 +257,8 @@
% Ultimately, we want to suppress the 'auto' feature of the
% axis labelling, so we manually set some values that normally
% don't need to be set.
-
- % if more then 4 regions, this will not work (well) ...
+
+ % if more then 4 regions, this will not work (well) ...
if ( plotdat.nregions > 2 )
ax1 = subplot(2,2,plotdat.region);
else
@@ -265,7 +269,8 @@
end
Stripes(plotdat.Xrange, plotdat.level_edges);
- set(ax1,'YDir', plotdat.YDir,'YTick',plotdat.level,'Layer','top')
+ set(ax1,'YTick',plotdat.level,'Layer','top')
+ set(ax1,'YDir', plotdat.YDir)
ylabel(plotdat.level_units)
%% draw the result of the experiment
@@ -301,7 +306,7 @@
h2 = line(nobs_poss,plotdat.level,'Color','b','Parent',ax2);
h3 = line(nobs_used,plotdat.level,'Color','b','Parent',ax2);
set(h2,'LineStyle','none','Marker','o');
- set(h3,'LineStyle','none','Marker','+');
+ set(h3,'LineStyle','none','Marker','+');
% use same Y ticks
set(ax2,'YTick', get(ax1,'YTick'), ...
@@ -354,6 +359,8 @@
function [y,ydims] = FindVerticalVars(x)
% Returns UNIQUE (i.e. base) vertical variable names
+% In this context, if the variable has a 'time' dimension
+% it cannot be a variable of interest.
if ( ~(isfield(x,'allvarnames') && isfield(x,'allvardims')))
error('Doh! no ''allvarnames'' and ''allvardims'' components')
end
@@ -363,8 +370,8 @@
basedims = struct([]);
for i = 1:length(x.allvarnames)
- indx = findstr('time',x.allvardims{i});
- if (isempty(indx))
+ dimnames = lower(x.allvardims{i});
+ if (isempty(strfind(dimnames,'time')))
j = j + 1;
basenames{j} = ReturnBase(x.allvarnames{i});
@@ -372,7 +379,7 @@
end
end
-[b,i,j] = unique(basenames);
+[~,i,j] = unique(basenames);
y = struct([]);
ydims = struct([]);
@@ -408,22 +415,23 @@
function s = ReturnBase(string1)
-inds = findstr('_guess',string1);
+s = [];
+inds = strfind(string1,'_guess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_analy',string1);
+inds = strfind(string1,'_analy');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPguess',string1);
+inds = strfind(string1,'_VPguess');
if (inds > 0 )
s = string1(1:inds-1);
end
-inds = findstr('_VPanaly',string1);
+inds = strfind(string1,'_VPanaly');
if (inds > 0 )
s = string1(1:inds-1);
end
@@ -436,7 +444,7 @@
%
% In this scope, y is bounded from below by 0.0
%
-% If the numbers are very small ...
+% If the numbers are very small ...
bob = [y.ges_copy(:); ...
y.anl_copy(:)];
@@ -449,7 +457,7 @@
ymin = min(glommed);
ymax = max(glommed);
- if ( ymax > 1.0 )
+ if ( ymax > 1.0 )
ymin = floor(min(glommed));
ymax = ceil(max(glommed));
end
@@ -457,7 +465,7 @@
if (ymin == 0 && ymax == 0)
ymax = 1;
end
-
+
if (ymin == ymax)
ymin = ymin - 0.1*ymin;
ymax = ymax + 0.1*ymax;
@@ -487,7 +495,7 @@
% FIXME:
% This really should be modified to add a percentage of the data
% range to provide space for the legend. Right now it is hardwired
-% to assume that we are plotting hPa, on a 'reverse' axis.
+% to assume that we are plotting hPa, on a 'reverse' axis.
hold on;
@@ -498,7 +506,7 @@
h = plot([min(x) max(x)],[min(edges) max(edges)]);
axlims = axis;
axlims(4) = max(edges);
-axlims(3) = -100; % This gives extra space for the legend.
+axlims(3) = axlims(3) - 0.1*(max(edges) - min(edges)); % extra space for the legend
axis(axlims)
xc = [ axlims(1) axlims(2) axlims(2) axlims(1) axlims(1) ];
@@ -517,7 +525,7 @@
%% This takes the existing X ticks from ax1 (presumed nice)
% and determines the matching labels for ax2 so we can keep
% at least one of the axes looking nice.
-
+
Dlimits = get(ax1,'XLim');
DXticks = get(ax1,'XTick');
nXticks = length(DXticks);
Modified: DART/branches/development/matlab/CheckModel.m
===================================================================
--- DART/branches/development/matlab/CheckModel.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/CheckModel.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -84,7 +84,7 @@
vars.fname = fname;
- case 'forced_lorenz_96'
+ case {'forced_lorenz_96'}
% This model has the state variables replicated, so there is a difference
% between num_state_vars and the length of the state variable.
@@ -124,7 +124,7 @@
vars.fname = fname;
- case 'lorenz_96_2scale'
+ case {'lorenz_96_2scale'}
num_X = dim_length(fname,'Xdim'); % # of X variables
Xdim = nc_varget(fname,'Xdim');
@@ -154,7 +154,7 @@
vars.fname = fname;
- case 'simple_advection'
+ case {'simple_advection'}
num_locs = dim_length(fname,'loc1d'); % # of X variables
loc1d = nc_varget(fname,'loc1d');
@@ -185,7 +185,7 @@
vars.vars = varnames;
vars.fname = fname;
- case 'wrf'
+ case {'wrf'}
% requires a 'domain' and 'bottom_top_d01' dimension.
% without both of these, it will fail in an ugly fashion.
@@ -221,13 +221,15 @@
'num_state_vars',num_vars, ...
'num_ens_members',num_copies, ...
'time_series_length',num_times, ...
- 'min_ens_mem',min(copy), ...
- 'max_ens_mem',max(copy) );
+ 'num_ens_members',ens_size, ...
+ 'ensemble_indices',ens_indices, ...
+ 'min_ens_mem',ens_indices(1), ...
+ 'max_ens_mem',ens_indices(ens_size) );
vars.vars = varnames;
vars.fname = fname;
- case {'cam','tiegcm','fms_bgrid','pe2lyr','mitgcm_ocean','pbl_1d','mpas_atm'}
+ case {'cam','tiegcm','fms_bgrid','pe2lyr','mitgcm_ocean','pbl_1d','mpas_atm','sqg'}
varnames = get_DARTvars(fname);
num_vars = length(varnames);
@@ -239,6 +241,8 @@
'num_copies',num_copies, ...
'num_ens_members',ens_size, ...
'ensemble_indices',ens_indices, ...
+ 'min_ens_mem',ens_indices(1), ...
+ 'max_ens_mem',ens_indices(ens_size), ...
'time',dates, ...
'time_series_length',num_times);
Modified: DART/branches/development/matlab/CheckModelCompatibility.m
===================================================================
--- DART/branches/development/matlab/CheckModelCompatibility.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/CheckModelCompatibility.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -252,55 +252,48 @@
switch lower(modelname)
- case 'wrf'
+ case {'wrf'}
dnum_lons = dim_length(fname, 'west_east_d01');
dnum_lats = dim_length(fname,'south_north_d01');
dnum_lvls = dim_length(fname, 'bottom_top_d01');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- case 'cam'
+ case {'pe2lyr','cam','sqg'}
dnum_lons = dim_length(fname,'lon');
dnum_lats = dim_length(fname,'lat');
dnum_lvls = dim_length(fname,'lev');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- case 'pe2lyr'
- dnum_lons = dim_length(fname,'lon');
- dnum_lats = dim_length(fname,'lat');
- dnum_lvls = dim_length(fname,'lev');
- x = 3;
- y = [dnum_lons dnum_lats dnum_lvls];
-
- case 'fms_bgrid'
+ case {'fms_bgrid'}
dnum_lons = dim_length(fname,'TmpI');
dnum_lats = dim_length(fname,'TmpJ');
dnum_lvls = dim_length(fname,'lev');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- case 'mitgcm_ocean'
+ case {'mitgcm_ocean'}
dnum_lons = dim_length(fname,'XG');
dnum_lats = dim_length(fname,'YG');
dnum_lvls = dim_length(fname,'ZG');
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
- case 'mpas_atm'
+ case {'mpas_atm'}
dnum_cells = dim_length(fname,'nCells');
dnum_lvls = dim_length(fname,'nVertLevels');
x = 2;
y = [dnum_cells dnum_lvls];
- case 'lorenz_96_2scale'
+ case {'lorenz_96_2scale'}
dnum_X = dim_length(fname,'Xdim');
dnum_Y = dim_length(fname,'Ydim');
x = 2;
y = [dnum_X dnum_Y];
- case 'simple_advection'
+ case {'simple_advection'}
y = dim_length(fname,'loc1d');
x = 1;
Modified: DART/branches/development/matlab/GetCamInfo.m
===================================================================
--- DART/branches/development/matlab/GetCamInfo.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/GetCamInfo.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -117,7 +117,8 @@
[level, lvlind] = GetLevel(pgvar,levels);
[ lat, latind] = GetLatitude(lat);
[ lon, lonind] = GetLongitude(lon);
- % [ lon, lonind] = GetCopies(pgvar,xxx);
+ copyindices = SetCopyID2(fname);
+ copy = length(copyindices);
pinfo.var_names = pgvar;
pinfo.truth_file = [];
@@ -129,8 +130,8 @@
pinfo.latindex = latind;
pinfo.longitude = lon;
pinfo.lonindex = lonind;
- pinfo.copies = 0;
- pinfo.copyindices = [];
+ pinfo.copies = copy;
+ pinfo.copyindices = copyindices;
if ( exist(pstruct.truth_file,'file') )
pinfo.truth_file = pstruct.truth_file;
Modified: DART/branches/development/matlab/GetNCindices.m
===================================================================
--- DART/branches/development/matlab/GetNCindices.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/GetNCindices.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -16,6 +16,9 @@
% pinfo.lonindex
% pinfo.stateindex
% pinfo.regionindex
+% pinfo.cellindex
+% pinfo.columnindex
+% pinfo.pftindex
%
% whichfile is a character string specifying which
% filename component of 'pinfo' will be used.
@@ -68,6 +71,8 @@
state1 = 0; stateN = -1;
region1 = 0; regionN = -1;
cell1 = 0; cellN = -1;
+column1 = 0; columnN = -1;
+pft1 = 0; pftN = -1;
if (isfield(pinfo,'timeindex'))
time1 = pinfo.timeindex - 1;
@@ -79,18 +84,40 @@
if (isfield(pinfo,'tcount'))
timeN = pinfo.tcount;
end
+
if (isfield(pinfo,'levelindex'))
level1 = pinfo.levelindex - 1;
levelN = 1;
end
+if (isfield(pinfo,'level1'))
+ level1 = pinfo.level1 - 1;
+end
+if (isfield(pinfo,'levelcount'))
+ levelN = pinfo.levelcount;
+end
+
if (isfield(pinfo,'latindex'))
lat1 = pinfo.latindex - 1;
latN = 1;
end
+if (isfield(pinfo,'lat1'))
+ lat1 = pinfo.lat1 - 1;
+end
+if (isfield(pinfo,'latcount'))
+ latN = pinfo.latcount;
+end
+
if (isfield(pinfo,'lonindex'))
lon1 = pinfo.lonindex - 1;
lonN = 1;
end
+if (isfield(pinfo,'lon1'))
+ lon1 = pinfo.lon1 - 1;
+end
+if (isfield(pinfo,'loncount'))
+ lonN = pinfo.loncount;
+end
+
if (isfield(pinfo,'stateindex'))
state1 = pinfo.stateindex - 1;
stateN = 1;
@@ -101,25 +128,62 @@
if (isfield(pinfo,'statecount'))
stateN = pinfo.statecount;
end
+
if (isfield(pinfo,'copyindex'))
copy1 = pinfo.copyindex - 1;
copyN = 1;
end
-if (isfield(pinfo,'copyindex1'))
- copy1 = pinfo.copyindex1 - 1;
+if (isfield(pinfo,'copy1'))
+ copy1 = pinfo.copy1 - 1;
end
if (isfield(pinfo,'copycount'))
copyN = pinfo.copycount;
end
+
if (isfield(pinfo,'regionindex'))
region1 = pinfo.regionindex - 1;
regionN = 1;
end
+if (isfield(pinfo,'region1'))
+ region1 = pinfo.region1 - 1;
+end
+if (isfield(pinfo,'regioncount'))
+ regionN = pinfo.regioncount;
+end
+
if (isfield(pinfo,'cellindex'))
cell1 = pinfo.cellindex - 1;
cellN = 1;
end
+if (isfield(pinfo,'cell1'))
+ cell1 = pinfo.cell1 - 1;
+end
+if (isfield(pinfo,'cellcount'))
+ cellN = pinfo.cellcount;
+end
+if (isfield(pinfo,'columnindex'))
+ column1 = pinfo.columnindex - 1;
+ columnN = 1;
+end
+if (isfield(pinfo,'column1'))
+ column1 = pinfo.column1 - 1;
+end
+if (isfield(pinfo,'columncount'))
+ columnN = pinfo.columncount;
+end
+
+if (isfield(pinfo,'pftindex'))
+ pft1 = pinfo.pftindex - 1;
+ pftN = 1;
+end
+if (isfield(pinfo,'pft1'))
+ pft1 = pinfo.pft1 - 1;
+end
+if (isfield(pinfo,'pftcount'))
+ pftN = pinfo.pftcount;
+end
+
% Determine shape of variable in question.
varinfo = nc_getvarinfo(fname,varname);
@@ -147,7 +211,7 @@
% So the XG coordinate dimension has 'cartesian_axis = X',
% for example.
- [len, status, value] = is_dimension_cartesian(fname, diminfo.Name);
+ [~, status, value] = is_dimension_cartesian(fname, diminfo.Name);
if (status > 0)
dimname = value;
@@ -164,7 +228,8 @@
case 'copy'
start(i) = copy1;
count(i) = copyN;
- case {'surf','unde','hlev','mlev','plev','heig','leve','bott','ilev','nver','levt','levs'}
+ case {'surf','unde','hlev','mlev','plev','heig','leve','bott', ...
+ 'ilev','nver','levt','levs'}
start(i) = level1;
count(i) = levelN;
case {'tmpj','sout'}
@@ -184,6 +249,9 @@
case 'ncel'
start(i) = cell1;
count(i) = cellN;
+ case 'colu'
+ start(i) = column1;
+ count(i) = columnN;
otherwise
fprintf('GetNCindices encountered unknown coordinate variable %s\n',dimname)
end
@@ -201,6 +269,9 @@
case {'lon','x'}
start(i) = lon1;
count(i) = lonN;
+ case 'pft'
+ start(i) = pft1;
+ count(i) = pftN;
otherwise
fprintf('GetNCindices encountered unknown coordinate variable %s\n',dimname)
end
Added: DART/branches/development/matlab/GetSqgInfo.m
===================================================================
--- DART/branches/development/matlab/GetSqgInfo.m (rev 0)
+++ DART/branches/development/matlab/GetSqgInfo.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -0,0 +1,353 @@
+function pinfo = GetSqgInfo(pstruct,fname,routine)
+%% GetSqgInfo prepares a structure of information needed by the subsequent "routine"
+% The information is gathered via rudimentary "input" routines.
+%
+% pinfo = GetSqgInfo(pinfo_in,fname,routine);
+%
+% pinfo_in Name of existing pinfo struct, e.g. output from CheckModelCompatibility
+% fname Name of the DART netcdf file
+% routine name of subsequent plot routine.
+
+%% 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
+% http://www.image.ucar.edu/DAReS/DART/DART_download
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
+
+pinfo = pstruct;
+model = nc_attget(fname,nc_global,'model');
+
+if strcmpi(model,'sqg') ~= 1
+ error('Not so fast, this is not a sqg model.')
+end
+
+copy = nc_varget(fname,'copy');
+times = nc_varget(fname,'time');
+levels = nc_varget(fname,'lev');
+lons = nc_varget(fname,'lon');
+lats = nc_varget(fname,'lat');
+
+switch lower(deblank(routine))
+
+ case {'plotbins','plotenserrspread','plotensmeantimeseries','plotenstimeseries'}
+
+ pgvar = GetVar(pinfo.vars);
+ [level, lvlind] = GetLevel(levels);
+ [lat , latind] = GetLatitude(lats);
+ [lon , lonind] = GetLongitude(lons);
+
+ pinfo.var = pgvar;
+ pinfo.level = level;
+ pinfo.levelindex = lvlind;
+ pinfo.longitude = lon;
+ pinfo.lonindex = lonind;
+ pinfo.latitude = lat;
+ pinfo.latindex = latind;
+
+
+ case 'plotcorrel'
+
+ disp('Getting information for the ''base'' variable.')
+ base_var = GetVar(pinfo.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
+ [base_lvl, base_lvlind] = GetLevel(levels);
+ [base_lat, base_latind] = GetLatitude(lats);
+ [base_lon, base_lonind] = GetLongitude(lons);
+
+ disp('Getting information for the ''comparison'' variable.')
+ comp_var = GetVar(pinfo.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel(levels, base_lvlind);
+
+ pinfo.base_var = base_var;
+ pinfo.comp_var = comp_var;
+ pinfo.base_time = base_time;
+ pinfo.base_tmeind = base_tmeind;
+ pinfo.base_lvl = base_lvl;
+ pinfo.base_lvlind = base_lvlind;
+ pinfo.base_lat = base_lat;
+ pinfo.base_latind = base_latind;
+ pinfo.base_lon = base_lon;
+ pinfo.base_lonind = base_lonind;
+ pinfo.comp_lvl = comp_lvl;
+ pinfo.comp_lvlind = comp_lvlind;
+
+
+ case 'plotvarvarcorrel'
+
+ disp('Getting information for the ''base'' variable.')
+ base_var = GetVar(pinfo.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
+ [base_lvl , base_lvlind] = GetLevel(levels);
+ [base_lat , base_latind] = GetLatitude(lats);
+ [base_lon , base_lonind] = GetLongitude(lons);
+
+ disp('Getting information for the ''comparison'' variable.')
+ comp_var = GetVar(pinfo.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel(levels, base_lvlind);
+ [comp_lat, comp_latind] = GetLatitude(lats, base_lat);
+ [comp_lon, comp_lonind] = GetLongitude(lons, base_lon);
+
+ pinfo.base_var = base_var;
+ pinfo.comp_var = comp_var;
+ pinfo.base_time = base_time;
+ pinfo.base_tmeind = base_tmeind;
+ pinfo.base_lvl = base_lvl;
+ pinfo.base_lvlind = base_lvlind;
+ pinfo.base_lat = base_lat;
+ pinfo.base_latind = base_latind;
+ pinfo.base_lon = base_lon;
+ pinfo.base_lonind = base_lonind;
+ pinfo.comp_lvl = comp_lvl;
+ pinfo.comp_lvlind = comp_lvlind;
+ pinfo.comp_lat = comp_lat;
+ pinfo.comp_latind = comp_latind;
+ pinfo.comp_lon = comp_lon;
+ pinfo.comp_lonind = comp_lonind;
+
+ case 'plotsawtooth'
+
+ pgvar = GetVar(pinfo.vars);
+ [level, lvlind] = GetLevel(levels);
+ [lat , latind] = GetLatitude(lats);
+ [lon , lonind] = GetLongitude(lons);
+ copyindices = SetCopyID2(fname);
+ copy = length(copyindices);
+
+ 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;
+ pinfo.latindex = latind;
+ pinfo.longitude = lon;
+ pinfo.lonindex = lonind;
+ pinfo.copies = copy;
+ pinfo.copyindices = copyindices;
+
+ if ( exist(pstruct.truth_file,'file') )
+ pinfo.truth_file = pstruct.truth_file;
+ end
+
+ case 'plotphasespace'
+
+ disp('Getting information for the ''X'' variable.')
+ var1 = GetVar(pinfo.vars);
+ [var1_lvl, var1_lvlind] = GetLevel(levels);
+ [var1_lat, var1_latind] = GetLatitude(lats);
+ [var1_lon, var1_lonind] = GetLongitude(lons);
+
+ disp('Getting information for the ''Y'' variable.')
+ var2 = GetVar(pinfo.vars, var1);
+ [var2_lvl, var2_lvlind] = GetLevel(levels, var1_lvlind);
+ [var2_lat, var2_latind] = GetLatitude(lats, var1_lat);
+ [var2_lon, var2_lonind] = GetLongitude(lons, var1_lon);
+
+ disp('Getting information for the ''Z'' variable.')
+ var3 = GetVar(pinfo.vars, var1);
+ [var3_lvl, var3_lvlind] = GetLevel(levels, var1_lvlind);
+ [var3_lat, var3_latind] = GetLatitude(lats, var1_lat);
+ [var3_lon, var3_lonind] = GetLongitude(lons, var1_lon);
+
+ % query for copy string ... ensemble member
+ [~, ens_mem_string] = GetCopyID(fname,1);
+ ens_mem = ens_mem_string{1}; % coerce to simple character string
+
+ % s1 = input('Input ensemble member metadata ID. <cr> for ''true state'' ','s');
+ % if isempty(s1), ens_mem = 'true state'; else ens_mem = s1; end
+
+ % query for line type
+ s1 = input('Input line type string. <cr> for ''k-'' ','s');
+ if isempty(s1), ltype = 'k-'; else ltype = s1; end
+
+ pinfo.var1name = var1;
+ pinfo.var2name = var2;
+ pinfo.var3name = var3;
+ pinfo.var1_lvl = var1_lvl;
+ pinfo.var1_lvlind = var1_lvlind;
+ pinfo.var1_lat = var1_lat;
+ pinfo.var1_latind = var1_latind;
+ pinfo.var1_lon = var1_lon;
+ pinfo.var1_lonind = var1_lonind;
+ pinfo.var2_lvl = var2_lvl;
+ pinfo.var2_lvlind = var2_lvlind;
+ pinfo.var2_lat = var2_lat;
+ pinfo.var2_latind = var2_latind;
+ pinfo.var2_lon = var2_lon;
+ pinfo.var2_lonind = var2_lonind;
+ pinfo.var3_lvl = var3_lvl;
+ pinfo.var3_lvlind = var3_lvlind;
+ pinfo.var3_lat = var3_lat;
+ pinfo.var3_latind = var3_latind;
+ pinfo.var3_lon = var3_lon;
+ pinfo.var3_lonind = var3_lonind;
+ pinfo.ens_mem = ens_mem;
+ pinfo.ltype = ltype;
+
+ case 'plottotalerr'
+
+ % Nothing required ... function defined in PlotTotalErr()
+
+ otherwise
+
+ error('plot type %s not supported yet for SQG model',routine)
+
+end
+
+
+
+function pgvar = GetVar(prognostic_vars, defvar)
+%----------------------------------------------------------------------
+if (nargin == 2), pgvar = defvar; else pgvar = 'theta'; end
+
+% If there is only one choice ... use it.
+if (length(prognostic_vars) == 1)
+ pgvar = prognostic_vars{1};
+ return
+end
+
+str = sprintf(' %s ',prognostic_vars{1});
+for i = 2:length(prognostic_vars),
+ str = sprintf(' %s %s ',str,prognostic_vars{i});
+end
+fprintf('Default variable is ''%s'', if this is OK, <cr>;\n',pgvar)
+fprintf('If not, please enter one of: %s\n',str)
+varstring = input('(no syntax required)\n','s');
+
+if ~isempty(varstring), pgvar = strtrim(varstring); end
+inds = strfind(pgvar,',');
+pgvar(inds) = '';
+
+
+
+function [time, timeind] = GetTime(times, deftime)
+%----------------------------------------------------------------------
+% Query for the time of interest.
+
+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), tindex = str2double(varstring); end
+
+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(levels, deflevel)
+%----------------------------------------------------------------------
+% level and lvlind will not be equal ... future expansion ...
+if (nargin == 3), lvlind = deflevel; else lvlind = 1; end
+
+fprintf('Default level (index) is %d, if this is OK, <cr>;\n',lvlind)
+fprintf('If not, enter a level between %d and %d, inclusive ...\n', ...
+ 1,length(levels))
+varstring = input('we''ll use the closest (no syntax required)\n','s');
+
+if ~isempty(varstring), lvlind = str2double(varstring); end
+
+level = levels(lvlind);
+
+
+
+function [lon, lonind] = GetLongitude(lons, deflon)
+%----------------------------------------------------------------------
+if (nargin == 2), lon = deflon; else lon = 255.0; end
+
+fprintf('Default longitude is %f, if this is OK, <cr>;\n',lon)
+fprintf('If not, enter a longitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lons),max(lons))
+varstring = input('(no syntax required)\n','s');
+
+if ~isempty(varstring), lon = str2double(varstring); end
+
+d = abs(lon - lons); % crude distance
+ind = find(min(d) == d); % multiple minima possible
+lonind = ind(1); % use the first one
+lon = lons(lonind);
+
+
+
+function [lat, latind] = GetLatitude(lats, deflat)
+%----------------------------------------------------------------------
+if (nargin == 2), lat = deflat; else lat = 40.0; end
+
+fprintf('Default latitude is %f, if this is OK, <cr>;\n',lat)
+fprintf('If not, enter a latitude between %.2f and %.2f, we use the closest.\n', ...
+ min(lats),max(lats))
+varstring = input('(no syntax required)\n','s');
+
+if ~isempty(varstring), lat = str2double(varstring); end
+
+d = abs(lat - lats); % crude distance
+ind = find(min(d) == d); % multiple minima possible
+latind = ind(1); % use the first one
+lat = lats(latind);
+
+
+
+
+function [copyid, copystrings] = GetCopyID(fname, numcopies)
+%----------------------------------------------------------------------
+%% GetCopyID queries for the copy indices in a specific netCDF file.
+
+if (exist(fname,'file') ~= 2), error('%s does not exist.',fname); end
+
+metadata = nc_varget(fname,'CopyMetaData'); % get all the metadata
+ncopies = size(metadata,1);
+
+if ( ncopies < 1 )
+ fprintf('%s has no valid ensemble members\n',fname)
+ 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)
+ error('netcdf file has no ensemble members.')
+end
+
+if (numcopies > ncopies)
+ error('Only have %d copies to work with, need %d',ncopies,numcopies)
+end
+
+for i = 1:ncopies
+ fprintf('ID %2d is %s\n',i,deblank(metadata(i,:)))
+end
+
+IDstring = input( sprintf('Enter %d IDs to plot.\n(no intervening syntax, please)\n',numcopies) ,'s');
+copyid = str2num(IDstring);
+
+if (length(copyid) ~= numcopies)
+ error('only entered %d, needed %d ... quitting',length(copyid),numcopies)
+end
+
+copystrings = cell(1,numcopies);
+for i=1:numcopies
+ copystrings{i} = metadata(copyid(i),:);
+end
+
Property changes on: DART/branches/development/matlab/GetSqgInfo.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/branches/development/matlab/PlotBins.m
===================================================================
--- DART/branches/development/matlab/PlotBins.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/PlotBins.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -61,7 +61,7 @@
'tindex1',pinfo.truth_time(1), 'tcount',pinfo.truth_time(2));
ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.var, ...
'stateindex',ivar, ...
- 'copyindex1',pinfo.ensemble_indices(1), ...
+ 'copy1',pinfo.ensemble_indices(1), ...
'copycount',pinfo.num_ens_members, ...
'tindex1',pinfo.diagn_time(1), 'tcount',pinfo.diagn_time(2));
@@ -92,7 +92,7 @@
'tindex1',pinfo.truth_time(1), 'tcount',pinfo.truth_time(2));
ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.var, ...
'stateindex',ivar, ...
- 'copyindex1',pinfo.ensemble_indices(1), ...
+ 'copy1',pinfo.ensemble_indices(1), ...
'copycount',pinfo.num_ens_members, ...
'tindex1',pinfo.diagn_time(1), 'tcount',pinfo.diagn_time(2));
@@ -122,7 +122,7 @@
'tindex1',pinfo.truth_time(1), 'tcount',pinfo.truth_time(2));
ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.var, ...
'stateindex',ivar, ...
- 'copyindex1',pinfo.ensemble_indices(1), ...
+ 'copy1',pinfo.ensemble_indices(1), ...
'copycount',pinfo.num_ens_members, ...
'tindex1',pinfo.diagn_time(1), 'tcount',pinfo.diagn_time(2));
@@ -141,7 +141,7 @@
axis tight
end
- case {'fms_bgrid','pe2lyr','mitgcm_ocean','cam','wrf','mpas_atm'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','cam','wrf','mpas_atm','mpas_ocn','sqg'}
% It is intended that all 3D models have all the required information
% set in the corresponding Get<model>Info.m script.
@@ -154,7 +154,7 @@
'tindex1',pinfo.truth_time(1), 'tcount',pinfo.truth_time(2));
ens = get_hyperslab('fname',pinfo.diagn_file, 'varname',pinfo.var, ...
'levelindex',pinfo.levelindex, ...
- 'copyindex1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
+ 'copy1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
'lonindex',pinfo.lonindex, 'latindex',pinfo.latindex, ...
'tindex1',pinfo.diagn_time(1), 'tcount',pinfo.diagn_time(2));
Modified: DART/branches/development/matlab/PlotCorrel.m
===================================================================
--- DART/branches/development/matlab/PlotCorrel.m 2013-05-09 21:28:46 UTC (rev 6124)
+++ DART/branches/development/matlab/PlotCorrel.m 2013-05-10 15:30:38 UTC (rev 6125)
@@ -55,12 +55,12 @@
% Get 'standard' ensemble series
base = get_hyperslab('fname',pinfo.fname, 'varname',pinfo.base_var, ...
- 'copyindex1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
+ 'copy1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
'stateindex',pinfo.base_var_index);
% Get (potentially large) state.
state_var = get_hyperslab('fname',pinfo.fname, 'varname',pinfo.base_var, ...
- 'copyindex1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
+ 'copy1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
'state1',pinfo.min_state_var,'statecount',pinfo.num_state_vars);
% It is efficient to preallocate correl storage ...
@@ -75,7 +75,7 @@
disp('Please be patient ... this usually takes a bit ...')
clf;
- [cs,h] = contour(pinfo.time, 1:pinfo.num_state_vars, correl, contourlevels);
+ [~,h] = contour(pinfo.time, 1:pinfo.num_state_vars, correl, contourlevels);
% clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
set(gca,'Clim',[-1 1])
hold on; % highlight the reference state variable and time
@@ -92,7 +92,7 @@
colorbar
- case 'fms_bgrid'
+ case {'fms_bgrid'}
% We are going to correlate one var/time/lvl/lat/lon with
% all other lats/lons for a var/time/lvl
@@ -115,12 +115,12 @@
nxny = nx*ny;
base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex',pinfo.base_lvlind, ...
'latindex', pinfo.base_latind, 'lonindex', pinfo.base_lonind );
bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
comp_ens = reshape(bob,[pinfo.num_ens_members, nxny]);
@@ -152,7 +152,7 @@
axis image
colorbar;
- case 'wrf'
+ case {'wrf'}
% We are going to correlate one var/time/lvl/lat/lon with
% all other lats/lons for a var/time/lvl
@@ -188,7 +188,7 @@
% Get the actual goods ... and perform the correlation
base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex',pinfo.base_lvlind, ...
'latindex', pinfo.base_latind, 'lonindex', pinfo.base_lonind );
@@ -199,7 +199,7 @@
end
bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
if (std(bob(:)) == 0.0)
@@ -243,7 +243,7 @@
axis image
colorbar;
- case 'mitgcm_ocean'
+ case {'mitgcm_ocean'}
% We are going to correlate one var/time/lvl/lat/lon with
% all other lats/lons for a var/time/lvl
@@ -271,12 +271,12 @@
nxny = nx*ny;
base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex',pinfo.base_lvlind, ...
'latindex', pinfo.base_latind, 'lonindex', pinfo.base_lonind );
bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
- 'copyindex1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
comp_ens = reshape(bob,[pinfo.num_ens_members, nxny]);
@@ -308,7 +308,7 @@
axis image
colorbar;
- case {'pe2lyr','cam'}
+ case {'pe2lyr','cam','sqg'}
% We are going to correlate one var/time/lvl/lat/lon with
% all other lats/lons for a var/time/lvl
@@ -322,12 +322,12 @@
nxny = nx*ny;
base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list