--Introducing python packages / modules often used for parsing astronomical data --For convenience of interest, I think that the handling of 3D (space: x, space: y, radial velocity: v) data is biased.

astropy

web : http://www.astropy.org/
document : http://docs.astropy.org/en/stable/
github : https://github.com/astropy/astropy

Astronomical data analysis package. Originally developed as a separate package, it was integrated around 2012. The impression is that the processing specific to astronomical analysis is generally covered, from the handling of data to the operation of unit systems and coordinate systems (the package is too large and I do not know everything).

Modules often used in astropy

astropy.io.fits

A module that reads and writes FITS files. Complete migration from pyfits (no further use of pyfits is recommended). FITS file is a file format with a high degree of freedom, created for the purpose of storing astronomical data. Normally, 2D or 3D image data is saved with coordinate information, but any data structure can be saved as a table. For details on FITS files, refer to "FITS User's Guide".

`ipython`


import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
hdu = fits[0]
header = hdu.header
data = hdu.data

#Check header information
print(repr(header))

naxis = header['NAXIS']   #Can be treated like a dictionary
naxis1 = header.get('NAXIS1')

#Check the data structure
type(data)
>>> numpy.ndarray

data.ndim
>>> 2

data.shape
>>> (480, 720)

astropy.wcs

WCS (World Coordinate System) is a specification that describes how the data recorded in FITS corresponds to the coordinates on the celestial sphere. For more information, see http://fits.gsfc.nasa.gov/fits_wcs.html. Based on information such as the data reference index (CRPIX), the physical coordinate value of the reference index (CRVAL), the physical coordinate correspondence value of the data grid (CDELT), and the projection method on the celestial sphere (CTYPE), the index value is calculated. Convert to physical coordinate values.

`ipython`


import astropy.wcs
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')

wcs = astropy.wcs.WCS(fits[0].header)
print(wcs)
>>> ...  #Information is output

wcs.wcs_world2pix([[210, -19]], 0)   # (l=210, b=-19)Get index of
>>> array([[ 359.00009815,  120.00221429]])

wcs.wcs_pix2world([[150, 60]], 0)   # (x=150, y=60)Get the physical coordinates of
>>> array([[ 213.48334194,  -20.0000389 ]])

astropy.coordinates A module that can convert and calculate the coordinates of the sky (Right ascension, declination coordinates, Milky Way galactic coordinates, etc.).

`ipython`


import numpy
import astropy.coordinates
import astropy.wcs
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
header = fits[0].header
data = fits[0].data

maxind = numpy.unravel_index(numpy.nanargmax(data), data.shape)
wcs = astropy.wcs.WCS(fits[0].header)
maxworld = wcs.wcs_pix2world([maxind[::-1]], 0) 
        # wcs_The order of the index axes required by pix2world is x, y(, z)
        # fits[0].The order of the axes of data is(z,) y, x
        #The two directions are opposite, so maxind[::-1]Inverted with
galx, galy = maxworld[0]
# >>> (208.99999963578261, -19.383371004831144)

#Created in galactic coordinates
coord = astropy.coordinates.SkyCoord(galx, galy, frame='galactic', unit='deg')
# >>> <SkyCoord (Galactic): (l, b) in deg
#         (208.99999964, -19.383371)>

#Check the constellations
coord.get_constellation()
>>> u'Orion'

#Right ascension declination(epoch=J2000)Output with
coord.fk5 
>>> <SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
        (83.81461251, -5.38035224)>

#Format and output RA and Declination
coord.fk5.to_string('hmsdms')
>>> u'05h35m15.507s -05d22m49.268s'

astropy.units A module that handles units used in astronomical fields.

`ipython`


import astropy.units

# 1: 5km/Let's calculate the time it takes for the molecular cloud core moving with s to move 5pc
v = 5 * astropy.units.kilometer / astropy.units.second
# >>> <Quantity 5.0 km / s>

d = 5 * astropy.units.pc
# >>> <Quantity 5.0 pc>

(d/v).to('megayear')
>>> <Quantity 0.9777922216731284 Myr>


# 2: 10^-6 Mo/Mass accretion rate of yr, kg/Try to sec
(1e-6 * astropy.units.solMass / astropy.units.year).to('kg/second')
>>> <Quantity 6.303077547088497e+16 kg / s>

APLpy

web : https://aplpy.github.io/
document : http://aplpy.readthedocs.org/en/stable/
github : https://github.com/aplpy/aplpy

FITS data visualization package. It can output quality images that can be used for publishing papers.

Example of use

`ipython`


import matplotlib.pyplot
import astropy.io.fits
import aplpy

