Resizing NVSS FITS format file in Python and operating on this file in astropy
<p>this questions is probably mainly for more or less advances astronomers.</p>

<p>Do you know how to transform NVSS fits file to the fits with only 2 (NOT 4!) axes? Or how to deal with the file which has 4 axis and generates the following errors in Python when I'm trying to overlapp nvss countours on optical DSS data, using astropy and other "astro" libraries for Python? (code below)</p>

<p>I've tried to do it and when there is thess 4 axes for NVSS FITS, there is error messages and warnings:</p>

<p>WARNING: FITSFixedWarning: The WCS transformation has more axes (4) than the image it is associated with (2) [astropy.wcs.wcs]
WARNING: FITSFixedWarning: 'datfix' made the change 'Invalid parameter value: invalid date '19970331''. [astropy.wcs.wcs]
<a href="http://stackoverflow.com/questions/33107224/re-sizing-a-fits-image-in-python">re-sizing a fits image in python</a></p>

<p>WARNING: FITSFixedWarning: 'datfix' made the change 'Invalid parameter value: invalid date '19970331''. [astropy.wcs.wcs]
Traceback (most recent call last):
  File "p.py", line 118, in 
    cont2 = ax[Header2].contour(opt.data, [-8,-2,2,4], colors="r", linewidth = 10, zorder = 2)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py", line 195, in contour
    return self._contour("contour", *XYCL, **kwargs)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py", line 167, in _contour
    ny, nx = C.shape
ValueError: too many values to unpack</p>

<p>I've also tried to use the FITS_tools/match_images.py to resize the NVSS FITS first to the normal 2 axis size of the DSS file, and then to use the corrected file instead of the original one, but it only gives me an error:</p>

<p>Traceback (most recent call last):
  File "p.py", line 64, in 
    im1,im2 = FITS_tools.match_fits(to_be_projected,reference_fits)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/match_images.py", line 105, in match_fits
    image1_projected = project_to_header(fitsfile1, header, **kwargs)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/match_images.py", line 64, in project_to_header
    **kwargs)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", line 49, in hcongrid
    grid1 = get_pixel_mapping(header1, header2)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", line 128, in get_pixel_mapping
    csys2 = _ctype_to_csys(wcs2.wcs)
  File "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", line 169, in _ctype_to_csys
    raise NotImplementedError("Non-fk4/fk5 equinoxes are not allowed")
NotImplementedError: Non-fk4/fk5 equinoxes are not allowed</p>

<p>I have no idea what to do. There is no similar problem with the FIRST.FITS files. Imsize in Python also tells me the NVSS is the only one who is 4 dimensional for example (1, 1, 250, 250). So it could not be overlaied properley. Do you have ANY idea? Please help me and I can donate your projects in revange. I spent few days trying to solve it and it's still not working, but I need it desperately.</p>

<p>CODE</p>

<pre><code># Import matplotlib modules

import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid1.axes_divider import make_axes_locatable
from matplotlib.axes import Axes
import matplotlib.cm as cm
from matplotlib.patches import Ellipse
import linecache
import FITS_tools
# Import numpy and scipy for filtering
import scipy.ndimage as nd
import numpy as np
import pyfits
import matplotlib.pyplot
import pylab

#Import astronomical libraries
from astropy.io import fits
import astropy.units as u
#from astroquery.ned import Ned
import pywcsgrid2

# Read and prepare the data

file1=open('/home/ela/file')
count=len(open('file', 'rU').readlines())
print count

for i in xrange(count):
  wiersz=file1.readline()
  title=str(wiersz)
  print title
  title2=title.strip("\n")
  print title2

  path = '/home/ela/'
  fitstitle = path+title2+'_DSS.FITS'
  fitstitle2 = path+title2+'_FIRST.FITS'
  fitstitle3 = path+title2+'_NVSS.FITS'
  datafile = path+title2 
  outtitle = path+title2+'.png'
  print outtitle
  print datafile

  nvss = fits.open(fitstitle)[0]
  first = fits.open(fitstitle2)[0]
  opt = fits.open(fitstitle3)[0]
  try:
     fsock = fits.open(fitstitle3)      #(2)   
  except IOError:                     
     print "Plik nie istnieje"
  print "Ta linia zawsze zostanie wypisana"  #(3)

  opt.verify('fix')
  first.verify('fix')
  nvss.verify('fix')

  Header = nvss.header
  Header2 = first.header
  Header3 = opt.header

  to_be_projected = path+title2+'_NVSS.FITS'
  reference_fits  = path+title2+'_DSS.FITS'
  im1,im2 = FITS_tools.match_fits(to_be_projected,reference_fits)

  print(opt.shape)
  print(first.shape)
  print(nvss.shape)
  print(im2.shape)


