# Basic Plotting in Python

This is an example of how to make a simple plot in python, using data stored in a .csv file.
We will:

•   Load the 2 columns of data from the file into a (numpy) array
•   Plot the data with pyplot.plot
•   Tweak some plot settings to make it pretty
•   Save the plot to a file, view the plot in a window, or both

You will need to have installed on your machine:

• python (I’m running 2.7 in this tutorial)
• numpy (I’m using 1.5.1 here)
• matplotlib (I’m using 1.1.0 here)

You can run this entire script (see the condensed version at the end of this post) from the command line by setting it to executable and typing:
./plottingexample.py
or, you can run the whole script in ipython (type ipython at the command line) by typing
execfile(‘plottingexample.py’)
or, you can run each command individually in ipython

This script is not optimized for resources or speed – it is simply meant to be an easy-to-follow introduction to simple plotting.  If you haven’t seen it already, the official Matplotlib example gallery is an invaluable resource for plotting in python.  I would also refer the reader to the overview of plotting with Matplotlib by Nicolas P. Rougier here.

First, we will need to load some python modules…
– numpy is the numerical module that allows you to do Matlab- and IDL-esque numerical operations
– pyplot (part of matplotlib) is the plotting library

import numpy
from matplotlib import pyplot as plt


Here, we will load the data into python.  The example csv file is separated by commas, and strings (the titles) are in double-quotes.  This particular file has headers in the first row.  There are several ways to load a csv file into python; this is simply one easy way.

Normally, we might use numpy.loadtxt to do this, but our file has a funny format: the + sign in the scientific notation.  numpy.genfromtxt is a more robust version of loadtxt which can interpret this correctly and deal with missing cells.  This creates a 2D array of the data columns.

data=numpy.genfromtxt('SampleData.csv',skiprows=1,delimiter=',')



If we really wanted to be pedantic with our separation, we could create separate arrays for each column:

luminosity=data[:,1]
mass=data[:,0]


We can explicitly make a log(luminosity) array with numpy: (note that numpy.log is the natural log)

logL=numpy.log10(luminosity)



Plot the data values.  Syntax is pyplot.plot(xarray,yarray, other kwargs).  There is also a plt.scatter command, but we can just set the linewidth to 0.  Note that we can do computations on the fly within pyplot!

plt.plot(numpy.log10(luminosity),mass,'*',linewidth=0)



And now to set the labels. You can use Latex inline equation syntax.

plt.title('M vs. L')
plt.xlabel('log(Luminosity) in some units...')
plt.ylabel(r'log(M$_{H_{2}}$) in some units...')



Plot simply creates a plot object.  To view it, we need to either show or savefig.  Uncomment the one you want.  Note that pyplot will recognize the filetype (png,eps,pdf, etc…) you give it in savefig and save accordingly!

plt.savefig('plottest1.png',dpi=100)



Then clear for the next example:

plt.clf()


…And this is what we get:

### ==== SECOND EXAMPLE ====

Now let’s make two subplots:
a.) Large red half-transparent squares with blue edges of all the data points
b.) Thin green diamonds for real points, downward pointing arrows for the upper limits
We first need to create a figure, then add the subplots

For this one, we’ll use the second csv file which has a column describing whether a mass value is an upper limit
–> genfromtxt by default expects floats, but will read strings if we set the expected datatype (dtype) to None

mass=numpy.genfromtxt('SampleData2.csv',skiprows=1,delimiter=',',usecols=0)
logL=numpy.log10(numpy.genfromtxt('SampleData2.csv',skiprows=1,delimiter=',',usecols=2))



An aside: The slickest way to use masks is with a numpy masked array.  See my tutorial on NaNs and masked arrays here.  For example:

realmask=numpy.ma.masked_where(massmask=='"Yes"',massmask) # Remember our mask column is '"Yes"' and '"No"'



If you want only the unmasked values, call with mass_limits.compressed(), etc… [end aside] Here we create separate lists for mass & luminosity of the upper limit values only and regular values only, using a standard for loop.  This method requires extra steps below.

limits=[]; regular=[]
for i in range(0,len(mass)):
else: regular.append([mass[i],logL[i]])



Convert the lists to numpy arrays and transpose to put the data in order for plotting

limits=numpy.array(limits).transpose()
regular=numpy.array(regular).transpose()



Plotting.
Create the figure to which you will add subplots:

fig1=plt.figure(1)



sub1=fig1.add_subplot(2,1,1)
# 2,1,1 means 2 rows, 1 column, 1st plot



Now use the plotting commands just like before.

plt.plot(logL,mass,linewidth=0,marker='s',markersize=9,markerfacecolor='r',alpha=.5,markeredgecolor='b')
plt.title('Subplot 1 - Everything Together',fontsize=12)
plt.xlabel('log(Luminosity) ',fontsize=12)
plt.ylabel(r'log(M$_{H_{2}}$) ',fontsize=12)



Make the second subplot

sub2=fig1.add_subplot(2,1,2)
# 2 rows, 1 column, 2nd plot



Plot the upper limits with pyplot.errorbar(xlocs,ylocs,xerr_length,yerr_length, other kwargs)
fmt=None means no marker for the point, lolims turns the bars into upper limit arrows, capsize=arrowhead, elinewidth=the width of the arrow shaft, mew=markeredgewidth= width of the bar

plt.errorbar(limits[1],limits[0],xerr=None,yerr=.15,fmt=None,ecolor='k', \
lolims=True,capsize=3,elinewidth=1,mew=0,label='Upper Lims')



Then plot as usual for the regular points

plt.plot(regular[1],regular[0],linewidth=0,marker='d',markerfacecolor='g',alpha=.5,label='Exact Mass')
plt.title('Subplot 2 - Separated',fontsize=12)
plt.xlabel('log(Luminosity) ',fontsize=12)
plt.ylabel(r'log(M$_{H_{2}}$) ',fontsize=12)



Change the limits, for kicks

plt.xlim(18,23)
plt.ylim(6,10)



Let’s add a legend, using the labels we defined in our plot calls. Properties defined in a dictionary.

plt.legend(loc='upper left',prop={'size':8},numpoints=1) #More args: ncol=1,loc=1,fancybox=True ...


fig1.suptitle('Mass vs. Luminosity',fontsize=16)



Tidy up the subplot spacing – numbers are fractions of subplot sizes.

plt.subplots_adjust(left=.2,right=.8,bottom=.15,top=.85,hspace=.5)



Make the tick labels smaller. This can also control putting ticks on top, etc…  You can also use plt.ticklabel_format() to change the notation style

plt.tick_params(axis='both',labelsize=10)



Change the font of all the text and labels to standard serif

plt.rc('font',family='serif')


Finally, show or save.

#plt.show()
plt.savefig('plottest2.png')

plt.clf()



Here is the resulting plot:

To recap, here is a list of some of the commands I used to plot and a quick explanation of some of the most useful keyword arguments (with example values). See the links for the official man page that explains all the keywords! The general page is
http://matplotlib.org/api/pyplot_api.html#module-matplotlib.pyplot

