Write and open fits file with Astropy


In [8]:
import numpy as np
from astropy.io import fits 
data = np.ones((10, 20, 30, 5))

# hdul = fits.HDUList()
# hdul.append(fits.PrimaryHDU())

# for img in data:
#     hdul.append(fits.ImageHDU(data=img))

hdu = fits.PrimaryHDU(data)
hdulist = fits.HDUList([hdu])

hdulist.writeto('output-1.fits')

In [9]:
from astropy.io import fits
fitsdata = fits.open('./output-1.fits')
fitsdata.info()


Filename: ./output-1.fits
No.    Name         Type      Cards   Dimensions   Format
0    PRIMARY     PrimaryHDU       8   (5, 30, 20, 10)   float64   

Convert galactic coor to cartesian coor from Tom Rice's paper

The code does not work now, d0 here is not clear yet!


In [ ]:
# read table data
from astropy.io import fits
import numpy as np
from math import pi, cos, sin

hdulist = fits.open('BGPS_3D_catalogue_HCO+.fit')
table_data = hdulist[1].data
# table_data.columns to radian
b = table_data.field('_Glat')*pi/180.
l = table_data.field('_Glon')*pi/180.


# use algorithm in Tom Rice catalogue paper to convert galatic coor to cartesian coor
x = d0*cos(l)*cos(b)
y = d0*sin(l)*cos(b)
z = d0*sin(b)

Extract vtk object from fits file


In [1]:
import vtk
from vtk.util import numpy_support
from mayavi import mlab
from astropy.io import fits

import numpy as np

hdulist = fits.open('l1448_13co.fits')
NumPy_data = hdulist[0].data

NumPy_data_shape = NumPy_data.shape
VTK_data = numpy_support.numpy_to_vtk(num_array=NumPy_data.ravel(), deep=True, array_type=vtk.VTK_FLOAT)
scalarfield = mlab.pipeline.scalar_field(NumPy_data)
print('done', scalarfield)
%tb
# export to file not found


An exception has occurred, use %tb to see the full traceback.

SystemExit: This program needs access to the screen.
Please run with a Framework build of python, and only when you are
logged in on the main display of your Mac.
/Users/penny/anaconda/lib/python2.7/site-packages/IPython/core/interactiveshell.py:2889: UserWarning: To exit: use 'exit', 'quit', or Ctrl-D.
  warn("To exit: use 'exit', 'quit', or Ctrl-D.", stacklevel=1)

In [ ]:
# http://www.vtk.org/Wiki/VTK/Writing_VTK_files_using_python

from pyevtk.hl import gridToVTK

from astropy.io import fits
import numpy as np

hdulist = fits.open('l1448_13co.fits')
NumPy_data = hdulist[0].data

noSlices = 5
juliaStacked = numpy.dstack([julia]*noSlices)

x = numpy.arange(0, w+1)
y = numpy.arange(0, h+1)
z = numpy.arange(0, noSlices+1)

gridToVTK("./julia", x, y, z, cellData = {'julia': juliaStacked})