+27833657551     admin@starscientist.org
    +27833657551     admin@starscientist.org

We will reproduce the below plot using aplpy. We will plot channel maps of the Holmberg I datacube (channel 11 to channel 22, counting from zero). The data can be obtained from https://www2.mpia-hd.mpg.de/THINGS/Data_files/HO_I_NA_CUBE_THINGS.FITS. To get it, simply type the following in your terminal:

wget https://www2.mpia-hd.mpg.de/THINGS/Data_files/HO_I_NA_CUBE_THINGS.FITS
plot_aplpy

Below is the python script to reproduce the above plot.

from astropy.io import fits
import matplotlib.pyplot as plt
import aplpy
import numpy
import os
import matplotlib as mpl
from matplotlib.backends.backend_pdf import PdfPages

font = 'Arial'  # 
mpl.rcParams['font.sans-serif'] = font
mpl.rc('mathtext', fontset='custom', it=font + ':italic')
mpl.rc('font', size=13)

fig = plt.figure(figsize=(7, 15))
############

cmap = plt.cm.cubehelix
majlen = 8.7
minlen = 5
pad = 16
dx = 0.35
dy = 0.15
xspacing = 20 # arcmin, avoid clustering of xticklabels
ncol = 3 # number of column
nrow = 4 # number of rows
np = ncol*nrow
jx = list(numpy.arange(0, ncol))*np
jy = numpy.repeat(numpy.arange(0, nrow),ncol)
x_len = 0.1
y_len = 0.1
v_left = list(numpy.arange(0, np-1, ncol))
v_right = list(numpy.arange(ncol-1, np-1, ncol))
h_bot = list(numpy.arange(np-ncol, np))

f = "HO_I_NA_CUBE_THINGS.FITS" # datacubes whose channel maps are to be plotted
hdu = fits.open(f)[0]
img = numpy.squeeze(hdu.data) # numpy.squeeze removes unnecessary dimensions
hdr = hdu.header
# put these otherwise aplpy.add_beam won't work
hdr["BMAJ"] = 4.0710E-03 # minor axis of beam
hdr["BMIN"] = 3.5366E-03 # major axis of beam
hdr["BPA"] = -41.60 #  Beam position angle
##
hdr["NAXIS"] = 2
delheader4 = ["NAXIS4",
             "CTYPE4",
             "CRVAL4",
             "CDELT4",
             "CRPIX4",
             "CROTA4"]
delheader3 = ["NAXIS3",
              "CTYPE3",
              "CRVAL3",
              "CDELT3",
              "CRPIX3",
              "CROTA3",
              "CUNIT3"]
allhead = delheader4 + delheader3
#Remove unnecessary keys in headers
for p in allhead:
    if p in hdr.keys():
        #print(p)
        del hdr[p]

t = 11
vmin = -0.007 * 1000 # converting to mJy
vmax = 0.0241 * 1000
xspacing = 20
axlist = list()
for j in range(0, np):
    print(j)
    xj = x_len + dx * jx[j]     
    yj = y_len - dy * jy[j]
    fits.writeto(f[:-5]+"_CH"+str(t)+".FITS", img[t,:,:] * 1000, hdr,  overwrite=True)
    fj = aplpy.FITSFigure(f[:-5]+"_CH"+str(t)+".FITS", figure=fig, subplot=[xj, yj, dx, dy])
    fj.add_beam()
    fj.beam.set_color('white')
    fj.show_colorscale(aspect='auto', cmap=cmap, vmin=vmin, vmax=vmax)
    fj.ticks.set_xspacing(xspacing/60)
    t = t + 1 
    fj.beam.show(facecolor="None", edgecolor='black')
    if j not in (v_left + v_right + h_bot): 
      fj.tick_labels.hide_y()
      fj.tick_labels.hide_x()
      fj.axis_labels.hide_y()
    if j not in list(numpy.arange(0, np, ncol)):
      fj.axis_labels.hide_y()
      fj.tick_labels.hide_y()
    if j in range(0, np-ncol):
      fj.tick_labels.hide_x()
      fj.axis_labels.hide_x()
    axj = plt.gca()
    axlist.append(axj)
    axj.tick_params(pad=pad)
    axj.tick_params(direction='in', length=majlen)
    axj.tick_params(which='minor', length=minlen)
sm = plt.cm.ScalarMappable(cmap=cmap, norm=plt.Normalize(vmin=vmin, vmax=vmax))
sm._A = []
cax = plt.axes([x_len+dx*ncol, yj, 0.03, dy*nrow])
cb = fig.colorbar(sm, ax=axlist, cax=cax, shrink=0.75)
cb.ax.set_ylabel("Flux "+r"$(\mathrm{mJy~beam^{-1}})$", labelpad=13)
cb.ax.tick_params(direction='in', left='on', pad=15)
pdf = PdfPages('plot_aplpy.pdf')
plt.savefig(pdf, format='pdf', bbox_inches='tight', dpi=400)
pdf.close()
pdf = None
os.system("xdg-open plot_aplpy.pdf") # Change to open on a mac, instead of xdg-open