pyplot.plot(x_array,y_array,*kwargs)

• marker=’s’ –> The marker to use for data points.
• color=’red’ –> The color to use for plotting.  Accepts standard python color names, 0.0-1.0 for grayscale, or hex.
• linewidth=0.5 –> The thickness of the line between points.  Set to 0 for no line.
• markersize=10 –> The point size of the marker.  The shorthand is ms=…
• markerfacecolor=’b’ –> The color of the marker’s body.  The shorthand is mfc=…
• markeredgecolor=’0.5′ –> The color of the marker’s outline.  The shorthand is mec=…
• alpha=0.5 –> Alpha sets the transparency level

pyplot.legend(*kwargs)

• Note that there are some subtleties with plt.legend(), depending on if you call it as an object or not.  (i.e., leg=plt.legend()…) and if you created your plot instances as objects.
• loc=3 –> The location where the legend is drawn
• numpoints=1 –> The number of merkers to show in the legend
• handlelength=0.5 –> The length of line to draw (in points) in the legend
• ncol=4 –> The number of columns to use in the legend
• fancybox=True –> Use a box with rounded edges
• prop={‘size’:8} –> This sets the font properties, such as size

pyplot.title(‘MyTitle’,*kwargs)   Also see pyplot.suptitle()

• loc=’left’ –> Justify in center, left, or right
• x=0.2, y=0.95 –> Manually set x,y coordinates (in figure fraction)

• Left, right, top, bottom, wspace, hspace

pyplot.savefig(‘MyPlotTitle.pdf’,*kwargs)

• Note that MyPlotTitle can be a full path, so you can save in any directory you like.
• Also note that the pyplot figures out file savetype from the extension you put in the title (png, pdf, eps…)
• dpi=300 –> The resolution of the output image in dots per inch
• facecolor=’0.8′ –> Sets the color of the plot border
• transparent=True –> Set the axes patch and figure patch to be transparent
• bbox_inches=’tight –> Sets the the bounding box size. ‘Tight’ forces a crop to cut out whitespace.

### ==== EVERYTHING CONDENSED ====

Here is a condensed version of the commands.  You can save this as plottingexample.py and run from the command line.

#!/usr/bin/python

import numpy
import csv
from matplotlib import pyplot as plt

data=numpy.genfromtxt('SampleData.csv',skiprows=1,delimiter=',')
luminosity=data[:,1]; logL=numpy.log10(luminosity)
mass=data[:,0]
plt.plot(numpy.log10(luminosity),mass,'*',linewidth=0)
plt.title('M vs. L')
plt.xlabel('log(Luminosity) in some units...')
plt.ylabel(r'log(M$_{H_{2}}$) in some units...')
plt.show()
#plt.savefig('plottest1.png')
plt.clf()

mass=numpy.genfromtxt('SampleData2.csv',skiprows=1,delimiter=',',usecols=0)
logL=numpy.log10(numpy.genfromtxt('SampleData2.csv',skiprows=1,delimiter=',',usecols=2))
limits=[]; regular=[]
for i in range(0,len(mass)):
else: regular.append([mass[i],logL[i]])

limits=numpy.array(limits).transpose()
regular=numpy.array(regular).transpose()

fig1=plt.figure(1)
plt.plot(logL,mass,linewidth=0,marker='s',markersize=9,markerfacecolor='r',alpha=.5,markeredgecolor='b')
plt.title('Subplot 1 - Everything Together',fontsize=12)
plt.xlabel('log(Luminosity) ',fontsize=12)
plt.ylabel(r'log(M$_{H_{2}}$) ',fontsize=12)

plt.errorbar(limits[1],limits[0],xerr=None,yerr=.15,fmt=None,ecolor='k',lolims=True,capsize=3,elinewidth=1,mew=0,label='Upper Lims')
plt.plot(regular[1],regular[0],linewidth=0,marker='d',markerfacecolor='g',alpha=.5,label='Exact Mass')
plt.title('Subplot 2 - Separated',fontsize=12)
plt.xlabel('log(Luminosity) ',fontsize=12)
plt.ylabel(r'log(M$_{H_{2}}$) ',fontsize=12)

plt.xlim(18,23); plt.ylim(6,10)
plt.legend(loc='upper left',prop={'size':8},numpoints=1)
fig1.suptitle('Mass vs. Luminosity',fontsize=16)
plt.tick_params(axis='both',labelsize=10)
plt.rc('font',family='serif')

plt.show()
#plt.savefig('plottest2.png')
plt.clf()



### Sample csv files used in this tutorial

"log M(H2) [Msun]","Luminosity [W/Hz]"
8.02,3.34E+19
7.94,3.81E+19
7.41,4.00E+20
8,3.91E+19
7.27,3.56E+20
7.55,1.76E+19
7.46,2.86E+20
7.63,2.69E+19
6.96,4.94E+18
7.19,3.17E+20
7.35,1.23E+19
7.73,9.42E+19
7.62,
7.28,1.72E+19
7.78,1.35E+19
7.36,1.66E+19
7.4,9.88E+18
7.66,1.83E+19
8.21,4.57E+20
7.77,5.82E+19
7.78,3.70E+19
7.91,5.81E+20
7.48,1.20E+19
7.53,
6.89,1.61E+19
7.5,1.74E+20
7.06,3.96E+21
8.13,8.32E+20
7.98,2.01E+20
7.22,2.98E+20
7.2,1.71E+20
7.74,
7.39,2.17E+20
7.11,1.62E+19
7.45,1.00E+22
7.12,1.11E+19
7.39,9.38E+18
7.25,1.21E+19
8.72,5.92E+20
7.58,9.38E+18
8.28,2.55E+20
7.31,1.23E+20
7.76,7.83E+21
7.79,2.58E+19
7.94,2.63E+19
8.33,9.80E+19
7.78,9.16E+19
7.85,3.30E+19
7.72,7.45E+20
7.6,2.40E+19
7.6,1.33E+19
7.82,2.08E+19
7.92,2.68E+19
7.67,1.92E+19
7.57,1.14E+19
7.68,8.57E+19
7.64,4.02E+19
8.47,4.77E+20
7.9,2.38E+19
7.98,2.66E+20
8,4.94E+20
7.58,1.14E+19
7.62,6.45E+19
7.83,5.69E+19
7.49,1.53E+19
7.41,1.40E+21
7.39,1.55E+19
6.96,4.31E+18
8.79,1.89E+21
8.03,1.10E+21
8.53,1.81E+21
7.9,4.27E+19
7.98,3.38E+19
7.92,3.60E+19
8,5.97E+19
7.65,7.36E+19
7.81,1.62E+19
7.62,1.98E+19
8.32,1.13E+21
8.58,1.34E+21
7.8,2.43E+19
8.77,3.52E+19
7.75,6.98E+19
7.87,2.18E+20
7.89,3.08E+19
7.52,1.64E+19
7.79,6.59E+19
7.54,7.42E+19
7.87,2.67E+19
9.19,5.14E+21
7.48,2.31E+19
8.65,1.83E+21
7.68,9.61E+19
7.61,1.71E+19
7.44,7.32E+18
7.12,1.43E+19
7.52,3.68E+19
7.47,1.24E+19
8.33,3.68E+20
7.91,5.20E+19

