[Dart-dev] [4834] DART/trunk/matlab: Full support for WRF for all routines that do not require a True_State.nc

nancy at ucar.edu nancy at ucar.edu
Thu Mar 31 17:12:11 MDT 2011


Revision: 4834
Author:   thoar
Date:     2011-03-31 17:12:11 -0600 (Thu, 31 Mar 2011)
Log Message:
-----------
Full support for WRF for all routines that do not require a True_State.nc
Those routines are untested until I get someone to perform a perfect model experiment.
The filename is now annotated in the title on a separate line.

The parse.m routine breaks a line into words - and then creates a cell array.

Modified Paths:
--------------
    DART/trunk/matlab/CheckModelCompatibility.m
    DART/trunk/matlab/GetBgridInfo.m
    DART/trunk/matlab/GetWRFInfo.m
    DART/trunk/matlab/PlotCorrel.m
    DART/trunk/matlab/PlotEnsErrSpread.m
    DART/trunk/matlab/PlotEnsMeanTimeSeries.m
    DART/trunk/matlab/PlotJeffCorrel.m
    DART/trunk/matlab/PlotPhaseSpace.m
    DART/trunk/matlab/PlotSawtooth.m
    DART/trunk/matlab/PlotTotalErr.m
    DART/trunk/matlab/SetCopyID2.m
    DART/trunk/matlab/get_copy_index.m
    DART/trunk/matlab/plot_correl.m
    DART/trunk/matlab/plot_ens_err_spread.m
    DART/trunk/matlab/plot_ens_mean_time_series.m
    DART/trunk/matlab/plot_ens_time_series.m
    DART/trunk/matlab/plot_jeff_correl.m
    DART/trunk/matlab/plot_phase_space.m
    DART/trunk/matlab/plot_sawtooth.m
    DART/trunk/matlab/plot_total_err.m
    DART/trunk/matlab/plot_var_var_correl.m

Added Paths:
-----------
    DART/trunk/matlab/parse.m

-------------- next part --------------
Modified: DART/trunk/matlab/CheckModelCompatibility.m
===================================================================
--- DART/trunk/matlab/CheckModelCompatibility.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/CheckModelCompatibility.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -22,8 +22,6 @@
 if (nargin == 1)      % better be a pinfo struct with at least these fields
   file1 = arg1.truth_file;  % string
   file2 = arg1.diagn_file;  % string
-  time1 = arg1.truth_time;  % [a,b] array
-  time2 = arg1.diagn_time;  % [a,b] array
   pinfo_out = arg1;
 elseif (nargin == 2)  % pair of filenames
   file1 = arg1;             % truth_file
@@ -38,10 +36,10 @@
 if ( exist(file2,'file') ~= 2 ), error('(file2) %s does not exist.',file2); end
 
 % set this up for later
-pinfo_out = setfield(pinfo_out, 'truth_time', [-1,-1]);
-pinfo_out = setfield(pinfo_out, 'diagn_time', [-1,-1]);
+pinfo_out.truth_time = [-1,-1];
+pinfo_out.diagn_time = [-1,-1];
 
-% Get some information from the file1
+%% Get some information from the file1
 tmodel  = nc_attget(file1,nc_global,'model');
 
 if (isempty(tmodel)) 
@@ -50,15 +48,18 @@
 
 tnum_copies = dim_length(file1,'copy');
 tnum_times  = dim_length(file1,'time');
-ttimes      =  nc_varget(file1,'time');
+times       = nc_varget( file1,'time');
+timeunits   = nc_attget( file1,'time','units');
+timebase    = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin  = datenum(timebase(1),timebase(2),timebase(3));
+ttimes      = times + timeorigin;
 
 [tnum_vars,tdims] = ModelDimension(file1,tmodel);
 if ( tnum_vars <= 0 )
    error('Unable to determine resolution of %s.',file1)
 end
 
-
-% Get some information from the file2
+%% Get some information from the file2
 dmodel  = nc_attget(file1,nc_global,'model');
 
 if (isempty(dmodel)) 
@@ -67,7 +68,11 @@
 
 dnum_copies = dim_length(file2,'copy');
 dnum_times  = dim_length(file2,'time');
-dtimes      =  nc_varget(file2,'time');
+times       = nc_varget( file2,'time');
+timeunits   = nc_attget( file2,'time','units');
+timebase    = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin  = datenum(timebase(1),timebase(2),timebase(3));
+dtimes      = times + timeorigin;
 
 [dnum_vars,ddims] = ModelDimension(file2,dmodel);
 if ( dnum_vars <= 0 )
@@ -80,7 +85,7 @@
    fprintf('%s has model %s\n',file2,dmodel)
    error('no No NO ... models must be the same')
 end
-pinfo_out = setfield(pinfo_out, 'model', tmodel);
+pinfo_out.model = tmodel;
 
 if (prod(tnum_vars) ~= prod(dnum_vars))
    fprintf('%s has %d state variables\n',file1,prod(tnum_vars))
@@ -223,12 +228,19 @@
 
 
 function [x,y] = ModelDimension(fname,modelname)
-
+% Check the base geometry of the grid
 x = 0;
 y = NaN;
 
 switch lower(modelname)
 
+   case 'wrf'
+      diminfo = nc_getdiminfo(fname,  'west_east_d01'); dnum_lons = diminfo.Length;
+      diminfo = nc_getdiminfo(fname,'south_north_d01'); dnum_lats = diminfo.Length;
+      diminfo = nc_getdiminfo(fname, 'bottom_top_d01'); dnum_lvls = diminfo.Length;
+         x = 3;
+         y = [dnum_lons dnum_lats dnum_lvls];
+
    case 'cam'
       dnum_lons = dim_length(fname,'lon');
       dnum_lats = dim_length(fname,'lat');

