Skip to content
Snippets Groups Projects
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)
    )