"log M(H2) [Msun]","Limit?","Luminosity [W/Hz]"
8.02,"Yes",33368520915374800000
7.94,"Yes",38084926326744200000
7.41,"Yes",400085178787008000000
8,"Yes",39093499194531000000
7.27,"No",356189034772789000000
7.55,"Yes",17607223238288400000
7.46,"Yes",286360603640652000000
7.63,"Yes",26903793690921700000
6.96,"Yes",4940972388221380000
7.19,"Yes",316985913740269000000
7.35,"Yes",12305895234497800000
7.73,"Yes",94207316869183400000
7.62,"Yes",
7.28,"Yes",17158987504017800000
7.78,"Yes",13459820262483600000
7.36,"No",16581684400708200000
7.4,"Yes",9876452632728350000
7.66,"Yes",18319504517883800000
8.21,"No",457109684170750000000
7.77,"Yes",58175985682106300000
7.78,"Yes",37032397255629900000
7.91,"Yes",580913054903842000000
7.48,"Yes",11957640390437600000
7.53,"Yes",
6.89,"Yes",16119768544855700000
7.5,"Yes",174471237022577000000
7.06,"Yes",3962390085961550000000
8.13,"No",831966530019003000000
7.98,"Yes",200540918322003000000
7.22,"No",298216099008093000000
7.2,"Yes",171430526010093000000
7.74,"Yes",
7.39,"No",217288431492744000000
7.11,"Yes",16244996383071900000
7.45,"Yes",10007558570290300000000
7.12,"Yes",11062127409777000000
7.39,"Yes",9376584789089050000
7.25,"Yes",12138109293741500000
8.72,"No",592494207532420000000
7.58,"Yes",9381812413953460000
8.28,"No",255214352066382000000
7.31,"No",123266901644388000000
7.76,"Yes",7834806106838980000000
7.79,"Yes",25784508237840400000
7.94,"Yes",26342837758954400000
8.33,"No",97954524158250200000
7.78,"Yes",91574782751217700000
7.85,"Yes",32966512351829900000
7.72,"Yes",744943872339070000000
7.6,"Yes",23984795989100900000
7.6,"Yes",13340264557425200000
7.82,"Yes",20838405136609500000
7.92,"Yes",26812741409761600000
7.67,"Yes",19167571244596500000
7.57,"Yes",11410282942283200000
7.68,"Yes",85709161169739700000
7.64,"Yes",40159401069740800000
8.47,"No",477030827216338000000
7.9,"Yes",23798511990555400000
7.98,"Yes",265559967454768000000
8,"Yes",493856772637804000000
7.58,"Yes",11432674116179900000
7.62,"Yes",64528933270989600000
7.83,"No",56854874842154900000
7.49,"Yes",15295148643743600000
7.41,"Yes",1396303374284150000000
7.39,"Yes",15549742541566500000
6.96,"Yes",4311608709239170000
8.79,"No",1891600184373810000000
8.03,"Yes",1096533532393990000000
8.53,"No",1808214825980810000000
7.9,"Yes",42676953645405700000
7.98,"Yes",33773625486922300000
7.92,"Yes",36021130653426900000
8,"No",59720274241014000000
7.65,"Yes",73563528146885700000
7.81,"Yes",16179737317458500000
7.62,"Yes",19779722768637700000
8.32,"No",1129890630260920000000
8.58,"No",1339581028953160000000
7.8,"Yes",24344112210928900000
8.77,"No",35223633965183700000
7.75,"Yes",69758692223945700000
7.87,"Yes",217710841075595000000
7.89,"Yes",30761806253667100000
7.52,"Yes",16390578158495200000
7.79,"Yes",65929145637791600000
7.54,"Yes",74233946864774700000
7.87,"Yes",26738071215568200000
9.19,"No",5137401048219190000000
7.48,"Yes",23142567658709700000
8.65,"No",1828953187691050000000
7.68,"Yes",96082657044091600000
7.61,"Yes",17061710594939600000
7.44,"Yes",7322727758242710000
7.12,"Yes",14288481341147900000
7.52,"No",36800735135954300000
7.47,"Yes",12386569692225800000
8.33,"No",368338671117939000000
7.91,"Yes",52012168051107600000


# Curve Fitting

This post goes through the basics of fitting functions to data sets.  There are many tricks to optimizing fitting foutines for speed, but I will only present a general method here in the interest of keeping things simple and easy to follow.  We start out with some generated data, and finish by fitting curves to a real data set.

We will also explore the possibilities of giving certain parameters limits in our fits.  Python has some great tools to help you with all your fitting needs!

For this tutorial, I am using:

•  python v2.7
•  numpy v1.5.1
•  matplotlib v1.1.0
•  scipy v0.8.0
•  lmfit v0.4 – for constraining parameters
•  pyfits 2.3.1 – for I/O of the .fits file

The basic method I will outline below:

• Load or define the data array that you want to fit
• Define a function for your desired equation
• Fit with scipy.optimize.leastsq  (or with lmfit.minimize, if you want to constrain certain parameters)

Import the modules we will use

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import leastsq
from lmfit import minimize, Parameters
import pyfits


Let’s start with a simple example: we’ll generate some data with random noise, and fit a quadratic curve. Then we’ll plot the data, the fit, and residuals.

Creating the example data.  X-range from -3 to 10 in increments of 0.2.  Our y-data will start out as x2

xarray1=np.arange(-3,10,.2)
yarray1=xarray1**2


Now add gaussian random noise to each element.

yarray1+=np.random.randn(yarray1.shape[0])*2


There is a cheap and easy way to perform fits with little setup if you are only interested in polynomials.  scipy.polyfit() allows you to give an x-array and a y-array as input, with the order of polynomial to fit, then spits out an array of the fit coefficients with little fuss.  To fit a quadratic to our data generated above, for example:

from scipy import polyfit

fitcoeffs=polyfit(xarray1,yarray1,2)

print fitcoeffs
# --> Returns array (<coeff. order 2>, <coeff. order 1>, <coeff. order 0>)


If we want to fit an arbitrary expression, though, we must define a python function which will compute our desired equation.  Our function will take two arrays as inputs: the list of coefficiencts (the parameters that we will fit – I will label them “p”), and the array of x-axis data points. Later we will define the initial guesses for the parameters, called “p0”.Our function will take the form

def quadratic(p,x):
y_out=p[0]*(x**2)+p[1]*x+p[2]
return y_out


There is a simpler way to do this for short/simple expressions: python’s lambda functions are essentially one-line function definitions.

