[Dart-dev] [3868] DART/trunk/observations/utilities/threed_sphere: Has landmasses at the surface of the planet, more faithful axis limits.
nancy at ucar.edu
nancy at ucar.edu
Fri May 8 12:05:13 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090508/89864b4a/attachment.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-08 13:48:56 UTC (rev 3867)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf.m 2009-05-08 18:05:13 UTC (rev 3868)
@@ -1,22 +1,22 @@
-function obsstruct = plot_obs_netcdf(fname, ObsTypeString, region, ObsString, ...
+function h = plot_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];
-% ObsString = 'NCEP BUFR observation';
+% CopyString = 'NCEP BUFR observation';
% QCString = 'DART quality control';
% maxQC = 2;
% verbose = 1; % anything > 0 == 'true'
%
-% bob = plot_obs_netcdf(fname, ObsTypeString, region, ObsString, QCString, maxQC, verbose);
+% bob = plot_obs_netcdf(fname, ObsTypeString, region, CopyString, QCString, maxQC, verbose);
% record the user input
obsstruct.fname = fname;
obsstruct.ObsTypeString = ObsTypeString;
obsstruct.region = region;
-obsstruct.ObsString = ObsString;
+obsstruct.CopyString = CopyString;
obsstruct.QCString = QCString;
obsstruct.maxQC = maxQC;
obsstruct.verbose = verbose;
@@ -30,7 +30,7 @@
t = nc_varget(fname,'time');
obs_type = nc_varget(fname,'obs_type');
-z = nc_varget(fname,'which_vert');
+z_type = nc_varget(fname,'which_vert');
loc = nc_varget(fname,'location');
obs = nc_varget(fname,'observations');
@@ -69,7 +69,7 @@
% Find observations of the correct type.
-mytypeind = get_copy_index(fname, ObsString);
+mytypeind = get_copy_index(fname, CopyString);
myind = strmatch(ObsTypeString,ObsTypeStrings);
inds = find(obs_type == myind);
mylocs = loc(inds,:);
@@ -90,19 +90,29 @@
obsstruct.lats = mylocs(inds,2);
obsstruct.z = mylocs(inds,3);
obsstruct.obs = myobs(inds);
+obsstruct.Ztyp = z_type(inds);
-if ( ~ isempty(myqc))
-obsstruct.qc = myqc(inds);
+if (isempty(myqc))
+ obsstruct.qc = [];
+else
+ obsstruct.qc = myqc(inds);
end
-% It might be great to have a histogram of the observations with particular QC
-% values
-
% subset based on qc value
-if ( (~ isempty(myqc)) & ( ~ isempty(maxQC)) )
+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))
@@ -110,11 +120,19 @@
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));
@@ -122,24 +140,167 @@
zmin = min(obsstruct.z);
zmax = max(obsstruct.z);
-y1 = 36;
-yN = 10*y1;
-x1 = min(obsstruct.obs);
-xN = max(obsstruct.obs);
+y1 = 36;
+yN = 10*y1;
+x1 = min(obsstruct.obs);
+xN = max(obsstruct.obs);
slope = (yN-y1)/(xN-x1);
-b = y1 - slope*x1;
-
+b = y1 - slope*x1;
scalarray = obsstruct.obs*slope + b;
-h = scatter(obsstruct.lons,obsstruct.lats,scalarray,obsstruct.obs,'d','filled');
-axis image
-% axis([xmin xmax -Inf Inf])
-worldmap;
-colorbar;
+scatter3(obsstruct.lons, obsstruct.lats, obsstruct.z, ...
+ scalarray, obsstruct.obs, 'd', 'filled');
+axis([xmin xmax ymin ymax zmin zmax])
+
str1 = sprintf('%s level (%.2f - %.2f)',ObsTypeString,zmin,zmax);
str2 = sprintf('(%d locations)',length(obsstruct.obs));
str3 = sprintf('%s - %s',timestring(1,:),timestring(2,:));
-title( {str1, str3}, 'Interpreter','none','FontSize',18);
-xlabel(str2,'FontSize',16)
+title( {str1, str3, str2}, 'Interpreter','none','FontSize',18);
+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;
+h = colorbar;
+set(get(h,'YLabel'),'String',ObsTypeString,'Interpreter','none')
+
+%-------------------------------------------------------------------------------
+% Create graphic of spatial distribution of 'bad' observations & their QC value.
+%-------------------------------------------------------------------------------
+
+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([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
+
+text(0.0,0.0,s{1})
+text(0.0,0.5,s{2})
+
+
+
+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;
+
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-08 13:48:56 UTC (rev 3867)
+++ DART/trunk/observations/utilities/threed_sphere/plot_obs_netcdf_diffs.m 2009-05-08 18:05:13 UTC (rev 3868)
@@ -1,16 +1,16 @@
function obsstruct = plot_obs_netcdf_diffs(fname, ObsTypeString, region, ...
- ObsString1, ObsString2, QCString, maxQC, verbose)
+ CopyString1, CopyString2, QCString, maxQC, verbose)
%
% fname = 'obs_sequence_001.nc';
% ObsTypeString = 'RADIOSONDE_U_WIND_COMPONENT';
% region = [0 360 -90 90 -Inf Inf];
-% ObsString1 = 'NCEP BUFR observation';
-% ObsString2 = 'prior ensemble mean';
+% 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, ObsString1, ObsString2, ...
+% bob = plot_obs_netcdf_diffs(fname, ObsTypeString, region, CopyString1, CopyString2, ...
% QCString, maxQC, verbose);
% record the user input
@@ -18,8 +18,8 @@
obsstruct.fname = fname;
obsstruct.ObsTypeString = ObsTypeString;
obsstruct.region = region;
-obsstruct.ObsString1 = ObsString1;
-obsstruct.ObsString2 = ObsString2;
+obsstruct.CopyString1 = CopyString1;
+obsstruct.CopyString2 = CopyString2;
obsstruct.QCString = QCString;
obsstruct.maxQC = maxQC;
obsstruct.verbose = verbose;
@@ -33,7 +33,7 @@
t = nc_varget(fname,'time');
obs_type = nc_varget(fname,'obs_type');
-z = nc_varget(fname,'which_vert');
+z_type = nc_varget(fname,'which_vert');
loc = nc_varget(fname,'location');
obs = nc_varget(fname,'observations');
@@ -72,8 +72,8 @@
% Find observations of the correct types.
-mytype1 = get_copy_index(fname, ObsString1);
-mytype2 = get_copy_index(fname, ObsString2);
+mytype1 = get_copy_index(fname, CopyString1);
+mytype2 = get_copy_index(fname, CopyString2);
myind = strmatch(ObsTypeString,ObsTypeStrings);
inds = find(obs_type == myind);
@@ -99,19 +99,29 @@
obsstruct.lats = mylocs(inds,2);
obsstruct.z = mylocs(inds,3);
obsstruct.obs = myobs(inds);
+obsstruct.Ztyp = z_type(inds);
-if ( ~ isempty(myqc))
-obsstruct.qc = myqc(inds);
+if (isempty(myqc))
+ obsstruct.qc = [];
+else
+ obsstruct.qc = myqc(inds);
end
-% It might be great to have a histogram of the observations with particular QC
-% values
-
% subset based on qc value
-if ( (~ isempty(myqc)) & ( ~ isempty(maxQC)) )
+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))
@@ -119,11 +129,19 @@
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));
@@ -131,24 +149,169 @@
zmin = min(obsstruct.z);
zmax = max(obsstruct.z);
-y1 = 24;
-yN = 10*y1;
-x1 = min(obsstruct.obs);
-xN = max(obsstruct.obs);
+y1 = 36;
+yN = 10*y1;
+x1 = min(obsstruct.obs);
+xN = max(obsstruct.obs);
slope = (yN-y1)/(xN-x1);
-b = y1 - slope*x1;
+b = y1 - slope*x1;
scalarray = obsstruct.obs*slope + b;
+scalearray = 128 * ones(size(obsstruct.obs));
-h = scatter(obsstruct.lons,obsstruct.lats,128,obsstruct.obs,'d','filled');
-axis image
-% axis([xmin xmax -Inf Inf])
-worldmap;
-colorbar;
+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)',ObsString2,ObsString1,length(obsstruct.obs));
+str2 = sprintf('%s - %s (%d locations)',CopyString2,CopyString1,length(obsstruct.obs));
str3 = sprintf('%s - %s',timestring(1,:),timestring(2,:));
-title( {str1, str3}, 'Interpreter','none','FontSize',18);
-xlabel(str2,'FontSize',16)
+title( {str1, str3, str2}, 'Interpreter','none','FontSize',18);
+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;
+h = colorbar;
+set(get(h,'YLabel'),'String',ObsTypeString,'Interpreter','none')
+
+%-------------------------------------------------------------------------------
+% Create graphic of spatial distribution of 'bad' observations & their QC value.
+%-------------------------------------------------------------------------------
+
+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
+
+text(0.0,0.0,s{1})
+text(0.0,0.5,s{2})
+
+
+
+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;
+
More information about the Dart-dev
mailing list