[Dart-dev] [4107] DART/trunk/diagnostics/matlab: uses only snctools functions
nancy at ucar.edu
nancy at ucar.edu
Fri Oct 16 10:33:56 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20091016/2695beda/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/PlotObsLocs.m
===================================================================
--- DART/trunk/diagnostics/matlab/PlotObsLocs.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/PlotObsLocs.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -208,8 +208,8 @@
if (~isequal(arg_epochlist,[]))
epochlist = arg_epochlist;
else
- epochtime = getnc(arg_ncfname,'time');
- epochbnds = getnc(arg_ncfname,'time_bounds');
+ epochtime = nc_varget(arg_ncfname,'time');
+ epochbnds = nc_varget(arg_ncfname,'time_bounds');
epochlist = 1:length(epochtime);
end
@@ -248,7 +248,7 @@
% in the 'Observation_Kind()' array. (this file is produced by
% running obs_diag -- full name is ObsDiagAtts.m)
- Observation_Kind = getnc(arg_ncfname,'ObservationTypes');
+ Observation_Kind = nc_varget(arg_ncfname,'ObservationTypes');
% Use a different marker and color for each observation type
% (12*6) = 72 combinations at the moment. Could extend these
Modified: DART/trunk/diagnostics/matlab/get_varnames.m
===================================================================
--- DART/trunk/diagnostics/matlab/get_varnames.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/get_varnames.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -1,11 +1,11 @@
-function bob = get_varnames(fid)
+function bob = get_varnames(fname)
% get_varnames returns JUST the 'atmospheric' variable names in the netCDF file - ie - not the coordinate variables.
%
% the result is a cell array of strings ... must use {} notation to address elements.
%
% EXAMPLE:
-% fid = netcdf('obs_seq.final.nc','nowrite');
-% varnames = get_varnames(fid);
+% fname = 'obs_seq.final.nc';
+% varnames = get_varnames(fname);
% varnames{:}
% nvars = length(varnames);
% disp(sprintf('first atmospheric variable (of %d) is %s',nvars,varnames{1}))
@@ -21,41 +21,26 @@
% $Revision$
% $Date$
-% Use the netcdf primitives 'var','coord' to get complete list
-% of (all) variables and the coordinate variables, respectively.
+fileinfo = nc_info(fname);
+Nvarnames = length(fileinfo.Dataset);
-variables = var(fid);
-coordvars = coord(fid);
-
-atmosvarinds = VariableIndices(coordvars,variables,fid);
-
-% coerce just the names into a cell array
-
-for i = 1:length(atmosvarinds)
- bob{i} = name(variables{atmosvarinds(i)});
-end
-
-
-
-function inds = VariableIndices(coordvars,variables,f)
-% returns the indices of the atmospheric variables.
-
inds = [];
-for i = 1:length(variables)
+for i = 1:Nvarnames
- varname = name(variables{i});
+ varname = fileinfo.Dataset(i).Name;
isatmosvar = 1;
-
- for j = 1:length(coordvars)
- if (strcmp( varname , name(coordvars{j}))), isatmosvar = 0; end
- end
-
+
+ % Reject the obvious coordinate variables and some that are
+ % specific to DART
+
+ if ( nc_iscoordvar(fname,varname)), isatmosvar = 0; end
if (strcmp( varname , 'time_bounds')), isatmosvar = 0; end
if (strcmp( varname , 'region_names')), isatmosvar = 0; end
if (strcmp( varname , 'CopyMetaData')), isatmosvar = 0; end
if (strcmp( varname , 'ObservationTypes')), isatmosvar = 0; end
-
+
+ % keep track of the 'good' variables
if (isatmosvar > 0)
inds = [inds i];
end
@@ -64,3 +49,9 @@
if (isempty(inds))
error('No atmospheric variables in %s',name(f))
end
+
+% coerce just the names into a cell array
+
+for i = 1:length(inds)
+ bob{i} = fileinfo.Dataset(inds(i)).Name;
+end
Modified: DART/trunk/diagnostics/matlab/get_varsNdims.m
===================================================================
--- DART/trunk/diagnostics/matlab/get_varsNdims.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/get_varsNdims.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -1,34 +1,44 @@
-function [y, ydims] = get_varsNdims(fid);
+function [y, ydims] = get_varsNdims(fname);
% Get the dimension (strings) for each atmospheric variable.
-% [y, ydims] = get_vars_dims(fid);
+% [y, ydims] = get_vars_dims(fname);
%
-% fid a netcdf file id (result of fid = netcdf(fname);
+% fname a netcdf file name
%
% y a cell array of variable names
-% ydims a cell array of the (appended) dimension names
+% ydims a cell array of the concatenated dimension names
%
% EXAMPLE:
%
-% fid = netcdf('obs_seq.final.nc','nowrite');
-% [y, ydims] = get_varsNdims(fid);
-% nvars = length(y);
-% disp(sprintf('variable %s has named dimensions %s',y{1},ydims{1}))
+% fname = 'obs_seq.final.nc';
+% [y, ydims] = get_varsNdims(fname);
+%
+% >> plotdat.allvarnames{20}
+%
+% AIRCRAFT_U_WIND_COMPONENT_guess
+%
+% >> plotdat.allvardims{20}
+% region plevel copy time
-ALLvarnames = get_varnames(fid);
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2007, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+ALLvarnames = get_varnames(fname);
Nvarnames = length(ALLvarnames);
-ix = 0;
-
for i = 1:Nvarnames
- ncvobj = ncvar(ALLvarnames{i},fid);
- dims = dim(ncvobj); % cell array ...
- dimnames = [];
- for j = 1:length(dims)
- dimnames = [dimnames ' ' name(dims{j})];
- end
+ varname = ALLvarnames{i};
+ varinfo = nc_getvarinfo(fname,varname);
- y{i} = ALLvarnames{i};
- ydims{i} = dimnames;
+ y{i} = varname;
+ ydims{i} = sprintf('%s ',varinfo.Dimension{:});
end
Modified: DART/trunk/diagnostics/matlab/nc_var_dims.m
===================================================================
--- DART/trunk/diagnostics/matlab/nc_var_dims.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/nc_var_dims.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -1,26 +1,32 @@
-function vdims = nc_var_dims(fid,varname);
+function vdims = nc_var_dims(ncfname,varname);
% Get the names of the coordinate variables for each
% of the dimensions of the variable.
%
-% vdims = nc_var_dims(fid,varname);
+% vdims = nc_var_dims(ncfname,varname);
%
-% fid a netcdf file id (result of fid = netcdf(fname);
+% ncfname file name of a netcdf file
% varname a variable names
% vdims a cell array of the coordinate variables
%
% EXAMPLE:
%
-% fid = netcdf('obs_seq.final.nc','nowrite');
+% ncfname = 'obs_seq.final.nc';
% varname = 'RADIOSONDE_TEMPERATURE_guess';
-% vdims = nc_var_dims(fid,varname);
+% vdims = nc_var_dims(ncfname,varname);
% for i = 1:length(vdims)
% disp(sprintf('variable %s dimension %d is %s',varname,i,vdims{i}))
% end
-ncvobj = ncvar(varname,fid);
-dims = dim(ncvobj);
-vdims = [];
+% Data Assimilation Research Testbed -- DART
+% Copyright 2004-2009, Data Assimilation Research Section
+% University Corporation for Atmospheric Research
+% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
-for j = 1:length(dims)
- vdims{j} = name(dims{j});
-end
+varinfo = nc_getvarinfo(ncfname,varname);
+vdims = varinfo.Dimension;
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -37,40 +37,37 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-f = netcdf(fname,'nowrite');
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.binseparation = f.bin_separation(:);
-plotdat.binwidth = f.bin_width(:);
-time_to_skip = f.time_to_skip(:);
-plotdat.rat_cri = f.rat_cri(:);
-plotdat.input_qc_threshold = f.input_qc_threshold(:);
-plotdat.lonlim1 = f.lonlim1(:);
-plotdat.lonlim2 = f.lonlim2(:);
-plotdat.latlim1 = f.latlim1(:);
-plotdat.latlim2 = f.latlim2(:);
-plotdat.biasconv = f.bias_convention(:);
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-plotdat.bincenters = f{'time', 1}(:);
-plotdat.binedges = f{'time_bounds', 1}(:);
-plotdat.mlevel = f{'mlevel', 1}(:);
-plotdat.plevel = f{'plevel', 1}(:);
-plotdat.plevel_edges = f{'plevel_edges',1}(:);
-plotdat.hlevel = f{'hlevel', 1}(:);
-plotdat.hlevel_edges = f{'hlevel_edges',1}(:);
-plotdat.nregions = length(f{'region'});
-plotdat.region_names = char(f{'region_names',1}(:));
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names')
+plotdat.region_names = deblank(region_names);
-if (plotdat.nregions == 1)
- plotdat.region_names = deblank(plotdat.region_names');
-end
-
% Coordinate between time types and dates
-timeunits = f{'time'}.units(:);
-calendar = f{'time'}.calendar(:);
+calendar = nc_attget(fname,'time','calendar');
+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));
skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
@@ -86,10 +83,10 @@
plotdat.xlabel = sprintf('bias (%s) and %s',plotdat.biasconv,copystring);
plotdat.linewidth = 2.0;
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(f);
+[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
-plotdat.nvars = length(plotdat.varnames);
+plotdat.nvars = length(plotdat.varnames);
plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.biasindex = get_copy_index(fname,'bias');
@@ -117,8 +114,8 @@
% get appropriate vertical coordinate variable
- guessdims = nc_var_dims(f, plotdat.guessvar);
- analydims = nc_var_dims(f, plotdat.analyvar);
+ guessdims = nc_var_dims(fname, plotdat.guessvar);
+ analydims = nc_var_dims(fname, plotdat.analyvar);
if ( findstr('surface',guessdims{2}) > 0 )
disp(sprintf('%s is a surface field.',plotdat.guessvar))
@@ -126,15 +123,15 @@
elseif ( findstr('undef',guessdims{2}) > 0 )
disp(sprintf('%s has no vertical definition.',plotdat.guessvar))
error('Cannot display this field this way.')
- else
- plotdat.level_org = f{guessdims{2}}(:);
- plotdat.level_units = f{guessdims{2}}.units(:);
- plotdat.nlevels = length(f{guessdims{2}});
- edgename = sprintf('%s_edges',guessdims{2});
- plotdat.level_edges = f{edgename}(:);
- plotdat.Yrange = [min(plotdat.level_edges) max(plotdat.level_edges)];
end
+ [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
+ plotdat.level_org = level_org;
+ plotdat.level_units = level_units;
+ plotdat.nlevels = nlevels;
+ plotdat.level_edges = level_edges;
+ plotdat.Yrange = Yrange;
+
% Matlab likes strictly ASCENDING order for things to be plotted,
% then you can impose the direction. The data is stored in the original
% order, so the sort indices are saved to reorder the data.
@@ -150,12 +147,9 @@
level_edges = sort(plotdat.level_edges);
plotdat.level_edges = level_edges;
- guess = getnc(fname, plotdat.guessvar,-1,-1,-1,-1,-1,-1,0);
- analy = getnc(fname, plotdat.analyvar,-1,-1,-1,-1,-1,-1,0);
-
- % guess2 = f{plotdat.guessvar,1}(:); Does not reflect _FillValue
- % analy2 = f{plotdat.analyvar,1}(:); Does not reflect _FillValue
-
+ guess = nc_varget(fname, plotdat.guessvar);
+ analy = nc_varget(fname, plotdat.analyvar);
+
% Check for one region ... if the last dimension is a singleton
% dimension, it is auto-squeezed - which is bad.
% We want these things to be 3D.
@@ -386,6 +380,29 @@
+function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname);
+% Find the vertical dimension and harvest some info
+
+varinfo = nc_getvarinfo(fname,varname);
+leveldim = [];
+
+for i = 1:length(varinfo.Dimension)
+ inds = strfind(varinfo.Dimension{i},'level');
+ if ( ~ isempty(inds)), leveldim = i; end
+end
+
+if ( isempty(leveldim) )
+ error('There is no level information for %s in %s',varname,fname)
+end
+
+level_org = nc_varget(fname,varinfo.Dimension{leveldim});
+level_units = nc_attget(fname,varinfo.Dimension{leveldim},'units');
+nlevels = varinfo.Size(leveldim);
+edgename = sprintf('%s_edges',varinfo.Dimension{leveldim});
+level_edges = nc_varget(fname, edgename);
+Yrange = [min(level_edges) max(level_edges)];
+
+
function s = ReturnBase(string1)
inds = findstr('_guess',string1);
if (inds > 0 )
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -37,36 +37,36 @@
plotdat.fname = fname;
plotdat.copystring = copystring;
-plotdat.bincenters = getnc(fname,'time');
-plotdat.binedges = getnc(fname,'time_bounds');
-plotdat.mlevel = getnc(fname,'mlevel');
-plotdat.plevel = getnc(fname,'plevel');
-plotdat.plevel_edges = getnc(fname,'plevel_edges');
-plotdat.hlevel = getnc(fname,'hlevel');
-plotdat.hlevel_edges = getnc(fname,'hlevel_edges');
-plotdat.nregions = length(getnc(fname,'region'));
-plotdat.region_names = getnc(fname,'region_names');
+plotdat.bincenters = nc_varget(fname,'time');
+plotdat.binedges = nc_varget(fname,'time_bounds');
+plotdat.mlevel = nc_varget(fname,'mlevel');
+plotdat.plevel = nc_varget(fname,'plevel');
+plotdat.plevel_edges = nc_varget(fname,'plevel_edges');
+plotdat.hlevel = nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = nc_varget(fname,'hlevel_edges');
+plotdat.ncopies = length(nc_varget(fname,'copy'));
+plotdat.nregions = length(nc_varget(fname,'region'));
+plotdat.region_names = nc_varget(fname,'region_names');
if (plotdat.nregions == 1)
plotdat.region_names = deblank(plotdat.region_names');
end
-f = netcdf(fname,'nowrite');
-plotdat.binseparation = f.bin_separation(:);
-plotdat.binwidth = f.bin_width(:);
-time_to_skip = f.time_to_skip(:);
-plotdat.rat_cri = f.rat_cri(:);
-plotdat.input_qc_threshold = f.input_qc_threshold(:);
-plotdat.lonlim1 = f.lonlim1(:);
-plotdat.lonlim2 = f.lonlim2(:);
-plotdat.latlim1 = f.latlim1(:);
-plotdat.latlim2 = f.latlim2(:);
-plotdat.biasconv = f.bias_convention(:);
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
-timeunits = f{'time'}.units(:);
-calendar = f{'time'}.calendar(:);
+calendar = nc_attget(fname,'time','calendar');
+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));
skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
@@ -81,11 +81,10 @@
plotdat.linewidth = 2.0;
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(f);
+[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindTemporalVars(plotdat);
+plotdat.nvars = length(plotdat.varnames);
-plotdat.nvars = length(plotdat.varnames);
-
plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
plotdat.Nusedindex = get_copy_index(fname,'Nused');
@@ -111,23 +110,32 @@
% get appropriate vertical coordinate variable
- guessdims = nc_var_dims(f, plotdat.guessvar);
- analydims = nc_var_dims(f, plotdat.analyvar);
+ guessdims = nc_var_dims(fname, plotdat.guessvar);
+ analydims = nc_var_dims(fname, plotdat.analyvar);
if ( findstr('surface',guessdims{3}) > 0 )
- plotdat.level = 1;
+ plotdat.level = 1;
plotdat.level_units = 'surface';
elseif ( findstr('undef',guessdims{3}) > 0 )
- plotdat.level = 1;
+ plotdat.level = 1;
plotdat.level_units = 'undefined';
else
- plotdat.level = getnc(fname, guessdims{3});
- plotdat.level_units = f{guessdims{3}}.units(:);
+ plotdat.level = nc_varget(fname, guessdims{3});
+ plotdat.level_units = nc_attget(fname, guessdims{3}, 'units');
end
+ plotdat.nlevels = length(plotdat.level);
- guess = getnc(fname, plotdat.guessvar,-1,-1,-1,-1,-1,-1,0);
- analy = getnc(fname, plotdat.analyvar,-1,-1,-1,-1,-1,-1,0);
-
+ % Here is the tricky part. Singleton dimensions are auto-squeezed ...
+ % single levels, single regions ...
+
+ guess_raw = nc_varget(fname, plotdat.guessvar);
+ guess = reshape(guess_raw, plotdat.Nbins, plotdat.ncopies, ...
+ plotdat.nlevels, plotdat.nregions);
+
+ analy_raw = nc_varget(fname, plotdat.analyvar);
+ analy = reshape(analy_raw, plotdat.Nbins, plotdat.ncopies, ...
+ plotdat.nlevels, plotdat.nregions);
+
% check to see if there is anything to plot
nposs = sum(guess(:,plotdat.Npossindex,:,:));
@@ -136,11 +144,8 @@
continue
end
- % here is where you have to check for one region ... the last
- % dimension (if it is a singleton dimension) is auto-squeezed ...
+ for ilevel = 1:plotdat.nlevels
- for ilevel = 1:length(plotdat.level)
-
plotdat.ges_copy = guess(:,plotdat.copyindex, ilevel,:);
plotdat.anl_copy = analy(:,plotdat.copyindex, ilevel,:);
@@ -235,7 +240,7 @@
switch lower(plotdat.copystring)
case 'bias'
% plot a zero-bias line
- h4 = line(t,0*t, 'Color','r','Parent',ax1);
+ h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
set(h4,'LineWidth',1.5,'LineSTyle',':')
plotdat.ylabel = sprintf('%s (%s)',plotdat.copystring,plotdat.biasconv);
otherwise
Modified: DART/trunk/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_profile.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/plot_profile.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -37,40 +37,37 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-f = netcdf(fname,'nowrite');
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.binseparation = f.bin_separation(:);
-plotdat.binwidth = f.bin_width(:);
-time_to_skip = f.time_to_skip(:);
-plotdat.rat_cri = f.rat_cri(:);
-plotdat.input_qc_threshold = f.input_qc_threshold(:);
-plotdat.lonlim1 = f.lonlim1(:);
-plotdat.lonlim2 = f.lonlim2(:);
-plotdat.latlim1 = f.latlim1(:);
-plotdat.latlim2 = f.latlim2(:);
-plotdat.biasconv = f.bias_convention(:);
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-plotdat.bincenters = f{'time', 1}(:);
-plotdat.binedges = f{'time_bounds', 1}(:);
-plotdat.mlevel = f{'mlevel', 1}(:);
-plotdat.plevel = f{'plevel', 1}(:);
-plotdat.plevel_edges = f{'plevel_edges',1}(:);
-plotdat.hlevel = f{'hlevel', 1}(:);
-plotdat.hlevel_edges = f{'hlevel_edges',1}(:);
-plotdat.nregions = length(f{'region'});
-plotdat.region_names = char(f{'region_names',1}(:));
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names')
+plotdat.region_names = deblank(region_names);
-if (plotdat.nregions == 1)
- plotdat.region_names = deblank(plotdat.region_names');
-end
-
% Coordinate between time types and dates
-timeunits = f{'time'}.units(:);
-calendar = f{'time'}.calendar(:);
+calendar = nc_attget(fname,'time','calendar');
+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));
skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
@@ -86,10 +83,10 @@
plotdat.xlabel = sprintf('%s',copystring);
plotdat.linewidth = 2.0;
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(f);
+[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
-plotdat.nvars = length(plotdat.varnames);
+plotdat.nvars = length(plotdat.varnames);
plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.Npossindex = get_copy_index(fname,'Nposs');
@@ -116,8 +113,8 @@
% get appropriate vertical coordinate variable
- guessdims = nc_var_dims(f, plotdat.guessvar);
- analydims = nc_var_dims(f, plotdat.analyvar);
+ guessdims = nc_var_dims(fname, plotdat.guessvar);
+ analydims = nc_var_dims(fname, plotdat.analyvar);
if ( findstr('surface',guessdims{2}) > 0 )
disp(sprintf('%s is a surface field.',plotdat.guessvar))
@@ -125,15 +122,15 @@
elseif ( findstr('undef',guessdims{2}) > 0 )
disp(sprintf('%s has no vertical definition.',plotdat.guessvar))
error('Cannot display this field this way.')
- else
- plotdat.level_org = f{guessdims{2}}(:);
- plotdat.level_units = f{guessdims{2}}.units(:);
- plotdat.nlevels = length(f{guessdims{2}});
- edgename = sprintf('%s_edges',guessdims{2});
- plotdat.level_edges = f{edgename}(:);
- plotdat.Yrange = [min(plotdat.level_edges) max(plotdat.level_edges)];
end
+ [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
+ plotdat.level_org = level_org;
+ plotdat.level_units = level_units;
+ plotdat.nlevels = nlevels;
+ plotdat.level_edges = level_edges;
+ plotdat.Yrange = Yrange;
+
% Matlab likes strictly ASCENDING order for things to be plotted,
% then you can impose the direction. The data is stored in the original
% order, so the sort indices are saved to reorder the data.
@@ -149,12 +146,9 @@
level_edges = sort(plotdat.level_edges);
plotdat.level_edges = level_edges;
- guess = getnc(fname, plotdat.guessvar,-1,-1,-1,-1,-1,-1,0);
- analy = getnc(fname, plotdat.analyvar,-1,-1,-1,-1,-1,-1,0);
-
- % guess2 = f{plotdat.guessvar,1}(:); Does not reflect _FillValue
- % analy2 = f{plotdat.analyvar,1}(:); Does not reflect _FillValue
-
+ guess = nc_varget(fname, plotdat.guessvar);
+ analy = nc_varget(fname, plotdat.analyvar);
+
% Check for one region ... if the last dimension is a singleton
% dimension, it is auto-squeezed - which is bad.
% We want these things to be 3D.
@@ -276,9 +270,6 @@
otherwise
end
-% disp(plotdat.myregion)
-% [plotdat.level CG' CA' ca' cg' plotdat.level_org]
-
set(gca,'YTick',plotdat.level,'Ylim',plotdat.Yrange)
ylabel(plotdat.level_units)
@@ -380,6 +371,29 @@
+function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname);
+% Find the vertical dimension and harvest some info
+
+varinfo = nc_getvarinfo(fname,varname);
+leveldim = [];
+
+for i = 1:length(varinfo.Dimension)
+ inds = strfind(varinfo.Dimension{i},'level');
+ if ( ~ isempty(inds)), leveldim = i; end
+end
+
+if ( isempty(leveldim) )
+ error('There is no level information for %s in %s',varname,fname)
+end
+
+level_org = nc_varget(fname,varinfo.Dimension{leveldim});
+level_units = nc_attget(fname,varinfo.Dimension{leveldim},'units');
+nlevels = varinfo.Size(leveldim);
+edgename = sprintf('%s_edges',varinfo.Dimension{leveldim});
+level_edges = nc_varget(fname, edgename);
+Yrange = [min(level_edges) max(level_edges)];
+
+
function s = ReturnBase(string1)
inds = findstr('_guess',string1);
if (inds > 0 )
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -38,36 +38,36 @@
plotdat.fname = fname;
plotdat.copystring = copystring;
-plotdat.bincenters = getnc(fname,'time');
-plotdat.binedges = getnc(fname,'time_bounds');
-plotdat.mlevel = getnc(fname,'mlevel');
-plotdat.plevel = getnc(fname,'plevel');
-plotdat.plevel_edges = getnc(fname,'plevel_edges');
-plotdat.hlevel = getnc(fname,'hlevel');
-plotdat.hlevel_edges = getnc(fname,'hlevel_edges');
-plotdat.nregions = length(getnc(fname,'region'));
-plotdat.region_names = getnc(fname,'region_names');
+plotdat.bincenters = nc_varget(fname,'time');
+plotdat.binedges = nc_varget(fname,'time_bounds');
+plotdat.mlevel = nc_varget(fname,'mlevel');
+plotdat.plevel = nc_varget(fname,'plevel');
+plotdat.plevel_edges = nc_varget(fname,'plevel_edges');
+plotdat.hlevel = nc_varget(fname,'hlevel');
+plotdat.hlevel_edges = nc_varget(fname,'hlevel_edges');
+plotdat.ncopies = length(nc_varget(fname,'copy'));
+plotdat.nregions = length(nc_varget(fname,'region'));
+plotdat.region_names = nc_varget(fname,'region_names');
if (plotdat.nregions == 1)
plotdat.region_names = deblank(plotdat.region_names');
end
-f = netcdf(fname,'nowrite');
-plotdat.binseparation = f.bin_separation(:);
-plotdat.binwidth = f.bin_width(:);
-time_to_skip = f.time_to_skip(:);
-plotdat.rat_cri = f.rat_cri(:);
-plotdat.input_qc_threshold = f.input_qc_threshold(:);
-plotdat.lonlim1 = f.lonlim1(:);
-plotdat.lonlim2 = f.lonlim2(:);
-plotdat.latlim1 = f.latlim1(:);
-plotdat.latlim2 = f.latlim2(:);
-plotdat.biasconv = f.bias_convention(:);
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
% Coordinate between time types and dates
-timeunits = f{'time'}.units(:);
-calendar = f{'time'}.calendar(:);
+calendar = nc_attget(fname,'time','calendar');
+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));
skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
@@ -83,11 +83,10 @@
plotdat.ylabel = sprintf('rmse and %s',copystring);
plotdat.linewidth = 2.0;
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(f);
+[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindTemporalVars(plotdat);
+plotdat.nvars = length(plotdat.varnames);
-plotdat.nvars = length(plotdat.varnames);
-
plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.rmseindex = get_copy_index(fname,'rmse');
plotdat.Npossindex = get_copy_index(fname,'Nposs');
@@ -114,23 +113,32 @@
% get appropriate vertical coordinate variable
- guessdims = nc_var_dims(f, plotdat.guessvar);
- analydims = nc_var_dims(f, plotdat.analyvar);
+ guessdims = nc_var_dims(fname, plotdat.guessvar);
+ analydims = nc_var_dims(fname, plotdat.analyvar);
if ( findstr('surface',guessdims{3}) > 0 )
- plotdat.level = 1;
+ plotdat.level = 1;
plotdat.level_units = 'surface';
elseif ( findstr('undef',guessdims{3}) > 0 )
- plotdat.level = 1;
+ plotdat.level = 1;
plotdat.level_units = 'undefined';
else
- plotdat.level = getnc(fname, guessdims{3});
- plotdat.level_units = f{guessdims{3}}.units(:);
+ plotdat.level = nc_varget(fname, guessdims{3});
+ plotdat.level_units = nc_attget(fname, guessdims{3}, 'units');
end
+ plotdat.nlevels = length(plotdat.level);
- guess = getnc(fname, plotdat.guessvar,-1,-1,-1,-1,-1,-1,0);
- analy = getnc(fname, plotdat.analyvar,-1,-1,-1,-1,-1,-1,0);
-
+ % Here is the tricky part. Singleton dimensions are auto-squeezed ...
+ % single levels, single regions ...
+
+ guess_raw = nc_varget(fname, plotdat.guessvar);
+ guess = reshape(guess_raw, plotdat.Nbins, plotdat.ncopies, ...
+ plotdat.nlevels, plotdat.nregions);
+
+ analy_raw = nc_varget(fname, plotdat.analyvar);
+ analy = reshape(analy_raw, plotdat.Nbins, plotdat.ncopies, ...
+ plotdat.nlevels, plotdat.nregions);
+
% check to see if there is anything to plot
nposs = sum(guess(:,plotdat.Npossindex,:,:));
@@ -139,15 +147,12 @@
continue
end
- % here is where you have to check for one region ... the last
- % dimension (if it is a singleton dimension) is auto-squeezed ...
+ for ilevel = 1:plotdat.nlevels
- for ilevel = 1:length(plotdat.level)
-
plotdat.ges_copy = guess(:,plotdat.copyindex, ilevel,:);
plotdat.anl_copy = analy(:,plotdat.copyindex, ilevel,:);
- plotdat.ges_rmse = guess(:,plotdat.rmseindex, ilevel,:);
- plotdat.anl_rmse = analy(:,plotdat.rmseindex, ilevel,:);
+ plotdat.ges_rmse = guess(:,plotdat.rmseindex, ilevel,:);
+ plotdat.anl_rmse = analy(:,plotdat.rmseindex, ilevel,:);
plotdat.ges_Nposs = guess(:,plotdat.Npossindex, ilevel,:);
plotdat.anl_Nposs = analy(:,plotdat.Npossindex, ilevel,:);
@@ -251,7 +256,7 @@
switch lower(plotdat.copystring)
case 'bias'
% plot a zero-bias line
- h4 = line(t,0*t, 'Color','r','Parent',ax1);
+ h4 = line(axlims(1:2),[0 0], 'Color','r','Parent',ax1);
set(h4,'LineWidth',1.5,'LineSTyle',':')
plotdat.ylabel = sprintf('%s (%s) and rmse',plotdat.copystring,plotdat.biasconv);
otherwise
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2009-10-12 23:31:04 UTC (rev 4106)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2009-10-16 16:33:55 UTC (rev 4107)
@@ -37,40 +37,37 @@
% Harvest plotting info/metadata from netcdf file.
-plotdat.fname = fname;
-plotdat.copystring = copystring;
+plotdat.fname = fname;
+plotdat.copystring = copystring;
-f = netcdf(fname,'nowrite');
+plotdat.binseparation = nc_attget(fname, nc_global, 'bin_separation');
+plotdat.binwidth = nc_attget(fname, nc_global, 'bin_width');
+time_to_skip = nc_attget(fname, nc_global, 'time_to_skip');
+plotdat.rat_cri = nc_attget(fname, nc_global, 'rat_cri');
+plotdat.input_qc_threshold = nc_attget(fname, nc_global, 'input_qc_threshold');
+plotdat.lonlim1 = nc_attget(fname, nc_global, 'lonlim1');
+plotdat.lonlim2 = nc_attget(fname, nc_global, 'lonlim2');
+plotdat.latlim1 = nc_attget(fname, nc_global, 'latlim1');
+plotdat.latlim2 = nc_attget(fname, nc_global, 'latlim2');
+plotdat.biasconv = nc_attget(fname, nc_global, 'bias_convention');
-plotdat.binseparation = f.bin_separation(:);
-plotdat.binwidth = f.bin_width(:);
-time_to_skip = f.time_to_skip(:);
-plotdat.rat_cri = f.rat_cri(:);
-plotdat.input_qc_threshold = f.input_qc_threshold(:);
-plotdat.lonlim1 = f.lonlim1(:);
-plotdat.lonlim2 = f.lonlim2(:);
-plotdat.latlim1 = f.latlim1(:);
-plotdat.latlim2 = f.latlim2(:);
-plotdat.biasconv = f.bias_convention(:);
+plotdat.bincenters = nc_varget(fname, 'time');
+plotdat.binedges = nc_varget(fname, 'time_bounds');
+plotdat.mlevel = nc_varget(fname, 'mlevel');
+plotdat.plevel = nc_varget(fname, 'plevel');
+plotdat.plevel_edges = nc_varget(fname, 'plevel_edges');
+plotdat.hlevel = nc_varget(fname, 'hlevel');
+plotdat.hlevel_edges = nc_varget(fname, 'hlevel_edges');
-plotdat.bincenters = f{'time', 1}(:);
-plotdat.binedges = f{'time_bounds', 1}(:);
-plotdat.mlevel = f{'mlevel', 1}(:);
-plotdat.plevel = f{'plevel', 1}(:);
-plotdat.plevel_edges = f{'plevel_edges',1}(:);
-plotdat.hlevel = f{'hlevel', 1}(:);
-plotdat.hlevel_edges = f{'hlevel_edges',1}(:);
-plotdat.nregions = length(f{'region'});
-plotdat.region_names = char(f{'region_names',1}(:));
+diminfo = nc_getdiminfo(fname,'region');
+plotdat.nregions = diminfo.Length;
+region_names = nc_varget(fname,'region_names')
+plotdat.region_names = deblank(region_names);
-if (plotdat.nregions == 1)
- plotdat.region_names = deblank(plotdat.region_names');
-end
-
% Coordinate between time types and dates
-timeunits = f{'time'}.units(:);
-calendar = f{'time'}.calendar(:);
+calendar = nc_attget(fname,'time','calendar');
+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));
skip_seconds = time_to_skip(4)*3600 + time_to_skip(5)*60 + time_to_skip(6);
@@ -86,10 +83,10 @@
plotdat.xlabel = sprintf('rmse and %s',copystring);
plotdat.linewidth = 2.0;
-[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(f);
+[plotdat.allvarnames, plotdat.allvardims] = get_varsNdims(fname);
[plotdat.varnames, plotdat.vardims] = FindVerticalVars(plotdat);
-plotdat.nvars = length(plotdat.varnames);
+plotdat.nvars = length(plotdat.varnames);
plotdat.copyindex = get_copy_index(fname,copystring);
plotdat.rmseindex = get_copy_index(fname,'rmse');
@@ -117,8 +114,8 @@
% get appropriate vertical coordinate variable
- guessdims = nc_var_dims(f, plotdat.guessvar);
- analydims = nc_var_dims(f, plotdat.analyvar);
+ guessdims = nc_var_dims(fname, plotdat.guessvar);
+ analydims = nc_var_dims(fname, plotdat.analyvar);
if ( findstr('surface',guessdims{2}) > 0 )
disp(sprintf('%s is a surface field.',plotdat.guessvar))
@@ -126,15 +123,15 @@
elseif ( findstr('undef',guessdims{2}) > 0 )
disp(sprintf('%s has no vertical definition.',plotdat.guessvar))
error('Cannot display this field this way.')
- else
- plotdat.level_org = f{guessdims{2}}(:);
- plotdat.level_units = f{guessdims{2}}.units(:);
- plotdat.nlevels = length(f{guessdims{2}});
- edgename = sprintf('%s_edges',guessdims{2});
- plotdat.level_edges = f{edgename}(:);
- plotdat.Yrange = [min(plotdat.level_edges) max(plotdat.level_edges)];
end
+ [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname, plotdat.guessvar);
+ plotdat.level_org = level_org;
+ plotdat.level_units = level_units;
+ plotdat.nlevels = nlevels;
+ plotdat.level_edges = level_edges;
+ plotdat.Yrange = Yrange;
+
% Matlab likes strictly ASCENDING order for things to be plotted,
% then you can impose the direction. The data is stored in the original
% order, so the sort indices are saved to reorder the data.
@@ -150,12 +147,9 @@
level_edges = sort(plotdat.level_edges);
plotdat.level_edges = level_edges;
- guess = getnc(fname, plotdat.guessvar,-1,-1,-1,-1,-1,-1,0);
- analy = getnc(fname, plotdat.analyvar,-1,-1,-1,-1,-1,-1,0);
-
- % guess2 = f{plotdat.guessvar,1}(:); Does not reflect _FillValue
- % analy2 = f{plotdat.analyvar,1}(:); Does not reflect _FillValue
-
+ guess = nc_varget(fname, plotdat.guessvar);
+ analy = nc_varget(fname, plotdat.analyvar);
+
% Check for one region ... if the last dimension is a singleton
% dimension, it is auto-squeezed - which is bad.
% We want these things to be 3D.
@@ -285,9 +279,6 @@
set(gca,'YDir', plotdat.YDir)
hold on; plot([0 0],plotdat.Yrange,'k-')
-% disp(plotdat.myregion)
-% [plotdat.level CG' CA' MG' MA' ma' mg' ca' cg' plotdat.level_org]
-
set(gca,'YTick',plotdat.level,'Ylim',plotdat.Yrange)
ylabel(plotdat.level_units)
@@ -389,6 +380,29 @@
+function [level_org level_units nlevels level_edges Yrange] = FindVerticalInfo(fname,varname);
+% Find the vertical dimension and harvest some info
+
+varinfo = nc_getvarinfo(fname,varname);
+leveldim = [];
+
+for i = 1:length(varinfo.Dimension)
+ inds = strfind(varinfo.Dimension{i},'level');
+ if ( ~ isempty(inds)), leveldim = i; end
+end
+
+if ( isempty(leveldim) )
+ error('There is no level information for %s in %s',varname,fname)
+end
+
+level_org = nc_varget(fname,varinfo.Dimension{leveldim});
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list