[Dart-dev] [3787] DART/trunk/models/wrf/matlab/plot_wind_vectors.m: New, improved ...
nancy at ucar.edu
nancy at ucar.edu
Thu Mar 5 11:04:27 MST 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090305/0008fb67/attachment.html
-------------- next part --------------
Modified: DART/trunk/models/wrf/matlab/plot_wind_vectors.m
===================================================================
--- DART/trunk/models/wrf/matlab/plot_wind_vectors.m 2009-03-03 22:51:14 UTC (rev 3786)
+++ DART/trunk/models/wrf/matlab/plot_wind_vectors.m 2009-03-05 18:04:27 UTC (rev 3787)
@@ -1,13 +1,31 @@
-function obs = plot_wind_vectors(fname,ncname,platform,levelindex)
-% plot_wind_vectors(fname,ncname,platform,ilevel)
+function obs = plot_wind_vectors( fname, ncname, platform, varargin )
+% plot_wind_vectors creates maps with wind vector images overlain.
%
+% The required arguments are:
+% fname ... the file name containing the wind observation pairs
+% ncname ... the output netCDF file from obs_diag.f90 (for some metadata)
+% platform ... a string to represent the observation platform.
+%
+% The optional arguments are:
+% levels ... specifies the vertical area of interest. If not present,
+% all vertical levels are used.
+% region ... specifies that horizontal area of interest. If not present,
+% all available observations are used.
+% scalefactor ... provides control over the plotted size of the
+% wind vectors. A smaller number results in a
+% bigger wind vector. If not present, a value of 10.0 is
+% used.
+%
% fname = 'wind_vectors.006.dat';
% ncname = 'obs_diag_output.nc';
-% platform = 'SAT';
-% ilevel = 1; % surface == 1?
-% ilevel = -1; % all levels == -1
+% platform = 'RADIOSONDE';
+% levels = [1020 500];
+% region = [0 360 0 90]; %
+% scalefactor = 5; % reference arrow magnitude
%
-% obs = plot_wind_vectors(fname,ncname,platform,ilevel);
+% obs = plot_wind_vectors(fname, ncname, platform, ...
+% 'levels', levels, 'region', region, 'scalefactor', scalefactor);
+%
% Data Assimilation Research Testbed -- DART
% Copyright 2004-2007, Data Assimilation Research Section
@@ -20,32 +38,52 @@
% $Revision$
% $Date$
-scalefactor = 10.0;
+narg = nargin;
+if narg == 3
+
+ levels = [];
+ region = [];
+ scalefactor = 10;
+
+else
+
+ [levels, region, scalefactor] = parseinput(varargin{:});
+
+end
+
+data.filename = fname;
+data.ncname = ncname;
+data.platform = platform;
+data.levels = levels;
+data.region = region;
+data.scalefactor = scalefactor;
+
f = netcdf(ncname,'nowrite');
-platforms = f{'ObservationTypes'}(:);
-timeunits = f{'time'}.units(:);
+data.platforms = f{'ObservationTypes'}(:);
+data.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));
-obs = Load_Parse(fname, platforms, platform, levelindex );
+timebase = sscanf(data.timeunits,'%*s%*s%d%*c%d%*c%d'); % YYYY MM DD
+data.timeorigin = datenum(timebase(1),timebase(2),timebase(3));
+obs = Load_Parse(data);
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));
+ h = text(0.5,0.5,sprintf('%s has no %s data',data.filename,data.platform));
set(h,'Interpreter','none')
return
end
-obs.times = obs.times + timeorigin;
+obs.times = obs.times + data.timeorigin;
t1 = datestr(min(obs.times));
t2 = datestr(max(obs.times));
clf;
-axlims = DrawBackground(obs,scalefactor);
+axlims = DrawBackground( obs );
+
goodUV = find( (obs.Uqc < 1) & (obs.Vqc < 1));
baadUV = find( (obs.Uqc > 1) & (obs.Vqc > 1));
goodU = find( (obs.Uqc < 1) & (obs.Vqc > 0));
@@ -55,78 +93,88 @@
legstr = {};
if ~ isempty(goodUV)
- hgood = obs2plot(obs,goodUV,[0 0 0],scalefactor);
+ hgood = obs2plot(obs, goodUV, [0 0 0] );
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];
+ hbaadUV = obs2plot(obs, baadUV, [1 0 0] );
+ legh = [legh; hbaadUV];
legstr{length(legstr)+1} = sprintf('%d badUbadV',length(baadUV));
end
if ~ isempty(goodU)
- hgoodU = obs2plot(obs,goodU,[0 1 0],scalefactor);
+ hgoodU = obs2plot(obs, goodU, [0 1 0] );
legh = [legh; hgoodU];
legstr{length(legstr)+1} = sprintf('%d goodUbadV',length(goodU));
end
if ~ isempty(goodV)
- hgoodV = obs2plot(obs,goodV,[0 0 1],scalefactor);
+ hgoodV = obs2plot(obs, goodV, [0 0 1] );
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)});
+ sprintf('levels %s ',obs.levelstring)});
set(h,'FontSize',18)
-h = xlabel(fname); set(h,'Interpreter','none');
+h = xlabel(data.filename); set(h,'Interpreter','none');
legend(legh,legstr,'Location','NorthWestOutside')
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) ];
+function axlims = DrawBackground( obs )
+%======================================================================
+
+if isempty(obs.region)
+ % Figure out bounds of the data
+ axlims = [ min(obs.lon) max(obs.lon) min(obs.lat) max(obs.lat) ] ;
+else
+ axlims = obs.region;
+end
+
% It is nice to have a little padding around the perimeter
dx = 0.05 * (axlims(2) - axlims(1));
dy = 0.05 * (axlims(4) - axlims(3));
-% axis(axlims + [-dx dx -dy dy])
-axis([115 145 5 35])
+axis(axlims + [-dx dx -dy dy])
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;
+% Plot a wind vector pair of known length for reference
+tx = axlims(1)+dx;
+ty = axlims(4)-dy;
-U = 10/xfactor; V = 0/xfactor;
+U = 10/obs.scalefactor;
+V = 0/obs.scalefactor;
h = quiver(tx, ty, U, V, 0.0 ,'LineWidth',4.0,'Color','k');
set(h,'LineWidth',3.0)
-U = 0/xfactor; V = 10/xfactor;
+U = 0/obs.scalefactor; V = 10/obs.scalefactor;
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');
+h = text(tx, ty-0.1*dy,'10 m/s','VerticalAlignment','top');
-function h1 = obs2plot(obs, mask, colspec, xfactor)
-%----------------------------------------------------------------------
+function h1 = obs2plot(obs, mask, colspec )
+%======================================================================
+% Actually plots the wind vectors. to defeat Matlab's automatic
+% scaling, we use a scalefactor of 0.0 and manually scale the
+% data from the user input (or default).
+
lon = obs.lon(mask);
lat = obs.lat(mask);
-U = obs.U(mask)/xfactor;
-V = obs.V(mask)/xfactor;
+U = obs.U(mask)/obs.scalefactor;
+V = obs.V(mask)/obs.scalefactor;
h1 = quiver(lon, lat, U, V, 0.0);
h2 = plot(lon, lat, '.','MarkerSize',4);
set(h1,'Color',colspec)
@@ -134,22 +182,39 @@
-function newobs = Load_Parse(fname, platforms, platform, levelindex)
-%----------------------------------------------------------------------
+function obs = Load_Parse( data )
+%======================================================================
% Makes no attempt to find/replace/identify MISSING values
+%
+% data.filename incoming filename with the data
+% data.ncname incoming filename with the metadata
+% data.platform the observation platform of interest
+% data.levels the top/bottom levels of interest
+% data.region the region of interest [xmin xmax ymin ymax]
+% data.scalefactor the reference wind vector magnitude
+% data.platforms the observation platforms in the incoming file
+% data.timeunits the units string e.g. 'days since yyyy-mm-dd'
+% data.timeorigin
-obsmat = load(fname);
+obsmat = load(data.filename);
+% Find the types of data in the obs_sequence file for this epoch.
+
platformIDs = unique(obsmat(:,1));
uid = floor(platformIDs/100);
vid = platformIDs - uid*100;
-Ustrings = platforms(uid,:);
-Vstrings = platforms(vid,:);
+Ustrings = data.platforms(uid,:);
+Vstrings = data.platforms(vid,:);
nplatforms = length(uid);
pid = [];
-newobs = [];
+obs = [];
+levelstring = [];
+regionstring = [];
+% This block divines the platform string and companion obs_kinds
+% from the netCDF file metadata - the only reason we need the netCDF file.
+
for i = 1:nplatforms
uindex = findstr(Ustrings(i,:),'_U_WIND_COMPONENT') - 1;
vindex = findstr(Vstrings(i,:),'_V_WIND_COMPONENT') - 1;
@@ -161,77 +226,207 @@
error('U and V wind component platforms do not match')
end
- if strcmp(lower(Ubase),lower(platform))
+ % Determine what numeric pid corresponds to the desired platform string
+ if strcmp(lower(Ubase),lower(data.platform))
pid = platformIDs(i);
end
+ % echo a little informational statement about how many obs of
+ % this type there are in this assimilation period.
inds = find(obsmat(:,1) == platformIDs(i));
nobs = length(inds);
- disp(sprintf('%6d %14s observations in %s (%4d)',nobs,Ubase,fname,platformIDs(i)))
+ disp(sprintf('%6d %14s observations in %s (%4d)',nobs,Ubase,data.filename,platformIDs(i)))
end
+% This block extracts just the desired observations based on platform.
+
if isempty(pid)
- disp(sprintf('no %s observations in %s',platform,fname))
+ disp(sprintf('no %s observations in %s', data.platform, data.filename))
return
end
inds = find(obsmat(:,1) == pid);
if isempty(inds)
- disp(sprintf('no observations of type %d in %s',pid,fname))
+ disp(sprintf('no %s observations (type %d) in %s',data.platform, pid, data.filename))
+ return
end
-obs.platform = obsmat(inds, 1);
-obs.day = obsmat(inds, 2);
-obs.seconds = obsmat(inds, 3);
-obs.lon = obsmat(inds, 4);
-obs.lat = obsmat(inds, 5);
-obs.level = obsmat(inds, 6);
-obs.levind = obsmat(inds, 7);
-obs.Uqc = obsmat(inds, 8);
-obs.Vqc = obsmat(inds, 9);
-obs.U = obsmat(inds,10);
-obs.V = obsmat(inds,11);
-obs.Upr = obsmat(inds,12);
-obs.Vpr = obsmat(inds,13);
-obs.Upo = obsmat(inds,14);
-obs.Vpo = obsmat(inds,15);
-obs.times = obs.day + obs.seconds/86400;
+platform = obsmat(inds, 1);
+day = obsmat(inds, 2);
+seconds = obsmat(inds, 3);
+lon = obsmat(inds, 4);
+lat = obsmat(inds, 5);
+level = obsmat(inds, 6);
+levind = obsmat(inds, 7);
+Uqc = obsmat(inds, 8);
+Vqc = obsmat(inds, 9);
+U = obsmat(inds,10);
+V = obsmat(inds,11);
+Upr = obsmat(inds,12);
+Vpr = obsmat(inds,13);
+Upo = obsmat(inds,14);
+Vpo = obsmat(inds,15);
+times = day + seconds/86400;
-% Sort by level
+%--------------------------------------------------
+% Subset the levels of interest
+%--------------------------------------------------
-if ( levelindex > 0 )
- levelinds = find(obs.levind == levelindex);
+if ( isempty(data.levels) )
+ % then we want all levels, do nothing ...
+ level1 = min(level);
+ levelN = max(level);
+
+ levelstring = sprintf('all (%.2f to %.2f)',level1,levelN);
else
- levelinds = 1:length(inds);
+ level1 = min(data.levels);
+ levelN = max(data.levels);
+
+ levelstring = sprintf('%.2f to %.2f',level1,levelN);
+
+ inds = find ((level >= level1) & (level <= levelN));
+
+ if (length(inds) == 0)
+ disp(sprintf('no %s observations in %s', data.platform, levelstring))
+ return
+ end
+
+ platform = platform(inds);
+ day = day(inds);
+ seconds = seconds(inds);
+ lon = lon(inds);
+ lat = lat(inds);
+ level = level(inds);
+ levind = levind(inds);
+ Uqc = Uqc(inds);
+ Vqc = Vqc(inds);
+ U = U(inds);
+ V = V(inds);
+ Upr = Upr(inds);
+ Vpr = Vpr(inds);
+ Upo = Upo(inds);
+ Vpo = Vpo(inds);
+ times = times(inds);
end
+%--------------------------------------------------
+% Subset the region of interest
+% for the moment, we are not supporting wrapping at Greenwich.
+%--------------------------------------------------
-if isempty(levelinds)
- error('no observations of type %d at level %d',pid,levelindex)
+if isempty(data.region)
+ % then we want the entire dataset, do nothing ...
+ regionstring = 'global';
else
- disp(sprintf('%d observations of type %d at level %d',length(levelinds),pid,levelindex))
+ lon1 = data.region(1);
+ lonN = data.region(2);
+ lat1 = data.region(3);
+ latN = data.region(4);
+
+ regionstring = sprintf('(%.2f -> %.2f, %.2f -> %.2f)',lon1,lonN,lat1,latN);
+
+ inds = find ((lon >= lon1) & (lon <= lonN) & ...
+ (lat >= lat1) & (lat <= latN));
+
+ if (length(inds) == 0)
+ disp(sprintf('no %s observations in %s', data.platform, regionstring))
+ return
+ end
+
+ platform = platform(inds);
+ day = day(inds);
+ seconds = seconds(inds);
+ lon = lon(inds);
+ lat = lat(inds);
+ level = level(inds);
+ levind = levind(inds);
+ Uqc = Uqc(inds);
+ Vqc = Vqc(inds);
+ U = U(inds);
+ V = V(inds);
+ Upr = Upr(inds);
+ Vpr = Vpr(inds);
+ Upo = Upo(inds);
+ Vpo = Vpo(inds);
+ times = times(inds);
end
-clear newobs
+%--------------------------------------------------
+% Insert NaN's for missing values
+%--------------------------------------------------
-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);
+U( U < -900) = NaN;
+V( V < -900) = NaN;
+Upr( Upr < -900) = NaN;
+Vpr( Vpr < -900) = NaN;
+Upo( Upo < -900) = NaN;
+Vpo( Vpo < -900) = NaN;
-clear obsmat obs
+%--------------------------------------------------
+% Create the output structure.
+%--------------------------------------------------
+
+obs.platform = platform;
+obs.day = day;
+obs.seconds = seconds;
+obs.lon = lon;
+obs.lat = lat;
+obs.level = level;
+obs.levind = levind;
+obs.Uqc = Uqc;
+obs.Vqc = Vqc;
+obs.U = U;
+obs.V = V;
+obs.Upr = Upr;
+obs.Vpr = Vpr;
+obs.Upo = Upo;
+obs.Vpo = Vpo;
+obs.times = times;
+obs.levels = data.levels;
+obs.region = data.region;
+obs.scalefactor = data.scalefactor;
+obs.levelstring = levelstring;
+
+
+function [levels, region, scalefactor] = parseinput(varargin)
+%======================================================================
+% try to parse the input pairs ... which must be pairs
+
+if (mod(length(varargin),2) ~= 0)
+ error('Wrong number (%d) of optional arguments. Must be parameter/value pairs: ''levels'',[1000 500]',length(varargin))
+end
+
+npairs = length(varargin)/2;
+
+levels = [];
+region = [];
+scalefactor = [];
+
+for i = 1:2:length(varargin)
+ switch lower(varargin{i})
+ case 'levels'
+ levels = varargin{i+1};
+ case 'region'
+ region = varargin{i+1};
+ case 'scalefactor'
+ scalefactor = varargin{i+1};
+ otherwise
+ error('Unknown parameter %s',varargin{i})
+ end
+end
+
+% Make sure the levels array has a top/bottom
+if ~ isempty(levels)
+end
+
+% Make sure the geographic array makes sense
+if ~ isempty(region)
+end
+
+% Make sure the scalefactor makes sense
+if ~ isempty(scalefactor)
+end
+
More information about the Dart-dev
mailing list