[IceCube FOV Plot] Script to plot location of IceCube 160427A #icecube #plot #160427A

README.md

Here is the plot produced by this script:

icecube-160427A-plot.py

#!/usr/local/bin/python
# encoding: utf-8
"""
icecube-160427A-plot.py
=======================
:Summary:
    Script to plot location of IceCube 160427A

:Author:
    David Young

:Date Created:
    August 23, 2016

Usage:
    icecube-160427A-plot

    -h, --help            show this help message
    -v, --version         show version
    -s, --settings        the settings file

..todo::

    @review: when complete pull all general functions and classes into dryxPython
"""
## USER DEFINED VARIABLES - ONLY PLACE THIS SCRIPT NEEDS MODIFIED ##


################# GLOBAL IMPORTS ####################
import sys
import os
os.environ['TERM'] = 'vt100'
import readline
import glob
import docopt
# SUPPRESS MATPLOTLIB WARNINGS
import warnings
warnings.filterwarnings("ignore")
import matplotlib as mpl
from fundamentals import tools, times
from astropy import wcs as awcs
import numpy as np
from matplotlib.pyplot import savefig
from crowdedText import adjust_text


def main(arguments=None):
    """
    The main function used when ``icecube-160427A-plot.py`` is run as a single script from the cl, or when installed as a cl command
    """

    # setup the command-line util settings
    su = tools(
        arguments=arguments,
        docString=__doc__,
        logLevel="DEBUG",
        options_first=False,
        projectName="lv",
        tunnel=False
    )
    arguments, settings, log, dbConn = su.setup()

    plot(
        log=log
    )

    return


