Redimensionner le fichier de format NVSS FITS en Python et opérer sur ce fichier dans astropy

Question

cette question est probablement principalement pour les astronomes avance plus ou moins. Savez-vous comment transformer le fichier NVSS pour qu'il s'adapte à 2 axes (NON 4!)? Ou comment traiter le fichier qui a 4 axes et génère les erreurs suivantes en Python lorsque j'essaye de superposer des comptes nvss sur des données optiques DSS, en utilisant astropy et d'autres bibliothèques "astro" pour Python? (Code ci-dessous)

J'ai essayé de le faire et quand il y a Thess 4 axes pour NVSS FITS, il y a des messages d'erreur et les avertissements:

AVERTISSEMENT: FITSFixedWarning: La transformation WCS a plusieurs axes (4) que l'image à laquelle il est associé (2) [astropy.wcs.wcs] AVERTISSEMENT: FITSFixedWarning: 'datfix' a fait la modification 'Valeur de paramètre non valide: date invalide' 19970331 ''. [Astropy.wcs.wcs] re-sizing a fits image in python

AVERTISSEMENT: FITSFixedWarning: 'datfix' a fait le changement 'valeur de paramètre non valide: date non valide '19970331''. [astropy.wcs.wcs] Traceback (dernier appel en dernier): Fichier "p.py", ligne 118, en cont2 = ax [En-tête2] .contour (opt.data, [-8, -2,2 , 4], colors = "r", linewidth = 10, zorder = 2) Fichier "/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py", ligne 195, dans le contour return self._contour ("contour", * XYCL, ** kwargs) Fichier "/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py", ligne 167 , dans _contour ny, nx = C.shape ValueError: trop de valeurs pour déballer

J'ai aussi essayé d'utiliser les FITS_tools/match_images.py pour redimensionner la NVSS FITS premier à l'axe 2 normale SIZ e du fichier DSS, puis d'utiliser le fichier corrigé au lieu de l'original, mais il me donne seulement une erreur:

Traceback (dernier appel dernier): Fichier "p.py", ligne 64, dans 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", la ligne 105, en match_fits image1_projected = project_to_header (fitsfile1, header, ** kwargs) Fichier "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/match_images.py", ligne 64, dans project_to_header ** kwargs) Fichier "/ home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py ", ligne 49, dans hcongrid grid1 = get_pixel_mapping (en-tête1, en-tête2) Fichier "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", ligne 128, dans get_pixel_mapping csys2 = _ctype_to_csys (wcs2.wcs) Fichier "/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py", ligne 169, dans _ctype_to_csys raise NotImplementedError ("Les équinoxes non-fk4/fk5 ne sont pas autorisées") NotImplementedError : Les équinoxes non-fk4/fk5 ne sont pas autorisés

Je n'ai aucune idée de ce qu'il faut faire. Il n'y a pas de problème similaire avec les fichiers FIRST.FITS. Imsize en Python me dit aussi que le NVSS est le seul qui soit en 4 dimensions par exemple (1, 1, 250, 250). Donc, il ne pouvait pas être recouvert correctement. Avez-vous une idée? Aidez-moi s'il vous plaît et je peux faire un don de vos projets dans Revange. J'ai passé quelques jours à essayer de le résoudre et ça ne marche toujours pas, mais j'en ai désespérément besoin.

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

Answer 1

précisiez usage général: il est une question sur la façon de traiter les fichiers FITS avec des axes dégénérés, qui sont généralement produites par les CASA data reduction program et d'autres outils de réduction des données radio; les axes dégénérés sont fréquence/longueur d'onde et stokes. Certains des outils affiliés astropy savent comment les gérer (par exemple, aplpy), mais beaucoup ne le font pas.

La réponse la plus simple consiste simplement à utiliser squeeze pour supprimer les axes dégénérés dans les données. Toutefois, si vous souhaitez conserver les métadonnées lors de cette opération, il y a un peu plus impliqué:

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

Cette approche vous donnera un fichier avec un en-tête céleste valide (RA/Dec), mais il perd tout l'autre information d'en-tête.

Si vous souhaitez conserver l'autre information d'en-tête, vous pouvez utiliser l'outil FITS_toolsflatten_header au lieu d'utiliser les opérations ci-dessus WCS:

new_header = FITS_tools.strip_headers.flatten_header(header)