Resizing NVSS FITS format file in Python and operating on this file in astropy

Question

this questions is probably mainly for more or less advances astronomers.

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)

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

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]
https://stackoverflow.com/questions/33107224/re-sizing-a-fits-image-in-python

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 <module>
    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

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:

Traceback (most recent call last):
  File "p.py", line 64, in <module>
    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

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.

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")

This is poorly framed Q. Please look up :-http://stackoverflow.com/help/how-to-ask — Dravidian, Sep 10 '16 at 17:41
I'm not exactly sure what you're trying to do, but if you're trying to overplot images, have a look at [reproject](http://reproject.readthedocs.io/en/stable/). As for the WCS, you can instantiate an [`astropy.wcs.WCS`](http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS) object with limited axes by using the `naxis` keyword to indicate which (2) axes to use; or use `WCS(header).celestial`. — , Sep 10 '16 at 22:55

score 2 · Answer 1 · answered Sep 11 '16 at 15:10

To clarify for general use: this is a question about how to handle FITS files with degenerate axes, which are commonly produced by the CASA data reduction program 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., aplpy), but many do not.

The simplest answer is to just use squeeze 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:

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')

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

If you want to preserve the other header information, you can use the FITS_tools tool flatten_header instead of using the WCS operations above:

new_header = FITS_tools.strip_headers.flatten_header(header)

Resizing NVSS FITS format file in Python and operating on this file in astropy

1 Answers1