quadratic = lambda p,x: p[0]*(x**2)+p[1]*x+p[2]


We will also need to define a function of the difference between the fitted function with the current iteration’s parameters and the original data.

quadraticerr = lambda p,x,y: quadratic(p,x)-y


Define our initial guesses of the quadratic’s coefficients (I’m making bad guesses here on purpose, just to show that it still works)

p0=[2,2,2]


Now for the actual fitting.  We will tell scipy.optimize.leastsq to minimize the DIFFERENCE between the original data and iterated set of data given by the defined equation with the current coefficients.  So we tell it to perform a least-squares fit on “quadraticerr”, and give it the x- and y-data arrays as arguments.

fitout=leastsq(quadraticerr,p0[:],args=(xarray1,yarray1))
paramsout=fitout[0]
#covar=out[1] #This is the covariance matrix output
print('Fitted Parameters:\na = %.2f , b = %.2f , c = %.2f' % (paramsout[0],paramsout[1],paramsout[2]))


Now let’s plot our results. I will create a smoother curve for the fit on the fly within pyplot.  I will also add a text box with the fit parameters and plot the residuals just below the data.

plt.rc('font',family='serif')
fig1=plt.figure(1)
#xstart, ystart, xwidth, yheight --> units are fraction of the image from bottom left


Plot the data in frame1

xsmooth=np.linspace(xarray1[0],xarray1[-1])
plt.plot(xarray1,yarray1,'r.')
frame1.set_xticklabels([]) #We will plot the residuals below, so no x-ticks on this plot
plt.ylabel('y-data')
plt.grid(True)
frame1.annotate('$y$ = %.2f$\cdot x^2$+%.2f$\cdot x$+%.2f'%(paramsout[0],paramsout[1],paramsout[2]), \
xy=(.05,.95),xycoords='axes fraction',ha="left",va="top",bbox=dict(boxstyle="round", fc='1'))


If we want to keep the plot tidy:

from matplotlib.ticker import MaxNLocator
plt.gca().yaxis.set_major_locator(MaxNLocator(prune='lower')) #Removes lowest ytick label


Now plot the residuals in frame2

frame2=fig1.add_axes((.1,.1,.8,.2))

plt.ylabel('Residuals')
plt.grid(True)

#plt.show()

plt.clf()


Here is the result:

### ==== CONSTRAINING PARAMETERS ====

