0

I have looked to find a solution for my issue, but couldn't. I have a FITS data cube, and I need to crop it by PyFITS. When I do it by my script, finally I'll have a 2-D FITS image! The first dimension is energy, and the second and thirds are longitude and latitude, respectively.

My script is as below:

#!/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])
        )
Dom
  • 1,687
  • 6
  • 27
  • 37
Soheila Abd
  • 11
  • 1
  • What exactly is the problem? This looks fine, more or less. A few tiny notes are that `pyfits` is deprecated and you should use `astropy.io.fits` instead (it won't require any code changes except for the import). And the `wrtieto` method already has the functionality built-in to overwrite an existing file (the ` clobber` option). You don't need extra code for that. But again, did you have a question? – Iguananaut Apr 02 '16 at 17:33
  • Hi Iguananaut, Thank you for your reply. The problem is that I need to crop a 3-D data cube but when I run the code, the results is a 2-D fits image! I don't know how can I fix it. I replaced pyrites by astropy.io.fits, but still the same problem. Thank you again for your help – Soheila Abd Apr 03 '16 at 02:47
  • So you still want it to have 3 dimensions, but the first dimension of size 1? – Iguananaut Apr 03 '16 at 09:35

2 Answers2

0

If I understand correctly you want to slice a 3D array but keep the third dimensions (even if it's just size 1).

This is a question about Numpy arrays. When you have an N-dimensional numpy array, passing a scalar index for one dimension returns an array of dimension N-1, sliced along the axis you indexed. For example:

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

You can also slice along a different axis like:

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

If you want, for whatever reason, to return an N-dimensional array instead of an N-1 dimensional array, the easiest way is to explicitly request a slice of size 1, instead of using scalar index. For example:

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

I'll give the same advice I've given on other questions like this: this isn't really a question about PyFITS beyond the fact that your data came out of a FITS file. PyFITS, like most scientific Python libraries, returns data as numpy arrays. These are the main data structure used for numeric data in most scientific Python applications, so learning some basics of numpy are sort of a prerequisite, for better or worse, to doing data analysis in Python. If you've ever used MATLAB, NumPy arrays are similar to arrays in MATLAB. You can start with my short tutorial, but there are others (and probably better ones too :) github.com/embray/notebooks/blob/master/numpy.ipynb

Iguananaut
  • 21,810
  • 5
  • 50
  • 63
  • Thank you very much for your detailed and useful answer. That link helped me a lot to find the answer :) yes I needed to know more about manupulating header of FITS file. It seems there are many ways to resize the 2-D or 3-D FITS files such as np.resize. I changed the code as below and it worked ;) – Soheila Abd Apr 04 '16 at 10:57
  • Regarding the header please heed the advice that it is not necessary to manually edit header keywords concerning image dimensions. If this answer solved your problem please accept the solution. – Iguananaut Apr 04 '16 at 21:00
0
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')
Soheila Abd
  • 11
  • 1
  • 1
    While this code may answer the question, providing additional context regarding _why_ and/or _how_ it answers the question would significantly improve its long-term value. Please [edit] your answer to add some explanation. – Toby Speight Apr 04 '16 at 12:09