Astro: Python modules / functions often used for analysis

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

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

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

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/

Recommended Posts

Astro: Python modules / functions often used for analysis
XPath Basics (3) -Functions often used for XPath
Python for Data Analysis Chapter 4
Python for Data Analysis Chapter 2
Keyword arguments for Python functions
Python for Data Analysis Chapter 3
[Updated from time to time] Python memos often used for data analysis [N division, etc.]
What is Python? What is it used for?
Azure Functions: Try Durable Functions for Python
Import modules that are often used when starting the python interpreter
Installation summary often used for AI projects
Python visualization tool for data analysis work
Python functions
Modules of frequently used functions in Python (such as reading external files)
++ and-cannot be used for increment / decrement in python
Functions that can be used in for statements
[CovsirPhy] COVID-19 Python Package for Data Analysis: Data loading
Julia Quick Note [22] Calling Python functions and Python modules
Techniques often used in python short coding (Notepad)
[Python] f strings should be used for embedding strings
Summary of frequently used Python arrays (for myself)
Code often used in Python / Django apps [prefectures]
2016-10-30 else for Python3> for:
python [for myself]
Data analysis python
#Python basics (functions)
[Beginner] Python functions
Python Easy-to-use functions
Python basics: functions
Data analysis for improving POG 1 ~ Web scraping with Python ~
Automatically generate a polarity dictionary used for sentiment analysis
[For beginners] How to study Python3 data analysis exam
A collection of code often used in personal Python
Python 3.4 Create Windows7-64bit environment (for financial time series analysis)
3. Natural language processing with Python 4-1. Analysis for words with KWIC
[Python] Privately created & used small-scale functions (file operations, etc.)
List of Python code used in big data analysis
[CovsirPhy] COVID-19 Python package for data analysis: SIR-F model
[CovsirPhy] COVID-19 Python package for data analysis: S-R trend analysis
[CovsirPhy] COVID-19 Python Package for Data Analysis: SIR model
[CovsirPhy] COVID-19 Python Package for Data Analysis: Parameter estimation
A collection of Excel operations often used in Python
[Python for Hikari-] Chapter 06-03 Functions (arguments and return value 2)
Settings for testing C ++ 11 Python modules on Travis CI
Python3> slice copy / slice notation> used in for statements, etc.