Now we will set limits for what values the coefficients can be fit to. This is particularly useful if you know (or want to impose) on the physics of the phenomenon represented by the data. scipy.optimize.minimize() has this functionality, and I will probably add an example using this in the future.  Another useful module is lmfit, which I will use from here on out.  Here we will fit a power law to data for the M-Sigma relation (data from Gultekin et al. 2009 – http://arxiv.org/abs/0903.4897 )

dat=np.genfromtxt('M-SigmaData.csv',skiprows=2)
M=dat[:,0]      #The y-data
Sigma=dat[:,1]  #The x-data


The general (simple) form for a power law is .  We can do this simply enough directly, or we can be clever and note that we could fit a line for .

lmfit requires us to define the coefficients as Parameters() objects. We then define their initial guesses, whether or not they can vary, and what their min/max values can be.  Finally, when we define our equation functions, we must call the parameters with p0[‘par_x’].value :

params=Parameters()

power=lambda p,x: params['alpha'].value+params['beta'].value*x
powererr=lambda p,x,y: power(p,x)-y


Perform the fit with lmfit.minimize.  This output gives us a whole suite of info, such as reduced chi-square values.  With lmfit, the parameters will be updated within the Parameters() object

fitout2=minimize(powererr,params,args=(np.log10(Sigma),np.log10(M)))

print('Fitted Parameters:\nalpha = %.2f , beta = %.2f' % (params['alpha'].value,params['beta'].value))

xsmooth2=np.linspace(np.min(Sigma),np.max(Sigma))
fit2=10**params['alpha'].value*xsmooth2**params['beta'].value


Plot the data and the fit

fig2=plt.figure(2)

plt.plot(Sigma,M,'b*',linewidth=0)
plt.plot(xsmooth2,fit2,'g')
plt.yscale('log')
plt.xscale('log')
plt.xlim(np.min(Sigma),np.max(Sigma))
plt.xlabel('$\sigma$ (km s$^{-1}$)')
plt.ylabel('M (M$_{\odot}$)')
plt.title('Power Law Fit Example')
sp1.annotate('log($y$) = %.2f+%.2f$\cdot$log($x$)\n$y$ = %.2f$\cdot x^{%.2f}$' \

%(params['alpha'].value,params['beta'].value,10**params['alpha'].value,params['beta'].value),xy=(.05,.95), \
xycoords='axes fraction',ha="left", va="top",bbox=dict(boxstyle="round", fc='1'))

plt.savefig('power_fit_example.png')
#plt.show()

plt.clf()


Here is the result:

Now is a good time to mention that lmfit.minimize() gives you some freebie goodness-of-fit diagnostics:

 nfev number of function evaluations success boolean (True/False) for whether fit succeeded. errorbars boolean (True/False) for whether uncertainties were estimated. message message about fit success. ier integer error value from scipy.optimize.leastsq lmdif_message message from scipy.optimize.leastsq nvarys number of variables in fit ndata number of data points: nfree degrees of freedom in fit: residual residual array (return of func(): chisqr chi-square: redchi reduced chi-square:

Call with, e.g., output.chisqr or output.residual

Again, see http://cars9.uchicago.edu/software/python/lmfit/ for details.

### ==== FITTING GAUSSIANS TO DATA FROM A .FITS FILE ====

Next, let’s read in data from a .fits file and perform a fit.  This particular  data cube is a neutral hydrogen VLA observation of dwarf galaxy WLM.  The cube  can be downloaded (~430MB) from  https://science.nrao.edu/science/surveys/littlethings/data/wlm.html  –> I have reproduced the ‘stackspec’ and ‘vels’ arrays at the end of this tutorial. Specifically, I will sum up all of the x-y data to make one big stacked z-array, or ‘Super Profile’.  Then I will fit various gaussian functions (Gaus + Hermite  polynomials and double gaussian) and plot them with their residuals.  Of course,  in real life you would want to ‘Shuffle’ the spectra and account for other factors, but this is just an example.

I am using the default Marquardt-Levenberg algorithm.  Note that fitting results will depend quite a bit on what you give as initial guesses – ML finds LOCAL  extrema quite well, but it doesn’t necessarily find the global extrema.  In short, do your best to provide a good first guess to the fit parameters.

Load in the data. We don’t care about the 4th axis, so cut it out. Also get the header and specify the velocity parameters.


cube=pyfits.getdata('WLM_NA_ICL001.FITS')[0,:,:,:]

cdelt3=cubehdr['CDELT3']/1000.; crval3=cubehdr['CRVAL3']/1000.; crpix3=cubehdr['CRPIX3'];
minvel=crval3+(-crpix3+1)*cdelt3; maxvel=crval3+(cube.shape[0]-crpix3)*cdelt3
chanwidth=abs(cdelt3)

#minvel=20.452084230999986; maxvel=-252.452084231
#cdelt3=-2.574567627; chanwidth=2.574567627


Stack all the spectra

stackspec=np.sum(np.sum(cube,axis=2),axis=1)

vels=np.arange(minvel,maxvel+int(cdelt3),cdelt3)


The velocity array in the cube goes from positive to negative, so let’s reverse it to make the fitting go smoother.

vels=vels[::-1]; stackspec=stackspec[::-1]


Alternatively, to skip using pyfits and load the velocity and spectrum data from the .csv file provided below:

    #dat=np.genfromtxt('GausFitData.csv',skiprows=1)
#vels=dat[:,0]
#stackspec=dat[:,1]
#chanwidth=2.574567627


Next, we define our fit functions:

Gaussian plus Hermite Polynomials H3/H4   (labeled below as ~_gh)
, where

p_gh=Parameters()

def gaussfunc_gh(paramsin,x):
amp=paramsin['amp'].value
center=paramsin['center'].value
sig=paramsin['sig'].value
c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
skew=paramsin['skew'].value
kurt=paramsin['kurt'].value

gaustot_gh=amp*np.exp(-.5*((x-center)/sig)**2)*(1+skew*(c1*((x-center)/sig)+c3*((x-center)/sig)**3)+kurt*(c5+c2*((x-center)/sig)**2+c4*((x-center)/sig)**4))
return gaustot_gh


Double Gaussian   (labeled below as ~_2g)

Bounds–> amp: 10% of max to max   center: velocity range   disp: channel width to velocity range

p_2g=Parameters()

def gaussfunc_2g(paramsin,x):
amp1=paramsin['amp1'].value; amp2=paramsin['amp2'].value;
center1=paramsin['center1'].value; center2=paramsin['center2'].value;
sig1=paramsin['sig1'].value; sig2=paramsin['sig2'].value;

gaus1=amp1*np.exp(-.5*((x-center1)/sig1)**2)
gaus2=amp2*np.exp(-.5*((x-center2)/sig2)**2)
gaustot_2g=(gaus1+gaus2)
return gaustot_2g


And now the functions that compute the difference between the fit iteration and data. In addition, define a function for a simple single gaussian.

gausserr_gh = lambda p,x,y: gaussfunc_gh(p,x)-y
gausserr_2g = lambda p,x,y: gaussfunc_2g(p,x)-y
gausssingle = lambda a,c,sig,x: a*np.exp(-.5*((x-c)/sig)**2)


We will minimize with lmfit again, in order to keep limits on parameters

fitout_gh=minimize(gausserr_gh,p_gh,args=(vels,stackspec))
fitout_2g=minimize(gausserr_2g,p_2g,args=(vels,stackspec))


Let’s make some shorthand definitions to keep things tidy

pars_gh=[p_gh['amp'].value,p_gh['center'].value,p_gh['sig'].value,p_gh['skew'].value,p_gh['kurt'].value]
pars_2g=[p_2g['amp1'].value,p_2g['center1'].value,p_2g['sig1'].value,p_2g['amp2'].value,p_2g['center2'].value,p_2g['sig2'].value]


Aside:
If you want to use the simple fitter, without defining bounds, revert  to scipy.optimize.leastsq().  Note, however, that with complex functions such as two gaussians, you run the risk of getting nonphysical results, for example one gaussian being negative.  That said, here is the code:

pars1=[np.max(stackspec),vels[50],2*abs(cdelt3),0,0]
pars2=[np.max(stackspec),vels[50],2*abs(cdelt3),np.max(stackspec),vels[50],2*cdelt3]

def gaussfunc_gh(p,x):
c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
return p[0]*np.exp(-.5*((x-p[1])/p[2])**2)*(1+p[3]*(c1*((x-p[1])/p[2])+c3*((x-p[1])/p[2])**3)+p[4]*(c5+c2*((x-p[1])/p[2])**2+c4*((x-p[1])/p[2])**4))

def gaussfunc_2g(p,x):
gaus1=p[0]*np.exp(-.5*((x-p[1])/p[2])**2)
gaus2=p[3]*np.exp(-.5*((x-p[4])/p[5])**2)
return gaus1+gaus2

fitout_gh=leastsq(gausserr_gh,pars1[:],args=(vels,stackspec))
pars_gh=fitout_gh[0]
fitout_2g=leastsq(gausserr_2g,pars2[:],args=(vels,stackspec))
pars_2g=fitout_2g[0]

print('Fitted Parameters (Gaus+Hermite):\nAmp = %.2f , Center = %.2f , Disp = %.2f\nSkew = %.2f , Kurt = %.2f' % (pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]))
print('Fitted Parameters (Double Gaussian):\nAmp1 = %.2f , Center1 = %.2f , Sig1 = %.2f\nAmp2 = %.2f , Center2 = %.2f , Sig2 = %.2f' % (pars_2g[0],pars_2g[1],pars_2g[2],pars_2g[3],pars_2g[4],pars_2g[5]))


[End aside]

fit_gh=gaussfunc_gh(pars_gh,vels)
fit_2g=gaussfunc_2g(pars_2g,vels)
resid_gh=fit_gh-stackspec
resid_2g=fit_2g-stackspec


Print the fitted parameters

print('Fitted Parameters (Gaus+Hermite):\nAmp = %.2f , Center = %.2f , Disp = %.2f\nSkew = %.2f , Kurt = %.2f' \
%(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]))

print('Fitted Parameters (Double Gaussian):\nAmp1 = %.2f , Center1 = %.2f , Sig1 = %.2f\nAmp2 = %.2f , Center2 = %.2f , Sig2 = %.2f' \
%(pars_2g[0],pars_2g[1],pars_2g[2],pars_2g[3],pars_2g[4],pars_2g[5]))


Now let’s plot the results

fig3=plt.figure(3)
#xstart, ystart, xwidth, yheight --> units are fraction of the image from bottom left


Plot the data in frame1

plt.plot(vels,stackspec,'k.')
pgh,=plt.plot(vels,fit_gh,'b')
p2g,=plt.plot(vels,fit_2g,'r')
p2ga,=plt.plot(vels,gausssingle(pars_2g[0],pars_2g[1],pars_2g[2],vels),'-.',color='orange')
p2gb,=plt.plot(vels,gausssingle(pars_2g[3],pars_2g[4],pars_2g[5],vels),'-.',color='green')
f1.set_xticklabels([]) #We will plot the residuals below, so no x-ticks on this plot
plt.title('Multiple Gaussian Fit Example')
plt.ylabel('Amplitude (Some Units)')
f1.legend([pgh,p2g,p2ga,p2gb],['Gaus-Hermite','2-Gaus','Comp. 1','Comp2'],prop={'size':10},loc='center left')

from matplotlib.ticker import MaxNLocator
plt.gca().yaxis.set_major_locator(MaxNLocator(prune='lower')) #Removes lowest ytick label


