[Dart-dev] [9718] DART/trunk/diagnostics/matlab: Now works with all versions of Matlab - the continental overlays
nancy at ucar.edu
nancy at ucar.edu
Thu Feb 4 17:05:51 MST 2016
Revision: 9718
Author: thoar
Date: 2016-02-04 17:05:51 -0700 (Thu, 04 Feb 2016)
Log Message:
-----------
Now works with all versions of Matlab - the continental overlays
were not being handled correctly with recent versions of Matlab.
Extended to color-code the scatterplot with colors based on observation value.
time labelling on axes is managed better.
Creates 3 figure windows with one new axes ... observation value vs. QC
All figures use larger, more readable font.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/link_obs.m
DART/trunk/diagnostics/matlab/linked_observations.m
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/link_obs.m
===================================================================
--- DART/trunk/diagnostics/matlab/link_obs.m 2016-02-05 00:02:23 UTC (rev 9717)
+++ DART/trunk/diagnostics/matlab/link_obs.m 2016-02-05 00:05:51 UTC (rev 9718)
@@ -1,24 +1,27 @@
function link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
-%% link_obs generates the 'brushable' observation plots.
+%% link_obs generates the 'brushable' observation plots. Be sure to click on the painbrush icon and graphically select some observations.
%
-% Three figures will be generated.
+% link_obs explores the observations in a netCDF file created by 'obs_seq_to_netcdf'.
+% All the data in the graphics are linked and are 'brushable'. Click on
+% the paintbrush icon and use the mouse to select observations in any
+% graphic. That same observation will be highlighted in ALL the graphics.
%
-% Figure 1 will have a 3D geographic scatterplot.
+% Three figures will be generated.
%
-% Figure 2 has multiple axes.
-% The bottom axes has a plot of the observation
-% value vs. the QC value.
-% The next axes provides information about the original
-% observation index in the observation sequence file.
-% The next axes provides the ability to select observations
-% by time - useful if multiple observation sequence files are
-% contained in the single input netCDF file.
-% The final (top) axes plots the QC value as a function of time.
+% Figure 1 will have a 3D geographic scatterplot. Click on the rotate
+% icon and drag the graphic around for the best view angle.
%
-% Figure 3 has a 2D scatterplot of (typically) the prior mean
-% observation vs. the original observation. Both of
-% these can be changed however - the allowable set is defined
-% by the CopyMetaData variable in the netCDF file.
+% Figure 2 has multiple axes.
+% * The bottom axes has a plot of the observation 'key'
+% (some of the linked list information in the original file).
+% * The middle axes provides information about the observation
+% density as a function of time.
+% * The final (top) axes plots the QC value as a function of time.
+%
+% Figure 3 has a 2D scatterplot of (typically) the prior mean observation
+% vs. the original observation. Both of these can be changed however -
+% the allowable settings are defined in the CopyMetaData variable
+% in the netCDF file.
%
% link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
%
@@ -27,10 +30,12 @@
% ObsCopyString - the COPY specifying the raw observation ( -v CopyMetaData )
% CopyString - the COPY specifying the copy to compare to the raw obs
% QCString - character string - one of (ncdump -v QCMetaData *.nc)
+% region - geographic extent of interest
%
% EXAMPLE 1:
% fname = '/ptmp/thoar/POP/CAM/POP8/obs_epoch_001.nc';
-% ObsTypeString = 'APB_TEMPERATURE';
+% fname = '/Users/thoar/svn/DART/trunk/models/POP/work/obs_epoch_001.nc';
+% ObsTypeString = 'FLOAT_TEMPERATURE';
% ObsCopyString = 'WOD observation';
% CopyString = 'prior ensemble mean';
% QCString = 'DART quality control';
@@ -39,9 +44,9 @@
% link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
%
% EXAMPLE 2:
-% fname = 'obs_epoch_001.nc';
-% ObsTypeString = 'RADIOSONDE_TEMPERATURE';
-% ObsCopyString = 'NCEP BUFR observation';
+% fname = 'obs_epoch_002.nc';
+% ObsTypeString = 'ACARS_U_WIND_COMPONENT';
+% ObsCopyString = 'observation';
% CopyString = 'prior ensemble mean';
% QCString = 'DART quality control';
% region = [220 300 20 60 -Inf Inf];
@@ -90,6 +95,8 @@
obs.ObsCopyString = obs.CopyString;
obs.CopyString = copy.CopyString;
+obs.region(5) = min(obs.z); % use observation Z to specify vertical region
+obs.region(6) = max(obs.z);
%% Now pack the data in the same fashion as the cell array of column labels.
@@ -129,7 +136,7 @@
end
%% create the linked plots
-linked_observations(obs)
+linked_observations(obs);
% <next few lines under version control, do not edit>
Modified: DART/trunk/diagnostics/matlab/linked_observations.m
===================================================================
--- DART/trunk/diagnostics/matlab/linked_observations.m 2016-02-05 00:02:23 UTC (rev 9717)
+++ DART/trunk/diagnostics/matlab/linked_observations.m 2016-02-05 00:05:51 UTC (rev 9718)
@@ -38,32 +38,31 @@
global obsmat
-% Create figure
-%figure1 = figure('XVisual',...
-% '0x24 (TrueColor, depth 24, RGB mask 0xff0000 0xff00 0x00ff)',...
-% 'Renderer','OpenGL');
+%% ------------------------------------------------------------------------
+% Create figure and axes for 3D scatterplot
+% -------------------------------------------------------------------------
+
figure1 = figure(1); clf(figure1);
-%% Create axes for 3D scatterplot
-% should figure out how to query the vertical coordinate to determine
-% direction of the Z axis ...
-axes0 = axes('Parent',figure1,'OuterPosition',[0 0 1 0.90],'FontSize',18);
-view(axes0,[-37.5 30]);
-grid(axes0,'on');
-hold(axes0,'all');
+fig1ax1 = axes('Parent',figure1,'OuterPosition',[0 0 1 0.90]);
+view(fig1ax1,[-37.5 30]);
+grid(fig1ax1,'on');
xstring = sprintf('obsmat(:,%d)',obs.lonindex);
ystring = sprintf('obsmat(:,%d)',obs.latindex);
zstring = sprintf('obsmat(:,%d)',obs.zindex );
+symbolsize = set_symbol_size(numel(obs.obs));
+
scatter3(obsmat(:,obs.lonindex), obsmat(:,obs.latindex), obsmat(:,obs.zindex), ...
- 'Parent',axes0,'DisplayName','observation locations', ...
+ symbolsize, obsmat(:,obs.obsindex), 'filled', ...
+ 'Parent',fig1ax1, ...
+ 'DisplayName','observation locations', ...
'XDataSource',xstring, ...
'YDataSource',ystring, ...
'ZDataSource',zstring);
-axlims = axis;
-axis([ obs.region(1:4) axlims(5) axlims(6)]);
+set(fig1ax1,'FontSize',18);
myworldmap(obs);
@@ -72,101 +71,126 @@
zlabel(obs.colnames{obs.zindex});
h = title({obs.ObsTypeString, ...
+ sprintf('"%s"',obs.ObsCopyString), ...
sprintf('%s ---> %s',obs.timestring(1,:),obs.timestring(2,:)) });
set(h,'Interpreter','none')
linkdata on
-%% Create figure for ancillary plots
+%% ------------------------------------------------------------------------
+% Create figure for ancillary plots
+% -------------------------------------------------------------------------
figure2 = figure(2); clf(figure2); orient tall; wysiwyg
%% Create axes for time VS. QC
-% FIXME ... choose a format for datetick based on the timespan
-axes4 = axes('Parent',figure2,'OuterPosition',[0 0.80 1 0.175],'FontSize',14);
-set(axes4,'XAxisLocation','top')
-box(axes4,'on');
-hold(axes4,'all');
+fig2ax1 = axes('Parent',figure2,'Position',[0.13 0.71 0.78 0.22]);
+set(fig2ax1,'XAxisLocation','bottom')
+box(fig2ax1,'on');
+
xstring = sprintf('obsmat(:,%d)',obs.timeindex);
ystring = sprintf('obsmat(:,%d)',obs.qcindex);
-scatter(obsmat(:,obs.timeindex),obsmat(:,obs.qcindex),'Parent',axes4, ...
- 'DisplayName','time vs qc', ...
- 'XDataSource',xstring, ...
- 'YDataSource',ystring);
-datetick(axes4,'x',6);
+scatter(obsmat(:,obs.timeindex),obsmat(:,obs.qcindex), ...
+ 'Parent', fig2ax1, ...
+ 'DisplayName', 'time vs qc', ...
+ 'XDataSource', xstring, ...
+ 'YDataSource', ystring);
+
+set(fig2ax1,'FontSize',14);
+
+ax = axis;
+ax(3:4) = [0 8]; %min/max range of DART QC values
+axis(ax)
+datetick(fig2ax1,'x')
ylabel(obs.colnames{obs.qcindex});
+grid(fig2ax1,'on');
-%% Create axes for observation index VS. time
-% FIXME ... choose a format for datetick based on the timespan
-%axes3 = axes('Parent',figure2,'OuterPosition',[0 0.575 1 0.175]);
-axes3 = axes('Parent',figure2,'OuterPosition',[0 0.600 1 0.2],'FontSize',14);
-box(axes3,'on');
-hold(axes3,'all');
+%% Create axes for observation count/density VS. time
+fig2ax2 = axes('Parent',figure2,'Position',[0.13 0.41 0.78 0.22]);
+box(fig2ax2,'on');
+
xstring = sprintf('obsmat(:,%d)',obs.timeindex);
ystring = sprintf('obsmat(:,%d)',obs.indindex);
-scatter(obsmat(:,obs.timeindex),obsmat(:,obs.indindex),'Parent',axes3, ...
+scatter(obsmat(:,obs.timeindex),obsmat(:,obs.indindex),'Parent',fig2ax2, ...
'DisplayName','time vs key', ...
'XDataSource',xstring, ...
'YDataSource',ystring);
+
+set(fig2ax2,'FontSize',14);
+
ylabel('obs count');
-datetick(axes3,'x',6);
+xlabel('time');
+datetick(fig2ax2,'x');
%% Create axes for observation index VS. linked list key
-%axes2 = axes('Parent',figure2,'OuterPosition',[0.0 0.375 1 0.2]);
-axes2 = axes('Parent',figure2,'OuterPosition',[0.0 0.35 1 0.25],'FontSize',14);
-box(axes2,'on');
-hold(axes2,'all');
+fig2ax3 = axes('Parent',figure2,'Position',[0.13 0.11 0.78 0.22]);
+box(fig2ax3,'on');
+
xstring = sprintf('obsmat(:,%d)',obs.indindex);
ystring = sprintf('obsmat(:,%d)',obs.keyindex);
-scatter(obsmat(:,obs.indindex),obsmat(:,obs.keyindex),'Parent',axes2, ...
+scatter(obsmat(:,obs.indindex),obsmat(:,obs.keyindex),'Parent',fig2ax3, ...
'DisplayName','count vs key', ...
'XDataSource',xstring, ...
'YDataSource',ystring);
+
+set(fig2ax3,'FontSize',14);
+
xlabel('obs count');
ylabel('key');
-%% Create axes for ObsVal vs. QC scatterplot
-%axes1 = axes('Parent',figure2,'Position',[0.05 0.05 0.6 0.25],'FontSize',14);
-%box(axes1,'on');
-%hold(axes1,'all');
-%
-%xstring = sprintf('obsmat(:,%d)',obs.obsindex);
-%ystring = sprintf('obsmat(:,%d)',obs.qcindex);
-%scatter(obsmat(:,obs.obsindex),obsmat(:,obs.qcindex),'Parent',axes1, ...
-% 'DisplayName','obs vs qc', ...
-% 'XDataSource',xstring, ...
-% 'YDataSource',ystring);
-%xlabel( obs.colnames{obs.obsindex});
-%h = title( obs.ObsTypeString);
-%set(h,'Interpreter','none');
+refreshdata
+linkdata on
+
+%% ------------------------------------------------------------------------
+% Create figure for 2D scatterplots plots
+% -------------------------------------------------------------------------
+
+figure3 = figure(3); clf(figure3); orient tall; wysiwyg
+
+%% Create axes for ObsVal vs. DART QC scatterplot
+fig3ax1 = axes('Parent',figure3,'Position',[0.15 0.675 0.7 0.25]);
+box(fig3ax1,'on');
+
+xstring = sprintf('obsmat(:,%d)',obs.obsindex);
+ystring = sprintf('obsmat(:,%d)',obs.qcindex);
+
+h1 = scatter(obsmat(:,obs.obsindex),obsmat(:,obs.qcindex), ...
+ 'Parent',fig3ax1, ...
+ 'DisplayName','obs vs qc', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+
+set(fig3ax1,'FontSize',14);
+
+xlabel(obs.colnames{obs.obsindex});
+ylabel(obs.colnames{obs.qcindex});
+h = title(obs.ObsTypeString);
+set(h,'Interpreter','none');
+axis([-Inf Inf 0 8])
+grid(fig3ax1,'on');
+
fprintf('QC summary follows:\n')
LabelQC(obs.colnames{obs.qcindex}, obs.qc)
-%refreshdata
-%linkdata on
-
%% Create axes for observation vs ensemble
% This figure is most useful when all the 'bad' obs have been
% replaced by Matlab's NAN so as not to blow the scale.
% The idea is - if both copies 'match', they line up on the diagonal.
-%figure3 = figure(3); clf(figure3);
-
-% axes5 = axes('Parent',figure3,'OuterPosition',[0 0 1 0.95],'FontSize',18);
-axes5 = axes('Parent',figure2,'Position',[0.20 0.05 0.6 0.25],'FontSize',14);
-grid(axes5,'on');
-hold(axes5,'all');
-
+fig3ax2 = axes('Parent',figure3,'Position',[0.15 0.05 0.7 0.5]);
xstring = sprintf('obsmat(:,%d)',obs.obsindex);
ystring = sprintf('obsmat(:,%d)',obs.copyindex);
-scatter(obsmat(:,obs.obsindex), obsmat(:,obs.copyindex), ...
- 'Parent',axes5,'DisplayName','copy1 v copy2', ...
- 'XDataSource',xstring, ...
- 'YDataSource',ystring);
+h2 = scatter(obsmat(:,obs.obsindex), obsmat(:,obs.copyindex), ...
+ 'Parent',fig3ax2, ...
+ 'DisplayName','copy1 v copy2', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+set(fig3ax2,'FontSize',14);
+
xlabel(obs.colnames{obs.obsindex});
ylabel(obs.colnames{obs.copyindex});
h = title({obs.ObsTypeString, ...
@@ -175,15 +199,19 @@
axlims = [min(axis) max(axis) min(axis) max(axis)];
axis(axlims)
-plot(axes5,[min(axis) max(axis)],[min(axis) max(axis)],'k-')
+line([min(axis) max(axis)],[min(axis) max(axis)],'LineWidth',1.5,'Color','k')
+grid(fig3ax2,'on');
-%% thats it folks
+if (sum(isfinite(get(h2,'YData')))) == 0
+ Print_Empty_Banner(obs.colnames{obs.copyindex});
+end
refreshdata
linkdata on
-function s = LabelQC(QCString, qcarray)
+
+function LabelQC(QCString, qcarray)
%% Create legend for (DART) QC values.
%
% 0 observation assimilated
@@ -218,9 +246,6 @@
qccount(i), dartqc_strings{qcvals(i)+1});
end
- % set(gca,'YTick',qcvals,'YAxisLocation','right')
- % set(gca,'YTickLabel',char(s{:}),'FontSize',12)
-
otherwise,
qcvals = unique(qcarray);
@@ -233,13 +258,13 @@
-function h = myworldmap(obs)
+function myworldmap(obs)
%%--------------------------------------------------------------------------
% GET THE ELEVATION DATA AND SET UP THE ASSOCIATED COORDINATE DATA
%---------------------------------------------------------------------------
-load topo; % GET Matlab-native [180x360] ELEVATION DATASET
+topo = load('topo'); % GET Matlab-native [180x360] ELEVATION DATASET
lats = -89.5:89.5; % CREATE LAT ARRAY FOR TOPO MATRIX
lons = 0.5:359.5; % CREATE LON ARRAY FOR TOPO MATRIX
nlon = length(lons);
@@ -250,11 +275,12 @@
% If we didn't explicitly tell it, make a guess.
%---------------------------------------------------------------------------
-ax = axis;
+axis(obs.region);
+ax = obs.region;
-if (ax(1) < -2)
+if (ax(1) < 0)
lons = lons - 180.0;
- topo = [ topo(:,nlon/2+1:nlon) topo(:,1:nlon/2) ];
+ topo.topo = [ topo.topo(:,nlon/2+1:nlon) topo.topo(:,1:nlon/2) ];
end
%%--------------------------------------------------------------------------
@@ -271,50 +297,92 @@
if (isempty(lat_ind1)), lat_ind1 = 1; end;
if (isempty(lat_ind2)), lat_ind2 = nlat; end;
-elev = topo(lat_ind1:lat_ind2,lon_ind1:lon_ind2);
+elev = topo.topo(lat_ind1:lat_ind2,lon_ind1:lon_ind2);
x = lons(lon_ind1:lon_ind2);
y = lats(lat_ind1:lat_ind2);
%%--------------------------------------------------------------------------
+% Augment the colormap and the CLim so that the lowest color index can be
+% forced to a light gray without compromising the data range.
+
+bob = colormap;
+ncolors = length(bob);
+bob(1,:) = 0.7; % set lowest color to be light gray.
+colormap(bob);
+
+cmin = min(obs.obs);
+cmax = max(obs.obs);
+dz = linspace(cmin,cmax,ncolors-1); % desired dynamic range
+dclim = dz(2) - dz(1);
+newcmin = cmin - dclim;
+clim = [newcmin cmax]; % add extra bin at bottom that no data will use.
+set(gca,'CLim',clim)
+colorbar
+
+%%--------------------------------------------------------------------------
% Contour the "subset" - and give the whole thing an appropriate zlevel
% so the continents are either at the top of the plot (for depth), or
% the bottom (for just about everything else.
-% There are differences between 6.5 and 7.0 that make changing the colors
-% of the filled contours a real pain.
-%---------------------------------------------------------------------------
orgholdstate = ishold;
hold on;
+cmd = '[~,h] = contourf(x,y,elev+newcmin,[newcmin newcmin],''k-'');';
+
switch lower(obs.Zunits)
- case 'depth'
- set(gca,'Zdir','reverse')
- zlevel = min(ax(5:6));
- case 'pressure'
- zlevel = max(ax(5:6));
- set(gca,'Zdir','reverse')
- otherwise
- set(gca,'Zdir','normal')
- zlevel = min(ax(5:6));
+ case 'surface'
+ set(gca,'Zdir','normal')
+ zlevel = ax(5);
+ case 'depth'
+ set(gca,'Zdir','reverse')
+ zlevel = ax(5); % minimum depth
+ cmd = '[~,h] = contour(x,y,elev+newcmin,[newcmin newcmin],''k-'');';
+ case 'pressure'
+ set(gca,'Zdir','reverse')
+ zlevel = ax(6); % maximum pressure
+ otherwise
+ set(gca,'Zdir','normal')
+ zlevel = ax(5); % minimum height
end
-fcolor = [0.7 0.7 0.7]; % light grey
+eval(cmd)
-[c,h] = contourf(x,y,elev,[0.0 0.0],'k-');
+t1 = hgtransform('Parent',gca);
+set(h,'Parent',t1);
+m = makehgtform('translate',[0 0 zlevel]);
+set(t1,'Matrix',m)
-h_patch = get(h, 'Children');
+if (orgholdstate == 0), hold off; end;
-for i = 1:numel(h_patch)
- y = get(h_patch(i), 'YData');
- s = size(y);
- set(h_patch(i), 'ZData', zlevel*ones(s),'FaceColor',fcolor);
- set(h_patch(i),'AlphaDataMapping','none','FaceVertexAlphaData',0.3)
- set(h_patch(i),'FaceAlpha',0.3)
+
+
+function y = set_symbol_size(xx)
+% seems like a symbol size of 200 is good for a handful of points
+% seems like a symbol size of 10 is good for about 5000 points (or more)
+% linearly scale between the two ... sucks
+
+if xx > 5000
+ y = 10;
+elseif xx < 1
+ y = 200;
+else
+ x = [ 1 300 500 900 2400 3700 5000];
+ z = [200 150 100 60 20 12 10];
+ cs = spline(x,[0 z 0]);
+ y = round(ppval(cs,xx));
end
-if (orgholdstate == 0), hold off; end;
+function Print_Empty_Banner(mystring)
+
+ax = axis;
+dx = mean(ax(1:2));
+dy = mean(ax(3:4));
+text(dx,dy,sprintf('No "%s" \n with a meaningful DART QC', mystring), ...
+ 'HorizontalAlignment','center', ...
+ 'FontSize',20)
+
% <next few lines under version control, do not edit>
% $URL$
% $Revision$
More information about the Dart-dev
mailing list