Modified: DART/trunk/matlab/GetBgridInfo.m
===================================================================
--- DART/trunk/matlab/GetBgridInfo.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/GetBgridInfo.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -173,9 +173,14 @@
       [var3_lat, var3_latind] = GetLatitude( var3, TmpJ, VelJ, var1_lat);
       [var3_lon, var3_lonind] = GetLongitude(var3, TmpI, VelI, var1_lon);
 
-      % query for ensemble member
-      s1 = input('Input ensemble member metadata STRING. <cr> for ''true state''  ','s');
-      if isempty(s1), ens_mem = 'true state'; else ens_mem = s1; end
+      % query for ensemble member string
+      metadata   = nc_varget(fname,'CopyMetaData');
+      [N,M]      = size(metadata);
+      cell_array = mat2cell(metadata, ones(1,N), M);
+      ens_mem    = strtrim(cell_array{1});
+      str1 = sprintf('Input ensemble member metadata STRING. <cr> for ''%s''   ',ens_mem);
+      s1   = input(str1,'s');
+      if ~ isempty(s1), ens_mem = s1; end
 
       % query for line type
       s1 = input('Input line type string. <cr> for ''k-''  ','s');

Modified: DART/trunk/matlab/GetWRFInfo.m
===================================================================
--- DART/trunk/matlab/GetWRFInfo.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/GetWRFInfo.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -166,19 +166,17 @@
    case 'plotsawtooth'
 
        pgvar          = GetVarString(pinfo_in.vars);
-      [level, lvlind] = GetLevel(fname,pgvar);   % Determine level and index
-      [lat  , latind] = GetLatitude( pgvar,TmpJ,VelJ);
-      [lon  , lonind] = GetLongitude(pgvar,TmpI,VelI);
-      %[copy , lonind] = GetCopies(pgvar,copy);
-      copyindices     = SetCopyID(fname);
+      [level, lvlind] = GetLevel(pinfo_in.prior_file, pgvar);
+      [lat, lon, latind, lonind] = GetLatLon(pinfo_in.prior_file, pgvar);
+      [copyindices, copymetadata]= SetCopyID2(pinfo_in.prior_file);
       copy            = length(copyindices);
 
       pinfo.model          = model;
+      pinfo.truth_file     = pinfo_in.truth_file;
+      pinfo.prior_file     = pinfo_in.prior_file;
+      pinfo.posterior_file = pinfo_in.posterior_file;
       pinfo.times          = dates;
       pinfo.var_names      = pgvar;
-     %pinfo.truth_file     = [];
-     %pinfo.prior_file     = pinfo.prior_file;
-     %pinfo.posterior_file = pinfo.posterior_file;
       pinfo.level          = level;
       pinfo.levelindex     = lvlind;
       pinfo.latitude       = lat;
@@ -188,29 +186,44 @@
       pinfo.copies         = copy;
       pinfo.copyindices    = copyindices;
 
+      % Hui has a habit of cutting out the rest of the ensemble members
+      % from the posterior, but leaving them in the prior. Thanks.
+      % So now I have to figure out if the posterior and prior copy metadata match.
+
+      for i = 1:copy,
+         copyi = get_copy_index(pinfo_in.posterior_file,copymetadata{i}); 
+         pstruct.postcopyindices = copyi;
+      end
+
+      if ( any(pstruct.postcopyindices < 1) )
+         error('The prior copy does not exist in the posterior file.')
+      end
+
    case 'plotphasespace'
 
       disp('Getting information for the ''X'' variable.')
        var1                   = GetVarString(pinfo_in.vars);
       [var1_lvl, var1_lvlind] = GetLevel(fname,var1);
-      [var1_lat, var1_latind] = GetLatitude( var1, TmpJ, VelJ);
-      [var1_lon, var1_lonind] = GetLongitude(var1, TmpI, VelI);
+      [var1_lat, var1_lon, var1_latind, var1_lonind] = GetLatLon(fname, var1);
 
       disp('Getting information for the ''Y'' variable.')
        var2                   = GetVarString(pinfo_in.vars,        var1    );
       [var2_lvl, var2_lvlind] = GetLevel(fname,var2,var1_lvl);
-      [var2_lat, var2_latind] = GetLatitude( var2, TmpJ, VelJ, var1_lat);
-      [var2_lon, var2_lonind] = GetLongitude(var2, TmpI, VelI, var1_lon);
+      [var2_lat, var2_lon, var2_latind, var2_lonind] = GetLatLon(fname, var2);
 
       disp('Getting information for the ''Z'' variable.')
        var3                   = GetVarString(pinfo_in.vars,        var1    );
       [var3_lvl, var3_lvlind] = GetLevel(fname,var3,var1_lvl);
-      [var3_lat, var3_latind] = GetLatitude( var3, TmpJ, VelJ, var1_lat);
-      [var3_lon, var3_lonind] = GetLongitude(var3, TmpI, VelI, var1_lon);
+      [var3_lat, var3_lon, var3_latind, var3_lonind] = GetLatLon(fname, var3);
 