Let’s also put in some text boxes with the fit results

f1.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f' \
%(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]),xy=(.05,.95), \
xycoords='axes fraction',ha="left", va="top", \
bbox=dict(boxstyle="round", fc='1'),fontsize=10)
f1.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f' \
%(pars_2g[0],pars_2g[3],pars_2g[1],pars_2g[4],pars_2g[2],pars_2g[5]),xy=(.95,.95), \
xycoords='axes fraction',ha="right", va="top", \
bbox=dict(boxstyle="round", fc='1'),fontsize=10)


Now plot the residuals in frame2

f2=fig3.add_axes((.1,.1,.8,.2))

resgh,res2g,=plt.plot(vels,resid_gh,'k--',vels,resid_2g,'k')

plt.ylabel('Residuals')
plt.xlabel('Velocity (km s$^{-1}$)')
f2.legend([resgh,res2g],['Gaus-Hermite','2-Gaus'],numpoints=4,prop={'size':9},loc='upper left')

plt.savefig('gaussian_fit_example.png',dpi=100)
#plt.show()

plt.clf()


Finally, here is the result:

### ==== EVERYTHING CONDENSED ====

Here is a condensed version of the commands.

Condensed Code
#!/usr/bin/python

import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import leastsq
from lmfit import minimize, Parameters
import pyfits

### Example 1: Quadratic fitting with scipy.optimize.leastsq() ###
xarray1=np.arange(-3,10,.2)
yarray1=xarray1**2
yarray1+=np.random.randn(yarray1.shape[0])*2

p0=[2,2,2] #Initial guesses

paramsout=fitout[0]
print('Fitted Parameters:\na = %.2f , b = %.2f , c = %.2f' % (paramsout[0],paramsout[1],paramsout[2]))

#Plotting
plt.rc('font',family='serif')
fig1=plt.figure(1)
frame1=fig1.add_axes((.1,.3,.8,.6)) #xstart, ystart, xwidth, yheight --> units are fraction of the image from bottom left

xsmooth=np.linspace(xarray1[0],xarray1[-1])
plt.plot(xarray1,yarray1,'r.')
frame1.set_xticklabels([]) #We will plot the residuals below, so no x-ticks on this plot
plt.ylabel('y-data')
plt.grid(True)
frame1.annotate('$y$ = %.2f$\cdot x^2$+%.2f$\cdot x$+%.2f' \
%(paramsout[0],paramsout[1],paramsout[2]),xy=(.05,.95), \
xycoords='axes fraction',ha="left", va="top", \
bbox=dict(boxstyle="round", fc='1'))

from matplotlib.ticker import MaxNLocator
plt.gca().yaxis.set_major_locator(MaxNLocator(prune='lower')) #Removes lowest ytick label

plt.ylabel('Residuals')
plt.grid(True)

#plt.show()
plt.clf()

### Example 2: power law with lmfit.minimize() ###

M=dat[:,0]
Sigma=dat[:,1]

params=Parameters()

power=lambda p,x: params['alpha'].value+params['beta'].value*x
powererr=lambda p,x,y: power(p,x)-y

fitout2=minimize(powererr,params,args=(np.log10(Sigma),np.log10(M)))
print('Fitted Parameters:\nalpha = %.2f , beta = %.2f' % (params['alpha'].value,params['beta'].value))

xsmooth2=np.linspace(np.min(Sigma),np.max(Sigma))
fit2=10**params['alpha'].value*xsmooth2**params['beta'].value

# Plot the data and the fit
fig2=plt.figure(2)
plt.plot(Sigma,M,'b*',linewidth=0)
plt.plot(xsmooth2,fit2,'g')
plt.yscale('log')
plt.xscale('log')
plt.xlim(np.min(Sigma),np.max(Sigma))
plt.xlabel('$\sigma$ (km s$^{-1}$)')
plt.ylabel('M (M$_{\odot}$)')
plt.title('Power Law Fit Example')
sp1.annotate('log($y$) = %.2f+%.2f$\cdot$log($x$)\n$y$ = %.2f$\cdot x^{%.2f}$' \
%(params['alpha'].value,params['beta'].value,10**params['alpha'].value,params['beta'].value),xy=(.05,.95), \
xycoords='axes fraction',ha="left", va="top", \
bbox=dict(boxstyle="round", fc='1'))

plt.savefig('power_fit_example.png')
#plt.show()
plt.clf()

### Example 3: multiple Gaussians from data cube with parameter limits ###

#cube=pyfits.getdata('WLM_NA_ICL001.FITS')[0,:,:,:]

#cdelt3=cubehdr['CDELT3']/1000.; crval3=cubehdr['CRVAL3']/1000.; crpix3=cubehdr['CRPIX3'];
#minvel=crval3+(-crpix3+1)*cdelt3; maxvel=crval3+(cube.shape[0]-crpix3)*cdelt3
#chanwidth=abs(cdelt3)

#stackspec=np.sum(np.sum(cube,axis=2),axis=1) # Stack all the spectra
#vels=np.arange(minvel,maxvel+int(cdelt3),cdelt3)
#vels=vels[::-1]; stackspec=stackspec[::-1] #Reversing the arrays to go from - to +

dat=np.genfromtxt('GausFitData.csv',skiprows=1)
vels=dat[:,0]
stackspec=dat[:,1]
chanwidth=2.574567627

# Gaussian plus Hermite Polynomials H3/H4   (labeled below as ..._gh)
p_gh=Parameters()

def gaussfunc_gh(paramsin,x):
amp=paramsin['amp'].value
center=paramsin['center'].value
sig=paramsin['sig'].value
c1=-np.sqrt(3); c2=-np.sqrt(6); c3=2/np.sqrt(3); c4=np.sqrt(6)/3; c5=np.sqrt(6)/4
skew=p_gh['skew'].value
kurt=p_gh['kurt'].value
return amp*np.exp(-.5*((x-center)/sig)**2)*(1+skew*(c1*((x-center)/sig)+c3*((x-center)/sig)**3)+kurt*(c5+c2*((x-center)/sig)**2+c4*((x-center)/sig)**4))

# Double Gaussian   (labeled below as ..._2g)
p_2g=Parameters()

def gaussfunc_2g(paramsin,x):
amp1=paramsin['amp1'].value; amp2=paramsin['amp2'].value;
center1=paramsin['center1'].value; center2=paramsin['center2'].value;
sig1=paramsin['sig1'].value; sig2=paramsin['sig2'].value;
gaus1=amp1*np.exp(-.5*((x-center1)/sig1)**2)
gaus2=amp2*np.exp(-.5*((x-center2)/sig2)**2)
return (gaus1+gaus2)

gausserr_gh = lambda p,x,y: gaussfunc_gh(p,x)-y
gausserr_2g = lambda p,x,y: gaussfunc_2g(p,x)-y
gausssingle = lambda a,c,sig,x: a*np.exp(-.5*((x-c)/sig)**2)

