[Dart-dev] [4241] DART/trunk/observations/utilities/threed_sphere: Renamed these files into DART/diagnostics/matlab/
nancy at ucar.edu
nancy at ucar.edu
Thu Jan 28 16:48:52 MST 2010
Revision: 4241
Author: thoar
Date: 2010-01-28 16:48:52 -0700 (Thu, 28 Jan 2010)
Log Message:
-----------
Renamed these files into DART/diagnostics/matlab/
Removed Paths:
-------------
DART/trunk/observations/utilities/threed_sphere/linked_observations.m
DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m
DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m
DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m
-------------- next part --------------
Deleted: DART/trunk/observations/utilities/threed_sphere/linked_observations.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/linked_observations.m 2010-01-28 23:47:24 UTC (rev 4240)
+++ DART/trunk/observations/utilities/threed_sphere/linked_observations.m 2010-01-28 23:48:52 UTC (rev 4241)
@@ -1,162 +0,0 @@
-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
Deleted: DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m 2010-01-28 23:47:24 UTC (rev 4240)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m 2010-01-28 23:48:52 UTC (rev 4241)
@@ -1,350 +0,0 @@
-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
Deleted: DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m 2010-01-28 23:47:24 UTC (rev 4240)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m 2010-01-28 23:48:52 UTC (rev 4241)
@@ -1,345 +0,0 @@
-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;
-
Deleted: DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m
===================================================================
--- DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m 2010-01-28 23:47:24 UTC (rev 4240)
+++ DART/trunk/observations/utilities/threed_sphere/read_obs_netcdf.m 2010-01-28 23:48:52 UTC (rev 4241)
@@ -1,187 +0,0 @@
-function obsstruct = read_obs_netcdf(fname, ObsTypeString, region, CopyString, ...
- QCString, maxQC, verbose)
-%% read_obs_netcdf reads in the netcdf flavor observation sequence file
-% and returns a subsetted structure.
-%
-% fname = 'obs_sequence_001.nc';
-% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT'; % or 'ALL' ...
-% 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);
-%
-% The return variable 'obs' is a structure. As an example ...
-%
-% 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
-% timestring: [2x20 char]
-% lons: [2343x1 double]
-% lats: [2343x1 double]
-% z: [2343x1 double]
-% obs: [2343x1 double]
-% Ztyp: [2343x1 double]
-% numbadqc: 993
-% qc: [2343x1 double]
-% badobs: [1x1 struct]
-
-%% 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
-
-%% 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);
-
- fprintf('N = %6d %s obs (type %3d) between levels %.2f and %.2f\n', ...
- 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.
-% If 'ALL' is requested ... do not subset.
-
-mytypeind = get_copy_index(fname, CopyString);
-
-switch lower(ObsTypeString)
- case 'all'
- inds = 1:size(obs,1);
-
- otherwise % subset the desired observation type
- myind = strmatch(ObsTypeString, ObsTypeStrings);
- if ( isempty(myind) )
- error('no %s observations ... stopping',obsstruct.ObsTypeString)
- end
- inds = find(obs_type == myind);
-end
-
-mylocs = loc(inds,:);
-myobs = obs(inds,mytypeind);
-
-%% Find desired QC values of those observations
-
-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.qc = [];
-obsstruct.numbadqc = 0;
-
-if ~ isempty(myqc)
- 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
-
- fprintf('Removing %d obs with a %s value greater than %f\n', ...
- 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
-
-if ( exist('badobs','var') )
- obsstruct.badobs = badobs;
-end
-
More information about the Dart-dev
mailing list