-      % query for ensemble member
-      s1 = input('Input ensemble member metadata STRING. <cr> for ''true state''  ','s');
-      if isempty(s1), ens_mem = 'true state'; else ens_mem = s1; end
+      % query for ensemble member string
+      metadata   = nc_varget(fname,'CopyMetaData');
+      [N,M]      = size(metadata);
+      cell_array = mat2cell(metadata, ones(1,N), M);
+      ens_mem    = strtrim(cell_array{1});
+      str1 = sprintf('Input ensemble member metadata STRING. <cr> for ''%s''   ',ens_mem);
+      s1   = input(str1,'s');
+      if ~ isempty(s1), ens_mem = s1; end
 
       % query for line type
       s1 = input('Input line type string. <cr> for ''k-''  ','s');

Modified: DART/trunk/matlab/PlotCorrel.m
===================================================================
--- DART/trunk/matlab/PlotCorrel.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotCorrel.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -86,11 +86,11 @@
       clf;
       
       contour(correl,-1:0.2:1);
-      s1 = sprintf('%s Correlation of variable %s index %d, T = %d of %s', ...
-               model, pinfo.base_var, base_var_index, base_time, pinfo.fname);
+      s1 = sprintf('%s Correlation of variable %s index %d, T = %d', ...
+               model, pinfo.base_var, base_var_index, base_time);
       s2 = sprintf('against all variables, all times, %d ensemble members', ...
                size(state_var,2)); 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel('time (timestep #)')
       ylabel('state variable (index)')
       set(gca,'YTick',1:num_vars)
@@ -143,13 +143,13 @@
       contour(lons,lats,correl,-1:0.2:1); hold on;
       plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
                  'MarkerSize',12,'MarkerFaceColor','k');
-      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
            model, pinfo.base_var, pinfo.base_lvl, ...
-             pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+             pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
 
       s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
                pinfo.comp_var, pinfo.comp_lvl, nmembers); 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
       ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
       worldmap;
@@ -290,13 +290,13 @@
       contour(lons,lats,correl,-1:0.2:1); hold on;
       plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
                  'MarkerSize',12,'MarkerFaceColor','k');
-      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
            model, pinfo.base_var, pinfo.base_lvl, ...
-             pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+             pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
 
       s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
                pinfo.comp_var, pinfo.comp_lvl, nmembers); 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
       ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
       worldmap;
@@ -337,13 +337,13 @@
       contour(lons,lats,correl,-1:0.2:1); hold on;
       plot(pinfo.base_lon, pinfo.base_lat, 'pk', ...
                  'MarkerSize',12,'MarkerFaceColor','k');
-      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f of %s', ...
+      s1 = sprintf('%s Correlation of ''%s'', level %d, (%.2f,%.2f) T = %f', ...
            model, pinfo.base_var, pinfo.base_lvl, ...
-             pinfo.base_lat, pinfo.base_lon, pinfo.base_time, pinfo.fname);
+             pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
 
       s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
                pinfo.comp_var, pinfo.comp_lvl, nmembers); 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
       ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
       worldmap;

Modified: DART/trunk/matlab/PlotEnsErrSpread.m
===================================================================
--- DART/trunk/matlab/PlotEnsErrSpread.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotEnsErrSpread.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -107,9 +107,8 @@
             subplot(3, 1, j);
                plot(times,err, 'b', ...
                     times,ens_spread(:, ivar), 'r');
-               s1 = sprintf('%s model Var %d Ensemble Error Spread for %s', ...
-                            tmodel, ivar, pinfo.diagn_file);
-               title(s1,'interpreter','none','fontweight','bold')
+               s1 = sprintf('%s model Var %d Ensemble Error Spread', tmodel, ivar);
+               title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
                legend(string1,string2,0)
                legend boxoff
                xlabel(sprintf('model time (%d timesteps)',tnum_times))
@@ -139,9 +138,8 @@
             subplot(length(pinfo.var_inds), 1, iplot);
                plot(times,err, 'b', ...
                     times,ens_spread(:, ivar), 'r');
-               s1 = sprintf('%s model Var %d Ensemble Error Spread for %s', ...
-                            tmodel, ivar, pinfo.diagn_file);
-               title(s1,'interpreter','none','fontweight','bold')
+               s1 = sprintf('%s model Var %d Ensemble Error Spread', tmodel, ivar);
+               title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
                legend(string1,string2,0)
                legend boxoff
                xlabel(sprintf('model time (%d timesteps)',tnum_times))
@@ -170,15 +168,11 @@
          string2 = ['time-mean Ensemble Spread = ' num2str(spreadTotal)];
 
          plot(times,err, 'b', times,ens_spread, 'r');
-         s1 = sprintf('%s model Var %s Ensemble Error Spread for %s', ...
-                            tmodel, pinfo.var, pinfo.diagn_file);
-         title(s1,'interpreter','none','fontweight','bold');
 
-      title({ ...
-        sprintf('Ensemble Mean Error, Ensemble Spread  %s ''%s'' for %s', ...
-                tmodel, pinfo.var, pinfo.diagn_file), ...
-        sprintf('level %d lat %.2f lon %.2f',pinfo.level, pinfo.latitude, ...
-                 pinfo.longitude)}, 'interpreter','none','fontweight','bold');
+         s1 = sprintf('Ensemble Mean Error, Ensemble Spread %s ''%s''',tmodel,pinfo.var);
+         s2 = sprintf('level %d lat %.2f lon %.2f', ...
+                       pinfo.level, pinfo.latitude, pinfo.longitude);
+         title({s1, s2, pinfo.diagn_file},'interpreter','none','fontweight','bold');
 
          legend(string1,string2,0);
          legend boxoff

Modified: DART/trunk/matlab/PlotEnsMeanTimeSeries.m
===================================================================
--- DART/trunk/matlab/PlotEnsMeanTimeSeries.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotEnsMeanTimeSeries.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -131,8 +131,8 @@
                legend('Ensemble Mean','True State',0);
             end
 