#When downloading a FITS file from the web
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
fig = aplpy.FITSFigure(fits)

#If the FITS file is stored locally
# fig = aplpy.FITSFigure('Orion.CO1221.Osaka.beam204.mom0.fits.gz')

fig.show_colorscale()
fig.add_colorbar()
# fig.save('Orion.CO1221.Osaka.beam204.mom0.png')  #When saving an image

matplotlib.pyplot.show()

numpy

web : http://www.numpy.org/

A module that provides the data structure and operations that form the basis for numerical calculations.

Functions often used in numpy

numpy.sum, .nansum Calculate the total value of numpy.ndarray or list. By specifying the axis argument, you can specify the axis to be integrated.

`ipython`


import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')

#Unobserved area(voxel)Often contains nan, so use nansum
numpy.nansum(fits[0].data)
>>> 1630958.2    #The unit is K km/s pix^2


#3D data(cube)From 2D data(Integral strength)Also used when issuing
# !!!!Download capacity note(1.8GB) !!!!
# fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.cube.fits.gz')
fits = astropy.io.fits.open('Orion.CO1221.Osaka.beam204.cube.fits.gz')

fits[0].data.shape
>>> (2521, 480, 720)    #z axis(speed),y axis(side),x axis(Vertical) 

mom0 = numpy.nansum(fits[0].data, axis=0)  #Sum the 0 axes
mom0.shape
>>> (480, 720)

numpy.min, .max, .nanmin, .nanmax Get the minimum and maximum of numpy.ndarray or list.

`ipython`


import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')

numpy.nanmin(fits[0].data)
>>> -0.61183316  # K km/s

numpy.nanmax(fits[0].data)
>>> 431.83881  # K km/s

numpy.argmin, .argmax, .nanargmin, .nanargmax Get the index that holds the minimum and maximum values of numpy.ndarray or list.

`ipython`


import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')
data = fits[0].data

maxind = numpy.nanargmax(data)
# >>> 70259

maxcoord = numpy.unravel_index(maxind, data.shape)
# >>> (97, 419)

data[maxcoord] == numpy.nanmax(data)
>>> True

numpy.meshgrid Create a two-dimensional array using two one-dimensional arrays. Use when you want a map of coordinate information.

`ipython`


import numpy
import astropy.io.fits
fits = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.mom0.fits.gz')

def get_axis(header, axis=1):
    crval = header.get('CRVAL%1d'%(axis))
    crpix = header.get('CRPIX%1d'%(axis))
    cdelt = header.get('CDELT%1d'%(axis))
    num = header.get('NAXIS%1d'%(axis))
    return crval + (numpy.arange(num) - crpix + 1) * cdelt

x1d = get_axis(fits[0].header, 1)
y1d = get_axis(fits[0].header, 2)

X2d, Y2d = numpy.meshgrid(x1d, y1d)
X2d.shape
>>> (480, 720)


#For 3d data(Since meshgrid cannot be used, do as follows)
# !!!!Download capacity note(1.8GB) !!!!
# fits_cube = astropy.io.fits.open('http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/data/Orion.CO1221.Osaka.beam204.cube.fits.gz')
fits_cube = astropy.io.fits.open('Orion.CO1221.Osaka.beam204.cube.fits.gz')

x1d_ = get_axis(fits_cube[0].header, 1)
y1d_ = get_axis(fits_cube[0].header, 2)
v1d_ = get_axis(fits_cube[0].header, 3)

X3d = x1d_[None, None, :]*1 + y1d_[None, :, None]*0 + v1d_[:, None, None]*0
Y3d = x1d_[None, None, :]*0 + y1d_[None, :, None]*1 + v1d_[:, None, None]*0
V3d = x1d_[None, None, :]*0 + y1d_[None, :, None]*0 + v1d_[:, None, None]*1
X3d.shape
>>> (2521, 480, 720)
#Numpy instead of None.newaxis can also be used(Is that recommended?)

numpy.polyfit, .polyval

A polynomial fitting and a function that restores a curve from fitting parameters.

scipy

Functions often used in scipy

scipy.signal.convolve scipy.ndimage.map_coordinates scipy.interpolate.interp1d

matplotlib ... (planned to be added) ...

About data

http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/

Astro: Python modules / functions often used for analysis

Modules often used in astropy

ipython

ipython

ipython

ipython

Example of use

ipython

Functions often used in numpy

ipython

ipython

ipython

ipython

Functions often used in scipy

About data

`ipython`

`ipython`

`ipython`

`ipython`

`ipython`

`ipython`

`ipython`

`ipython`

`ipython`