[Dart-dev] [4356] DART/trunk:
nancy at ucar.edu
nancy at ucar.edu
Wed Apr 21 17:00:05 MDT 2010
Revision: 4356
Author: thoar
Date: 2010-04-21 17:00:04 -0600 (Wed, 21 Apr 2010)
Log Message:
-----------
The scripts now appropriately put the land mask at the top or bottom,
depending on the the 'which_vert' portion of the location AND some
knowledge of the origin of the observation. At present observations
of oceanic quantities (i.e. observations with a 'WOD observation' copy)
get plotted with a Z label of 'depth' even though their which_vert
is 'height' ... pretty tricksy. Atmospheric quantities are handled
correctly, as far as I have been able to test.
The figure with multiple panels is simpler, and the bottom graphic
is generally a plot of the observation value vs. the ensemble mean estimate.
If these things are alike, they line up on the diagonal.
The summary of how many obs got rejected, etc. is ONLY echoed to
the command window.
Modified Paths:
--------------
DART/trunk/diagnostics/matlab/linked_observations.m
DART/trunk/models/wrf/matlab/link_obs.m
Property Changed:
----------------
DART/trunk/models/wrf/matlab/link_obs.m
-------------- next part --------------
Modified: DART/trunk/diagnostics/matlab/linked_observations.m
===================================================================
--- DART/trunk/diagnostics/matlab/linked_observations.m 2010-04-21 21:53:29 UTC (rev 4355)
+++ DART/trunk/diagnostics/matlab/linked_observations.m 2010-04-21 23:00:04 UTC (rev 4356)
@@ -66,7 +66,7 @@
'YDataSource',ystring, ...
'ZDataSource',zstring);
-worldmap('light');
+myworldmap(obs);
xlabel(obs.colnames{obs.lonindex});
ylabel(obs.colnames{obs.latindex});
@@ -127,32 +127,34 @@
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');
+%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');
+fprintf('QC summary follows:\n')
LabelQC(obs.colnames{obs.qcindex}, obs.qc)
-refreshdata
-linkdata on
+%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);
+%figure3 = figure(3); clf(figure3);
-axes5 = axes('Parent',figure3,'OuterPosition',[0 0 1 0.95],'FontSize',18);
+% 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');
@@ -180,7 +182,7 @@
linkdata on
-function LabelQC(QCString, qcarray)
+function s = LabelQC(QCString, qcarray)
%% Create legend for (DART) QC values.
%
% 0 observation assimilated
@@ -215,10 +217,93 @@
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{:}),'FontSize',12)
+ % set(gca,'YTick',qcvals,'YAxisLocation','right')
+ % set(gca,'YTickLabel',char(s{:}),'FontSize',12)
otherwise,
str = sprintf('no way to interpret values of %s',strtrim(QCString));
text(0.0, 0.0, str)
end
+
+
+
+function h = myworldmap(obs)
+
+%%--------------------------------------------------------------------------
+% 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 = find(ax(1) <= lons, 1);
+lon_ind2 = find(ax(2) <= lons, 1);
+lat_ind1 = find(ax(3) <= lats, 1);
+lat_ind2 = find(ax(4) <= lats, 1);
+
+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" - 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;
+
+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));
+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;
+
+
Modified: DART/trunk/models/wrf/matlab/link_obs.m
===================================================================
--- DART/trunk/models/wrf/matlab/link_obs.m 2010-04-21 21:53:29 UTC (rev 4355)
+++ DART/trunk/models/wrf/matlab/link_obs.m 2010-04-21 23:00:04 UTC (rev 4356)
@@ -1,15 +1,60 @@
-function link_obs(fname,ObsTypeString)
+function link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
%% link_obs generates the 'brushable' observation plots.
%
+% Three figures will be generated. In order to make any sense of this
+%
+% Figure 1 will have a 3D geographic scatterplot.
+%
+% 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 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.
+%
+% link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
+%
+% fname - name of netcdf file that results from obs_seq_to_netcdf
+% ObsTypeString - the TYPE of observation (ncdump -v ObsTypesMetaData *.nc)
+% 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)
+%
% EXAMPLE 1:
-% fname = 'obs_sequence_013.nc';
+% fname = '/ptmp/thoar/POP/CAM/POP8/obs_epoch_001.nc';
+% ObsTypeString = 'APB_TEMPERATURE';
+% ObsCopyString = 'WOD observation';
+% CopyString = 'prior ensemble mean';
+% QCString = 'DART quality control';
+% region = [0 360 -90 90 -Inf Inf];
+% global obsmat;
+% link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
+%
+% EXAMPLE 2:
+% fname = 'obs_epoch_001.nc';
% ObsTypeString = 'RADIOSONDE_TEMPERATURE';
-% link_obs(fname,ObsTypeString)
+% ObsCopyString = 'NCEP BUFR observation';
+% CopyString = 'prior ensemble mean';
+% QCString = 'DART quality control';
+% region = [220 300 20 60 -Inf Inf];
+% global obsmat;
+% link_obs(fname, ObsTypeString, ObsCopyString, CopyString, QCString, region)
%
% EXAMPLE 3:
-% link_obs('obs_sequence_003.nc','RADIOSONDE_TEMPERATURE')
+% global obsmat;
+% region = [0 360 -90 90 -Inf Inf];
+% link_obs('obs_epoch_002.nc','RADIOSONDE_TEMPERATURE', 'observation', ...
+% 'prior ensemble member 3', 'DART quality control', region)
-%% 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
%
@@ -19,32 +64,65 @@
% $Revision$
% $Date$
-region = [0 360 -90 90 -Inf Inf];
-CopyString = 'observation';
-CopyString = 'NCEP BUFR observation';
-QCString = 'DART quality control';
-verbose = 1;
+if (exist(fname,'file') ~= 2)
+ error('%s does not exist.',fname)
+end
-obs = read_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, verbose);
+verbose = 1;
+obs = read_obs_netcdf(fname, ObsTypeString, region, ObsCopyString, QCString, verbose);
+copy = read_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, 0);
+
+if ( length(obs.lons) == 0 )
+ error('There are no %s observations in the region specified in %s', ObsTypeString, fname)
+end
+
+l1 = length(copy.obs);
+l2 = length(obs.obs);
+
+if (l1 ~= l2 )
+ error('there are %d %s and only %d %s',l1, CopyString, l2, ObsCopyString)
+end
+
+obs.ObsCopyString = obs.CopyString;
+obs.CopyString = copy.CopyString;
+
+%% Now pack the data in the same fashion as the cell array of column labels.
+
obs.lonindex = 1;
obs.latindex = 2;
obs.zindex = 3;
obs.obsindex = 4;
-obs.qcindex = 5;
-obs.keyindex = 6;
-obs.timeindex = 7;
-obs.indindex = 8;
+obs.copyindex = 5;
+obs.qcindex = 6;
+obs.keyindex = 7;
+obs.timeindex = 8;
+obs.indindex = 9;
global obsmat
-obsmat = zeros(length(obs.lons),5);
-obsmat(:,obs.lonindex ) = obs.lons;
-obsmat(:,obs.latindex ) = obs.lats;
-obsmat(:,obs.zindex ) = obs.z;
-obsmat(:,obs.obsindex ) = obs.obs;
-obsmat(:,obs.qcindex ) = obs.qc;
-obsmat(:,obs.keyindex ) = obs.keys;
-obsmat(:,obs.timeindex) = obs.time;
-obsmat(:,obs.indindex ) = [1:length(obs.time)];
+obsmat = zeros(length(obs.lons),9);
+obsmat(:,obs.lonindex ) = obs.lons; obs.colnames{obs.lonindex} = 'longitude';
+obsmat(:,obs.latindex ) = obs.lats; obs.colnames{obs.latindex} = 'latitude';
+obsmat(:,obs.zindex ) = obs.z ; obs.colnames{obs.zindex} = obs.Zunits;;
+obsmat(:,obs.obsindex ) = obs.obs ; obs.colnames{obs.obsindex} = ObsCopyString;
+obsmat(:,obs.copyindex) = copy.obs; obs.colnames{obs.copyindex} = CopyString;
+obsmat(:,obs.qcindex ) = obs.qc ; obs.colnames{obs.qcindex} = QCString;
+obsmat(:,obs.keyindex ) = obs.keys; obs.colnames{obs.keyindex} = 'obskey';
+obsmat(:,obs.timeindex) = obs.time; obs.colnames{obs.timeindex} = 'time';
+obsmat(:,obs.indindex ) = 1:length(obs.time); obs.colnames{obs.indindex} = 'index';
-linked_observations(obsmat,obs)
+%% Replace all ill-posed copies with
+% This should really check to see if the copy is a 'mean' or 'spread' ...
+
+iscalculated = ~isempty(findstr(lower(obs.colnames{obs.copyindex}),'mean')) | ...
+ ~isempty(findstr(lower(obs.colnames{obs.copyindex}),'spread'));
+
+if iscalculated
+ disp('replacing copies with [1 < QC flag < 5] with NaN')
+ allindices = obsmat(:,obs.qcindex);
+ badinds = ((allindices > 1) & (allindices < 5));
+ obsmat(badinds,obs.copyindex) = NaN;
+end
+
+%% create the linked plots
+linked_observations(obs)
Property changes on: DART/trunk/models/wrf/matlab/link_obs.m
___________________________________________________________________
Deleted: svn:executable
- *
More information about the Dart-dev
mailing list