fitout_gh=minimize(gausserr_gh,p_gh,args=(vels,stackspec))
fitout_2g=minimize(gausserr_2g,p_2g,args=(vels,stackspec))

pars_gh=[p_gh['amp'].value,p_gh['center'].value,p_gh['sig'].value,p_gh['skew'].value,p_gh['kurt'].value]
pars_2g=[p_2g['amp1'].value,p_2g['center1'].value,p_2g['sig1'].value,p_2g['amp2'].value,p_2g['center2'].value,p_2g['sig2'].value]
fit_gh=gaussfunc_gh(pars_gh,vels)
fit_2g=gaussfunc_2g(pars_2g,vels)
resid_gh=fit_gh-stackspec
resid_2g=fit_2g-stackspec

print('Fitted Parameters (Gaus+Hermite):\nAmp = %.2f , Center = %.2f , Disp = %.2f\nSkew = %.2f , Kurt = %.2f' \
%(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]))
print('Fitted Parameters (Double Gaussian):\nAmp1 = %.2f , Center1 = %.2f , Sig1 = %.2f\nAmp2 = %.2f , Center2 = %.2f , Sig2 = %.2f' \
%(pars_2g[0],pars_2g[1],pars_2g[2],pars_2g[3],pars_2g[4],pars_2g[5]))

#Plotting
fig3=plt.figure(3)
plt.plot(vels,stackspec,'k.')
pgh,=plt.plot(vels,fit_gh,'b')
p2g,=plt.plot(vels,fit_2g,'r')
p2ga,=plt.plot(vels,gausssingle(pars_2g[0],pars_2g[1],pars_2g[2],vels),'-.',color='orange')
p2gb,=plt.plot(vels,gausssingle(pars_2g[3],pars_2g[4],pars_2g[5],vels),'-.',color='green')
f1.set_xticklabels([]) #We will plot the residuals below, so no x-ticks on this plot
plt.title('Multiple Gaussian Fit Example')
plt.ylabel('Amplitude (Some Units)')
f1.legend([pgh,p2g,p2ga,p2gb],['Gaus-Hermite','2-Gaus','Comp. 1','Comp2'],prop={'size':10},loc='center left')

from matplotlib.ticker import MaxNLocator
plt.gca().yaxis.set_major_locator(MaxNLocator(prune='lower')) #Removes lowest ytick label

f1.annotate('Gauss-Hermite:\nAmp = %.2f\nCenter = %.2f\n$\sigma$ = %.2f\nH3 = %.2f\nH4 = %.2f'%(pars_gh[0],pars_gh[1],pars_gh[2],pars_gh[3],pars_gh[4]),xy=(.05,.95),xycoords='axes fraction',ha="left", va="top",bbox=dict(boxstyle="round", fc='1'),fontsize=10)
f1.annotate('Double Gaussian:\nAmp$_1$ = %.2f\nAmp$_2$ = %.2f\nCenter$_1$ = %.2f\nCenter$_2$ = %.2f\n$\sigma_1$ = %.2f\n$\sigma_2$ = %.2f'%(pars_2g[0],pars_2g[3],pars_2g[1],pars_2g[4],pars_2g[2],pars_2g[5]),xy=(.95,.95),xycoords='axes fraction',ha="right", va="top",bbox=dict(boxstyle="round", fc='1'),fontsize=10)

resgh,res2g,=plt.plot(vels,resid_gh,'k--',vels,resid_2g,'k')
plt.ylabel('Residuals')
plt.xlabel('Velocity (km s$^{-1}$)')
f2.legend([resgh,res2g],['Gaus-Hermite','2-Gaus'],numpoints=4,prop={'size':9},loc='upper left')
plt.savefig('gaussian_fit_example.png',dpi=100)
#plt.show()
plt.clf()



Here is the sample data for the Gaussian fitting, for those who opt not to download the data cube.  Formatted velocities on the left, stacked spectra values on the right.

-252.452084231    7.0471
-249.877516604    5.41091
-247.302948977    -0.0805578
-244.72838135    9.80676
-242.153813723    -2.64569
-239.579246096    3.44144
-237.004678469    -0.797753
-234.430110842    -2.29659
-231.855543215    -0.593432
-229.280975588    5.25571
-226.706407961    8.75715
-224.131840334    9.0627
-221.557272707    -6.08199
-218.98270508    -3.14464
-216.408137453    -2.99775
-213.833569826    -4.07002
-211.259002199    -2.21975
-208.684434572    -2.19285
-206.109866945    -3.57909
-203.535299318    -1.38921
-200.960731691    -2.03753
-198.386164064    -9.93365
-195.811596437    -1.14918
-193.23702881    3.11781
-190.662461183    -1.26111
-188.087893556    -7.38781
-185.513325929    2.30879
-182.938758302    0.976938
-180.364190675    -4.54513
-177.789623048    -2.34182
-175.215055421    0.835287
-172.640487794    0.578699
-170.065920167    7.0345
-167.49135254    9.30069
-164.916784913    13.3074
-162.342217286    21.4498
-159.767649659    39.6186
-157.193082032    61.6891
-154.618514405    71.2766
-152.043946778    84.5251
-149.469379151    92.3035
-146.894811524    103.288
-144.320243897    103.722
-141.74567627    105.536
-139.171108643    110.039
-136.596541016    119.542
-134.021973389    112.635
-131.447405762    111.95
-128.872838135    133.203
-126.298270508    135.486
-123.723702881    145.488
-121.149135254    156.818
-118.574567627    176.232
-116    190.496
-113.425432373    198.189
-110.850864746    195.464
-108.276297119    184.057
-105.701729492    170.15
-103.127161865    149.234
-100.552594238    118.705
-97.978026611    101.411
-95.403458984    92.2328
-92.828891357    73.3222
-90.25432373    67.2127
-87.679756103    41.8572
-85.105188476    37.7619
-82.530620849    20.5096
-79.956053222    13.5301
-77.381485595    15.9504
-74.806917968    13.5154
-72.232350341    5.63028
-69.657782714    3.38025
-67.083215087    -2.62706
-64.50864746    2.74693
-61.934079833    1.86001
-59.359512206    1.59868
-56.784944579    2.25844
-54.210376952    -0.369037
-51.635809325    0.334506
-49.061241698    -0.196296
-46.486674071    -1.56728
-43.912106444    -1.1896
-41.337538817    1.94682
-38.76297119    0.268496
-36.188403563    7.77246
-33.613835936    -5.04572
-31.039268309    -4.57877
-28.464700682    5.39095
-25.890133055    -0.00106338
-23.315565428    4.24355
-20.740997801    -2.42964
-18.166430174    -3.30451
-15.591862547    -0.803972
-13.01729492    -8.57761
-10.442727293    -14.549
-7.868159666    -16.3109
-5.293592039    10.3016
-2.719024412    16.9872
-0.144456785    14.6043
2.430110842    2.62755
5.004678469    -0.388005
7.579246096    -10.0996
10.153813723    -3.00706
12.72838135    2.11268
15.302948977    5.23242
17.877516604    4.40532
20.452084231    -1.07435