-            title(sprintf('%s Variable %d of %s',pinfo.model,ivar,pinfo.diagn_file), ...
-                     'interpreter','none','fontweight','bold')
+            s1 = sprintf('%s Variable %d',pinfo.model,ivar);
+            title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
             xlabel(sprintf('model time (%d timesteps)',num_times))
             legend boxoff
       end
@@ -178,8 +178,8 @@
                hold on; plot(times,truth,'b'); hold off;
                legend('Ensemble Mean','True State',0)
             end
-            title(sprintf('%s Variable %d of %s',pinfo.model,ivar,pinfo.diagn_file), ...
-                     'interpreter','none','fontweight','bold')
+            s1 = sprintf('%s Variable %d',pinfo.model,ivar);
+            title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
             xlabel(sprintf('model time (%d timesteps)',num_times))
             legend boxoff
       end
@@ -205,7 +205,7 @@
                             pinfo.model, pinfo.var, pinfo.diagn_file);
          s2 = sprintf('level %d lat %.2f lon %.2f', ...
                     pinfo.level, pinfo.latitude, pinfo.longitude);
-         title({s1,s2},'interpreter','none','fontweight','bold')
+         title({s1,s2,pinfo.diagn_file},'interpreter','none','fontweight','bold')
 
          if (have_truth)
             truth = GetCopy(pinfo.truth_file, truth_index,      pinfo , ...
@@ -244,8 +244,7 @@
          subplot(2,1,2)
             plot(times,ens_mean,'r','LineWidth',2);
             legend('Ensemble Mean', 0)
-            s1 = sprintf('%s model ''%s'' %s Ensemble Mean ', ...
-                 pinfo.model, pinfo.var, pinfo.diagn_file);
+            s1 = sprintf('%s model ''%s'' Ensemble Mean', pinfo.model, pinfo.var);
             s2 = sprintf('level index %d lat %.2f lon %.2f', ...
                        pinfo.levelindex, pinfo.latitude, pinfo.longitude);
 
@@ -254,13 +253,13 @@
                                   pinfo.truth_time(1), pinfo.truth_time(2)) ;
                hold on; plot(times,truth,'b','LineWidth',2); hold off;
                legend('Ensemble Mean','True State',0);
-               s1 = sprintf('%s model ''%s'' %s Truth and Ensemble Mean ', ...
-                               pinfo.model, pinfo.var, pinfo.diagn_file);
+               s1 = sprintf('%s model ''%s'' Truth and Ensemble Mean', ...
+                               pinfo.model, pinfo.var);
             end
 
             %plot(times,ens_mean,'r','LineWidth',2); %      again - on top
 
-            title({s1,s2},'interpreter','none','fontweight','bold')
+            title({s1,s2,pinfo.diagn_file},'interpreter','none','fontweight','bold')
             xlabel(sprintf('time (%s) %d timesteps',timeunits, num_times))
             ylabel(varunits)
             legend boxoff

Modified: DART/trunk/matlab/PlotJeffCorrel.m
===================================================================
--- DART/trunk/matlab/PlotJeffCorrel.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotJeffCorrel.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -46,7 +46,7 @@
 
 switch lower(model)
 
-   case {'fms_bgrid','pe2lyr','mitgcm_ocean'}
+   case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
 
       clf;
 
@@ -69,11 +69,11 @@
           model, pinfo.base_var, pinfo.base_time, pinfo.base_lvl, ...
           pinfo.base_lat, pinfo.base_lon);
 
-      s2 = sprintf('with ''%s'', lvl = %d, lat = %.2f, lon= %.2f, %d ensemble members -- %s', ...
+      s2 = sprintf('with ''%s'', lvl = %d, lat = %.2f, lon= %.2f, %d ensemble members', ...
           pinfo.comp_var, pinfo.comp_lvl, pinfo.comp_lat, pinfo.comp_lon, ...
-          nmembers, pinfo.fname); 
+          nmembers); 
 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel(sprintf('time (%s) %d timesteps',timeunits, num_times))
       ylabel('correlation')
       
@@ -83,7 +83,6 @@
       plot([pinfo.base_time pinfo.base_time],[ -1 1 ],'k:', ...
            [ax(1)         ax(2)],[  0 0 ],'k:')
 
-
    otherwise
 
       num_vars   = dim_length(pinfo.fname,'StateVariable');
@@ -123,8 +122,8 @@
       s1 = sprintf('%s Correlation of variable %s %d, T = %d, with variable %s %d', ...
                model, pinfo.base_var, pinfo.base_var_index, pinfo.base_time, ...
                       pinfo.state_var, pinfo.state_var_index);
-      s2 = sprintf('%d ensemble members -- %s', nmembers, pinfo.fname); 
-      title({s1,s2},'interpreter','none','fontweight','bold')
+      s2 = sprintf('%d ensemble members', nmembers); 
+      title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
       xlabel('time (timestep #)')
       ylabel('correlation')
       

Modified: DART/trunk/matlab/PlotPhaseSpace.m
===================================================================
--- DART/trunk/matlab/PlotPhaseSpace.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotPhaseSpace.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -163,11 +163,8 @@
       end
       legend boxoff
 
-   case {'fms_bgrid','pe2lyr','mitgcm_ocean'}
+   case {'fms_bgrid','pe2lyr','mitgcm_ocean','wrf'}
 
