Reputation: 19
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")
Upvotes: 0
Views: 1151
Reputation: 19215
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)
Upvotes: 2