--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
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).
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
FITS data visualization package. It can output quality images that can be used for publishing papers.
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
A module that provides the data structure and operations that form the basis for numerical calculations.
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
scipy.signal.convolve scipy.ndimage.map_coordinates scipy.interpolate.interp1d
matplotlib ... (planned to be added) ...
http://www.astro.s.osakafu-u.ac.jp/~nishimura/Orion/
Recommended Posts