[Dart-dev] [4831] DART/trunk/matlab: The correlation routine works, as long as you don 't pick a variable
nancy at ucar.edu
nancy at ucar.edu
Tue Mar 29 14:43:57 MDT 2011
Revision: 4831
Author: thoar
Date: 2011-03-29 14:43:57 -0600 (Tue, 29 Mar 2011)
Log Message:
-----------
The correlation routine works, as long as you don't pick a variable
that is zero for all the ensemble members. Must guard against that.
Modified Paths:
--------------
DART/trunk/matlab/CheckModel.m
DART/trunk/matlab/GetWRFInfo.m
DART/trunk/matlab/PlotCorrel.m
DART/trunk/matlab/plot_correl.m
-------------- next part --------------
Modified: DART/trunk/matlab/CheckModel.m
===================================================================
--- DART/trunk/matlab/CheckModel.m 2011-03-28 22:31:52 UTC (rev 4830)
+++ DART/trunk/matlab/CheckModel.m 2011-03-29 20:43:57 UTC (rev 4831)
@@ -52,6 +52,8 @@
'max_ens_mem',max(copy), ...
'def_state_vars',def_state_vars);
+ vars.fname = fname;
+
case {'lorenz_96', 'lorenz_04'}
num_vars = dim_length(fname,'StateVariable'); % determine # of state varbls
@@ -73,6 +75,8 @@
'max_ens_mem',max(copy), ...
'def_state_vars',def_state_vars);
+ vars.fname = fname;
+
case 'forced_lorenz_96'
% This model has the state variables replicated, so there is a difference
@@ -106,6 +110,8 @@
'def_state_vars',def_state_vars, ...
'def_force_vars',def_force_vars);
+ vars.fname = fname;
+
case 'lorenz_96_2scale'
num_X = dim_length(fname,'Xdim'); % # of X variables
@@ -131,6 +137,8 @@
'min_ens_mem', min(copy), 'max_ens_mem', max(copy), ...
'def_state_vars',def_X_inds);
+ vars.fname = fname;
+
case 'simple_advection'
num_locs = dim_length(fname,'loc1d'); % # of X variables
@@ -157,8 +165,10 @@
'max_state_var' ,num_locs, ...
'def_state_vars' ,def_inds, ...
'num_vars' ,length(varnames));
- vars.vars = varnames;
+ vars.vars = varnames;
+ vars.fname = fname;
+
case 'fms_bgrid'
% A more robust way would be to use the netcdf low-level ops:
@@ -184,7 +194,8 @@
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy));
- vars.vars = varnames;
+ vars.vars = varnames;
+ vars.fname = fname;
case 'cam'
@@ -204,9 +215,10 @@
'time_series_length',num_times, ...
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy) );
- % 'max_ens_mem',max(copy), ...
- % 'varnames',varnames);
+ vars.vars = varnames;
+ vars.fname = fname;
+
case 'pbl_1d'
% A more robust way would be to use the netcdf low-level ops:
@@ -228,6 +240,8 @@
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy));
+ vars.fname = fname;
+
case 'pe2lyr'
% Since this is a 3D model, only the most rudimentary information
@@ -245,7 +259,8 @@
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy) );
- vars.vars = varnames;
+ vars.vars = varnames;
+ vars.fname = fname;
case 'mitgcm_ocean'
@@ -266,7 +281,8 @@
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy) );
- vars.vars = varnames;
+ vars.vars = varnames;
+ vars.fname = fname;
case 'wrf'
@@ -289,7 +305,8 @@
'min_ens_mem',min(copy), ...
'max_ens_mem',max(copy));
- vars.vars = varnames;
+ vars.vars = varnames;
+ vars.fname = fname;
otherwise
Modified: DART/trunk/matlab/GetWRFInfo.m
===================================================================
--- DART/trunk/matlab/GetWRFInfo.m 2011-03-28 22:31:52 UTC (rev 4830)
+++ DART/trunk/matlab/GetWRFInfo.m 2011-03-29 20:43:57 UTC (rev 4831)
@@ -37,6 +37,14 @@
copy = nc_varget(fname,'copy');
times = nc_varget(fname,'time');
+% Coordinate between time types and dates
+
+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));
+dates = times + timeorigin;
+
+
dx = varget(fname, 'DX');
dy = varget(fname, 'DY');
truelat1 = varget(fname, 'TRUELAT1');
@@ -97,15 +105,14 @@
case 'plotcorrel'
disp('Getting information for the ''base'' variable.')
- base_var = GetVarString(pinfo_in.vars);
- [base_time, base_tmeind] = GetTime(times);
- [base_lvl, base_lvlind] = GetLevel(fname,base_var);
- [base_lat, base_latind] = GetLatitude( base_var,TmpJ,VelJ);
- [base_lon, base_lonind] = GetLongitude(base_var,TmpI,VelI);
+ base_var = GetVarString(pinfo_in.vars);
+ [base_time, base_tmeind] = GetTime(dates);
+ [base_lvl, base_lvlind] = GetLevel(fname, base_var);
+ [base_lat, base_lon, base_latind, base_lonind] = GetLatLon(fname, base_var);
disp('Getting information for the ''comparison'' variable.')
- comp_var = GetVarString(pinfo_in.vars, base_var);
- [comp_lvl, comp_lvlind] = GetLevel(fname,comp_var,base_lvl);
+ comp_var = GetVarString(pinfo_in.vars, base_var);
+ [comp_lvl, comp_lvlind] = GetLevel(fname, comp_var, base_lvl);
pinfo.model = model;
pinfo.fname = fname;
@@ -126,7 +133,7 @@
disp('Getting information for the ''base'' variable.')
base_var = GetVarString(pinfo_in.vars);
- [base_time, base_tmeind] = GetTime(times);
+ [base_time, base_tmeind] = GetTime(dates);
[base_lvl , base_lvlind] = GetLevel(fname,base_var);
[base_lat , base_latind] = GetLatitude( base_var,TmpJ,VelJ);
[base_lon , base_lonind] = GetLongitude(base_var,TmpI,VelI);
@@ -271,21 +278,37 @@
function [time, timeind] = GetTime(times, deftime)
%----------------------------------------------------------------------
-if (nargin == 2), time = deftime; else time = mean(times); end
+% Query for the time of interest.
-fprintf('Default time is %f, if this is OK, <cr>;\n',time)
-fprintf('If not, enter a time between %.4f and %.4f, we use the closest.\n', ...
- min(times),max(times))
+ntimes = length(times);
+
+% Determine a sensible default.
+if (nargin == 2),
+ time = 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), time = str2num(varstring); end
+if ~isempty(varstring), tindex = str2num(varstring); end
-d = abs(time - times); % crude distance
-ind = find(min(d) == d); % multiple minima possible
-timeind = ind(1); % use the first one
-time = times(timeind);
+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, deflevel)
%----------------------------------------------------------------------
%
Modified: DART/trunk/matlab/PlotCorrel.m
===================================================================
--- DART/trunk/matlab/PlotCorrel.m 2011-03-28 22:31:52 UTC (rev 4830)
+++ DART/trunk/matlab/PlotCorrel.m 2011-03-29 20:43:57 UTC (rev 4831)
@@ -22,7 +22,7 @@
% pinfo.base_time = 238; % ditto
% PlotCorrel(pinfo) % generates a plot
-%% 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
%
@@ -160,6 +160,81 @@
ax = get(h,'Position');
%set(h,'Position',[ax(1) ax(2) ax(3)/2 ax(4)]);
+ case 'wrf'
+
+ % We are going to correlate one var/time/lvl/lat/lon with
+ % all other lats/lons for a var/time/lvl
+
+ clf;
+
+ times = nc_varget(pinfo.fname, 'time');
+
+ % 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;
+
+ % 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');
+
+ inds = (lonmat < 0); % Convert to 0,360 to minimize dateline probs.
+ lonmat(inds) = lonmat(inds) + 360.0;
+
+ % Get the actual goods ... and perform the correlation
+
+ base_mem = Get1Ens( pinfo.fname, pinfo.base_var, pinfo.base_tmeind, ...
+ pinfo.base_lvlind, pinfo.base_latind, pinfo.base_lonind );
+
+ comp_ens = GetEnsLevel( pinfo.fname, pinfo.comp_var, ...
+ pinfo.base_tmeind, pinfo.comp_lvlind);
+
+ nmembers = size(comp_ens,1);
+
+ corr = zeros(nxny,1);
+
+ for i = 1:nxny,
+ x = corrcoef(base_mem, comp_ens(:, i));
+ corr(i) = x(1, 2);
+ end
+
+ correl = reshape(corr,[ny nx]);
+
+ % Plot it up ...
+
+ [cs,h] = contour(lonmat,latmat,correl,-1:0.2:1);
+ clabel(cs,h,'FontSize',12,'Color','k','Rotation',0);
+ 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', ...
+ 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, nmembers);
+ title({s1,s2,pinfo.fname},'interpreter','none','fontweight','bold')
+ xlabel(sprintf('longitude (%s)',lonunits),'interpreter','none')
+ ylabel(sprintf('latitude (%s)',latunits),'interpreter','none')
+ worldmap;
+ axis image
+ h = colorbar;
+ ax = get(h,'Position');
+ %set(h,'Position',[ax(1) ax(2) ax(3)/2 ax(4)]);
+
case 'mitgcm_ocean'
% We are going to correlate one var/time/lvl/lat/lon with
@@ -289,9 +364,9 @@
% find which are actual ensemble members
metadata = nc_varget(fname,'CopyMetaData'); % get all the metadata
-copyindices = strmatch('ensemble member',metadata); % find all 'member's
+copyindices = strmatch('ensemble member',metadata); % find all 'member's
-if ( isempty(copyindices) )
+if ( length(copyindices) == 1 && copyindices == 1 )
fprintf('%s has no valid ensemble members\n',fname)
disp('To be a valid ensemble member, the CopyMetaData for the member')
disp('must start with the character string ''ensemble member''')
@@ -325,7 +400,7 @@
% find which are actual ensemble members
metadata = nc_varget(fname,'CopyMetaData'); % get all the metadata
-copyindices = strmatch('ensemble member',metadata); % find all 'member's
+copyindices = strmatch('ensemble member',metadata); % find all 'member's
if ( isempty(copyindices) )
fprintf('%s has no valid ensemble members\n',fname)
Modified: DART/trunk/matlab/plot_correl.m
===================================================================
--- DART/trunk/matlab/plot_correl.m 2011-03-28 22:31:52 UTC (rev 4830)
+++ DART/trunk/matlab/plot_correl.m 2011-03-29 20:43:57 UTC (rev 4831)
@@ -99,6 +99,10 @@
pinfo = GetCamInfo(pinfo, diagn_file, 'PlotCorrel');
+ case 'wrf'
+
+ pinfo = GetWRFInfo(pinfo, diagn_file, 'PlotCorrel');
+
case 'pe2lyr'
pinfo = GetPe2lyrInfo(pinfo, diagn_file, 'PlotCorrel');
More information about the Dart-dev
mailing list