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

from astropy.io import fits
fitsdata = fits.open('./output-1.fits')
fitsdata.info()

"""
Explanation: Write and open fits file with Astropy
End of explanation
"""


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

"""
Explanation: Convert galactic coor to cartesian coor from Tom Rice's paper
The code does not work now, d0 here is not clear yet!
End of explanation
"""


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

# 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})

"""
Explanation: Extract vtk object from fits file
End of explanation
"""


