0

我一直在寻找解决我的问题的方法，但找不到。我有一个 FITS 数据立方体，我需要通过 PyFITS 对其进行裁剪。当我用我的脚本来做这件事时，最后我会得到一个 2-D FITS 图像！第一个维度是能量，第二个和第三个维度分别是经度和纬度。

我的脚本如下：

#!/usr/bin/env python
import pyfits
import os
import sys


def CropFitsFile( src, dst, xs, xe, ys, ye):
    fh = pyfits.open(src)
    for eng in range(0,2):
        img = fh[0].data[eng,ys:ye,xs:xe]
        header = fh[0].header
        newfh=pyfits.PrimaryHDU(data=img,header=header)
        if os.path.exists(dst):
            os.remove(dst)
        newfh.writeto(dst)


if __name__ == "__main__":
    CropFitsFile(
        src=sys.argv[1],
        dst=sys.argv[2],
        xs=int(sys.argv[3]),
        xe=int(sys.argv[4]),
        ys=int(sys.argv[5]),
        ye=int(sys.argv[6])
        )
4

2 回答 2

0

如果我理解正确，您想要切片 3D 数组但保留第三维（即使它只是大小 1）。

这是一个关于 Numpy 数组的问题。当您有一个 N 维 numpy 数组时，传递一个维度的标量索引会返回一个维度为 N-1 的数组，沿您索引的轴切片。例如：

>>> arr = np.arange(27).reshape(3, 3, 3)
>>> arr
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]],

       [[18, 19, 20],
        [21, 22, 23],
        [24, 25, 26]]])
>>> arr[0]
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> arr[1]
array([[ 9, 10, 11],
       [12, 13, 14],
       [15, 16, 17]])

您还可以沿不同的轴切片，例如：

>>> arr[:,0,:]
array([[ 0,  1,  2],
       [ 9, 10, 11],
       [18, 19, 20]])

如果出于某种原因想要返回 N 维数组而不是 N-1 维数组，最简单的方法是显式请求大小为 1 的切片，而不是使用标量索引。例如：

>>> arr[0:1]
array([[[0, 1, 2],
        [3, 4, 5],
        [6, 7, 8]]])

我将给出与其他类似问题相同的建议：除了您的数据来自 FITS 文件这一事实之外，这并不是关于 PyFITS 的真正问题。PyFITS 与大多数科学 Python 库一样，将数据作为 numpy 数组返回。这些是大多数科学 Python 应用程序中用于数值数据的主要数据结构，因此无论好坏，学习一些 numpy 基础知识都是在 Python 中进行数据分析的先决条件。如果您曾经使用过 MATLAB，那么 NumPy 数组类似于 MATLAB 中的数组。你可以从我的简短教程开始，但还有其他的（可能还有更好的 :) github.com/embray/notebooks/blob/master/numpy.ipynb

于 2016-04-03T09:55:59.370 回答
-1 
from astropy.io import fits

Ccube = fits.open('Cha_binned_ccube.fits', mode='update')
Ccube.info()
Ccube[0].shape
Ccube[0].data = Ccube[0].data[0:3,0:181,0:402]
Ccube[0].writeto('Cha_binned_ccube_resize.fits')
于 2016-04-04T10:59:17.740 回答