def plot(
        log,
        fileFormats=["pdf"],
        folderName="",
        projection="tan",
        outputDirectory=False):
    """
    """
    log.info('starting the ``plot`` method')

    import healpy as hp
    from matplotlib.font_manager import FontProperties
    font = FontProperties()
    import matplotlib.pyplot as plt

    # HEALPY REQUIRES RA, DEC IN RADIANS AND AS TWO SEPERATE ARRAYS
    import math
    pi = (4 * math.atan(1.0))
    DEG_TO_RAD_FACTOR = pi / 180.0
    RAD_TO_DEG_FACTOR = 180.0 / pi

    # VARIABLES
    font.set_family("Arial")
    pixelSizeDeg = 0.066667
    unit = "likelihood"
    cmap = "YlOrRd"
    colorBar = False

    # PARAMETERS
    raCentre = 240.57
    decCentre = 9.34
    centralCoordinate = [raCentre, decCentre]
    plotParameters = {}
    plotParameters["raRange"] = 1.5
    plotParameters["decRange"] = 1.5
    plotTitle = "IceCube-160427A-QSOs-PS16cgx-Location"

    # INITIALISE FIGURE
    fig = plt.figure()

    # pixelSizeDeg = 1.0 / 60.

    # CREATE A NEW WCS OBJECT
    w = awcs.WCS(naxis=2)
    # SET THE REQUIRED PIXEL SIZE
    w.wcs.cdelt = np.array([pixelSizeDeg, pixelSizeDeg])
    # WORLD COORDINATES AT REFERENCE PIXEL
    w.wcs.crval = centralCoordinate

    projectionDict = {
        "mollweide": "MOL",
        "aitoff": "AIT",
        "hammer": "AIT",
        "lambert": "ZEA",
        "polar": "TAN",
        "rectilinear": "MER"
    }

    if projection in ["mollweide", "aitoff", "hammer", "lambert", "polar", "rectilinear"]:
        # MAP VISULISATION RATIO IS ALWAYS 1/2
        xRange = 2000
        yRange = xRange / 2.

        # SET THE REFERENCE PIXEL TO THE CENTRE PIXEL
        w.wcs.crpix = [xRange / 2., yRange / 2.]

        # FULL-SKY MAP SO PLOT FULL RA AND DEC RANGES
        # DEC FROM 180 to 0
        theta = np.linspace(np.pi, 0, yRange)
        latitude = np.radians(np.linspace(-90, 90, yRange))
        # RA FROM -180 to +180
        phi = np.linspace(-np.pi, np.pi, xRange)
        longitude = np.radians(np.linspace(-180, 180, xRange))
        X, Y = np.meshgrid(longitude, latitude)

        # PROJECT THE MAP TO A RECTANGULAR MATRIX xRange X yRange
        PHI, THETA = np.meshgrid(phi, theta)
        healpixIds = hp.ang2pix(nside, THETA, PHI)
        probs = aMap[healpixIds]

        # healpixIds = np.reshape(healpixIds, (1, -1))[0]

        # CTYPE FOR THE FITS HEADER
        thisctype = projectionDict[projection]
        w.wcs.ctype = ["RA---%(thisctype)s" %
                       locals(), "DEC--%(thisctype)s" % locals()]

        stampProb = np.sum(aMap)
        print "Probability for the plot stamp is %(stampProb)s" % locals()

        # MATPLOTLIB IS DOING THE PROJECTION
        # ax = fig.add_subplot(111, projection=projection)
        # USE WCS AS THE PROJECTION

        if projection == "mollweide":
            ax = plt.axes(projection=ccrs.Mollweide())
            ax.set_extent([80, 170, -45, 30])
        else:
            ax = fig.add_axes([0.15, 0.1, 0.8, 0.8], projection=projection)

        print ax.get_extent()

        # RASTERIZED MAKES THE MAP BITMAP WHILE THE LABELS REMAIN VECTORIAL
        # FLIP LONGITUDE TO THE ASTRO CONVENTION
        ax.contourf(longitude[
            ::-1], latitude, probs,
            transform=ccrs.PlateCarree(),
            cmap=cmap)

        # ax.coastlines()
        ax.set_global()
        plt.show()

        sys.exit(0)

        image = ax.pcolormesh(longitude[
            ::-1], latitude, probs, rasterized=True, cmap=cmap)

        plt.show()

        # GRATICULE
        if projection not in ["polar", "rectilinear", "mollweide"]:
            ax.set_longitude_grid(60)
            ax.set_latitude_grid(45)
            ax.xaxis.set_major_formatter(ThetaFormatterShiftPi(60))
            ax.set_longitude_grid_ends(90)

        # CONTOURS - NEED TO ADD THE CUMMULATIVE PROBABILITY
        i = np.flipud(np.argsort(aMap))
        cumsum = np.cumsum(aMap[i])
        cls = np.empty_like(aMap)
        cls[i] = cumsum * 99.99999999 * stampProb

        # EXTRACT CONTOUR VALUES AT HEALPIX INDICES
        contours = []
        contours[:] = [cls[i] for i in healpixIds]
        # contours = np.reshape(np.array(contours), (yRange, xRange))

        CS = ax.contour(longitude[::-1], latitude,
                        contours, linewidths=.5, alpha=0.7, zorder=2, origin='image')

        if projection not in ["lambert", "aitoff"]:
            CS.set_alpha(0.5)
            CS.clabel(fontsize=10, inline=True,
                      fmt='%2.1f', fontproperties=font, alpha=0.0)

        # COLORBAR
        if colorBar:
            cb = fig.colorbar(image, orientation='horizontal',
                              shrink=.6, pad=0.05, ticks=[0, 1])
            cb.ax.xaxis.set_label_text("likelihood")
            cb.ax.xaxis.labelpad = -8
            # WORKAROUND FOR ISSUE WITH VIEWERS, SEE COLORBAR DOCSTRING
            cb.solids.set_edgecolor("face")

        ax.tick_params(axis='x', labelsize=16)
        ax.tick_params(axis='y', labelsize=16)
        # lon.set_ticks_position('bt')
        # lon.set_ticklabel_position('b')
        # lon.set_ticklabel(size=20)
        # lat.set_ticklabel(size=20)
        # lon.set_axislabel_position('b')
        # lat.set_ticks_position('lr')
        # lat.set_ticklabel_position('l')
        # lat.set_axislabel_position('l')

        # # REMOVE TICK LABELS
        # ax.xaxis.set_ticklabels([])
        # ax.yaxis.set_ticklabels([])
        # # REMOVE GRID
        # ax.xaxis.set_ticks([])
        # ax.yaxis.set_ticks([])

        # REMOVE WHITE SPACE AROUND FIGURE
        spacing = 0.01
        plt.subplots_adjust(bottom=spacing, top=1 - spacing,
                            left=spacing, right=1 - spacing)

        plt.grid(True)

    elif projection in ["tan"]:

        # UNPACK THE PLOT PARAMETERS
        raRange = plotParameters["raRange"]
        decRange = plotParameters["decRange"]

        raMax = centralCoordinate[0] + raRange / 2.
        raMin = centralCoordinate[0] - raRange / 2.
        decMax = centralCoordinate[1] + decRange / 2.
        decMin = centralCoordinate[1] - decRange / 2.

        # DETERMINE THE PIXEL GRID X,Y RANGES
        xRange = int(raRange / pixelSizeDeg)
        yRange = int(decRange / pixelSizeDeg)

        # SET THE REFERENCE PIXEL TO THE CENTRE PIXEL
        w.wcs.crpix = [xRange / 2., yRange / 2.]

        # USE THE "GNOMONIC" PROJECTION ("COORDINATESYS---PROJECTION")
        w.wcs.ctype = ["RA---TAN", "DEC--TAN"]

        # CREATE A PIXEL GRID - 2 ARRAYS OF X, Y
        columns = []
        px = np.tile(np.arange(0, xRange), yRange)
        py = np.repeat(np.arange(0, yRange), xRange)

        # CONVERT THE PIXELS TO WORLD COORDINATES
        wr, wd = w.wcs_pix2world(px, py, 1)

        # MAKE SURE RA IS +VE
        nr = []
        nr[:] = [r if r > 0 else r + 360. for r in wr]
        wr = np.array(nr)

        # CREATE THE FITS HEADER WITH WCS
        header = w.to_header()

        # CREATE THE FITS FILE

        # GRAB THE WCS FROM HEADER GENERATED EARLIER
        from wcsaxes import datasets, WCS
        from astropy.wcs import WCS
        from wcsaxes import WCSAxes

        wcs = WCS(header)
        # USE WCS AS THE PROJECTION
        ax = WCSAxes(fig, [0.15, 0.1, 0.8, 0.8], wcs=wcs)
        # note that the axes have to be explicitly added to the figure
        ax = fig.add_axes(ax)

        # RESET THE AXES TO THE FRAME OF THE FITS FILE
        ax.set_xlim(-0.5, xRange - 0.5)
        ax.set_ylim(-0.5, yRange - 0.5)

        # THE COORDINATES USED IN THE PLOT CAN BE ACCESSED USING THE COORDS
        # ATTRIBUTE (NOT X AND Y)
        lon = ax.coords[0]
        lat = ax.coords[1]

        lon.set_axislabel('RA (deg)', minpad=0.87,
                          size=20)
        lat.set_axislabel('DEC (deg)', minpad=0.87,
                          size=20)

        lon.set_major_formatter('d.d')
        lat.set_major_formatter('d.d')

        # THE SEPARATORS FOR ANGULAR COORDINATE TICK LABELS CAN ALSO BE SET BY
        # SPECIFYING A STRING
        lat.set_separator(':-s')
        # SET THE APPROXIMATE NUMBER OF TICKS, WITH COLOR & PREVENT OVERLAPPING
        # TICK LABELS FROM BEING DISPLAYED.
        lon.set_ticks(number=4, color='#657b83',
                      exclude_overlapping=True, size=10)
        lat.set_ticks(number=10, color='#657b83',
                      exclude_overlapping=True, size=10)

        # MINOR TICKS NOT SHOWN BY DEFAULT
        lon.display_minor_ticks(True)
        lat.display_minor_ticks(True)
        lat.set_minor_frequency(2)

        # CUSTOMISE TICK POSITIONS (l, b, r, t == left, bottom, right, or
        # top)
        lon.set_ticks_position('bt')
        lon.set_ticklabel_position('b')
        lon.set_ticklabel(size=20)
        lat.set_ticklabel(size=20)
        lon.set_axislabel_position('b')
        lat.set_ticks_position('lr')
        lat.set_ticklabel_position('l')
        lat.set_axislabel_position('l')

        # HIDE AXES
        # lon.set_ticklabel_position('')
        # lat.set_ticklabel_position('')
        # lon.set_axislabel('', minpad=0.5, fontsize=12)
        # lat.set_axislabel('', minpad=0.5, fontsize=12)

        # ADD A GRID
        ax.coords.grid(color='#657b83', alpha=0.5, linestyle='dashed')
        plt.gca().invert_xaxis()

    else:
        log.error(
            'please give a valid projection. The projection given was `%(projection)s`.' % locals())

    header = w.to_header()
    # CREATE THE FITS FILE

    subTitle = ""

    raDeg = raCentre
    decDeg = decCentre

    # PLOT ICECUBE FOOTPRINT
    from matplotlib.patches import Circle
    from matplotlib.patches import Ellipse
    height = 36 * 2 / 60.
    width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)

    # MULTIPLE CIRCLES
    if projection in ["tan"]:
        circ = Ellipse(
            (raDeg, decDeg), width=width, height=height, alpha=0.2, color='#268bd2', fill=True, transform=ax.get_transform('fk5'), zorder=3)
    else:
        if raDeg > 180.:
            raDeg = raDeg - 360.
        circ = Ellipse(
            (-raDeg * DEG_TO_RAD_FACTOR, decDeg * DEG_TO_RAD_FACTOR), width=width * DEG_TO_RAD_FACTOR, height=height * DEG_TO_RAD_FACTOR, alpha=0.2, color='#268bd2', fill=True, zorder=3)

    ax.add_patch(circ)

    # # LEGEND FOR PS1
    # if len(ps1Pointings):
    #     raDeg = 90.
    #     decDeg = 10.
    #     height = 3.5

    #     width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)
    #     if projection in ["tan"]:
    #         circ = Ellipse(
    #             (raDeg, decDeg), width=width, height=height, alpha=0.2, color='#859900', fill=True, transform=ax.get_transform('fk5'), zorder=3)
    #         ax.text(
    #             raDeg - 3,
    #             decDeg - 1.,
    #             "PS1",
    #             fontsize=14,
    #             zorder=4,
    #             family='monospace',
    #             transform=ax.get_transform('fk5')
    #         )
    #     else:
    #         height = 6.
    #         width = height / math.cos(decDeg * DEG_TO_RAD_FACTOR)
    #         raDeg = 285.
    #         decDeg = 38.
    #         if raDeg > 180.:
    #             raDeg = raDeg - 360.
    #         circ = Ellipse(
    #             (-raDeg * DEG_TO_RAD_FACTOR, decDeg * DEG_TO_RAD_FACTOR), width=width * DEG_TO_RAD_FACTOR, height=height * DEG_TO_RAD_FACTOR, alpha=0.2, color='#859900', fill=True, zorder=3)
    #         ax.text(
    #             (-raDeg + 8.) * DEG_TO_RAD_FACTOR,
    #             (decDeg - 3.) * DEG_TO_RAD_FACTOR,
    #             "PS1",
    #             fontsize=14,
    #             zorder=4,
    #             family='monospace'
    #         )
    #     ax.add_patch(circ)

    # ADD DATA POINTS FOR TRANSIENTS
    names = []
    ra = []
    dec = []
    raRad = []
    decRad = []
    texts = []
    ps1QSOs = []

    ps1SN = [
        {"ps1_designation": "PS16cgv",
            "ra_psf": "15 59 53.84", "dec_psf": "+09 11 08.4"},
        {"ps1_designation": "PS16cgw",
            "ra_psf": "16 00 39.69", "dec_psf": "+09 39 21.2"},
        {"ps1_designation": "PS16cfu",
            "ra_psf": "16 01 15.66", "dec_psf": "+09 25 04.7"},
        {"ps1_designation": "PS16cfz",
            "ra_psf": "16 02 11.96", "dec_psf": "+09 54 07.9"},
        {"ps1_designation": "PS16cgb",
            "ra_psf": "16 02 19.12", "dec_psf": "+09 34 50.1"},
        {"ps1_designation": "PS16cgi", "ra_psf": "16 03 27.72", "dec_psf": "+08 54 05.2"}
    ]

    ps1BigSN = [{"ps1_designation": "PS16cgx",
                 "ra_psf": "16 01 18.60", "dec_psf": "+09 51 53.1"}]

    # ps1QSOs = [
    #     {"ps1_designation": ""}
    # ]
    from astrocalc.coords import unit_conversion
    converter = unit_conversion(
        log=log
    )
    for trans in ps1SN:
        # if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
        #     continue

        name = trans["ps1_designation"]

        names.append(name)
        decDeg = converter.dec_sexegesimal_to_decimal(
            dec=trans["dec_psf"]
        )
        raDeg = converter.ra_sexegesimal_to_decimal(
            ra=trans["ra_psf"]
        )

        ra.append(raDeg)
        dec.append(decDeg)
        raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
        decRad.append(decDeg * DEG_TO_RAD_FACTOR)

    if len(ra) > 0:
        # MULTIPLE CIRCLES
        if projection in ["tan"]:
            ax.scatter(
                x=np.array(ra),
                y=np.array(dec),
                transform=ax.get_transform('fk5'),
                s=15,
                c='black',
                edgecolor='black',
                alpha=1,
                zorder=4
            )
            xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
            print xx, yy
            # ADD TRANSIENT LABELS
            for r, d, n in zip(xx, yy, names):
                # if n == "PS16cgv":
                #     r = r + 0.5
                texts.append(ax.text(
                    r,
                    d,
                    n,
                    fontsize=16,
                    zorder=4,
                    family='monospace'
                ))
        else:
            ax.scatter(
                x=np.array(raRad),
                y=np.array(decRad),
                s=15,
                c='#dc322f',
                edgecolor='#dc322f',
                alpha=1,
                zorder=4
            )

    ps1QSOs = [
        {"ps1_designation": "PS16cfw",
            "ra_psf": "16:01:39.10", "dec_psf": "+09:21:59.5"},
        {"ps1_designation": "PS16cfx",
            "ra_psf": "16:01:49.65", "dec_psf": "+08:50:29.7"},
        {"ps1_designation": "PS16cfy",
            "ra_psf": "16:02:05.20", "dec_psf": "+09:46:30.4"},
        {"ps1_designation": "PS16cgg",
            "ra_psf": "16:03:09.19", "dec_psf": "+09:20:19.5"},
        {"ps1_designation": "PS16cgh",
            "ra_psf": "16:03:17.92", "dec_psf": "+09:00:38.0"},
        {"ps1_designation": "PS16cgm",
            "ra_psf": "15:59:47.90", "dec_psf": "+09:10:22.4"},
        {"ps1_designation": "PS16cgn",
            "ra_psf": "15:59:54.76", "dec_psf": "+09:14:26.3"},
        {"ps1_designation": "PS16cgo",
            "ra_psf": "16:01:06.27", "dec_psf": "+08:46:05.8"}
    ]

    # ADD DATA POINTS FOR TRANSIENTS
    names = []
    ra = []
    dec = []
    raRad = []
    decRad = []

    for trans in ps1QSOs:
        # if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
        #     continue

        name = trans["ps1_designation"]

        names.append(name)
        decDeg = converter.dec_sexegesimal_to_decimal(
            dec=trans["dec_psf"]
        )
        raDeg = converter.ra_sexegesimal_to_decimal(
            ra=trans["ra_psf"]
        )

        ra.append(raDeg)
        dec.append(decDeg)
        raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
        decRad.append(decDeg * DEG_TO_RAD_FACTOR)

    if len(ra) > 0:
        # MULTIPLE CIRCLES
        if projection in ["tan"]:
            ax.scatter(
                x=np.array(ra),
                y=np.array(dec),
                transform=ax.get_transform('fk5'),
                s=15,
                c='#268bd2',
                edgecolor='#268bd2',
                alpha=1,
                zorder=4,
                marker="D"
            )

            xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
            print xx, yy
            # ADD TRANSIENT LABELS
            for r, d, n in zip(xx, yy, names):
                # if n in ["PS16cgn", "PS16cgm"]:
                #     r = r + 0.8
                texts.append(ax.text(
                    r,
                    d,
                    n,
                    color='#268bd2',
                    fontsize=16,
                    zorder=4,
                    family='monospace'
                ))

        else:
            ax.scatter(
                x=np.array(raRad),
                y=np.array(decRad),
                s=15,
                c='#268bd2',
                edgecolor='#268bd2',
                alpha=1,
                zorder=4
            )

     # ADD DATA POINTS FOR TRANSIENTS
    names = []
    ra = []
    dec = []
    raRad = []
    decRad = []

    for trans in ps1BigSN:
        # if trans["ps1_designation"] in ["PS15dpg", "PS15dpp", "PS15dpq", "PS15don", "PS15dpa", "PS15dom"]:
        #     continue

        name = trans["ps1_designation"]

        names.append(name)
        decDeg = converter.dec_sexegesimal_to_decimal(
            dec=trans["dec_psf"]
        )
        raDeg = converter.ra_sexegesimal_to_decimal(
            ra=trans["ra_psf"]
        )

        ra.append(raDeg)
        dec.append(decDeg)
        raRad.append(-raDeg * DEG_TO_RAD_FACTOR)
        decRad.append(decDeg * DEG_TO_RAD_FACTOR)

    if len(ra) > 0:
        # MULTIPLE CIRCLES
        if projection in ["tan"]:

            ax.scatter(
                x=np.array(ra),
                y=np.array(dec),
                transform=ax.get_transform('fk5'),
                s=50,
                facecolor='none',
                edgecolor='#dc322f',
                alpha=1,
                zorder=4
            )
            ax.scatter(
                x=np.array(ra),
                y=np.array(dec),
                transform=ax.get_transform('fk5'),
                s=15,
                c='#dc322f',
                edgecolor='#dc322f',
                alpha=1,
                zorder=4
            )
            xx, yy = w.wcs_world2pix(np.array(ra), np.array(dec), 0)
            # ADD TRANSIENT LABELS
            for r, d, n in zip(xx, yy, names):

                # if n == "PS16cgv":
                #     r = r + 0.5
                texts.append(ax.text(
                    r,
                    d,
                    n,
                    color='#dc322f',
                    fontsize=16,
                    zorder=4,
                    family='monospace'
                ))
        else:
            ax.scatter(
                x=np.array(raRad),
                y=np.array(decRad),
                s=15,
                c='#dc322f',
                edgecolor='#dc322f',
                alpha=1,
                zorder=4
            )

    from random import shuffle
    shuffle(texts)
    adjust_text(
        xx,
        yy,
        texts,
        expand_text=(1.0, 1.0),
        expand_points=(1.0, 1.0),
        va='center',
        ha='right',
        autoalign=False,
        force_text=1.,
        force_points=1.,
        lim=100,
        precision=0,
        only_move={},
        text_from_text=True,
        text_from_points=True,
        save_steps=False,
        save_prefix='',
        save_format='png',
        add_step_numbers=True,
        min_arrow_sep=3.,
        draggable=True,
        arrowprops=dict(arrowstyle="->", facecolor='black', lw=1.2,
                        patchB=None, shrinkB=0, connectionstyle="arc3,rad=0.1", zorder=3, alpha=0.5),
        fontsize=16,
        family='monospace',
        repel_from_axes=True
    )

    # TIME-RANGE LABEL
    fig = plt.gcf()

    fig.set_size_inches(8.0, 8.0)
    fWidth, fHeight = fig.get_size_inches()
    print fWidth, fHeight
    if projection == "tan":
        plt.text(
            xRange * 0.25,
            # xRange * 0.95,
            yRange * 0.93,
            "",
            fontsize=16,
            zorder=4,
            color="#dc322f",
            fontproperties=font
        )
        plt.text(
            xRange * 0.1,
            # xRange * 0.95,
            yRange * 0.93,
            "",
            fontsize=20,
            zorder=4,
            color="black",
            fontproperties=font
        )
    else:
        plt.text(
            fWidth * 0.7,
            # xRange * 0.95,
            fHeight * 0.95,
            timeRangeLabel,
            fontsize=16,
            zorder=4,
            color="#dc322f",
            fontproperties=font
        )
        plt.text(
            fWidth * 0.87,
            # xRange * 0.95,
            fHeight * 0.315,
            "",
            fontsize=20,
            zorder=4,
            color="black",
            fontproperties=font
        )

    # Recursively create missing directories
    plotDir = "./"
    plotTitle = plotTitle.replace(" ", "_").replace(
        "<", "lt").replace(">", "gt").replace(",", "").replace("\n", "_").replace("&", "and")
    figureName = """%(plotTitle)s""" % locals(
    )

    for f in fileFormats:
        figurePath = "%(plotDir)s/%(figureName)s_%(projection)s.%(f)s" % locals()
        savefig(figurePath, bbox_inches='tight', dpi=300)

    log.info('completed the ``plot`` method')
    return None

if __name__ == '__main__':
    main()