Изменение размера файла формата NVSS FITS в Python и работа с данным файлом в astropy

Question

Эти вопросы, вероятно, в основном для более или менее продвигают астрономов.

Знаете ли вы, как преобразовать файл NVSS в соответствие с только 2 (НЕ 4!) Осями? Или как обращаться с файлом, который имеет 4 оси и генерирует следующие ошибки в Python, когда я пытаюсь перекрывать nvss-счетчики на оптических данных DSS, используя астропию и другие «астрологические» библиотеки для Python? (Код ниже)

Я пытался сделать это, и когда есть Фессалоникийцам 4 осей для NVSS припадки, есть сообщения об ошибках и предупреждения:

ВНИМАНИЕ: FITSFixedWarning: Преобразование WCS имеет более осей (4) чем изображение, связанное с (2) [astropy.wcs.wcs] ПРЕДУПРЕЖДЕНИЕ: FITSFixedWarning: 'datfix' произнес изменение 'Недопустимое значение параметра: неверная дата' 19970331 ''. [Astropy.wcs.wcs] re-sizing a fits image in python

ВНИМАНИЕ: FITSFixedWarning: 'datfix' внес изменения 'недопустимое значение параметра: неверная дата' 19970331'. [astropy.wcs.wcs] Traceback (последний последний звонок): Файл «p.py», строка 118, в cont2 = ax [Header2] .contour (opt.data, [-8, -2,2 , 4], colors = "r", ширина линии = 10, zorder = 2) Файл "/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py", строка 195, в контуре return self._contour («контур», * XYCL, ** kwargs) Файл «/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py», строка 167 в _contour пу, пх = C.shape ValueError: слишком много значений для распаковки

Я также попытался использовать FITS_tools/match_images.py, чтобы изменить размер NVSS FITS первым к нормальной 2 оси СИЗ e файла DSS, а затем использовать исправленный файл вместо исходного, но это только дает мне ошибку:

Traceback (последний последний звонок): Файл «p.py», строка 64, в IM1, IM2 = FITS_tools.match_fits (to_be_projected, reference_fits) Файл "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/match_images.py", строка 105, в match_fits image1_projected = project_to_header (fitsfile1, header, ** kwargs) Файл «/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/match_images.py», строка 64, в project_to_header ** kwargs) Файл "/ home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py ", строка 49, в hcongrid grid1 = get_pixel_mapping (header1, header2) Файл «/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py», строка 128, в get_pixel_mapping csys2 = _ctype_to_csys (wcs2.wcs) Файл "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", строка 169, в _ctype_to_csys повышение NotImplementedError ("Non-FK4/FK5 равноденствия не допускаются") NotImplementedError : Неденежные знаки не-fk4/fk5 не допускаются

Я понятия не имею, что делать. Подобной проблемы с файлами FIRST.FITS нет. Imsize в Python также говорит мне, что 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") 

Answer 1

Для уточнения общего пользования: это вопрос о том, как обрабатывать FITS файлов с вырожденными осями, которые обычно получают с помощью CASA data reduction program и других средств обработки информации радио; вырожденные оси - это частота/длина волны и сток. Некоторые из афроамериканских инструментов знают, как их обрабатывать (например, 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/дек), но он потеряет все другая информация заголовка.

Если вы хотите сохранить другую информацию заголовка, вы можете использовать инструмент FITS_toolsflatten_header вместо использования операций WCS выше:

new_header = FITS_tools.strip_headers.flatten_header(header)