#We select the range we want to plot
  minmax_image = [np.average(nvss.data)-6.*np.std(nvss.data), np.average(nvss.data)+5.*np.std(nvss.data)]   #Min and max value for the image
  minmax_PM = [-500., 500.]

# PREPARE PYWCSGRID2 AXES AND FIGURE
  params = {'text.usetex': True,'font.family': 'serif', 'font.serif': 'Times New Roman'}
  plt.rcParams.update(params)

#INITIALIZE FIGURE
  fig = plt.figure(1)
  ax = pywcsgrid2.subplot(111, header=Header)

#CREATE COLORBAR AXIS
  divider = make_axes_locatable(ax)
  cax = divider.new_horizontal("5%", pad=0.1, axes_class=Axes)
  #fig.add_axes(cax)

#Configure axis
  ax.grid()                                                                              #Will plot a grid in the figure
#  ax.set_ticklabel_type("arcmin", center_pixel=[Header['CRPIX1'],Header['CRPIX2']])      #Coordinates centered at the galaxy
  ax.set_ticklabel_type("arcmin")      #Coordinates centered at the galaxy
  ax.set_display_coord_system("fk5")
#  ax.add_compass(loc=3)                                                                  #Add a compass at the bottom left of the image

#Plot the u filter image
  i = ax.imshow(nvss.data, origin="lower", interpolation="nearest", cmap=cm.bone_r, vmin= minmax_image[0], vmax = minmax_image[1], zorder = 0)

#Plot contours from the infrared image

  cont = ax[Header2].contour(nd.gaussian_filter(first.data,4),2 , colors="r", linewidth = 20, zorder = 2)
#  cont = ax[Header2].contour(first.data, [-2,0,2], colors="r", linewidth = 20, zorder = 1)

#  cont2 = ax[Header2].contour(opt.data, [-8,-2,2,4], colors="r", linewidth = 10, zorder = 2)

#Plot PN positions with color coded velocities
# Velocities = ax['fk5'].scatter(Close_to_M31_PNs['RA(deg)'], Close_to_M31_PNs['DEC(deg)'], c = Close_to_M31_PNs['Velocity'], s = 30, cmap=cm.RdBu, edgecolor = 'none',
# vmin = minmax_PM[0], vmax = minmax_PM[1], zorder = 2)


  f2=open(datafile)
  count2=len(open('f2', 'rU').readlines())
  print count2

  for i in xrange(count):
    xx=f2.readline()
   # print xx
    yy=f2.readline()
    xxx=float(xx)
    print xxx
    yyy=float(yy)
    print yyy

    Velocities = ax['fk5'].scatter(xxx, yyy ,c=40,   s = 200, marker='x', edgecolor = 'red', vmin = minmax_PM[0], vmax = minmax_PM[1], zorder = 1)

  it2 = ax.add_inner_title(title2, loc=1)

#  Plot the colorbar, with the v_los of the PN
#  cbar = plt.colorbar(Velocities, cax=cax)
#  cbar.set_label(r'$v_{los}[$m s$^{-1}]$')
#  set_label('4444')
  plt.show()
  plt.savefig(outtitle)
#plt.savefig("image1.png")
</code></pre>

<p>To clarify for general use: this is a question about how to handle FITS files with degenerate axes, which are commonly produced by the <a href="https://casa.nrao.edu/" rel="nofollow">CASA data reduction program</a> and other radio data reduction tools; the degenerate axes are frequency/wavelength and stokes.  Some of the astropy affiliated tools know how to handle these (e.g., <a href="https://aplpy.github.io/" rel="nofollow">aplpy</a>), but many do not.</p>

<p>The simplest answer is to just use <code>squeeze</code> to drop the degenerate axes in the data.  However, if you want to preserve the metadata when doing this, there's a little more involved:</p>

<pre><code>from astropy.io import fits
from astropy import wcs

fh = fits.open('file.fits')
data = fh[0].data.squeeze() # drops the size-1 axes
header = fh[0].header
mywcs = wcs.WCS(header).celestial
new_header = mywcs.to_header()
new_fh = fits.PrimaryHDU(data=data, header=new_header)
new_fh.writeto('new_file.fits')
</code></pre>

<p>That approach will give you a file with a valid celestial (RA/Dec) header, but it will lose all of the other header information.</p>

<p>If you want to preserve the other header information, you can use the <a href="http://fits-tools.readthedocs.io/en/latest/fits_tools.html" rel="nofollow"><code>FITS_tools</code></a> tool <a href="https://github.com/keflavich/FITS_tools/blob/master/FITS_tools/strip_headers.py" rel="nofollow"><code>flatten_header</code></a> instead of using the WCS operations above:</p>

<pre><code>new_header = FITS_tools.strip_headers.flatten_header(header)
</code></pre>

python astronomy astropy fits pyfits