[Dart-dev] [4240] DART/trunk/diagnostics/matlab: new terms of use - no other changes
nancy at ucar.edu
nancy at ucar.edu
Thu Jan 28 16:47:25 MST 2010
Revision: 4240
Author: thoar
Date: 2010-01-28 16:47:24 -0700 (Thu, 28 Jan 2010)
Log Message:
-----------
new terms of use - no other changes
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/PlotObsLocs.m
DART/trunk/diagnostics/matlab/PlotWindObs.m
DART/trunk/diagnostics/matlab/get_varnames.m
DART/trunk/diagnostics/matlab/get_varsNdims.m
DART/trunk/diagnostics/matlab/nc_var_dims.m
DART/trunk/diagnostics/matlab/obs_num_time.m
DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_evolution.m
DART/trunk/diagnostics/matlab/plot_observation_locations.m
DART/trunk/diagnostics/matlab/plot_profile.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
DART/trunk/diagnostics/matlab/plot_wind_vectors.m
Added Paths:
-----------
DART/trunk/diagnostics/matlab/linked_observations.m
DART/trunk/diagnostics/matlab/plot_obs_netcdf.m
DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m
DART/trunk/diagnostics/matlab/read_obs_netcdf.m
Property Changed:
----------------
DART/trunk/diagnostics/matlab/PlotWindObs.m
DART/trunk/diagnostics/matlab/plot_wind_vectors.m
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/PlotObsLocs.m
===================================================================
--- DART/trunk/diagnostics/matlab/PlotObsLocs.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/PlotObsLocs.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -96,10 +96,9 @@
% subsequent plots will overwrite earlier ones, so this input
% can be used to customize the name.
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Modified: DART/trunk/diagnostics/matlab/PlotWindObs.m
===================================================================
--- DART/trunk/diagnostics/matlab/PlotWindObs.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/PlotWindObs.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -1,9 +1,8 @@
-%
+% PlotWindObs ... a function in progress ... not for general use.
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Property changes on: DART/trunk/diagnostics/matlab/PlotWindObs.m
___________________________________________________________________
Modified: svn:keywords
- Date Rev Author URL Id
+ Date Revision Author HeadURL Id
Modified: DART/trunk/diagnostics/matlab/get_varnames.m
===================================================================
--- DART/trunk/diagnostics/matlab/get_varnames.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/get_varnames.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -10,10 +10,9 @@
% nvars = length(varnames);
% disp(sprintf('first atmospheric variable (of %d) is %s',nvars,varnames{1}))
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Modified: DART/trunk/diagnostics/matlab/get_varsNdims.m
===================================================================
--- DART/trunk/diagnostics/matlab/get_varsNdims.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/get_varsNdims.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -19,10 +19,9 @@
% >> plotdat.allvardims{20}
% region plevel copy time
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Copied: DART/trunk/diagnostics/matlab/linked_observations.m (from rev 4237, DART/trunk/observations/utilities/threed_sphere/linked_observations.m)
===================================================================
--- DART/trunk/diagnostics/matlab/linked_observations.m (rev 0)
+++ DART/trunk/diagnostics/matlab/linked_observations.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -0,0 +1,162 @@
+function linked_observations(obsmat,obs)
+% linked_observations(obs)
+%
+% obs is a structure with the following required components
+%
+% obs.lons longitudes of the observations
+% obs.lats latitudes of the observations
+% obs.z vertical level (depth) of the observations
+% obs.obs observation values
+% obs.qc observation DART QC code
+
+%% DART software - Copyright \xA9 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% Create figure
+%figure1 = figure('XVisual',...
+% '0x24 (TrueColor, depth 24, RGB mask 0xff0000 0xff00 0x00ff)',...
+% 'Renderer','OpenGL');
+figure1 = figure(1); clf(figure1);;
+
+%% Create axes for 3D plot
+axes0 = axes('Parent',figure1,'OuterPosition',[0 0 1 0.90],'FontSize',12);
+view(axes0,[-37.5 30]);
+grid(axes0,'on');
+hold(axes0,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.lonindex);
+ystring = sprintf('obsmat(:,%d)',obs.latindex);
+zstring = sprintf('obsmat(:,%d)',obs.zindex );
+
+h0 = scatter3(obsmat(:,obs.lonindex), obsmat(:,obs.latindex), obsmat(:,obs.zindex), ...
+ 'Parent',axes0,'DisplayName','observation locations', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring, ...
+ 'ZDataSource',zstring);
+
+worldmap('light');
+xlabel('longitude');
+ylabel('latitude');
+zlabel('depth');
+h = title({obs.ObsTypeString, ...
+ sprintf('%s ---> %s',obs.timestring(1,:),obs.timestring(2,:)) });
+set(h,'Interpreter','none')
+linkdata on
+
+%% Create figure for ancillary plots
+
+figure2 = figure(2); clf(figure2); orient tall; wysiwyg
+
+%% Create axes for time VS. QC
+axes4 = axes('Parent',figure2,'OuterPosition',[0 0.80 1 0.175]);
+set(axes4,'XAxisLocation','top')
+box(axes4,'on');
+hold(axes4,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.timeindex);
+ystring = sprintf('obsmat(:,%d)',obs.qcindex);
+h4 = scatter(obsmat(:,obs.timeindex),obsmat(:,obs.qcindex),'Parent',axes4, ...
+ 'DisplayName','time vs qc', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+datetick(axes4,'x',6);
+ylabel(obs.QCString);
+
+%% Create axes for observation index VS. time
+axes3 = axes('Parent',figure2,'OuterPosition',[0 0.575 1 0.175]);
+box(axes3,'on');
+hold(axes3,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.timeindex);
+ystring = sprintf('obsmat(:,%d)',obs.indindex);
+h3 = scatter(obsmat(:,obs.timeindex),obsmat(:,obs.indindex),'Parent',axes3, ...
+ 'DisplayName','time vs key', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+ylabel('obs count');
+datetick(axes3,'x',6);
+
+%% Create axes for observation index VS. linked list key
+axes2 = axes('Parent',figure2,'OuterPosition',[0.0 0.400 1 0.15]);
+box(axes2,'on');
+hold(axes2,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.indindex);
+ystring = sprintf('obsmat(:,%d)',obs.keyindex);
+h2 = scatter(obsmat(:,obs.indindex),obsmat(:,obs.keyindex),'Parent',axes2, ...
+ 'DisplayName','count vs key', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+xlabel('obs count');
+ylabel('key');
+
+%% Create axes for QC vs. ObsVal scatterplot
+axes1 = axes('Parent',figure2,'Position',[0.05 0.05 0.6 0.25]);
+box(axes1,'on');
+hold(axes1,'all');
+
+xstring = sprintf('obsmat(:,%d)',obs.obsindex);
+ystring = sprintf('obsmat(:,%d)',obs.qcindex);
+h1 = scatter(obsmat(:,obs.obsindex),obsmat(:,obs.qcindex),'Parent',axes1, ...
+ 'DisplayName','obs vs qc', ...
+ 'XDataSource',xstring, ...
+ 'YDataSource',ystring);
+xlabel(obs.CopyString);
+title(obs.QCString);
+
+
+LabelQC(obs.QCString, obs.qc)
+
+
+refreshdata
+linkdata on
+
+
+function LabelQC(QCString, qcarray)
+%% Create legend for (DART) QC values.
+%
+% 0 observation assimilated
+% 1 observation evaluated only
+% --- everything above this means the prior and posterior are OK
+% 2 assimilated, but the posterior forward operator failed
+% 3 Evaluated only, but the posterior forward operator failed
+% --- everything above this means only the prior is OK
+% 4 prior forward operator failed
+% 5 not used
+% 6 prior QC rejected
+% 7 outlier rejected
+
+dartqc_strings = { 'assimilated', ...
+ 'observation evaluated only', ...
+ 'assimilated, post fwd op failed', ...
+ 'eval only, post fwd op failed', ...
+ 'prior forward operator failed', ...
+ 'not used', ...
+ 'prior QC rejected', ...
+ 'outlier rejected', ...
+ 'reserved'};
+
+switch lower(strtrim(QCString))
+ case 'dart quality control',
+
+ qcvals = unique(qcarray);
+ qccount = zeros(size(qcvals));
+ for i = 1:length(qcvals)
+ qccount(i) = sum(qcarray == qcvals(i));
+ s{i} = sprintf('%d - %s %d obs',qcvals(i), dartqc_strings{qcvals(i)+1}, qccount(i));
+ end
+
+ set(gca,'YTick',qcvals,'YAxisLocation','right')
+ set(gca,'YTickLabel',char(s))
+
+ otherwise,
+ str = sprintf('no way to interpret values of %s',strtrim(QCString));
+ text(0.0, 0.0, str)
+end
Modified: DART/trunk/diagnostics/matlab/nc_var_dims.m
===================================================================
--- DART/trunk/diagnostics/matlab/nc_var_dims.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/nc_var_dims.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -17,10 +17,9 @@
% disp(sprintf('variable %s dimension %d is %s',varname,i,vdims{i}))
% end
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2009, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Modified: DART/trunk/diagnostics/matlab/obs_num_time.m
===================================================================
--- DART/trunk/diagnostics/matlab/obs_num_time.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/obs_num_time.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -17,11 +17,9 @@
%
% obs_num_time
-
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Modified: DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_bias_xxx_profile.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -18,11 +18,9 @@
% copystring = 'totalspread'; % 'copy' string == quantity of interest
% plotdat = plot_bias_xxx_profile(fname,copystring);
-
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL:
Modified: DART/trunk/diagnostics/matlab/plot_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_evolution.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_evolution.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -18,10 +18,9 @@
% copystring = 'bias'; % 'copy' string == quantity of interest
% plotdat = plot_evolution(fname,copystring);
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Copied: DART/trunk/diagnostics/matlab/plot_obs_netcdf.m (from rev 4237, DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m)
===================================================================
--- DART/trunk/diagnostics/matlab/plot_obs_netcdf.m (rev 0)
+++ DART/trunk/diagnostics/matlab/plot_obs_netcdf.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -0,0 +1,350 @@
+function obsstruct = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, ...
+ QCString, maxQC, verbose)
+%% plot_obs_netcdf will plot the locations and values of the observations in a DART netcdf file.
+% any observations with a QC value greater than 'maxgoodQC' will get
+% plotted on a separate figure ... color-coded to its QC value, not the
+% observation value.
+%
+%--------------------------------------------------
+% EXAMPLE 1: plotting just one type of observation
+%--------------------------------------------------
+% fname = 'obs_sequence_001.nc';
+% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT';
+% region = [0 360 -90 90 -Inf Inf];
+% CopyString = 'NCEP BUFR observation';
+% QCString = 'DART quality control';
+% maxgoodQC = 2;
+% verbose = 1; % anything > 0 == 'true'
+%
+% bob = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxgoodQC, verbose);
+%
+% view(0,90); % for a traditional '2D' plot
+%
+%--------------------------------------------------
+% EXAMPLE 2: plotting all the observation types
+%--------------------------------------------------
+% fname = 'obs_sequence_001.nc';
+% ObsTypeString = 'ALL';
+% region = [0 360 -90 90 -Inf Inf];
+% CopyString = 'WOD observation';
+% QCString = 'WOD QC';
+% maxgoodQC = 0;
+% verbose = 1; % anything > 0 == 'true'
+%
+% bob = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxgoodQC, verbose);
+
+%% DART software - Copyright \xA9 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
+
+%% Read the observation sequence
+
+obsstruct = read_obs_netcdf(fname, ObsTypeString, region, ...
+ CopyString, QCString, maxQC, verbose);
+
+%% Create graphic with area-weighted symbols for the good observations.
+% It has happened that there have been zero good observations in a file.
+
+xmin = min(region(1:2));
+xmax = max(region(1:2));
+ymin = min(region(3:4));
+ymax = max(region(3:4));
+zmin = min(obsstruct.z);
+zmax = max(obsstruct.z);
+
+pstruct.colorbarstring = obsstruct.ObsTypeString;
+pstruct.region = region;
+pstruct.str3 = sprintf('%s - %s',obsstruct.timestring(1,:),obsstruct.timestring(2,:));
+
+if ( length(obsstruct.obs) < 1 )
+ fprintf('There are no ''good'' observations to plot\n')
+else
+
+ figure(gcf+1); clf
+
+ % choose a symbol size based on the number of obs to plot.
+
+ if (length(obsstruct.obs) < 1000)
+ pstruct.scalearray = scaleme(obsstruct.obs, 36);
+ else
+ pstruct.scalearray = 50.0 * ones(size(obsstruct.obs));
+ end
+ pstruct.clim = [min(obsstruct.obs) max(obsstruct.obs)];
+ pstruct.str2 = sprintf('%s (%d locations)',obsstruct.CopyString,length(obsstruct.obs));
+
+ % If all the observations live on the same level ... make a 2D plot.
+
+ if ( zmin ~= zmax )
+
+ pstruct.axis = [xmin xmax ymin ymax zmin zmax];
+ pstruct.str1 = sprintf('%s level (%.2f - %.2f)',obsstruct.ObsTypeString,zmin,zmax);
+
+ plot_3D(obsstruct, pstruct)
+
+ else
+
+ pstruct.axis = [xmin xmax ymin ymax];
+ pstruct.str1 = sprintf('%s',obsstruct.ObsTypeString);
+
+ plot_2D(obsstruct, pstruct)
+
+ end
+end
+
+%% Create graphic of spatial distribution of 'bad' observations & their QC value.
+%
+% 0 observation assimilated
+% 1 observation evaluated only
+% --- everything above this means the prior and posterior are OK
+% 2 assimilated, but the posterior forward operator failed
+% 3 Evaluated only, but the posterior forward operator failed
+% --- everything above this means only the prior is OK
+% 4 prior forward operator failed
+% 5 not used
+% 6 prior QC rejected
+% 7 outlier rejected
+
+dartqc_strings = { ...
+ '''observation evaluated only''', ...
+ '''assimilated, but the posterior forward operator failed''', ...
+ '''evaluated only, but the posterior forward operator failed''',...
+ '''prior forward operator failed''',...
+ '''not used''',...
+ '''prior QC rejected''',...
+ '''outlier rejected''',...
+ '''reserved for future use'''};
+
+if (obsstruct.numbadqc > 0 ) % if there are bad observation to plot ... carry on.
+
+ figure(gcf+1); clf
+
+ subplot('position',[0.1 0.20 0.8 0.65])
+
+ zmin = min(obsstruct.badobs.z);
+ zmax = max(obsstruct.badobs.z);
+
+ pstruct.scalearray = 128 * ones(size(obsstruct.badobs.obs));
+ pstruct.colorbarstring = QCString;
+ pstruct.clim = [min(obsstruct.badobs.qc) max(obsstruct.badobs.qc)];
+ pstruct.str1 = sprintf('%s level (%.2f - %.2f)',obsstruct.ObsTypeString,zmin,zmax);
+ pstruct.str2 = sprintf('%s (%d bad observations)', ...
+ obsstruct.CopyString, ...
+ length(obsstruct.badobs.obs));
+
+ obsstruct.badobs.obs = obsstruct.badobs.qc; % plot QC values, not obs values
+ if ( zmin ~= zmax )
+
+ pstruct.axis = [xmin xmax ymin ymax zmin zmax];
+
+ plot_3D(obsstruct.badobs, pstruct)
+
+ else
+
+ pstruct.axis = [xmin xmax ymin ymax];
+
+ plot_2D(obsstruct.badobs, pstruct)
+
+ end
+
+ subplot('position',[0.1 0.05 0.8 0.10])
+ axis off
+
+ %% If the QC is a DART QC, we know how to interpret them.
+
+ switch lower(strtrim(QCString))
+ case 'dart quality control',
+
+ qcvals = unique(obsstruct.badobs.qc);
+ qccount = zeros(size(qcvals));
+ for i = 1:length(qcvals)
+ qccount(i) = sum(obsstruct.badobs.qc == qcvals(i));
+ s{i} = sprintf('%d obs with qc == %d %s',qccount(i),qcvals(i), ...
+ dartqc_strings{qcvals(i)});
+ end
+
+ dy = 1.0/length(s);
+ for i = 1:length(s)
+ text(0.0, (i-1)*dy ,s{i})
+ end
+
+ otherwise,
+ str = sprintf('no way to interpret values of %s',strtrim(QCString));
+ text(0.0, 0.0, str)
+ end
+end
+
+
+
+function h = myworldmap
+
+%%--------------------------------------------------------------------------
+% GET THE ELEVATION DATA AND SET UP THE ASSOCIATED COORDINATE DATA
+%---------------------------------------------------------------------------
+
+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);
+nlat = length(lats);
+
+%%--------------------------------------------------------------------------
+% IF WE NEED TO SWAP HEMISPHERES, DO SO NOW.
+% If we didn't explicitly tell it, make a guess.
+%---------------------------------------------------------------------------
+
+ax = axis;
+
+if (ax(1) < -2)
+ lons = lons - 180.0;
+ topo = [ topo(:,nlon/2+1:nlon) topo(:,1:nlon/2) ];
+end
+
+%%--------------------------------------------------------------------------
+% We need to determine the geographic subset of the elevation matrix.
+%---------------------------------------------------------------------------
+
+lon_ind1 = min(find(ax(1) <= lons));
+lon_ind2 = min(find(ax(2) <= lons));
+lat_ind1 = min(find(ax(3) <= lats));
+lat_ind2 = min(find(ax(4) <= lats));
+
+if (isempty(lon_ind1)), lon_ind1 = 1; end;
+if (isempty(lon_ind2)), lon_ind2 = nlon; end;
+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);
+x = lons(lon_ind1:lon_ind2);
+y = lats(lat_ind1:lat_ind2);
+
+%%--------------------------------------------------------------------------
+% Contour the "subset"
+% 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;
+
+switch get(gca,'ZDir')
+ case 'reverse'
+ zlevel = max(ax(5:6));
+ otherwise
+ zlevel = min(ax(5:6));
+end
+
+fcolor = [0.7 0.7 0.7]; % light grey
+
+[c,h] = contourf(x,y,elev,[0.0 0.0],'k-');
+
+h_patch = get(h, 'Children');
+
+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);
+end
+
+if (orgholdstate == 0), hold off; end;
+
+
+
+
+function s = scaleme(x,minsize)
+% scaleme returns a uniformly scaled array the same size as the input
+% array where the maximum is 10 times the minimum
+maxsize = 10*minsize;
+minx = min(x);
+maxx = max(x);
+slope = (maxsize-minsize)/(maxx-minx);
+b = minsize - slope*minx;
+
+s = x*slope + b;
+
+
+
+function plot_3D(obsstruct, pstruct)
+
+if (pstruct.clim(1) == pstruct.clim(2))
+ % If all the observations have the same value, setting the
+ % colorbar limits is a real pain. Fundamentally, I am
+ % forcing the plot symbols to be the lowest color of the
+ % colormap and setting the colorbar to have some more
+ % colors 'on top' - that are never used.
+ cmap = colormap;
+ h = plot3(obsstruct.lons, obsstruct.lats, obsstruct.z, 'bo');
+ set(h,'MarkerFaceColor',cmap(1,:),'MarkerEdgeColor',cmap(1,:))
+ set(gca,'Clim',[pstruct.clim(1) pstruct.clim(2)+1])
+ set(gca,'XGrid','on','YGrid','on','ZGrid','on')
+
+else
+ scatter3(obsstruct.lons, obsstruct.lats, obsstruct.z, ...
+ pstruct.scalearray, obsstruct.obs, 'd', 'filled');
+end
+
+clim = get(gca,'CLim');
+
+axis(pstruct.axis)
+
+title( {pstruct.str1, pstruct.str3, pstruct.str2}, 'Interpreter','none','FontSize',16);
+xlabel('longitude')
+ylabel('latitude')
+
+if (obsstruct.Ztyp(1) == -2) % VERTISUNDEF = -2
+ zlabel('unspecified')
+elseif (obsstruct.Ztyp(1) == -1) % VERTISSURFACE = -1
+ zlabel('surface')
+elseif (obsstruct.Ztyp(1) == 1) % VERTISLEVEL = 1
+ zlabel('level')
+elseif (obsstruct.Ztyp(1) == 2) % VERTISPRESSURE = 2
+ set(gca,'ZDir','reverse')
+ zlabel('pressure')
+elseif (obsstruct.Ztyp(1) == 3) % VERTISHEIGHT = 3
+ zlabel('height')
+end
+
+myworldmap;
+set(gca,'CLim',clim)
+h = colorbar;
+set(get(h,'YLabel'),'String',pstruct.colorbarstring,'Interpreter','none')
+
+
+
+
+function plot_2D(obsstruct, pstruct)
+
+axis(pstruct.axis); hold on; worldmap('light');
+
+if (pstruct.clim(1) == pstruct.clim(2))
+ cmap = colormap;
+ h = plot(obsstruct.lons, obsstruct.lats, 'bo');
+ set(h,'MarkerFaceColor',cmap(1,:),'MarkerEdgeColor',cmap(1,:))
+ set(gca,'Clim',[pstruct.clim(1) pstruct.clim(2)+1])
+ set(gca,'XGrid','on','YGrid','on')
+
+else
+
+ scatter(obsstruct.lons, obsstruct.lats, ...
+ pstruct.scalearray, obsstruct.obs, 'd', 'filled');
+end
+
+clim = get(gca,'CLim');
+
+title( {pstruct.str1, pstruct.str3, pstruct.str2}, 'Interpreter','none','FontSize',16);
+xlabel('longitude')
+ylabel('latitude')
+
+set(gca,'CLim',clim)
+h = colorbar;
+set(get(h,'YLabel'),'String',pstruct.colorbarstring,'Interpreter','none')
+hold off
Copied: DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m (from rev 4237, DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m)
===================================================================
--- DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m (rev 0)
+++ DART/trunk/diagnostics/matlab/plot_obs_netcdf_diffs.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -0,0 +1,345 @@
+function obsstruct = plot_obs_netcdf_diffs(fname, ObsTypeString, region, ...
+ CopyString1, CopyString2, QCString, maxQC, verbose)
+%
+% fname = 'obs_sequence_001.nc';
+% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT';
+% region = [0 360 -90 90 -Inf Inf];
+% CopyString1 = 'NCEP BUFR observation';
+% CopyString2 = 'prior ensemble mean';
+% QCString = 'DART quality control';
+% maxQC = 1;
+% verbose = 1; % anything > 0 == 'true'
+%
+% bob = plot_obs_netcdf_diffs(fname, ObsTypeString, region, CopyString1, CopyString2, ...
+% QCString, maxQC, verbose);
+
+% record the user input
+
+%% DART software - Copyright \xA9 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+obsstruct.fname = fname;
+obsstruct.ObsTypeString = ObsTypeString;
+obsstruct.region = region;
+obsstruct.CopyString1 = CopyString1;
+obsstruct.CopyString2 = CopyString2;
+obsstruct.QCString = QCString;
+obsstruct.maxQC = maxQC;
+obsstruct.verbose = verbose;
+
+% get going
+
+ObsTypes = nc_varget(fname,'ObsTypes');
+ObsTypeStrings = nc_varget(fname,'ObsTypesMetaData');
+CopyStrings = nc_varget(fname,'CopyMetaData');
+QCStrings = nc_varget(fname,'QCMetaData');
+
+t = nc_varget(fname,'time');
+obs_type = nc_varget(fname,'obs_type');
+z_type = nc_varget(fname,'which_vert');
+
+loc = nc_varget(fname,'location');
+obs = nc_varget(fname,'observations');
+qc = nc_varget(fname,'qc');
+
+my_types = unique(obs_type); % only ones in the file, actually.
+timeunits = nc_attget(fname,'time','units');
+timerange = nc_attget(fname,'time','valid_range');
+calendar = nc_attget(fname,'time','calendar');
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+timestring = datestr(timerange + timeorigin);
+
+% Echo summary if requested
+
+if ( verbose > 0 )
+ for i = 1:length(my_types)
+ obtype = my_types(i);
+ inds = find(obs_type == obtype);
+ myz = loc(inds,3);
+
+ disp(sprintf('N = %6d %s obs (type %3d) between levels %.2f and %.2f', ...
+ length(inds), ObsTypeStrings(obtype,:), obtype, ...
+ unique(min(myz)), unique(max(myz))))
+ end
+
+% uniquelevels = unique(loc(:,3));
+
+% for i = 1:length(uniquelevels)
+% mylevel = uniquelevels(i);
+% inds = find(loc(:,3) == mylevel);
+% disp(sprintf('level %2d %f has %d observations',i,mylevel,length(inds)))
+% end
+
+end
+
+% Find observations of the correct types.
+
+myind = strmatch(ObsTypeString,ObsTypeStrings);
+
+if ( isempty(myind) )
+ error('no %s observations ... stopping',obsstruct.ObsTypeString)
+end
+
+mytype1 = get_copy_index(fname, CopyString1);
+mytype2 = get_copy_index(fname, CopyString2);
+inds = find(obs_type == myind);
+mylocs = loc(inds,:);
+myobs1 = obs(inds,mytype1);
+myobs2 = obs(inds,mytype2);
+myobs = myobs2 - myobs1;
+
+if ~ isempty(QCString)
+ myQCind = get_qc_index(fname, QCString);
+ myqc = qc(inds,myQCind);
+else
+ myqc = [];
+end
+
+clear myobs1 myobs2 obs loc qc
+
+% geographic subset if needed
+
+inds = locations_in_region(mylocs,region);
+
+obsstruct.lons = mylocs(inds,1);
+obsstruct.lats = mylocs(inds,2);
+obsstruct.z = mylocs(inds,3);
+obsstruct.obs = myobs(inds);
+obsstruct.Ztyp = z_type(inds);
+obsstruct.numbadqc = 0;
+
+if (isempty(myqc))
+ obsstruct.qc = [];
+else
+ obsstruct.qc = myqc(inds);
+end
+
+% subset based on qc value
+
+if ( (~ isempty(myqc)) & (~ isempty(maxQC)) )
+
+ inds = find(obsstruct.qc > maxQC);
+
+ obsstruct.numbadqc = length(inds);
+
+ if (~isempty(inds))
+ badobs.lons = obsstruct.lons(inds);
+ badobs.lats = obsstruct.lats(inds);
+ badobs.Ztyp = obsstruct.Ztyp(inds);
+ badobs.z = obsstruct.z( inds);
+ badobs.obs = obsstruct.obs(inds);
+ badobs.qc = obsstruct.qc(inds);
+ end
+
+ disp(sprintf('Removing %d obs with a %s value greater than %f', ...
+ length(inds),QCString,maxQC))
+
+ inds = find(obsstruct.qc <= maxQC);
+
+ bob = obsstruct.lons(inds); obsstruct.lons = bob;
+ bob = obsstruct.lats(inds); obsstruct.lats = bob;
+ bob = obsstruct.Ztyp(inds); obsstruct.Ztyp = bob;
+ bob = obsstruct.z( inds); obsstruct.z = bob;
+ bob = obsstruct.obs( inds); obsstruct.obs = bob;
+ bob = obsstruct.qc( inds); obsstruct.qc = bob;
+
+end
+
+%-------------------------------------------------------------------------------
+% Create graphic with area-weighted symbols for the good observations.
+%-------------------------------------------------------------------------------
+
+figure(1); clf
+
+xmin = min(region(1:2));
+xmax = max(region(1:2));
+ymin = min(region(3:4));
+ymax = max(region(3:4));
+zmin = min(obsstruct.z);
+zmax = max(obsstruct.z);
+
+scalearray = scaleme(obsstruct.obs,36);
+scalearray = 128 * ones(size(obsstruct.obs));
+
+scatter3(obsstruct.lons, obsstruct.lats, obsstruct.z, ...
+ scalearray, obsstruct.obs,'d','filled');
+
+axis([xmin xmax ymin ymax zmin zmax])
+
+str1 = sprintf('%s level (%.2f - %.2f)',ObsTypeString,zmin,zmax);
+str2 = sprintf('%s - %s (%d locations)',CopyString2,CopyString1,length(obsstruct.obs));
+str3 = sprintf('%s - %s',timestring(1,:),timestring(2,:));
+
+title( {str1, str3, str2}, 'Interpreter','none','FontSize',16);
+xlabel('longitude')
+ylabel('latitude')
+
+if (obsstruct.Ztyp(1) == -2) % VERTISUNDEF = -2
+ zlabel('curious ... undefined')
+elseif (obsstruct.Ztyp(1) == -1) % VERTISSURFACE = -1
+ zlabel('surface')
+elseif (obsstruct.Ztyp(1) == 1) % VERTISLEVEL = 1
+ zlabel('level')
+elseif (obsstruct.Ztyp(1) == 2) % VERTISPRESSURE = 2
+ set(gca,'ZDir','reverse')
+ zlabel('pressure')
+elseif (obsstruct.Ztyp(1) == 3) % VERTISHEIGHT = 3
+ zlabel('height')
+end
+
+myworldmap;
+set(gca,'CLim',[min(obsstruct.obs) max(obsstruct.obs)])
+h = colorbar;
+set(get(h,'YLabel'),'String',ObsTypeString,'Interpreter','none')
+
+%-------------------------------------------------------------------------------
+% Create graphic of spatial distribution of 'bad' observations & their QC value.
+%-------------------------------------------------------------------------------
+
+if (obsstruct.numbadqc > 0 )
+
+ figure(2); clf
+
+ subplot('position',[0.1 0.20 0.8 0.65])
+ scalearray = 128 * ones(size(badobs.obs));
+
+ zmin = min(badobs.z);
+ zmax = max(badobs.z);
+
+ scatter3(badobs.lons, badobs.lats, badobs.z, scalearray, badobs.qc,'filled')
+
+ title( {str1, str3, 'Bad Observations'}, 'Interpreter','none','FontSize',16);
+ xlabel('longitude')
+ ylabel('latitude')
+
+ if (badobs.Ztyp(1) == -2) % VERTISUNDEF = -2
+ zlabel('curious ... undefined')
+ elseif (badobs.Ztyp(1) == -1) % VERTISSURFACE = -1
+ zlabel('surface')
+ elseif (badobs.Ztyp(1) == 1) % VERTISLEVEL = 1
+ zlabel('level')
+ elseif (badobs.Ztyp(1) == 2) % VERTISPRESSURE = 2
+ set(gca,'ZDir','reverse')
+ zlabel('pressure')
+ elseif (badobs.Ztyp(1) == 3) % VERTISHEIGHT = 3
+ zlabel('height')
+ end
+
+ axis([region(1) region(2) ymin ymax zmin zmax])
+
+ myworldmap;
+ set(gca,'CLim',[min(badobs.qc) max(badobs.qc)])
+ h = colorbar;
+ set(get(h,'YLabel'),'String',QCString,'Interpreter','none')
+
+ subplot('position',[0.1 0.05 0.8 0.10])
+ axis off
+
+ qcvals = unique(badobs.qc);
+ qccount = zeros(size(qcvals));
+ for i = 1:length(qcvals)
+ qccount(i) = sum(badobs.qc == qcvals(i));
+ s{i} = sprintf('%d obs with qc == %d',qccount(i),qcvals(i));
+ end
+
+ dy = 1.0/length(s);
+ for i = 1:length(s)
+ text(0.0, (i-1)*dy ,s{i})
+ end
+
+end
+
+function h = myworldmap
+
+%---------------------------------------------------------------------------
+% GET THE ELEVATION DATA AND SET UP THE ASSOCIATED COORDINATE DATA
+%---------------------------------------------------------------------------
+
+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);
+nlat = length(lats);
+
+%---------------------------------------------------------------------------
+% IF WE NEED TO SWAP HEMISPHERES, DO SO NOW.
+% If we didn't explicitly tell it, make a guess.
+%---------------------------------------------------------------------------
+
+ax = axis;
+
+if (ax(1) < -2)
+ lons = lons - 180.0;
+ topo = [ topo(:,nlon/2+1:nlon) topo(:,1:nlon/2) ];
+end
+
+%---------------------------------------------------------------------------
+% We need to determine the geographic subset of the elevation matrix.
+%---------------------------------------------------------------------------
+
+lon_ind1 = min(find(ax(1) <= lons));
+lon_ind2 = min(find(ax(2) <= lons));
+lat_ind1 = min(find(ax(3) <= lats));
+lat_ind2 = min(find(ax(4) <= lats));
+
+if (isempty(lon_ind1)) lon_ind1 = 1; end;
+if (isempty(lon_ind2)) lon_ind2 = nlon; end;
+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);
+x = lons(lon_ind1:lon_ind2);
+y = lats(lat_ind1:lat_ind2);
+
+%---------------------------------------------------------------------------
+% Contour the "subset"
+% There are differences between 6.5 and 7.0 that make changing the colors
+% of the filled contours a real pain. Providing both solutions.
+%---------------------------------------------------------------------------
+
+orgholdstate = ishold;
+hold on;
+
+switch get(gca,'ZDir')
+ case 'reverse'
+ zlevel = max(ax(5:6));
+ otherwise
+ zlevel = min(ax(5:6));
+end
+
+fcolor = [0.7 0.7 0.7]; % light grey
+
+[c,h] = contourf(x,y,elev,[0.0 0.0],'k-');
+
+new_level = 1000;
+
+h_patch = get(h, 'Children');
+
+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);
+end
+
+if (orgholdstate == 0) hold off; end;
+
+
+function s = scaleme(x,minsize)
+% scaleme returns a uniformly scaled array the same size as the input
+% array where the maximum is 10 times the minimum
+maxsize = 10*minsize;
+minx = min(x);
+maxx = max(x);
+slope = (maxsize-minsize)/(maxx-minx);
+b = minsize - slope*minx;
+
+s = x*slope + b;
+
Modified: DART/trunk/diagnostics/matlab/plot_observation_locations.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_observation_locations.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_observation_locations.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -11,10 +11,9 @@
% See the documentation for PlotObsLocs() -- it has a lot of arguments in the
% calling sequence.
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Modified: DART/trunk/diagnostics/matlab/plot_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_profile.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_profile.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -19,10 +19,9 @@
% plotdat = plot_profile(fname,copystring);
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL:
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_evolution.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -18,10 +18,9 @@
% copystring = 'totalspread'; % 'copy' string == quantity of interest
% plotdat = plot_rmse_xxx_evolution(fname,copystring);
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL:
Modified: DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_rmse_xxx_profile.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -18,11 +18,9 @@
% copystring = 'totalspread'; % 'copy' string == quantity of interest
% plotdat = plot_rmse_xxx_profile(fname,copystring);
-
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL:
Modified: DART/trunk/diagnostics/matlab/plot_wind_vectors.m
===================================================================
--- DART/trunk/diagnostics/matlab/plot_wind_vectors.m 2010-01-28 23:41:56 UTC (rev 4239)
+++ DART/trunk/diagnostics/matlab/plot_wind_vectors.m 2010-01-28 23:47:24 UTC (rev 4240)
@@ -27,10 +27,9 @@
% 'levels', levels, 'region', region, 'scalefactor', scalefactor);
%
-% Data Assimilation Research Testbed -- DART
-% Copyright 2004-2007, Data Assimilation Research Section
-% University Corporation for Atmospheric Research
-% Licensed under the GPL -- www.gpl.org/licenses/gpl.html
+%% DART software - Copyright \xA9 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
%
% <next few lines under version control, do not edit>
% $URL$
Property changes on: DART/trunk/diagnostics/matlab/plot_wind_vectors.m
___________________________________________________________________
Modified: svn:keywords
- Date Rev Author URL Id
+ Date Revision Author HeadURL Id
Copied: DART/trunk/diagnostics/matlab/read_obs_netcdf.m (from rev 4237, DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m)
===================================================================
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list