# Plotting Sky Maps with the Kapteyn Package

This script demonstrates simple map visualization of astronomical data with the Kapteyn package.  There are other excellent modules for displaying astronomical maps, such as pywcsgrid2, which I will cover in later tutorials.
The optical fits file (of M101) and examples used in this script come from the Kapteyn package website, http://www.astro.rug.nl/software/kapteyn/
The VLA fits file can be found on NED.

First, import the pertinent python modules.


import numpy
import pyfits
from matplotlib import pyplot as plt
from kapteyn import maputils



The simplest way to view a fits file is to load the data into an array with pyfits, then plot the image with pyplot.imshow

* Note I’m only doing the bare minimum to load.  Much more sophisticated options are available.


dat=pyfits.getdata('m101.fits')



numpy arrays are indexed as [ydata,xdata] for 2D, or [zdata,ydata,xdata] for 3D…  By default, plt.imshow tries to plot from the top left, but pyfits reads from the bottom left.  So we specify that our first data points should be plotted at the bottom of the image


plot1=plt.imshow(dat,origin='lower')

plt.show()

plt.clf() #Clearing the plotter for the next example.


The Kapteyn package has a great wealth of built-in routines for dealing with wcs coordinates, defining colorbars, annotations, etc., and many standard matplotlib commands will work, too.  The basic steps for displaying with kapteyn.maputils:

2.  Create an Annotatedimage()
3.  display the data as Image() or Contours()
5.  Annotatedimage().plot()
m101dat=maputils.FITSimage('m101.fits')


* Note: We can also use an external dataset and header with data=maputils.FITSimage(externaldata=…,externalheader=…) There are options for loading certain slices for 3D cubes, e.g.: data.set_imageaxes(axnr1=1,axnr2=2,slicepos=17)

Set up the figure and subplot

fig1=plt.figure(1)
fig1.suptitle('Mapping Demo Using the Kapteyn Package',fontsize=16)

#In a grid of 2 rows, 1 column, add the 1st subplot
plt.title('Subplot 1 - Simple Map',fontsize=12)


Here we create the plottable image with maputils.  Let’s give it the ‘jet’ colormap for fun.  Standard matplotlib.cm maps are understood

annim=m101dat.Annotatedimage(sub1,cmap='jet')
im1=annim.Image() #kwargs include interpolation commands, etc.


If we want to display the wcs coordinates, we create a Graticule() This will automatically generate coords from a valid header

grat1=annim.Graticule()


Add a colorbar.  Note the capital “C” in the call. Specifying a matplotlib Axes frame adjusts the size.

cbar1=annim.Colorbar(fontsize=8)
cbar1.set_label(label='Some Units',fontsize=12)


Note that we are not using pyplot.plot here, but rather the Kapteyn built-in routine

annim.plot()


### ==== OVERLAYS ====

Here we will overlay optical contours on a VLA 21cm background image.  This requires regridding the second image onto the first’s wcs coordinates.

sub2=fig1.add_subplot(2,1,2) #In a grid of 2 rows, 1 column, add the 2nd subplot
plt.title('Subplot 2 - Grayscale (VLA 21cm) \nwith Heat Contours (Optical)',fontsize=12)


Load up the data for the background (optical) and overlay (VLA)

bg=maputils.FITSimage('m101_VLA.fits')
ol=maputils.FITSimage('m101.fits')


Create the re-projection object

Reprojfits=ol.reproject_to(bg)


Create the background image the same way as before, this time in grayscale

bgim=bg.Annotatedimage(sub2,cmap='gist_yarg')
bgim.Image()
grat2=bgim.Graticule()
grat2.setp_gratline(visible=False)
#This turns off the grid lines, but keeps the ticks


Now we make an Annotatedimage for the overlay based on the bg

olim=bg.Annotatedimage(sub2,cmap='gist_heat',boxdat=Reprojfits.boxdat)


Let’s define the contour levels explicitly and create the Contours() object

contlevs=numpy.arange(3e3,15e3,2e3)
cont=olim.Contours(filled=False,levels=contlevs,alpha=0.3)


We could add a beam, a skyobject polygon, a ruler, etc. here…

Now plot both the background and overlay images

bgim.plot()
olim.plot()


Let’s add a text object to sub2 with the pyplot.subplot built-in text annotator: We will print ‘M101’, in a box anchored at x=.05,y=.95 (based on the axes of the subplot), aligned to the top left of that

sub2.text(.05,.95,'M101',ha='left',va='top',transform=sub2.transAxes,size=10,bbox=None)


Tweaking the spacing of the subplots

plt.subplots_adjust(hspace=.5)


Making the default font ‘serif’

plt.rc('font',family='serif')


Finally, view the end result

plt.savefig('mapexample.png',dpi=150)
plt.show()
plt.clf()


### ==== CONDENSED COMMANDS ====

This produces the same results, just without the comments from above.

Condensed Code

#!/usr/bin/python

import numpy
import pyfits
from matplotlib import pyplot as plt
from kapteyn import maputils

dat=pyfits.getdata('m101.fits')
plot1=plt.imshow(dat,origin='lower')
plt.show()
plt.clf()

#============================================================

m101dat=maputils.FITSimage('m101.fits')
fig1=plt.figure(1)
fig1.suptitle('Mapping Demo Using the Kapteyn Package',fontsize=16)
sub1=fig1.add_subplot(2,1,1) #In a grid of 2 rows, 1 column, add the 1st subplot
plt.title('Subplot 1 - Simple Map',fontsize=12)

annim=m101dat.Annotatedimage(sub1,cmap='jet')
im1=annim.Image() #kwargs include interpolation commands, etc.
grat1=annim.Graticule()
cbar1.set_label(label='Some Units',fontsize=12)

annim.plot()

sub2=fig1.add_subplot(2,1,2) #In a grid of 2 rows, 1 column, add the 2nd subplot
plt.title('Subplot 2 - Grayscale with Heat Contours',fontsize=12)

bg=maputils.FITSimage('m101.fits')
ol=maputils.FITSimage('m101.fits') #This should be a different file, but I'll use this again for simplicity
Reprojfits=ol.reproject_to(bg)

bgim=m101dat.Annotatedimage(sub2,cmap='gist_yarg')
bgim.Image()
grat2=bgim.Graticule()
grat2.setp_gratline(visible=False) #This turns off the grid lines, but keeps the ticks

olim=bg.Annotatedimage(sub2,cmap='gist_heat',boxdat=Reprojfits.boxdat)
contlevs=numpy.arange(3e3,15e3,2e3)
cont=olim.Contours(filled=False,levels=contlevs,alpha=0.3)
bgim.plot()
olim.plot()

sub2.text(.05,.95,'M101',ha='left',va='top',transform=sub2.transAxes,size=10,bbox=None)