[Dart-dev] [3885] DART/trunk/observations/utilities/threed_sphere: The second graphic of bad observation locations is only generated
nancy at ucar.edu
nancy at ucar.edu
Fri May 22 16:39:28 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090522/6e7911fd/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m 2009-05-22 19:39:40 UTC (rev 3884)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m 2009-05-22 22:39:27 UTC (rev 3885)
@@ -13,127 +13,22 @@
%
% view(0,90); % for a traditional '2D' plot
-% record the user input
+% 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
-obsstruct.fname = fname;
-obsstruct.ObsTypeString = ObsTypeString;
-obsstruct.region = region;
-obsstruct.CopyString = CopyString;
-obsstruct.QCString = QCString;
-obsstruct.maxQC = maxQC;
-obsstruct.verbose = verbose;
+% Read the observation sequence
-% get going
+obsstruct = read_obs_netcdf(fname, ObsTypeString, region, ...
+ CopyString, QCString, maxQC, verbose);
-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 type.
-
-myind = strmatch(ObsTypeString,ObsTypeStrings);
-
-if ( isempty(myind) )
- error('no %s observations ... stopping',obsstruct.ObsTypeString)
-end
-
-mytypeind = get_copy_index(fname, CopyString);
-inds = find(obs_type == myind);
-mylocs = loc(inds,:);
-myobs = obs(inds,mytypeind);
-
-if ~ isempty(QCString)
- myQCind = get_qc_index(fname, QCString);
- myqc = qc(inds,myQCind);
-else
- myqc = [];
-end
-
-% 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);
-
-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);
-
- 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.
%-------------------------------------------------------------------------------
@@ -156,7 +51,7 @@
str1 = sprintf('%s level (%.2f - %.2f)',ObsTypeString,zmin,zmax);
str2 = sprintf('%s (%d locations)',CopyString,length(obsstruct.obs));
-str3 = sprintf('%s - %s',timestring(1,:),timestring(2,:));
+str3 = sprintf('%s - %s',obsstruct.timestring(1,:),obsstruct.timestring(2,:));
title( {str1, str3, str2}, 'Interpreter','none','FontSize',18);
xlabel('longitude')
@@ -183,58 +78,64 @@
% Create graphic of spatial distribution of 'bad' observations & their QC value.
%-------------------------------------------------------------------------------
-figure(2); clf
+if (obsstruct.numbadqc > 0 )
-subplot('position',[0.1 0.20 0.8 0.65])
-scalearray = 128 * ones(size(badobs.obs));
+ 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')
+
+ str2 = sprintf('%s (%d bad observations)',CopyString,length(badobs.obs));
+
+ title( {str1, str3, str2}, 'Interpreter','none','FontSize',18);
+ 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([xmin xmax ymin ymax zmin zmax])
+
+ myworldmap;
+ 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
-zmin = min(badobs.z);
-zmax = max(badobs.z);
-
-scatter3(badobs.lons, badobs.lats, badobs.z, scalearray, badobs.qc,'filled')
-
-str2 = sprintf('%s (%d bad observations)',CopyString,length(badobs.obs));
-
-title( {str1, str3, str2}, 'Interpreter','none','FontSize',18);
-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([xmin xmax ymin ymax zmin zmax])
-myworldmap;
-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
-
-
-
function h = myworldmap
%---------------------------------------------------------------------------
Modified: DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m 2009-05-22 19:39:40 UTC (rev 3884)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m 2009-05-22 22:39:27 UTC (rev 3885)
@@ -15,6 +15,17 @@
% record the user input
+% 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
obsstruct.fname = fname;
obsstruct.ObsTypeString = ObsTypeString;
obsstruct.region = region;
@@ -104,6 +115,7 @@
obsstruct.z = mylocs(inds,3);
obsstruct.obs = myobs(inds);
obsstruct.Ztyp = z_type(inds);
+obsstruct.numbadqc = 0;
if (isempty(myqc))
obsstruct.qc = [];
@@ -116,6 +128,8 @@
if ( (~ isempty(myqc)) & (~ isempty(maxQC)) )
inds = find(obsstruct.qc > maxQC);
+
+ obsstruct.numbadqc = length(inds);
if (~isempty(inds))
badobs.lons = obsstruct.lons(inds);
@@ -190,56 +204,58 @@
% Create graphic of spatial distribution of 'bad' observations & their QC value.
%-------------------------------------------------------------------------------
-figure(2); clf
+if (obsstruct.numbadqc > 0 )
-subplot('position',[0.1 0.20 0.8 0.65])
-scalearray = 128 * ones(size(badobs.obs));
+ 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',18);
+ 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;
+ 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
-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',18);
-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;
-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
-
-
-
function h = myworldmap
%---------------------------------------------------------------------------
Added: DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m (rev 0)
+++ DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m 2009-05-22 22:39:27 UTC (rev 3885)
@@ -0,0 +1,151 @@
+function obsstruct = read_obs_netcdf(fname, ObsTypeString, region, CopyString, ...
+ QCString, maxQC, verbose)
+%
+% 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';
+% maxQC = 2;
+% verbose = 1; % anything > 0 == 'true'
+%
+% obs = read_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxQC, verbose);
+%
+
+% 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
+%
+% <next few lines under version control, do not edit>
+% $URL$
+% $Id$
+% $Revision$
+% $Date$
+
+% record the user input
+
+obsstruct.fname = fname;
+obsstruct.ObsTypeString = ObsTypeString;
+obsstruct.region = region;
+obsstruct.CopyString = CopyString;
+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);
+
+obsstruct.timestring = timestring;
+
+% 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 type.
+
+myind = strmatch(ObsTypeString,ObsTypeStrings);
+
+if ( isempty(myind) )
+ error('no %s observations ... stopping',obsstruct.ObsTypeString)
+end
+
+mytypeind = get_copy_index(fname, CopyString);
+inds = find(obs_type == myind);
+mylocs = loc(inds,:);
+myobs = obs(inds,mytypeind);
+
+if ~ isempty(QCString)
+ myQCind = get_qc_index(fname, QCString);
+ myqc = qc(inds,myQCind);
+else
+ myqc = [];
+end
+
+% 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
+
Property changes on: DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Revision Author HeadURL Id
Name: svn:eol-style
+ native
More information about the Dart-dev
mailing list