Plotting Sky Maps with the Kapteyn Package

mapexample


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:

  1.  Load the data into a maputils.FITSreader() object
  2.  Create an Annotatedimage()
  3.  display the data as Image() or Contours()
  4.  Add colorbar, graticule, etc…
  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)

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)

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)
      #...,frame=fig1.add_axes((0.65,0.58,.02,.32))
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()

mapexample



==== 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=annim.Colorbar(fontsize=8) #...,frame=fig1.add_axes((0.65,0.58,.02,.32))
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)
plt.subplots_adjust(hspace=.5)
plt.rc('font',family='serif')

plt.savefig('mapexample.png')
plt.show()
plt.clf()

8 thoughts on “Plotting Sky Maps with the Kapteyn Package

  1. I’m really impressed with your writing skills and also with the layout on your blog. Is this a paid theme or did you customize it yourself? Either way keep up the excellent quality writing, it is rare to see a great blog like this one nowadays..

  2. I am also using Kapteyn maputils, but I have encountered an issue. I wanted to change the tick labels of the colormap but could not! I can do that in matplotlib simply set_ticklabels(…), but I cold not find anything in Kapteyn maputils. Do you know whether we are able to change this or not?

    Cheers,
    Mustafa

    • Hello, I believe there are a few ways to do it, but basically kapteyn has its own methods for generating the labels, so I don’t think you can use the standard matplotlib commands to change ticklabels. One way that should work is to set the graticule’s tick labels with [AnnotatedImage()].Graticule().setp_ticklabel(fontsize=…) formatting. See https://www.astro.rug.nl/software/kapteyn/maputilstutorial.html#graticule-tick-labels for an example. Also, if you don’t actually want to plot the graticule lines, you can just turn them off with wcsgrat.setp_gratline(False)

Leave a Reply to manchesterafa.org Cancel reply

Your email address will not be published. Required fields are marked *