Skip to content
Snippets Groups Projects
Commit 402e6080 authored by Jae-Joon Lee's avatar Jae-Joon Lee
Browse files

[galactic_extinction] initial commit for lallerment2019

parent 349952d1
Branches master
No related merge requests found
This is small helper module for Galactic extinction data.
For now, it includes helper for 3d cube data from Lallement+2019.
# Install
It has not been fully packaged yet. So, you need to manually copy the python file (it is a single file for now). And point it with a data path.
The data can be downloaded from
[lallerment2019_cube_converted.fits](http://210.110.233.67:9005/various-data-jjl/lallermant2019/lallerment2019_cube_converted.fits)
from distutils.core import setup
setup(name='galactic_extinction',
version='0.1',
description='JJ\'s small python modules for astronomy',
author='Jae-Joon Lee',
author_email='lee.j.joon@gmail.com',
url='https://data.kasi.re.kr/gitlab/leejjoon/small-astronomy-modules',
packages=['galactic_extinction'],
package_dir={'galactic_extinction':'src'}
)
import numpy as np
import pandas as pd
import astropy.io.fits as pyfits
from scipy.interpolate import RegularGridInterpolator
from astropy.coordinates import SkyCoord, Galactic
def _load_data():
fn = "../STILISM_cube.fits.gz"
f = pyfits.open(fn)
d = f[1].data
nx = int((3000+3000)/5. + 1)
ny = int((3000+3000)/5. + 1)
nz = int((400 + 400)/5. + 1)
zmax = 400
X = d["X"].reshape((nz, ny, nx))
Y = d["Y"].reshape((nz, ny, nx))
Z = d["Z"].reshape((nz, ny, nx))
dA0_dD = d['dA0_dD'].reshape((nz, ny, nx))
return ((Z[:, 0, 0], Y[0, :, 0], X[0, 0, :]),
dA0_dD)
def _save_converted_cube():
(z, y, x), dA0_dD = _load_data()
hdu = pyfits.PrimaryHDU(dA0_dD)
hdu_x = pyfits.ImageHDU(x, name="X")
hdu_y = pyfits.ImageHDU(y, name="y")
hdu_z = pyfits.ImageHDU(z, name="z")
outname = "lallerment2019_cube_converted.fits"
pyfits.HDUList([hdu, hdu_x, hdu_y, hdu_z]).writeto(outname, overwrite=True)
def _load_converted_cube(converted_cube_name):
# outname = "lallerment2019_cube_converted.fits"
f = pyfits.open(converted_cube_name)
cube = f[0].data
x = f["X"].data
y = f["Y"].data
z = f["Z"].data
return (z, y, x), cube
class ExtinctionCube(object):
def __init__(self, converted_cube_name=None):
if converted_cube_name is None:
converted_cube_name = "lallerment2019_cube_converted.fits"
(z, y, x), cube = _load_converted_cube(converted_cube_name)
self.zmax = np.max(z)
self.knots = RegularGridInterpolator((z, y, x), cube)
def get_a0_lb(self, l, b, unit="deg"):
# c = SkyCoord(30, 1, unit='deg', frame=Galactic)
c = SkyCoord(l, b, unit=unit, frame=Galactic)
xyz = c.cartesian.xyz.value
r_pc = np.arange(3000)
xyz_lb = xyz * r_pc[:, np.newaxis]
zmsk = np.abs(xyz_lb[:, -1]) < self.zmax
a0 = self.knots(xyz_lb[zmsk, ::-1])
return r_pc[zmsk], a0
def get_A0_lb(self, l, b, unit="deg"):
r_pc, a0 = self.get_a0_lb(l, b, unit=unit)
A0 = np.add.accumulate(a0)
return r_pc, A0
if __name__ == "__main__":
l, b = 30, 1
ec = ExtinctionCube()
r_pc, A0 = ec.get_A0_lb(l, b)
print(r_pc, A0)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment