In [1]:

    

from dustmaps.config import config
config['data_dir'] = '/Users/penny/Works/bayestar_dustmap'

import dustmaps.bayestar
dustmaps.bayestar.fetch()


    

    




File exists. Not overwriting.

In [2]:

    

from __future__ import print_function
from astropy.coordinates import SkyCoord
from dustmaps.bayestar import BayestarQuery
import numpy as np
import astropy.units as units

bayestar = BayestarQuery(max_samples=1)

# get cube in specific region
def get_cube(l0, b0, delta=5., step=0.05):
    # set up a grid of coordinates to plot, centered on the Aquila South cloud
    l = np.arange(l0 - delta, l0 + delta, step)
    b = np.arange(b0 - delta, b0 + delta, step)
    l, b = np.meshgrid(l, b)
    coords = SkyCoord(l*units.deg, b*units.deg, frame='galactic')
    
    # TODO: why in matplotlib the value is squared?
#     np.sqrt(Av)[::,::-1]

    # used the coefficient from Table 6 of Schlafly & Finkbeiner (2011) to convert SFD and Bayestar reddenings to magnitudes of AV.
#     av_cube = 2.742*bayestar(coords, mode='median')
    cube = bayestar(coords, mode='median')
    
    print(cube.shape)
    return np.rot90(cube)


    

In [3]:

    

# convert cube from accumutative 
def subtract_cube(np_array):
    d = np_array.shape[2]
    sort_array = np.zeros(np_array.shape)
    print(list(reversed(range(d))))
    for i in list(reversed(range(d))):
        if i>0:
            sort_array[:, :, i] = np_array[:, :, i] - np_array[:, :, i-1]
        else:
            sort_array[:, :, i] = np_array[:, :, i]
    return sort_array


    

In [4]:

    

# resample cube to match spatial coordinate: coorelate distance with slice z index, add empty slices between
# not currently necessary
def resample_cube(cube):
    distances = np.array(bayestar.distances)
    unit = bayestar.distances.unit
    d_shape = int(distances[-1])
    sam_ebv_cube = np.zeros((cube.shape[0], cube.shape[1], d_shape+1))
    print(sam_ebv_cube.shape, cube.shape[2])
    for item in range(cube.shape[2]):
        sam_ebv_cube[:, :, int(distances[item])] += cube[:, :, item]

    # cut distance < 1kpc part
    sam_ebv_cube[:, :, 0] = 0
    print(sam_ebv_cube.shape, np.nanmax(sam_ebv_cube), np.nanmin(sam_ebv_cube))
    print(np.nanmax(sam_ebv_cube[:, :, 63]))
    return sam_ebv_cube


    

In [5]:

    

# centered on the Aquila South cloud distance < 1.25kpc data 
AS_cloud = get_cube(37., -17, 5., 0.05)
sub_AS_cube = subtract_cube(AS_cloud)

dist_mask=np.where(bayestar.distances.value*1000.0 <= 1250.0) 
# try cut out from 1kpc, not solve the 'straight line' problem

sub_AS_cube=sub_AS_cube[:,:,dist_mask[0]]
sub_AS_cube.shape


    

    




(200, 200, 31)
[30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]






    
Out[5]:





(200, 200, 13)

In [6]:

    

from IPython.core.interactiveshell import InteractiveShell
InteractiveShell.ast_node_interactivity = "all"

import ipyvolume
from astropy.io import fits

# display AS cloud
ipyvolume.quickvolshow(sub_AS_cube.T, stero=True, data_min=None, level=[0.28, 0.19, 0.19], lighting=True, ambient=0.67, width=512, height=256, stereo=True, opacity=0.06)


    

    




/Users/penny/anaconda/lib/python2.7/site-packages/ipyvolume/serialize.py:27: RuntimeWarning: invalid value encountered in divide
  gradient = gradient / np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)






    






 
 

In [7]:

    

# check subtract: no output means pass

for i in range(sub_AS_cube.shape[2]-1):
    if not (AS_cloud[0, 0, i+1] - AS_cloud[0, 0, i])==sub_AS_cube[0, 0, i+1]:
        print('subtract wront on %i'%i)


    

In [8]:

    

# for i in range(0,sub_cube.shape[2]):
#     plt.imshow(sub_cube[:,:,i])
#     plt.show()


    

In [8]:

    

gamma_ring = get_cube(195, -11, step=0.05)
sub_gamma_ring = subtract_cube(gamma_ring)


dist_mask=np.where(bayestar.distances.value*1000.0 <= 1000.0)
sub_gamma_ring=sub_gamma_ring[:,:,dist_mask[0]]
sub_gamma_ring.shape

ipyvolume.quickvolshow(sub_gamma_ring.T, stero=True, data_min=None, level=[0.28, 0.21, 0.19], lighting=True, ambient=0.67, width=512, height=256, stereo=True, opacity=0.06)


    

    




(200, 200, 31)
[30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]






    
Out[8]:





(200, 200, 13)






    






 
 

In [12]:

    

orion_ring = get_cube(205, -15, delta=15, step=0.1)
sub_orion_ring = subtract_cube(orion_ring)

dist_mask=np.where(bayestar.distances.value*1000.0 <= 1000.0)
sub_orion_ring=sub_orion_ring[:,:,dist_mask[0]]

ipyvolume.quickvolshow(sub_orion_ring.T, stero=True, data_min=None, level=[0.22, 0.16, 0.20], \
                       lighting=True, ambient=0.67, width=512, height=256, stereo=True, opacity=0.06)


    

    




(300, 300, 31)
[30, 29, 28, 27, 26, 25, 24, 23, 22, 21, 20, 19, 18, 17, 16, 15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0]






    






 
 

In [10]:

    

dist_mask=(bayestar.distances.value*1000.0 > 200.0) & (bayestar.distances.value*1000.0 < 1000.0)


    

In [12]:

    

dist_mask


    

    
Out[12]:





array([False, False, False, False, False, False,  True,  True,  True,
        True,  True,  True, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False], dtype=bool)

In [13]:

    

dist_mask=np.where((bayestar.distances.value*1000.0 > 200.0) & (bayestar.distances.value*1000.0 < 1000.0)==True)


    

In [14]:

    

dist_mask


    

    
Out[14]:





(array([ 6,  7,  8,  9, 10, 11]),)

In [ ]: