vima_demo.vima_dd_simple

An simple example to read & plot DD data.

""" An simple example to read & plot DD data.
"""
import sys
import datetime
import dateutil.parser

import matplotlib.pyplot as plt
import matplotlib.colors as pltcolors

from irfpy.vima import imaextra as vextra
from irfpy.vima import energy


import logging
logger = logging.getLogger('irfpy.vima_dd_simple')
vmax = 10000
vmin = 0.1

def set_logger_verbose():
    handler = logging.StreamHandler()
    handler.setLevel(logging.DEBUG)
    formatter = logging.Formatter('[irfpy] %(levelname)s:%(asctime)s:%(name)s:%(message)s', datefmt='%FT%T')
    handler.setFormatter(formatter)
    logging.getLogger('irfpy').addHandler(handler)
    logging.getLogger('irfpy').setLevel(logging.DEBUG)


def plot_extra(t0, t1, species):
    # Load imaextra file
    imaex = vextra.getImaExtra(t0, t1)
    logger.info(imaex.ndat)   # ndat is the number of data

    tlist = imaex.getobstime()     # obstime contains the observation time
    logger.info((tlist[0], tlist[-1]))

    if species == 'H':
        hp = imaex.getHpSpec()      # Proton counts. (A16, E96, P16, Time)
    elif species == 'O':
        hp = imaex.getHeavySpec()
    else:
        raise ValueError('Only H or O supported for species (given: {})'.format(species))
    logger.info(hp.shape)       # The shape is (A16, E96, P16, Time)

    # hp is 4D data. This will be converted to 2-D data, (E, Time), with max in azimuth, collapse polar to time.
    hpet = hp.max(axis=0)    # Maximum in azimuth direction.  Now (E96, P16, Time)
    import numpy as np
    hpet = np.reshape(hpet, (96, 16 * imaex.ndat), order='F')    # Fortran order is specified, then the axis=0 is kept

    # Time list should also be updated
    tlist2 = []
    for _t in tlist:
        for _i in range(16):
            tlist2.append(_t + _i * datetime.timedelta(seconds=12))
    tlist = tlist2

    # Energy table is obtained.
    etbl = imaex.data['ETableN']   # It is the index of energy table. 0 is for version 1, and 1 is for version 3
    logger.info(etbl[0])

    if etbl[0] == 0:
        ene = energy.get_default_table_v1(negative=False)
    else:
        ene = energy.get_default_table_v3(negative=False)

    # Plotting
    plt.pcolormesh(tlist, ene, hpet, norm=pltcolors.LogNorm(vmin=vmin, vmax=vmax, clip=True))
    plt.yscale('log')
    plt.ylim(10, 15000)


def main(t0, t1):
    plt.subplot(211)
    plot_extra(t0, t1, 'H')
    plt.subplot(212)
    plot_extra(t0, t1, 'O')

    plt.savefig('vima_dd_simple.png')

def mainexec():
    t0 = dateutil.parser.parse(sys.argv[1])
    t1 = dateutil.parser.parse(sys.argv[2])

    set_logger_verbose()

    main(t0, t1)

if __name__ == '__main__':
    mainexec()