[Dart-dev] [8536] DART/trunk/matlab: POP is now supported for all the meaningful state-space diagnostics.
nancy at ucar.edu
nancy at ucar.edu
Wed Sep 9 11:46:42 MDT 2015
Revision: 8536
Author: thoar
Date: 2015-09-09 11:46:42 -0600 (Wed, 09 Sep 2015)
Log Message:
-----------
POP is now supported for all the meaningful state-space diagnostics.
plot_total_err is a conundrum ... so I punted. I don't have a good
metric for 'total error' when the number of grid cells changes with depth,
when each variable has a different dynamic range, when each level ....
you get the picture ...
Modified Paths:
--------------
DART/trunk/matlab/CheckModel.m
DART/trunk/matlab/CheckModelCompatibility.m
DART/trunk/matlab/GetNCindices.m
DART/trunk/matlab/PlotBins.m
DART/trunk/matlab/PlotCorrel.m
DART/trunk/matlab/PlotEnsErrSpread.m
DART/trunk/matlab/PlotEnsMeanTimeSeries.m
DART/trunk/matlab/PlotEnsTimeSeries.m
DART/trunk/matlab/PlotJeffCorrel.m
DART/trunk/matlab/PlotPhaseSpace.m
DART/trunk/matlab/PlotSawtooth.m
DART/trunk/matlab/PlotTotalErr.m
DART/trunk/matlab/PlotVarVarCorrel.m
DART/trunk/matlab/plot_bins.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
DART/trunk/matlab/rank_hist.m
Added Paths:
-----------
DART/trunk/matlab/GetPOPInfo.m
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModel.m
===================================================================
--- DART/trunk/matlab/CheckModel.m 2015-09-09 17:44:26 UTC (rev 8535)
+++ DART/trunk/matlab/CheckModel.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -218,14 +218,12 @@
'num_ens_members',num_copies, ...
'time_series_length',num_times, ...
'num_ens_members',ens_size, ...
- 'ensemble_indices',ens_indices, ...
- 'min_ens_mem',ens_indices(1), ...
- 'max_ens_mem',ens_indices(ens_size) );
+ 'ensemble_indices',ens_indices);
vars.vars = varnames;
vars.fname = fname;
- case {'cam','tiegcm','fms_bgrid','pe2lyr','mitgcm_ocean','pbl_1d','mpas_atm','sqg'}
+ case {'cam','tiegcm','fms_bgrid','pe2lyr','mitgcm_ocean','pbl_1d','mpas_atm','sqg','pop'}
varnames = get_DARTvars(fname);
num_vars = length(varnames);
@@ -237,8 +235,6 @@
'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/trunk/matlab/CheckModelCompatibility.m
===================================================================
--- DART/trunk/matlab/CheckModelCompatibility.m 2015-09-09 17:44:26 UTC (rev 8535)
+++ DART/trunk/matlab/CheckModelCompatibility.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -276,6 +276,13 @@
x = 3;
y = [dnum_lons dnum_lats dnum_lvls];
+ case {'pop'}
+ dnum_lons = dim_length(fname,'i');
+ dnum_lats = dim_length(fname,'j');
+ dnum_lvls = dim_length(fname,'k');
+ x = 3;
+ y = [dnum_lons dnum_lats dnum_lvls];
+
case {'mpas_atm'}
dnum_cells = dim_length(fname,'nCells');
@@ -290,11 +297,11 @@
y = [dnum_X dnum_Y];
case {'simple_advection'}
- y = dim_length(fname,'loc1d');
+ y = dim_length(fname,'loc1d');
x = 1;
otherwise
- y = dim_length(fname,'StateVariable');
+ y = dim_length(fname,'StateVariable');
x = 1;
end
Modified: DART/trunk/matlab/GetNCindices.m
===================================================================
--- DART/trunk/matlab/GetNCindices.m 2015-09-09 17:44:26 UTC (rev 8535)
+++ DART/trunk/matlab/GetNCindices.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -304,13 +304,13 @@
case 't'
start(i) = time1;
count(i) = timeN;
- case {'lev','z'}
+ case {'lev','z','k'}
start(i) = level1;
count(i) = levelN;
- case {'lat','y','iy'}
+ case {'lat','y','iy','j'}
start(i) = lat1;
count(i) = latN;
- case {'lon','x','ix'}
+ case {'lon','x','ix','i'}
start(i) = lon1;
count(i) = lonN;
case 'pft'
Added: DART/trunk/matlab/GetPOPInfo.m
===================================================================
--- DART/trunk/matlab/GetPOPInfo.m (rev 0)
+++ DART/trunk/matlab/GetPOPInfo.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -0,0 +1,439 @@
+function pinfo = GetPOPInfo(pinfo_in,fname,routine)
+%% GetPOPInfo prepares a structure of information needed by the subsequent "routine"
+% The information is gathered via rudimentary "input" routines.
+%
+% pinfo = GetPOPInfo(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 - 2013 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
+%
+% DART $Id$
+
+if ( exist(fname,'file') ~= 2 ), error('%s does not exist.',fname); end
+
+pinfo = pinfo_in;
+
+if strcmpi(pinfo.model,'POP') ~= 1
+ error('Not so fast, this is not a POP model.')
+end
+
+%% Get the required information.
+
+global copy ULON ULAT TLON TLAT ZG ZC KMU KMT
+
+varexist(fname, {'copy','ULON','ULAT','TLON','TLAT','ZG','ZC','KMT','KMU'});
+
+copy = nc_varget(fname, 'copy');
+ULON = nc_varget(fname, 'ULON');
+ULAT = nc_varget(fname, 'ULAT');
+TLON = nc_varget(fname, 'TLON');
+TLAT = nc_varget(fname, 'TLAT');
+ZG = nc_varget(fname, 'ZG');
+ZC = nc_varget(fname, 'ZC');
+KMU = nc_varget(fname, 'KMU');
+KMT = nc_varget(fname, 'KMT');
+
+depthunits = nc_attget(fname,'ZG','units');
+
+
+%% The POP model has different number of possible levels for each location.
+% Consequently, it is required to query the lat/lon before the number of
+% levels.
+
+switch lower(deblank(routine))
+
+ case {'plotbins','plotenserrspread','plotensmeantimeseries','plotenstimeseries'}
+
+ pgvar = GetVarString(pinfo_in.vars); % set prognostic variable
+ [lat, lon, latind, lonind] = GetLatLon(fname, pgvar);
+ [level, lvlind] = GetLevel(fname,pgvar,lonind,latind);
+
+ pinfo.fname = fname;
+ 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 = GetVarString(pinfo_in.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
+ [base_lat, base_lon, base_latind, base_lonind] = GetLatLon(fname, base_var);
+ [base_lvl, base_lvlind] = GetLevel(fname, base_var, base_lonind, base_latind);
+
+ disp('Getting information for the ''comparison'' variable.')
+ comp_var = GetVarString(pinfo_in.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel(fname, base_var, base_lonind, base_latind, base_lvl);
+
+ pinfo.fname = fname;
+ 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.depthunits = depthunits;
+
+ case 'plotvarvarcorrel'
+
+ disp('Getting information for the ''base'' variable.')
+ base_var = GetVarString(pinfo_in.vars);
+ [base_time, base_tmeind] = GetTime(pinfo.time);
+ [base_lat, base_lon, base_latind, base_lonind] = GetLatLon(fname, base_var);
+ [base_lvl , base_lvlind] = GetLevel(fname, base_var, base_lonind, base_latind);
+
+ disp('Getting information for the ''comparison'' variable.')
+ comp_var = GetVarString(pinfo_in.vars, base_var);
+ [comp_lat, comp_lon, comp_latind, comp_lonind] = GetLatLon(fname, comp_var, ...
+ base_latind, base_lonind);
+ [comp_lvl, comp_lvlind] = GetLevel(fname, comp_var, comp_lonind, comp_latind, base_lvl);
+
+ pinfo.fname = fname;
+ 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;
+ pinfo.depthunits = depthunits;
+
+ case 'plotsawtooth'
+
+ pgvar = GetVarString(pinfo_in.vars);
+ [lat, lon, latind, lonind] = GetLatLon(pinfo_in.prior_file, pgvar);
+ [level, lvlind] = GetLevel(pinfo_in.prior_file, pgvar, lonind, latind);
+ [copyindices, copymetadata]= SetCopyID2(pinfo_in.prior_file);
+ copy = length(copyindices);
+
+ pinfo.truth_file = pinfo_in.truth_file;
+ pinfo.prior_file = pinfo_in.prior_file;
+ pinfo.posterior_file = pinfo_in.posterior_file;
+ pinfo.var_names = pgvar;
+ pinfo.level = level;
+ pinfo.levelindex = lvlind;
+ pinfo.latitude = lat;
+ pinfo.latindex = latind;
+ pinfo.longitude = lon;
+ pinfo.lonindex = lonind;
+ pinfo.copies = copy;
+ pinfo.copyindices = copyindices;
+ pinfo.copymetadata = copymetadata;
+
+ % XXX 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_lat, var1_lon, var1_latind, var1_lonind] = GetLatLon(fname, var1);
+ [var1_lvl, var1_lvlind] = GetLevel(fname,var1,var1_lonind,var1_latind);
+
+ disp('Getting information for the ''Y'' variable.')
+ var2 = GetVarString(pinfo_in.vars, var1 );
+ [var2_lat, var2_lon, var2_latind, var2_lonind] = GetLatLon(fname, var2, var1_latind, var1_lonind);
+ [var2_lvl, var2_lvlind] = GetLevel(fname,var2,var2_lonind,var2_latind,var1_lvl);
+
+ disp('Getting information for the ''Z'' variable.')
+ var3 = GetVarString(pinfo_in.vars, var1 );
+ [var3_lat, var3_lon, var3_latind, var3_lonind] = GetLatLon(fname, var3, var1_latind, var1_lonind);
+ [var3_lvl, var3_lvlind] = GetLevel(fname,var3,var3_lonind,var3_latind,var1_lvl);
+
+ % query for ensemble member string
+ metadata = nc_varget(fname,'CopyMetaData');
+ [N,M] = size(metadata);
+ if M == 1
+ cell_array{1} = metadata';
+ else
+ cell_array = mat2cell(metadata, ones(1,N), M);
+ end
+ 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');
+ if isempty(s1), ltype = 'k-'; else ltype = s1; end
+
+ pinfo.fname = fname;
+ 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;
+
+ otherwise
+
+end
+
+
+
+function pgvar = GetVarString(prognostic_vars, defvar)
+%----------------------------------------------------------------------
+pgvar = prognostic_vars{1};
+if (nargin == 2), pgvar = defvar; 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 = deblank(varstring); end
+
+
+
+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 = str2num(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(fname, pgvar, lonind, latind, deflevel)
+%----------------------------------------------------------------------
+%
+global ZC ZG KMU KMT
+
+if (nargin == 5), level = deflevel; else level = []; end
+
+varinfo = nc_getvarinfo(fname,pgvar);
+leveldim = find(strncmp('k',varinfo.Dimension,1) > 0);
+
+if (isempty(leveldim))
+
+ fprintf('''%s'' has no vertical level, imply using 1.\n',pgvar)
+ level = 1;
+ lvlind = 1;
+
+else
+
+ levelvar = varinfo.Dimension{leveldim};
+ dinfo = nc_getdiminfo(fname,levelvar);
+
+ levels = 1:dinfo.Length;
+
+ % must find the lowest level at this gridcell for this variable
+ velcomp = regexp(varinfo.Name,'VEL','ONCE');
+
+ if ( isempty(velcomp) )
+ levbottom = KMT(lonind,latind);
+ levels = ZC(1:levbottom);
+ else
+ levbottom = KMU(lonind,latind);
+ levels = ZG(1:levbottom);
+ end
+
+ if (isempty(level)), level = levels(1); end
+
+ fprintf('Default level is %0.2f, if this is OK, <cr>;\n',level)
+ fprintf('If not, enter a level between %0.2f and %0.2f, inclusive ...\n', ...
+ min(levels),max(levels))
+ varstring = input('we''ll use the closest (no syntax required)\n','s');
+
+ if ~isempty(varstring), level = str2num(varstring); end
+
+ d = abs(level - levels); % crude distance
+ ind = find(min(d) == d); % multiple minima possible
+ lvlind = ind(1); % use the first one
+ level = levels(lvlind);
+
+end
+
+
+
+function [lat, lon, latind, lonind] = GetLatLon(fname, pgvar, deflat, deflon)
+%----------------------------------------------------------------------
+% The POP domain is not a simple 1D lat 1D lon domain, so it makes more
+% sense to query about lat & lon indices rather than values.
+
+global ULON ULAT TLON TLAT
+
+varinfo = nc_getvarinfo(fname, pgvar);
+
+londim = find(strncmp('i',varinfo.Dimension,1) > 0);
+latdim = find(strncmp('j',varinfo.Dimension,1) > 0);
+
+if (isempty(latdim) || isempty(londim))
+ error('''%s'' does not have both lats and lons.\n',pgvar)
+end
+
+nlat = varinfo.Size(latdim);
+nlon = varinfo.Size(londim);
+
+% both U,V velocity components are on the ULON, ULAT coordinates.
+% The T,S components use TLON, TLAT
+
+velcomp = regexp(varinfo.Name,'VEL','ONCE');
+
+if ( isempty(velcomp) )
+ latmat = TLAT;
+ lonmat = TLON;
+else
+ latmat = ULAT;
+ lonmat = ULON;
+end
+
+% Determine a sensible default.
+if (nargin > 2)
+ latind = deflat;
+ lonind = deflon;
+else % in the middle of the domain.
+ latind = round(nlat/2);
+ lonind = round(nlon/2);
+end
+lat = latmat(latind,lonind);
+lon = lonmat(latind,lonind);
+
+% Ask if they like the default.
+
+fprintf('There are %d latitudes and %d longitudes for %s\n',nlat,nlon,pgvar)
+fprintf('Default lat/lon indices are %d,%d ... lat/lon(%f,%f)\n',latind,lonind,lat,lon)
+fprintf('If this is OK, <cr>; if not, enter grid indices lat/lon pair ...\n')
+varstring = input('we''ll use the closest (no syntax required)\n','s');
+
+if ~isempty(varstring)
+ nums = str2num(varstring);
+ if (length(nums) ~= 2)
+ error('Did not get two indices for the lat lon pair.')
+ end
+ latind = nums(1);
+ lonind = nums(2);
+end
+
+% Find the closest lat/lon location to the requested one.
+% not-so-Simple search over a matrix of lats/lons.
+
+% londiff = lonmat - lon;
+% latdiff = latmat - lat;
+% dist = sqrt(londiff.^2 + latdiff.^2);
+% [latind, lonind] = find(min(dist(:)) == dist);
+% lat = latmat(latind,lonind);
+% lon = lonmat(latind,lonind);
+
+% Much simpler search to find closest lat index.
+latinds = 1:nlat;
+loninds = 1:nlon;
+
+d = abs(latind - latinds); % crude distance
+ind = find(min(d) == d); % multiple minima possible
+latind = ind(1); % use the first one
+
+d = abs(lonind - loninds); % crude distance
+ind = find(min(d) == d); % multiple minima possible
+lonind = ind(1); % use the first one
+
+lat = latmat(latind,lonind);
+lon = lonmat(latind,lonind);
+
+
+
+function varexist(filename, varnames)
+%% We already know the file exists by this point.
+% Lets check to make sure that file contains all needed variables.
+% Fatally die if they do not exist.
+
+nvars = length(varnames);
+gotone = ones(1,nvars);
+
+for i = 1:nvars
+ gotone(i) = nc_isvar(filename,varnames{i});
+ if ( ~ gotone(i) )
+ fprintf('\n%s is not a variable in %s\n',varnames{i},filename)
+ end
+end
+
+if ~ all(gotone)
+ error('missing required variable ... exiting')
+end
+
+
+
+% <next few lines under version control, do not edit>
+% $URL$
+% $Revision$
+% $Date$
+
Property changes on: DART/trunk/matlab/GetPOPInfo.m
___________________________________________________________________
Added: svn:mime-type
+ text/plain
Added: svn:keywords
+ Date Rev Author HeadURL Id
Added: svn:eol-style
+ native
Modified: DART/trunk/matlab/PlotBins.m
===================================================================
--- DART/trunk/matlab/PlotBins.m 2015-09-09 17:44:26 UTC (rev 8535)
+++ DART/trunk/matlab/PlotBins.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -37,6 +37,8 @@
if isempty(pinfo.num_ens_members)
error('no ensemble members in %s, cannot create rank histogram.',pinfo.diagn_file)
+elseif (pinfo.num_ens_members < 2)
+ error('not ensemble members in %s, cannot create rank histogram.',pinfo.diagn_file)
end
% Get the state for the truth
@@ -137,7 +139,7 @@
axis tight
end
- case {'fms_bgrid','pe2lyr','mitgcm_ocean','cam','wrf','mpas_atm','mpas_ocn','sqg'}
+ case {'fms_bgrid','pe2lyr','mitgcm_ocean','cam','wrf','mpas_atm','mpas_ocn','sqg','pop'}
% It is intended that all 3D models have all the required information
% set in the corresponding Get<model>Info.m script.
Modified: DART/trunk/matlab/PlotCorrel.m
===================================================================
--- DART/trunk/matlab/PlotCorrel.m 2015-09-09 17:44:26 UTC (rev 8535)
+++ DART/trunk/matlab/PlotCorrel.m 2015-09-09 17:46:42 UTC (rev 8536)
@@ -34,500 +34,565 @@
switch(lower(pinfo.model))
- case {'9var','lorenz_63','lorenz_84','lorenz_96','lorenz_96_2scale', ...
- 'lorenz_04','forced_lorenz_96','ikeda','simple_advection'}
+ case {'9var','lorenz_63','lorenz_84','lorenz_96','lorenz_96_2scale', ...
+ 'lorenz_04','forced_lorenz_96','ikeda','simple_advection'}
- % The Base Variable Index must be a valid state variable
- if ( pinfo.base_var_index > pinfo.num_state_vars )
- fprintf('%s only has %d state variables\n', pinfo.fname, pinfo.num_state_vars)
- error('you wanted variable # %d ', pinfo.base_var_index)
- end
+ % The Base Variable Index must be a valid state variable
+ if ( pinfo.base_var_index > pinfo.num_state_vars )
+ fprintf('%s only has %d state variables\n', pinfo.fname, pinfo.num_state_vars)
+ error('you wanted variable # %d ', pinfo.base_var_index)
+ end
- % The Time must be within range also.
- if ( pinfo.base_time > pinfo.time_series_length )
- fprintf('%s only has %d output times\n', pinfo.fname, pinfo.time_series_length)
- error('you wanted time # %d ', pinfo.base_time)
- end
+ % The Time must be within range also.
+ if ( pinfo.base_time > pinfo.time_series_length )
+ fprintf('%s only has %d output times\n', pinfo.fname, pinfo.time_series_length)
+ error('you wanted time # %d ', pinfo.base_time)
+ end
- % Get 'standard' ensemble series
- base = get_hyperslab('fname',pinfo.fname, 'varname',pinfo.base_var, ...
- 'copy1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
- 'stateindex',pinfo.base_var_index);
+ % Get 'standard' ensemble series
+ base = get_hyperslab('fname',pinfo.fname, 'varname',pinfo.base_var, ...
+ '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, ...
- 'copy1',pinfo.ensemble_indices(1),'copycount',pinfo.num_ens_members, ...
- 'state1',pinfo.min_state_var,'statecount',pinfo.num_state_vars);
+ % Get (potentially large) state.
+ state_var = get_hyperslab('fname',pinfo.fname, 'varname',pinfo.base_var, ...
+ '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 ...
- correl = zeros(pinfo.num_state_vars,pinfo.time_series_length);
+ % It is efficient to preallocate correl storage ...
+ correl = zeros(pinfo.num_state_vars,pinfo.time_series_length);
- % Need to loop through all variables in the ensemble
- for i = 1:pinfo.num_state_vars,
- correl(i, :) = ens_correl(base, pinfo.base_time, state_var(:,:,i));
- end
+ % Need to loop through all variables in the ensemble
+ for i = 1:pinfo.num_state_vars,
+ correl(i, :) = ens_correl(base, pinfo.base_time, state_var(:,:,i));
+ end
- % Now for the plotting part ...
- disp('Please be patient ... this usually takes a bit ...')
- clf;
+ % Now for the plotting part ...
+ disp('Please be patient ... this usually takes a bit ...')
+ clf;
- [~,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
- plot(pinfo.time(pinfo.base_time), pinfo.base_var_index,'kh','MarkerSize',12,'MarkerFaceColor','k')
+ [~,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
+ plot(pinfo.time(pinfo.base_time), pinfo.base_var_index,'kh','MarkerSize',12,'MarkerFaceColor','k')
- s1 = sprintf('%s Correlation of variable %s index %d, timestep = %d', ...
- pinfo.model, pinfo.base_var, pinfo.base_var_index, pinfo.base_time);
- s2 = sprintf('against all variables, all times, %d ensemble members', ...
- size(state_var,2));
- title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
- xlabel(sprintf('model "days" (%d timesteps)',pinfo.time_series_length))
- ylabel('state variable (index)')
- set(gca,'YTick',1:pinfo.num_state_vars)
- colorbar
+ s1 = sprintf('%s Correlation of variable %s index %d, timestep = %d', ...
+ pinfo.model, pinfo.base_var, pinfo.base_var_index, pinfo.base_time);
+ s2 = sprintf('against all variables, all times, %d ensemble members', ...
+ size(state_var,2));
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
+ xlabel(sprintf('model "days" (%d timesteps)',pinfo.time_series_length))
+ ylabel('state variable (index)')
+ set(gca,'YTick',1:pinfo.num_state_vars)
+ 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
+ % We are going to correlate one var/time/lvl/lat/lon with
+ % all other lats/lons for a var/time/lvl
- clf;
+ clf;
- switch lower(pinfo.comp_var)
- case {'ps','t'}
- lats = nc_varget(pinfo.fname,'TmpJ'); ny = length(lats);
- lons = nc_varget(pinfo.fname,'TmpI'); nx = length(lons);
- latunits = nc_attget(pinfo.fname,'TmpJ','units');
- lonunits = nc_attget(pinfo.fname,'TmpI','units');
- otherwise
- lats = nc_varget(pinfo.fname,'VelJ'); ny = length(lats);
- lons = nc_varget(pinfo.fname,'VelI'); nx = length(lons);
- latunits = nc_attget(pinfo.fname,'VelJ','units');
- lonunits = nc_attget(pinfo.fname,'VelI','units');
- end
+ switch lower(pinfo.comp_var)
+ case {'ps','t'}
+ lats = nc_varget(pinfo.fname,'TmpJ'); ny = length(lats);
+ lons = nc_varget(pinfo.fname,'TmpI'); nx = length(lons);
+ latunits = nc_attget(pinfo.fname,'TmpJ','units');
+ lonunits = nc_attget(pinfo.fname,'TmpI','units');
+ otherwise
+ lats = nc_varget(pinfo.fname,'VelJ'); ny = length(lats);
+ lons = nc_varget(pinfo.fname,'VelI'); nx = length(lons);
+ latunits = nc_attget(pinfo.fname,'VelJ','units');
+ lonunits = nc_attget(pinfo.fname,'VelI','units');
+ end
- nxny = nx*ny;
+ nxny = nx*ny;
- base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
- '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 );
+ base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
+ '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, ...
- 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
- 'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
+ bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
+ '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]);
- corr = zeros(nxny,1);
+ comp_ens = reshape(bob,[pinfo.num_ens_members, nxny]);
+ corr = zeros(nxny,1);
- for i = 1:nxny,
- x = corrcoef(base_mem, comp_ens(:, i));
- corr(i) = x(1, 2);
- end
+ for i = 1:nxny,
+ x = corrcoef(base_mem, comp_ens(:, i));
+ corr(i) = x(1, 2);
+ end
- correl = reshape(corr,[ny nx]);
+ correl = reshape(corr,[ny nx]);
- [cs, h] = contour(lons,lats,correl,contourlevels);
- clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
- set(gca,'Clim',[-1 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', ...
- pinfo.model, pinfo.base_var, pinfo.base_lvl, ...
- pinfo.base_lat, pinfo.base_lon, pinfo.base_time);
+ [cs, h] = contour(lons,lats,correl,contourlevels);
+ clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
+ set(gca,'Clim',[-1 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', ...
+ pinfo.model, pinfo.base_var, pinfo.base_lvl, ...
+ 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, pinfo.num_ens_members);
- title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
- xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
- ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
- continents;
- axis image
- colorbar;
+ s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
+ pinfo.comp_var, pinfo.comp_lvl, pinfo.num_ens_members);
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
+ xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
+ ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
+ continents;
+ 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
+ % We are going to correlate one var/time/lvl/lat/lon with
+ % all other lats/lons for a var/time/lvl
- clf;
+ clf;
- % Get the plotting lat/lon for the comparison variable.
- % This is the variable that has a spatial extent.
+ % Get the plotting lat/lon for the comparison variable.
+ % This is the variable that has a spatial extent.
- varinfo = nc_getvarinfo(pinfo.fname, pinfo.comp_var);
- latdim = find(strncmp('south_north',varinfo.Dimension,length('south_north')) > 0);
- londim = find(strncmp( 'west_east',varinfo.Dimension,length( 'west_east')) > 0);
- ny = varinfo.Size(latdim);
- nx = varinfo.Size(londim);
- nxny = nx*ny;
+ varinfo = nc_getvarinfo(pinfo.fname, pinfo.comp_var);
+ latdim = find(strncmp('south_north',varinfo.Dimension,length('south_north')) > 0);
+ londim = find(strncmp( 'west_east',varinfo.Dimension,length( 'west_east')) > 0);
+ ny = varinfo.Size(latdim);
+ nx = varinfo.Size(londim);
+ nxny = nx*ny;
- % Each of the WRF variables has a 'coordinate' attribute signifying which
- % of the 6 possible lat/lon variables is appropriate.
+ % Each of the WRF variables has a 'coordinate' attribute signifying which
+ % of the 6 possible lat/lon variables is appropriate.
- coordinates{1} = sscanf(nc_attget(pinfo.fname,pinfo.comp_var,'coordinates'),'%s %*s');
- coordinates{2} = sscanf(nc_attget(pinfo.fname,pinfo.comp_var,'coordinates'),'%*s %s');
- latcoord = find(strncmp('XLAT',coordinates,length('XLAT')) > 0);
- loncoord = find(strncmp('XLON',coordinates,length('XLON')) > 0);
- latmat = nc_varget(pinfo.fname,coordinates{latcoord});
- lonmat = nc_varget(pinfo.fname,coordinates{loncoord});
- latunits = nc_attget(pinfo.fname,coordinates{latcoord},'units');
- lonunits = nc_attget(pinfo.fname,coordinates{latcoord},'units');
+ coordinates{1} = sscanf(nc_attget(pinfo.fname,pinfo.comp_var,'coordinates'),'%s %*s');
+ coordinates{2} = sscanf(nc_attget(pinfo.fname,pinfo.comp_var,'coordinates'),'%*s %s');
+ latcoord = find(strncmp('XLAT',coordinates,length('XLAT')) > 0);
+ loncoord = find(strncmp('XLON',coordinates,length('XLON')) > 0);
+ latmat = nc_varget(pinfo.fname,coordinates{latcoord});
+ lonmat = nc_varget(pinfo.fname,coordinates{loncoord});
+ latunits = nc_attget(pinfo.fname,coordinates{latcoord},'units');
+ lonunits = nc_attget(pinfo.fname,coordinates{latcoord},'units');
- inds = (lonmat < 0); % Convert to 0,360 to minimize dateline probs.
- lonmat(inds) = lonmat(inds) + 360.0;
- if (pinfo.base_lon < 0), pinfo.base_lon = pinfo.base_lon + 360.0; end
+ inds = (lonmat < 0); % Convert to 0,360 to minimize dateline probs.
+ lonmat(inds) = lonmat(inds) + 360.0;
+ if (pinfo.base_lon < 0), pinfo.base_lon = pinfo.base_lon + 360.0; end
- % Get the actual goods ... and perform the correlation
+ % Get the actual goods ... and perform the correlation
- base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
- '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 );
+ base_mem = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.base_var, ...
+ '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 );
- if (std(base_mem) == 0.0)
- warning('%s at level %d lat %d lon %d time %s is a constant\n',pinfo.base_var,...
- pinfo.base_lvlind,pinfo.base_latind,pinfo.base_lonind,datestr(pinfo.base_time))
- error('Cannot calculate correlation coefficient with a constant.')
- end
+ if (std(base_mem) == 0.0)
+ warning('%s at level %d lat %d lon %d time %s is a constant\n',pinfo.base_var,...
+ pinfo.base_lvlind,pinfo.base_latind,pinfo.base_lonind,datestr(pinfo.base_time))
+ error('Cannot calculate correlation coefficient with a constant.')
+ end
- bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
- 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
- 'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
+ bob = get_hyperslab('fname', pinfo.fname, 'varname', pinfo.comp_var, ...
+ 'copy1', pinfo.ensemble_indices(1), 'copycount',pinfo.num_ens_members, ...
+ 'timeindex', pinfo.base_tmeind, 'levelindex', pinfo.comp_lvlind);
- if (std(bob(:)) == 0.0)
- warning('%s at level %d time %s is a constant\n',pinfo.comp_var,...
- pinfo.comp_lvlind, datestr(pinfo.base_time))
- error('Cannot calculate correlation coefficient with a constant.')
- end
+ if (std(bob(:)) == 0.0)
+ warning('%s at level %d time %s is a constant\n',pinfo.comp_var,...
+ pinfo.comp_lvlind, datestr(pinfo.base_time))
+ error('Cannot calculate correlation coefficient with a constant.')
+ end
- comp_ens = reshape(bob,[pinfo.num_ens_members, nxny]);
- corr = zeros(nxny,1);
+ comp_ens = reshape(bob,[pinfo.num_ens_members, nxny]);
+ corr = zeros(nxny,1);
- % Really should check to see if each comp_ens is a constant value as
- % well - this is slow enough already.
+ % Really should check to see if each comp_ens is a constant value as
+ % well - this is slow enough already.
- fprintf('Performing correlations at %d locations ...\n',nxny)
- for i = 1:nxny,
- x = corrcoef(base_mem, comp_ens(:, i));
- corr(i) = x(1, 2);
- end
+ fprintf('Performing correlations at %d locations ...\n',nxny)
+ for i = 1:nxny,
+ x = corrcoef(base_mem, comp_ens(:, i));
+ corr(i) = x(1, 2);
+ end
- correl = reshape(corr,[ny nx]);
+ correl = reshape(corr,[ny nx]);
- % Plot it up ...
+ % Plot it up ...
- [cs,h] = contour(lonmat,latmat,correl,contourlevels);
- clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
- set(gca,'Clim',[-1 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 = %s', ...
- pinfo.model, pinfo.base_var, pinfo.base_lvl, ...
- pinfo.base_lat, pinfo.base_lon, datestr(pinfo.base_time));
+ [cs,h] = contour(lonmat,latmat,correl,contourlevels);
+ clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
+ set(gca,'Clim',[-1 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 = %s', ...
+ pinfo.model, pinfo.base_var, pinfo.base_lvl, ...
+ pinfo.base_lat, pinfo.base_lon, datestr(pinfo.base_time));
- s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
- pinfo.comp_var, pinfo.comp_lvl, pinfo.num_ens_members);
- title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
- xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
- ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
- continents;
- axis image
- colorbar;
+ s2 = sprintf('against ''%s'', entire level %d, same time, %d ensemble members', ...
+ pinfo.comp_var, pinfo.comp_lvl, pinfo.num_ens_members);
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
+ xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
+ ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
+ continents;
+ 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
+ % We are going to correlate one var/time/lvl/lat/lon with
+ % all other lats/lons for a var/time/lvl
- clf;
+ clf;
- switch lower(pinfo.comp_var)
- case {'u'}
- lats = nc_varget(pinfo.fname,'YC'); ny = length(lats);
- lons = nc_varget(pinfo.fname,'XG'); nx = length(lons);
- latunits = nc_attget(pinfo.fname,'YC','units');
- lonunits = nc_attget(pinfo.fname,'XG','units');
- case {'v'}
- lats = nc_varget(pinfo.fname,'YG'); ny = length(lats);
- lons = nc_varget(pinfo.fname,'XC'); nx = length(lons);
- latunits = nc_attget(pinfo.fname,'YG','units');
- lonunits = nc_attget(pinfo.fname,'XC','units');
- otherwise
- lats = nc_varget(pinfo.fname,'YC'); ny = length(lats);
- lons = nc_varget(pinfo.fname,'XC'); nx = length(lons);
- latunits = nc_attget(pinfo.fname,'YC','units');
- lonunits = nc_attget(pinfo.fname,'XC','units');
- end
+ switch lower(pinfo.comp_var)
+ case {'u'}
+ lats = nc_varget(pinfo.fname,'YC'); ny = length(lats);
+ lons = nc_varget(pinfo.fname,'XG'); nx = length(lons);
+ latunits = nc_attget(pinfo.fname,'YC','units');
+ lonunits = nc_attget(pinfo.fname,'XG','units');
+ case {'v'}
+ lats = nc_varget(pinfo.fname,'YG'); ny = length(lats);
+ lons = nc_varget(pinfo.fname,'XC'); nx = length(lons);
+ latunits = nc_attget(pinfo.fname,'YG','units');
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list