0

我正在尝试根据从Header. 但是，以我对 Python 的了解和astropy. 我的代码是这样的：

from astropy.io import fits
import numpy as np

Wise1 = fits.open('Image1.fits')
im1 = Wise1[0].data

im1 = np.where(im1 > *latitude1, 0, im1)

newhdu = fits.PrimaryHDU(im1)
newhdulist = fits.HDUList([newhdu])
newhdulist.writeto('1b1_Bg_Removed_2.fits')

这latitude1将是一个以度为单位的值，在从标题中调用后被识别。所以我需要完成两件事：

  1. 如何调用标题来识别银河纬度？
  2. 以这样的方式拼接数组，使其仅包含纬度范围的值，其他所有值都为 0。
4

1 回答 1

3

根据您展示的示例，我认为“拼接”的意思是“剪裁”或“裁剪”。

astropy.nddata有一个基于世界坐标系（即 lat/lon 或 ra/dec）剪切的例程

但是，在您处理的简单情况下，您只需要每个像素的坐标。通过制作 WCS 来做到这一点：

from astropy import wcs
w = wcs.WCS(Wise1[0].header)
xx,yy = np.indices(im.shape)
lon,lat = w.wcs_pix2world(xx,yy,0)

newim = im[lat > my_lowest_latitude]

但是，如果您想保留标题信息，最好使用剪切工具，因为您不必手动管理它。

from astropy.nddata import Cutout2D
from astropy import coordinates
from astropy import units as u

# example coordinate - you'll have to figure one out that's in your map
center = coordinates.SkyCoord(mylon*u.deg, mylat*u.deg, frame='fk5')

# then make an array cutout
co = nddata.Cutout2D(im, center, size=[0.1,0.2]*u.arcmin, wcs=w)

# create a new FITS HDU
hdu = fits.PrimaryHDU(data=co.data, header=co.wcs.to_header())

# write to disk
hdu.writeto('cropped_file.fits')

astropy 文档中有一个示例用例。

于 2015-11-04T07:17:22.410 回答