2016-09-10 65 views
0

這個問題可能主要是爲了或多或少地推動天文學家的進步。在Python中調整NVSS FITS格式文件的大小並在這個文件中進行運算

你知道如何將NVSS適合文件轉換爲只有2(不是4!)軸?或者當我試圖在光學DSS數據上覆蓋nvss計數時,如何使用Python的astropy和其他「astro」庫來處理具有4軸的文件,並在Python中產生以下錯誤? (下面的代碼)

我試圖做到這一點,當有帖4軸爲NVSS FITS,有錯誤和警告信息:

警告: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」,linewidth = 10,zorder = 2) 文件「/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py」,第195行,在輪廓 return self._contour(「contour」,* XYCL,** kwargs) 文件「/home/ela/anaconda2/lib/python2.7/site-packages/mpl_toolkits/axes_grid1/parasite_axes.py」,第167行在_contour NY中,nx = C.shape ValueError異常:值過多解壓

我還試圖使用FITS_tools/match_images.py來調整NVSS第一裝配到正常的2軸SIZ然後使用正確的文件,而不是原來的,但它只會給我一個錯誤:

回溯(最近呼叫最後): 文件「p.py」,第64行, in 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) get_pixel_mapping中的文件「/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py」,第128行 csys2 = _ctype_to_csys(wcs2.wcs) 文件「/home/ela/anaconda2/lib/python2.7/site-packages/FITS_tools/hcongrid.py」,第169行,在_ctype_to_csys中 raise NotImplementedError(「不允許非fk4/fk5分點」) NotImplementedError :非fk4/fk5分點不允許

我不知道該怎麼做。 FIRST.FITS文件沒有類似的問題。 Python中的Imsize也告訴我NVSS是唯一一個例如4維的(1,1,250,250)。所以它不能疊加在properley上。你有什麼主意嗎?請幫助我,我可以捐助你的項目。我花了幾天的時間試圖解決它,但它仍然無法正常工作,但我絕望地需要它。

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

這是誣陷,不良Q.請查閱:-http:// stackoverflow.com/help/how-to-ask – Dravidian

+0

我並不確定你想要做什麼,但是如果你想要重複繪製圖片,請查看[reproject](http:// reproject)。 readthedocs.io/en/stable/)。至於WCS,你可以實例化一個['astropy.wcs.WCS'](http://docs.astropy.org/en/stable/api/astropy.wcs.WCS.html#astropy.wcs.WCS)對象通過使用'naxis'關鍵字來指示哪些(2)軸使用;或者使用'WCS(header).celestial'。 – Evert

回答

2

爲了闡明一般用途:這是一個有關如何處理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_tools工具flatten_header代替以上WCS操作:

new_header = FITS_tools.strip_headers.flatten_header(header) 
相關問題