-      disp('PlotPhaseSpace')
-      pinfo
-
       ens_mem_id = get_copy_index(pinfo.fname, pinfo.ens_mem);   % errors out if no ens_mem 
       
       x = Get1Copy(pinfo.fname, pinfo.var1name, ens_mem_id,  ...
@@ -190,19 +187,22 @@
          % title(sprintf('%s ensemble member %d of %s',model,ens_mem_id,pinfo.fname), ...
          %    'interpreter','none','fontweight','bold')
          title('The Attractor','fontweight','bold')
-         xlabel(sprintf('%s %d %.2f %.2f', ...
-               pinfo.var1name, pinfo.var1_lvl, pinfo.var1_lat, pinfo.var1_lon))
-         ylabel(sprintf('%s %d %.2f %.2f', ...
-               pinfo.var2name, pinfo.var2_lvl, pinfo.var2_lat, pinfo.var2_lon))
-         zlabel(sprintf('%s %d %.2f %.2f', ...
-               pinfo.var3name, pinfo.var3_lvl, pinfo.var3_lat, pinfo.var3_lon))
+         hx = xlabel(sprintf('%s %d %.2f %.2f', ...
+               pinfo.var1name, pinfo.var1_lvl, pinfo.var1_lat, pinfo.var1_lon));
+         hy = ylabel(sprintf('%s %d %.2f %.2f', ...
+               pinfo.var2name, pinfo.var2_lvl, pinfo.var2_lat, pinfo.var2_lon));
+         hz = zlabel(sprintf('%s %d %.2f %.2f', ...
+               pinfo.var3name, pinfo.var3_lvl, pinfo.var3_lat, pinfo.var3_lon));
+         set(hx,'interpreter','none')
+         set(hy,'interpreter','none')
+         set(hz,'interpreter','none')
       
          s = sprintf('%s %s %s %s %s %s', pinfo.var1name, pinfo.var2name, pinfo.var3name, ...
                                               model, pinfo.fname, pinfo.ens_mem);
          h = legend(s,0);
          [legh, objh, outh, outm] = legend;
          set(objh(1),'interpreter','none')
-      
+
       else
          % Must add salient information to the legend.
          % legh     handle to the legend axes

Modified: DART/trunk/matlab/PlotSawtooth.m
===================================================================
--- DART/trunk/matlab/PlotSawtooth.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotSawtooth.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -40,7 +40,7 @@
 % pinfo.longitude  = 45.67;
 % PlotSawtooth( pinfo )
 
-%% DART software - Copyright \xA9 2004 - 2010 UCAR. This open source software is
+%% DART software - Copyright � 2004 - 2010 UCAR. This open source software is
 % provided by UCAR, "as is", without charge, subject to all terms of use at
 % http://www.image.ucar.edu/DAReS/DART/DART_download
 %
@@ -54,6 +54,7 @@
 pstruct = CheckModelCompatibility(pinfo.prior_file, pinfo.posterior_file);
 pinfo.prior_times     = pstruct.truth_time;
 pinfo.posterior_times = pstruct.diagn_time;
+pinfo.model           = pstruct.model;
 
 % Get some information from the truth_file, if it exists.
 if ( exist(pinfo.truth_file,'file') == 2 ) 
@@ -68,20 +69,17 @@
    truth = [];
 end
 
-% Get some information from the prior_file 
+%% Get some information from the prior_file 
+%  CheckModelCompatibility ensures that getting this from one is sufficient.
 prior.model     = nc_attget(pinfo.prior_file, nc_global, 'model');
 prior.num_times = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
-prior.times     = nc_varget(pinfo.prior_file, 'time', ...
+times           = nc_varget(pinfo.prior_file, 'time', ...
                         pinfo.prior_times(1)-1, prior.num_times);
+timeunits       = nc_attget(pinfo.prior_file,'time','units');
+timebase        = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin      = datenum(timebase(1),timebase(2),timebase(3));
+prior.times     = times + timeorigin;
 
-% Get some information from the posterior_file 
-post.model     = nc_attget(pinfo.posterior_file, nc_global, 'model');
-post.num_times = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
-post.times     = nc_varget(pinfo.posterior_file, 'time', ...
-                        pinfo.posterior_times(1)-1, post.num_times);
-post.timeunits = nc_attget(pinfo.posterior_file, 'time', 'units');
-post.calendar  = nc_attget(pinfo.posterior_file, 'time', 'calendar');
-
 % Get the indices for the true state and ensemble mean
 % The metadata is queried to determine which "copy" is appropriate.
 prior.ens_mean_index = get_copy_index(pinfo.prior_file,     'ensemble mean');
@@ -96,217 +94,215 @@
 pinfo.xax  = x(:);
 metadata   = nc_varget(pinfo.prior_file,'CopyMetaData'); % get all the metadata
 
-pinfo
-
 if isfield(pinfo,'var_inds')
    PlotGivenIndices( pinfo, prior, post, truth, metadata);
 elseif isfield(pinfo,'var_names')
-   PlotGivenVariable(pinfo, prior, post, truth, metadata);
+   PlotGivenVariable(pinfo, truth, metadata);
 end
 
-%------------------------------------------------------------------------------
-% Plot given an index into the state-space vector
-%------------------------------------------------------------------------------
 
+
 function PlotGivenIndices(pinfo, prior, post, truth, metadata)
 
-   nfigs = length(pinfo.var_inds);  % each variable gets its own figure
-   iplot = 0;
+%% Plot given an index into the state-space vector
 
-   for ivar = pinfo.var_inds,
+nfigs = length(pinfo.var_inds);  % each variable gets its own figure
+iplot = 0;
 
-      % Get the data from the netcdf files
+for ivar = pinfo.var_inds,
 
-      po_ens_mean = get_var_series(pinfo.posterior_file, pinfo.var, ...
-                    post.ens_mean_index, ivar, ...
-                    pinfo.posterior_times(1), pinfo.posterior_times(2));
-      pr_ens_mean = get_var_series(pinfo.prior_file, pinfo.var, ...
-                    prior.ens_mean_index, ivar, ...
-                    pinfo.prior_times(1), pinfo.prior_times(2));
+   % Get the data from the netcdf files
 
-      % Now we paste them together in a clever way to show 
-      % the effect of the assimilation
-      ens_mean(  1,:) = pr_ens_mean;
-      ens_mean(  2,:) = po_ens_mean;
-      a               = ens_mean(:);
+   po_ens_mean = get_var_series(pinfo.posterior_file, pinfo.var, ...
+                 post.ens_mean_index, ivar, ...
+                 pinfo.posterior_times(1), pinfo.posterior_times(2));
+   pr_ens_mean = get_var_series(pinfo.prior_file, pinfo.var, ...
+                 prior.ens_mean_index, ivar, ...
+                 pinfo.prior_times(1), pinfo.prior_times(2));
 
-      % Plot the true trajectory if it exists; the ens mean; annotate
+   % Now we paste them together in a clever way to show 
+   % the effect of the assimilation
+   ens_mean(  1,:) = pr_ens_mean;
+   ens_mean(  2,:) = po_ens_mean;
+   a               = ens_mean(:);
 
-      iplot = iplot + 1;
-      figure(iplot); clf; 
+   % Plot the true trajectory if it exists; the ens mean; annotate
 
-      if ( exist(pinfo.truth_file,'file') == 2 ) 
-         true_trajectory = get_var_series(pinfo.truth_file, pinfo.var, truth.truth_index, ivar);
+   iplot = iplot + 1;
+   figure(iplot); clf; 
+
+   if ( exist(pinfo.truth_file,'file') == 2 ) 
+      true_trajectory = get_var_series(pinfo.truth_file, pinfo.var, truth.truth_index, ivar);
 	 h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
 	 h = plot(pinfo.xax, a, 'k-','linewidth',2.0);
 	 legend('truth','ensemble mean')
-      else 
+   else 
 	 h = plot(pinfo.xax, a, 'k-','linewidth',2.0); hold on;
 	 legend('ensemble mean')
-      end
+   end
 
-      ylabel(sprintf('''%s'' index %d', pinfo.var, ivar ))
-      xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
-      title(sprintf('%s Trajectories',pinfo.model), ...  
-                     'interpreter','none','fontweight','bold')
+   ylabel(sprintf('''%s'' index %d', pinfo.var, ivar ))
+   xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
+   title(sprintf('%s Trajectories',pinfo.model), ...  
+                  'interpreter','none','fontweight','bold')
 
-      % Now check to see if we are overlaying any individual ensemble members.
-      % if pinfo.copyindices = [], nothing happens. 
+   % Now check to see if we are overlaying any individual ensemble members.
+   % if pinfo.copyindices = [], nothing happens. 
 
-      ens_colors = get(gca,'ColorOrder');   % trying to cycle through colors
-      ncolors = size(ens_colors,1) - 1;     % last one is black, already used.
+   ens_colors = get(gca,'ColorOrder');   % trying to cycle through colors
+   ncolors = size(ens_colors,1) - 1;     % last one is black, already used.
 
-      nmem = 0;
-      for imem = pinfo.copyindices
-         nmem = nmem + 1;
+   nmem = 0;
+   for imem = pinfo.copyindices
+      nmem = nmem + 1;
 
 	 %str1 = sprintf('ensemble member %d',imem);
 	 str1 = deblank(metadata(imem,:));
-         copy_index = get_copy_index(pinfo.prior_file, str1);
-         po_series  = get_var_series(pinfo.posterior_file, pinfo.var, ...
-                copy_index, ivar, pinfo.posterior_times(1), pinfo.posterior_times(2));
-         pr_series  = get_var_series(pinfo.prior_file,     pinfo.var, ...
-                copy_index, ivar, pinfo.prior_times(1), pinfo.prior_times(2));
+      copy_index = get_copy_index(pinfo.prior_file, str1);
+      po_series  = get_var_series(pinfo.posterior_file, pinfo.var, ...
+             copy_index, ivar, pinfo.posterior_times(1), pinfo.posterior_times(2));
+      pr_series  = get_var_series(pinfo.prior_file,     pinfo.var, ...
+             copy_index, ivar, pinfo.prior_times(1), pinfo.prior_times(2));
 
 	 ens_member(1,:) = pr_series;
 	 ens_member(2,:) = po_series;
 	 b               = ens_member(:);
 
-   	 hold on; 
+	 hold on; 
 	 memcolor = 1 + mod(nmem-1,ncolors); % cycles through colors [1,6]
 	 h = plot(pinfo.xax, b,'linewidth',0.5,'Color',ens_colors(memcolor,:));
 
 	 [legh, objh, outh, outm] = legend;
 	 nlines = length(outm);
-         outm{nlines+1} = str1;
+      outm{nlines+1} = str1;
 	 [legh, objh, outh, outm] = legend([outh; h],outm,0);
-      end
-      legend boxoff
-
    end
+   legend boxoff
 
-%------------------------------------------------------------------------------
-% Plot given an index into the state-space vector
-%------------------------------------------------------------------------------
+end
 
-function PlotGivenVariable(pinfo, prior, post, truth, metadata)
 
-   var_names = strread(pinfo.var_names,'%s','delimiter',' ');
 
-   nfigs = length(var_names);  % each variable gets its own figure
-   iplot = 0;
+function PlotGivenVariable(pinfo, truth, metadata)
+%% Plot given an index into the state-space vector
 
-   for ivar = 1:nfigs
+var_names = parse(pinfo.var_names);
 
-      iplot = iplot + 1;
-      figure(iplot); clf; 
+nfigs = length(var_names);  % each variable gets its own figure
+iplot = 0;
 
-      % Get the data from the netcdf files
+for ivar = 1:nfigs
 
-      vname  = var_names{ivar};
-            
-      [start, count] = GetNCindices(pinfo,'prior',vname);
+   iplot = iplot + 1;
+   figure(iplot); clf; 
 
-      ens_colors = get(gca,'ColorOrder');   % trying to cycle through colors
-      ncolors = size(ens_colors,1) - 1;     % last one is black, already used.
+   % Get the data from the netcdf files
 
-      for i = 1:pinfo.copies
+   vname  = var_names{ivar};
+         
+   [start, count] = GetNCindices(pinfo,'prior',vname);
 
-         imem = pinfo.copyindices(i);
+   ens_colors = get(gca,'ColorOrder');   % trying to cycle through colors
+   ncolors = size(ens_colors,1) - 1;     % last one is black, already used.
 
-         varinfo = nc_getvarinfo(pinfo.prior_file, vname);
+   for i = 1:pinfo.copies
 
-         % Get the Prior 
+      imem = pinfo.copyindices(i);
 
-         for j = 1:length(varinfo.Dimension)
-            switch( lower(varinfo.Dimension{j}))
-               case{'time'}
-                  start(j) = pinfo.prior_times(1) - 1;
-                  count(j) = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
-               case{'copy'}
-                  start(j) = imem - 1;
-                  count(j) = 1;
-               otherwise
-            end
+      varinfo = nc_getvarinfo(pinfo.prior_file, vname);
+
+      %% Get the Prior 
+
+      for j = 1:length(varinfo.Dimension)
+         switch( lower(varinfo.Dimension{j}))
+            case{'time'}
+               start(j) = pinfo.prior_times(1) - 1;
+               count(j) = pinfo.prior_times(2) - pinfo.prior_times(1) + 1;
+            case{'copy'}
+               start(j) = imem - 1;
+               count(j) = 1;
+            otherwise
          end
+      end
 
-         pr_series = nc_varget(pinfo.prior_file,     vname, start, count);
+      pr_series = nc_varget(pinfo.prior_file,     vname, start, count);
 
-         % Get the Posterior 
+      %% Get the Posterior 
 
-         for j = 1:length(varinfo.Dimension)
-            switch( lower(varinfo.Dimension{j}))
-               case{'time'}
-                  start(j) = pinfo.posterior_times(1) - 1;
-                  count(j) = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
-                  break
-               otherwise
-            end
+      for j = 1:length(varinfo.Dimension)
+         switch( lower(varinfo.Dimension{j}))
+            case{'time'}
+               start(j) = pinfo.posterior_times(1) - 1;
+               count(j) = pinfo.posterior_times(2) - pinfo.posterior_times(1) + 1;
+               break
+            otherwise
          end
+      end
 
-         po_series = nc_varget(pinfo.posterior_file, vname, start, count);
+      po_series = nc_varget(pinfo.posterior_file, vname, start, count);
 
-         % Paste Prior/Posterior into one series
+      %% Paste Prior/Posterior into one series
 
-	 ens_member(1,:) = pr_series;
-	 ens_member(2,:) = po_series;
-	 b               = ens_member(:);
+      ens_member(1,:) = pr_series;
+      ens_member(2,:) = po_series;
+      b               = ens_member(:);
 
-         % plot it
+      %% plot it
 
-   	 hold on; 
-	 memcolor = 1 + mod(i-1,ncolors); % cycles through colors [1,6]
+      hold on; 
+      memcolor = 1 + mod(i-1,ncolors); % cycles through colors [1,6]
 
-	 str1 = deblank(metadata(imem,:));
+      str1 = deblank(metadata(imem,:));
 
-         if (strcmp(str1,'ensemble mean') == 1)
+      if (strcmp(str1,'ensemble mean') == 1)
 	    h = plot(pinfo.xax, b,'k-','linewidth',2.0);
-         else
+      else
 	    h = plot(pinfo.xax, b,'linewidth',0.5,'Color',ens_colors(memcolor,:));
-         end
-
-	 [legh, objh, outh, outm] = legend;
-	 nlines = length(outm);
-         outm{nlines+1} = str1;
-	 [legh, objh, outh, outm] = legend([outh; h],outm,0);
       end
 
-      % Plot the true trajectory if it exists; maybe the ens mean
+      [legh, objh, outh, outm] = legend;
+      nlines                   = length(outm);
+      outm{nlines+1}           = str1;
+      [legh, objh, outh, outm] = legend([outh; h],outm,0);
+   end
 
-      if ( exist(pinfo.truth_file,'file') == 2 )
+   %% Plot the true trajectory if it exists
 
-         varinfo = nc_getvarinfo(pinfo.truth_file, vname);
+   if ( exist(pinfo.truth_file,'file') == 2 )
 
-         for j = 1:length(varinfo.Dimension)
-            switch( lower(varinfo.Dimension{j}))
-               case{'time'}
-                  start(j) =  0;
-                  count(j) = -1;
-               case{'copy'}
-                  start(j) = truth.truth_index - 1;
-                  count(j) = 1;
-               otherwise
-            end
+      varinfo = nc_getvarinfo(pinfo.truth_file, vname);
+
+      for j = 1:length(varinfo.Dimension)
+         switch( lower(varinfo.Dimension{j}))
+            case{'time'}
+               start(j) =  0;
+               count(j) = -1;
+            case{'copy'}
+               start(j) = truth.truth_index - 1;
+               count(j) = 1;
+            otherwise
          end
+      end
 
-         true_trajectory = nc_varget(pinfo.truth_file, vname, start, count);
+      true_trajectory = nc_varget(pinfo.truth_file, vname, start, count);
 
-	 h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
+      h = plot(truth.time, true_trajectory, 'k-','linewidth',1.0); hold on;
 
-	 [legh, objh, outh, outm] = legend;
-	 nlines = length(outm);
-         outm{nlines+1} = 'truth';
-	 [legh, objh, outh, outm] = legend([outh; h],outm,0);
-      end
+      [legh, objh, outh, outm] = legend;
+      nlines                   = length(outm);
+      outm{nlines+1}           = 'truth';
+      [legh, objh, outh, outm] = legend([outh; h],outm,0);
+   end
 
-      % annotate
+   % annotate
 
-      ylabel(sprintf('''%s''  (%.3f,%.3fE) level index %d', vname, ...
-                     pinfo.latitude, pinfo.longitude, pinfo.levelindex ))
-      xlabel(sprintf('model time (%d timesteps)', prior.num_times ))
-      title(sprintf('%s Trajectories',pinfo.model), ...  
-                     'interpreter','none','fontweight','bold')
-      legend boxoff
+   h = ylabel(sprintf('''%s''  (%.3f,%.3fE) level index %d', vname, ...
+                  pinfo.latitude, pinfo.longitude, pinfo.levelindex ));
+   set(h,'Interpreter','none')
+   datetick('x','yyyymmmdd HH:MM')
+   title(sprintf('%s Trajectories',pinfo.model), ...  
+                  'Interpreter','none','fontweight','bold')
+   legend boxoff
 
-   end
+end
 

Modified: DART/trunk/matlab/PlotTotalErr.m
===================================================================
--- DART/trunk/matlab/PlotTotalErr.m	2011-03-30 17:38:38 UTC (rev 4833)
+++ DART/trunk/matlab/PlotTotalErr.m	2011-03-31 23:12:11 UTC (rev 4834)
@@ -75,9 +75,8 @@
       plot(times,err, 'b', times,err_spread, 'r');
       legend(string1,string2,0)
       legend boxoff
-      title(sprintf('%s Total Error over all %d variables for %s',...
-                    model, num_vars, pinfo.diagn_file), ...
-            'interpreter','none','fontweight','bold')
+      s1 = sprintf('%s Total Error over all %d variables', model, num_vars);
+      title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
       xlabel(sprintf('model time (%d timesteps)',num_times))
       ylabel('Total Error')
 
@@ -116,9 +115,8 @@
       plot(times,err, 'b', times,err_spread, 'r');
       legend(string1,string2,0)
       legend boxoff
-      title(sprintf('%s Total Error over all %d variables for %s',...
-                    model, num_vars, pinfo.diagn_file), ...
-            'interpreter','none','fontweight','bold')
+      s1 = sprintf('%s Total Error over all %d variables',model,num_vars);
+      title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
       xlabel(sprintf('model time (%d timesteps)',num_times))
       ylabel('Total Error')
 
@@ -167,9 +165,8 @@
       plot(times,err, 'b', times,err_spread, 'r');
       legend(string1,string2,0)
       legend boxoff
-      title(sprintf('%s Total Error over statevars %d to %d for %s',...
-                    model, ind1, indN, pinfo.diagn_file), ...
-            'interpreter','none','fontweight','bold')
+      s1 = sprintf('%s Total Error over statevars %d to %d', model, ind1, indN);
+      title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
       xlabel(sprintf('model time (%d timesteps)',num_times))
       ylabel('Total Error')
 
@@ -198,9 +195,8 @@
       plot(times,err, 'b', times,err_spread, 'r');
       legend(string1,string2,0)
       legend boxoff
-      title(sprintf('%s Total Error over statevars %d to %d for %s',...
-                    model, ind1, indN, pinfo.diagn_file), ...
-            'interpreter','none','fontweight','bold')
+      s1 = sprintf('%s Total Error over statevars %d to %d', model, ind1, indN);
+      title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
       xlabel(sprintf('model time (%d timesteps)',num_times))
       ylabel('Total Error')
 
@@ -427,7 +423,7 @@
  
       varunits = nc_attget(pinfo.truth_file, pinfo.vars{ivar}, 'units');
 
-      s1 = sprintf('%s %s Ensemble Mean for %s', model,pinfo.vars{ivar},pinfo.diagn_file);
+      s1 = sprintf('%s %s Ensemble Mean', model,pinfo.vars{ivar});
 
       switch lower(pinfo.vars{ivar})
          case {'ps'}
@@ -451,7 +447,7 @@
       grid on;
       xlabel(sprintf('time (%s) %d timesteps',timeunits,num_times))
       ylabel(sprintf('global-area-weighted distance (%s)',varunits))
-      title(s1,'interpreter','none','fontweight','bold')
+      title({s1,pinfo.diagn_file},'interpreter','none','fontweight','bold')
 
 end
 
@@ -584,7 +580,7 @@
    figure(ivar); clf;
       varunits = nc_attget(pinfo.truth_file,pinfo.vars{ivar},'units');
 
-      s1 = sprintf('%s %s Ensemble Mean for %s', model,pinfo.vars{ivar},pinfo.diagn_file);

@@ Diff output truncated at 40000 characters. @@


More information about the Dart-dev mailing list