title: How could I convert a contour plot (matplotlib) to FITS format with header?
body: <p>I need to make a contour plot and overlie the contours on the image. I have used aplpy library to overlie the contours on an astronomical image.
  I have downloaded the 2MASS data in APlpy website(<a href="https://github.com/aplpy/aplpy-examples/tree/master/data" rel="nofollow">https://github.com/aplpy/aplpy-examples/tree/master/data</a>) and wrote the following piece of code:</p>
  
  <pre><code>import numpy as np
  import aplpy
  import atpy
  from pyavm import AVM
  import scipy
  import matplotlib.pyplot as plt
  import scipy.ndimage
  from matplotlib import cm
  import montage_wrapper
  import matplotlib.colors as colors
  from scipy import stats
  
  def get_contour_verts(cn):
      contours = []
      # for each contour line
      for cc in cn.collections:
          paths = []
          # for each separate section of the contour line
          for pp in cc.get_paths():
              xy = []
              # for each segment of that section
              for vv in pp.iter_segments():
                  xy.append(vv[0])
              paths.append(np.vstack(xy))
          contours.append(paths)
      return contours
  
  # Convert all images to common projection
  aplpy.make_rgb_cube(['2MASS_h.fits', '2MASS_j.fits', '2MASS_k.fits'], '2MASS_rgb.fits')
  
  # Customize it
  aplpy.make_rgb_image('2MASS_rgb.fits','2MASS_rgb_contour.png',embed_avm_tags=True)
  
  # Launch APLpy figure of 2D cube
  img = aplpy.FITSFigure('2MASS_rgb_contour.png')
  img.show_rgb()
  
  # Modify the tick labels for precision and format
  img.tick_labels.set_xformat('hhmm')
  img.tick_labels.set_yformat('ddmm')
  # Move the tick labels
  img.tick_labels.set_xposition('top')
  img.tick_labels.set_yposition('right')
  img.tick_labels.set_font(size='small', family='sans-serif', style='normal')
  
  data = scipy.loadtxt('sources.txt')
  m1=data[:,0]
  m2=data[:,1]
  xmin = m1.min()
  xmax = m1.max()
  ymin = m2.min()
  ymax = m2.max()
  #Gridding the data 
  
  X, Y = np.mgrid[xmin:xmax:100j, ymin:ymax:100j]
  positions = np.vstack([X.ravel(), Y.ravel()])
  values = np.vstack([m1, m2])
  kernel = stats.gaussian_kde(values)
  Z = np.reshape(kernel(positions).T, X.shape)
  fig = plt.figure()
  ax = fig.add_subplot(111)
  ax.imshow(np.rot90(Z), cmap=plt.cm.gist_earth_r, extent=[xmin, xmax, ymin, ymax])
  ax.plot(m1, m2, 'k.', markersize=2)
  ax.set_xlim([xmin, xmax])
  ax.set_ylim([ymin, ymax])
  
  CS = plt.contour(X, Y, Z)
  Contour_arrays=get_contour_verts(CS)
  
  #Adding contour plots
  for contours_at_level in Contour_arrays:
      radec = [img.pixel2world(*verts.T) for verts in contours_at_level]
      new_radec=[]
      for coor in radec:
          #new_radec.append(np.column_stack((coor[0], coor[1])))
          new_radec.append(np.vstack((coor[0], coor[1])))
      print new_radec
      img.show_lines(new_radec,color='red', zorder=100)
  
  img.save('tutorial.png')
  </code></pre>
  
  <p>It seems it does not still work at all.</p>