0

这个问题可能主要是针对或多或少进步的天文学家。

你知道如何将 NVSS 拟合文件转换为只有 2 个（不是 4 个！）轴的拟合吗？或者，当我尝试在光学 DSS 数据上重叠 nvss 计数时，如何处理具有 4 轴并在 Python 中生成以下错误的文件，使用 astropy 和其他用于 Python 的“astro”库？（下面的代码）

我试过这样做，当 NVSS FITS 有 4 个轴时，会出现错误消息和警告：

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

我还尝试使用 FITS_tools/match_images.py 首先将 NVSS FITS 调整为 DSS 文件的正常 2 轴大小，然后使用更正后的文件而不是原始文件，但这只会给我一个错误：

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

我不知道该怎么做。FIRST.FITS 文件没有类似的问题。Python 中的 Imsize 还告诉我 NVSS 是唯一一个是 4 维的，例如（1、1、250、250）。所以它不能被适当地覆盖。你有什么主意吗？请帮助我，我可以捐赠你的项目作为报复。我花了几天时间试图解决它，但它仍然无法正常工作，但我非常需要它。

代码

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

1 回答 1

2

澄清一下一般用途：这是一个关于如何处理退化轴的 FITS 文件的问题，这些文件通常由CASA 数据缩减程序和其他无线电数据缩减工具生成；简并轴是频率/波长和斯托克斯。一些 astropy 附属工具知道如何处理这些（例如，aplpy），但许多不知道。

最简单的答案是仅用于squeeze删除数据中的退化轴。但是，如果您想在执行此操作时保留元数据，则涉及更多内容：

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

该方法将为您提供具有有效天体 (RA/Dec) 标头的文件，但它会丢失所有其他标头信息。

如果要保留其他头信息，可以使用FITS_tools工具flatten_header代替上面的WCS操作：

new_header = FITS_tools.strip_headers.flatten_header(header)
于 2016-09-11T15:10:57.337 回答