[Dart-dev] [3779] DART/trunk/models/wrf/matlab: Modified the plot_wind_vectors function to plot at a fixed
nancy at ucar.edu
nancy at ucar.edu
Thu Feb 26 15:59:12 MST 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090226/a4da941a/attachment.html
-------------- next part --------------
Added: DART/trunk/models/wrf/matlab/PlotWindObs.m
===================================================================
--- DART/trunk/models/wrf/matlab/PlotWindObs.m (rev 0)
+++ DART/trunk/models/wrf/matlab/PlotWindObs.m 2009-02-26 22:59:12 UTC (rev 3779)
@@ -0,0 +1,27 @@
+%
+
+% 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$
+
+for periods = 2:37;
+
+ fname = sprintf('wind_vectors.%03d.dat',periods);
+ ncname = 'obs_diag_output.nc';
+ platform = 'SAT';
+ level = -1;
+
+ obs = plot_wind_vectors(fname,ncname,platform,level);
+
+ disp('Pausing, hit any key to continue ...')
+ pause
+
+end
+
Property changes on: DART/trunk/models/wrf/matlab/PlotWindObs.m
___________________________________________________________________
Name: svn:mime-type
+ text/x-matlab
Name: svn:keywords
+ Date Rev Author URL Id
Name: svn:eol-style
+ native
Modified: DART/trunk/models/wrf/matlab/plot_wind_vectors.m
===================================================================
--- DART/trunk/models/wrf/matlab/plot_wind_vectors.m 2009-02-26 00:21:36 UTC (rev 3778)
+++ DART/trunk/models/wrf/matlab/plot_wind_vectors.m 2009-02-26 22:59:12 UTC (rev 3779)
@@ -1,20 +1,13 @@
-function obs = plot_wind_vectors(fname,ncname,platform)
-% plot_wind_vectors(fname,ncname,platform)
+function obs = plot_wind_vectors(fname,ncname,platform,levelindex)
+% plot_wind_vectors(fname,ncname,platform,ilevel)
%
-% 102 RADIOSONDE
-% 2223 SAT
-% 4243 BUOY
-% 4748 SHIP
-% 5758 AIREP
-% 6768 PILOT
-% 8182 QKSWND
+% fname = 'wind_vectors.006.dat';
+% ncname = 'obs_diag_output.nc';
+% platform = 'SAT';
+% ilevel = 1; % surface == 1?
+% ilevel = -1; % all levels == -1
%
-% fname = 'wind_vectors.006.dat';
-% ncname = 'obs_diag_output.nc';
-% platform = 8182; % QKSWND
-% ilevel = 1; % surface == 1?
-%
-% plot_wind_vectors(fname,ncname,platform)
+% obs = plot_wind_vectors(fname,ncname,platform,ilevel);
% Data Assimilation Research Testbed -- DART
% Copyright 2004-2007, Data Assimilation Research Section
@@ -27,54 +20,122 @@
% $Revision$
% $Date$
+scalefactor = 10.0;
+
f = netcdf(ncname,'nowrite');
platforms = f{'ObservationTypes'}(:);
timeunits = f{'time'}.units(:);
close(f);
-timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
-timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+timebase = sscanf(timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+obs = Load_Parse(fname, platforms, platform, levelindex );
-obs = Load_Parse(fname, platforms, platform );
+if isempty(obs)
+ clf;
+ axis([0 1 0 1]); axis image
+ h = text(0.5,0.5,sprintf('%s has no %s data',fname,platform));
+ set(h,'Interpreter','none')
+ return
+end
-obs.times = obs.times + timeorigin;
-t1 = datestr(min(obs.times));
-t2 = datestr(max(obs.times));
+obs.times = obs.times + timeorigin;
+t1 = datestr(min(obs.times));
+t2 = datestr(max(obs.times));
clf;
-axlims = DrawBackground(obs);
-h1 = quiver(obs.lon, obs.lat, obs.U, obs.V);
-h2 = plot(obs.lon, obs.lat, '.');
-set(h2,'MarkerSize',8)
-hold off;
+axlims = DrawBackground(obs,scalefactor);
-h = title({sprintf('%s %s %s',t1,platform,t2),'all levels - including 99s bad UVs'});
+goodUV = find( (obs.Uqc < 1) & (obs.Vqc < 1));
+baadUV = find( (obs.Uqc > 1) & (obs.Vqc > 1));
+goodU = find( (obs.Uqc < 1) & (obs.Vqc > 0));
+goodV = find( (obs.Vqc < 1) & (obs.Uqc > 0));
+
+legh = [];
+legstr = {};
+
+if ~ isempty(goodUV)
+ hgood = obs2plot(obs,goodUV,[0 0 0],scalefactor);
+ legh = [legh; hgood];
+ legstr{length(legstr)+1} = sprintf('%d good',length(goodUV));
+end
+
+if ~ isempty(baadUV)
+ hbaadUV = obs2plot(obs,baadUV,[1 0 0],scalefactor);
+ legh = [legh; hbaadUV];
+ legstr{length(legstr)+1} = sprintf('%d badUbadV',length(baadUV));
+end
+
+if ~ isempty(goodU)
+ hgoodU = obs2plot(obs,goodU,[0 1 0],scalefactor);
+ legh = [legh; hgoodU];
+ legstr{length(legstr)+1} = sprintf('%d goodUbadV',length(goodU));
+end
+
+if ~ isempty(goodV)
+ hgoodV = obs2plot(obs,goodV,[0 0 1],scalefactor);
+ legh = [legh; hgoodV];
+ legstr{length(legstr)+1} = sprintf('%d badUgoodV',length(goodV));
+end
+
+h = title({sprintf('%s %s %s',t1,platform,t2), ...
+ sprintf('binlevel %d',levelindex)});
set(h,'FontSize',18)
-xlabel(fname)
+h = xlabel(fname); set(h,'Interpreter','none');
+legend(legh,legstr,'Location','NorthWestOutside')
-function axlims = DrawBackground(obs, platform)
+hold off;
+
+function axlims = DrawBackground(obs, xfactor)
%----------------------------------------------------------------------
% Figure out bounds of the data
axlims = [ min(obs.lon) max(obs.lon) min(obs.lat) max(obs.lat) ];
% It is nice to have a little padding around the perimeter
-dx = 0.025 * (axlims(2) - axlims(1));
-dy = 0.025 * (axlims(4) - axlims(3));
+dx = 0.05 * (axlims(2) - axlims(1));
+dy = 0.05 * (axlims(4) - axlims(3));
-axis(axlims + [-dx dx -dy dy])
+% axis(axlims + [-dx dx -dy dy])
+axis([115 145 5 35])
+axis image
+axlims = axis;
% It is nice to know where the land is
worldmap('light');
hold on;
+tx = axlims(1)+2;
+ty = axlims(4)-2;
+U = 10/xfactor; V = 0/xfactor;
+h = quiver(tx, ty, U, V, 0.0 ,'LineWidth',4.0,'Color','k');
+set(h,'LineWidth',3.0)
+U = 0/xfactor; V = 10/xfactor;
+h = quiver(tx, ty, U, V, 0.0 ,'LineWidth',4.0,'Color','k');
+set(h,'LineWidth',3.0)
+h = text(tx, ty,'10 m/s','VerticalAlignment','top');
-function obs = Load_Parse(fname, platforms, platform)
+
+
+function h1 = obs2plot(obs, mask, colspec, xfactor)
%----------------------------------------------------------------------
+lon = obs.lon(mask);
+lat = obs.lat(mask);
+U = obs.U(mask)/xfactor;
+V = obs.V(mask)/xfactor;
+h1 = quiver(lon, lat, U, V, 0.0);
+h2 = plot(lon, lat, '.','MarkerSize',4);
+set(h1,'Color',colspec)
+set(h2,'Color',colspec)
+
+
+
+function newobs = Load_Parse(fname, platforms, platform, levelindex)
+%----------------------------------------------------------------------
% Makes no attempt to find/replace/identify MISSING values
obsmat = load(fname);
@@ -87,6 +148,7 @@
nplatforms = length(uid);
pid = [];
+newobs = [];
for i = 1:nplatforms
uindex = findstr(Ustrings(i,:),'_U_WIND_COMPONENT') - 1;
@@ -99,22 +161,26 @@
error('U and V wind component platforms do not match')
end
- disp(sprintf('Available platforms in %s are %s',fname,Ubase))
-
if strcmp(lower(Ubase),lower(platform))
pid = platformIDs(i);
end
+ inds = find(obsmat(:,1) == platformIDs(i));
+ nobs = length(inds);
+
+ disp(sprintf('%6d %14s observations in %s (%4d)',nobs,Ubase,fname,platformIDs(i)))
+
end
if isempty(pid)
- error('no %s observations in %s',platform,fname)
+ disp(sprintf('no %s observations in %s',platform,fname))
+ return
end
inds = find(obsmat(:,1) == pid);
if isempty(inds)
- error('no observations of type %d in %s',pid,fname)
+ disp(sprintf('no observations of type %d in %s',pid,fname))
end
obs.platform = obsmat(inds, 1);
@@ -134,4 +200,38 @@
obs.Vpo = obsmat(inds,15);
obs.times = obs.day + obs.seconds/86400;
-clear obsmat
+% Sort by level
+
+if ( levelindex > 0 )
+ levelinds = find(obs.levind == levelindex);
+else
+ levelinds = 1:length(inds);
+end
+
+
+if isempty(levelinds)
+ error('no observations of type %d at level %d',pid,levelindex)
+else
+ disp(sprintf('%d observations of type %d at level %d',length(levelinds),pid,levelindex))
+end
+
+clear newobs
+
+newobs.platform = obs.platform(levelinds);
+newobs.day = obs.day( levelinds);
+newobs.seconds = obs.seconds( levelinds);
+newobs.lon = obs.lon( levelinds);
+newobs.lat = obs.lat( levelinds);
+newobs.level = obs.level( levelinds);
+newobs.levind = obs.levind( levelinds);
+newobs.Uqc = obs.Uqc( levelinds);
+newobs.Vqc = obs.Vqc( levelinds);
+newobs.U = obs.U( levelinds);
+newobs.V = obs.V( levelinds);
+newobs.Upr = obs.Upr( levelinds);
+newobs.Vpr = obs.Vpr( levelinds);
+newobs.Upo = obs.Upo( levelinds);
+newobs.Vpo = obs.Vpo( levelinds);
+newobs.times = obs.times( levelinds);
+
+clear obsmat obs
More information about the Dart-dev
mailing list