-
Greg Quinn authoredGreg Quinn authored
plot.py 5.53 KiB
import matplotlib
import numpy
import sys
from matplotlib.backends.backend_agg import FigureCanvasAgg
from matplotlib.figure import Figure
from mpl_toolkits.basemap import Basemap
from pyhdf.SD import SD
# helper function for creating plots
def plot_helper(out_file, data,
x_range,
title, x_label, y_label=None,
data_labels=None):
# plot setup; params determined via trial and error
fig = Figure((4, 2))
ax = fig.add_axes((0.15, 0.15, 0.8, 0.7))
# plot the given data
x = numpy.linspace(x_range[0], x_range[1], len(data[0]))
for y in data:
ax.plot(x, y)
# add labels and such
ax.set_xlim(x_range[0], x_range[1])
ax.set_title(title)
ax.set_xlabel(x_label)
if y_label:
ax.set_ylabel(y_label)
if data_labels:
ax.legend(data_labels)
# write out the image file
canvas = FigureCanvasAgg(fig)
canvas.print_figure(out_file)
# helper for ploting a single line based on the given SDS
def plot_sds(sds, band, out_file, **kwargs):
# just grab the data and x-axis limits and pass everything to
# plot_helper
data = (sds[band],)
x_range = (sds.min_bucket, sds.max_bucket)
plot_helper(out_file, data, x_range, **kwargs)
# helper function for creating map plots
def map_helper(out_file, data, title):
# imshow doesn't deal well with NaNs, so mask them off
data = numpy.ma.masked_array(data, mask=numpy.isnan(data))
# plot setup; params determined via trial and error
fig = Figure((5, 2.25))
ax = fig.add_axes((0.01, 0.01, 0.89, 0.88))
color_ax = fig.add_axes((0.91, 0.05, 0.02, 0.90))
# use basemap to plot the given data and draw geographical
# boundaries
bm = Basemap(projection='cyl', ax=ax)
data = bm.transform_scalar(data,
numpy.arange(-179.5, 180.0, 1.0),
numpy.arange(-89.5, 90.0, 1.0),
360, 180)
im = bm.imshow(data)
bm.drawmapboundary()
bm.drawcoastlines()
# add the given title and a colorbar
ax.set_title(title)
fig.colorbar(im, cax=color_ax)
# write out the image file
canvas = FigureCanvasAgg(fig)
canvas.print_figure(out_file)
# plot histograms of all BT samples for both AIRS and MODIS
def plot_bt_hist(sd, band, out_file):
# grab the data and x-axis limits and invoke plot_helper
airs_sds = sd.select('AIRS_BT_Hist')
modis_sds = sd.select('MODIS_BT_Hist')
data = (airs_sds[band], modis_sds[band])
x_range = (airs_sds.min_bucket, airs_sds.max_bucket)
plot_helper(out_file,
data,
x_range,
title='BT Histogram',
x_label='BT (K)',
data_labels=('AIRS', 'MODIS'))
# plot histogram of BT differences
def plot_bt_diff_hist(sd, band, out_file):
sds = sd.select('BT_Diff_Hist')
plot_sds(sds, band, out_file,
title='BT Difference Histogram',
x_label='AIRS-MODIS BT Diff (K)')
# plot mean BT difference as a function of BT
def plot_bt_diff_by_airs_bt(sd, band, out_file):
sds = sd.select('BT_Diff_By_AIRS_BT')
plot_sds(sds, band, out_file,
title='BT Difference By AIRS BT',
x_label='AIRS BT (K)',
y_label='Mean AIRS-MODIS BT Diff (K)')
# plot mean BT difference as a function of scan angle
def plot_bt_diff_by_scan_angle(sd, band, out_file):
sds = sd.select('BT_Diff_By_Scan_Angle')
plot_sds(sds, band, out_file,
title='BT Difference By Scan Angle',
x_label='Scan Angle (Degrees)',
y_label='Mean AIRS-MODIS BT Diff (K)')
# plot mean BT difference as a function of solar zenith angle
def plot_bt_diff_by_solar_zenith(sd, band, out_file):
sds = sd.select('BT_Diff_By_Solar_Zenith')
plot_sds(sds, band, out_file,
title='BT Difference By Solar Zenith Angle',
x_label='Solar Zenith Angle (Degrees)',
y_label='Mean AIRS-MODIS BT Diff (K)')
# create a map showing AIRS BTs
def plot_airs_bt_by_location(sd, band, out_file):
bt = sd.select('AIRS_BT_By_Location')[band]
map_helper(out_file, bt, 'Mean AIRS BT (K)')
# create a map showing BT differences
def plot_bt_diff_by_location(sd, band, out_file):
diff = sd.select('BT_Diff_By_Location')[band]
map_helper(out_file, diff, 'Mean AIRS-MODIS BT Diff (K)')
# check command-line arguments
if len(sys.argv) != 3:
print 'usage: python %s <aggregate_file> <output_dir>' % sys.argv[0]
sys.exit(1)
aggr_file = sys.argv[1]
output_dir = sys.argv[2]
# the default size for plot titles and labels is too big
matplotlib.rc('font', size=8)
# open the aggregate file and produce plots for all bands
sd = SD(aggr_file)
for band in range(16):
plot_bt_hist(
sd,
band,
'%s/bt_hist_%02d.png' % (output_dir, band)
)
plot_bt_diff_hist(
sd,
band,
'%s/bt_diff_hist_%02d.png' % (output_dir, band)
)
plot_bt_diff_by_airs_bt(
sd,
band,
'%s/bt_diff_by_airs_bt_%02d.png' % (output_dir, band)
)
plot_bt_diff_by_scan_angle(
sd,
band,
'%s/bt_diff_by_scan_angle_%02d.png' % (output_dir, band)
)
plot_bt_diff_by_solar_zenith(
sd,
band,
'%s/bt_diff_by_solar_zenith_%02d.png' % (output_dir, band)
)
plot_airs_bt_by_location(
sd,
band,
'%s/airs_bt_by_location_%02d.png' % (output_dir, band)
)
plot_bt_diff_by_location(
sd,
band,
'%s/bt_diff_by_location_%02d.png' % (output_dir, band)
)