In [21]:

    

from astropy.io import fits

hdulist = fits.open('l1448_13co.fits')
cube = hdulist[0].data
# And to write the cube to disk as 8-bit integers in a binary file (as expected by my loader on Unity's side):
cube_8bit = ((cube-cube.min())/(cube.max()-cube.min())*255).astype('uint8')
cube_8bit


    

In [32]:

    

# cube size needs to be the same on width, height, and depth - NOT necessary

# unity has non-power of two

import numpy as np
# npad is a tuple of (n_before, n_after) for each dimension
b = np.pad(cube_8bit, pad_width=((0, 11), (0,23), (0,23)), mode='constant', constant_values=0)
b.transpose().tofile('l1448co.raw')
b.shape


    

    
Out[32]:





(64, 128, 128)

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 [19]:

    

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, deltal=5., deltab=5., step=0.05):
    # set up a grid of coordinates to plot, centered on the Aquila South cloud
    l = np.arange(l0 - deltal, l0 + deltal, step)
    b = np.arange(b0 - deltab, b0 + deltab, 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 [65]:

    

orion_ring = get_cube(205, -15, deltal=18, deltab=10, step=0.1)
sub_orion_ring = subtract_cube(orion_ring)


# dist_mask=np.where(bayestar.distances.value*1000.0 <= 1000.0)
dist_mask=np.where((bayestar.distances.value*1000.0 > 200.0) & (bayestar.distances.value*1000.0 < 1000.0)==True)
sub_orion_ring=sub_orion_ring[:,:,dist_mask[0]]
sub_orion_ring.shape

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

import ipyvolume
from astropy.io import fits
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)


    

    




(200, 360, 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[65]:





(360, 200, 6)






    






 
 

In [69]:

    

# And to write the cube to disk as 8-bit integers in a binary file (as expected by my loader on Unity's side):
# TODO: a better way to normalize display, the median of array is only 0.004, while the max() is 3.82
# ring_cube_8bit = ((sub_orion_ring-sub_orion_ring.min())/(sub_orion_ring.max()-sub_orion_ring.min())*255).astype('uint8')
sub_orion_ring[sub_orion_ring<sub_orion_ring.mean()] = 0.
sub_orion_ring.max()
ring_cube_8bit = ((sub_orion_ring-sub_orion_ring.min())/(sub_orion_ring.max()-sub_orion_ring.min())*255).astype('uint8')
print(sub_orion_ring.min(), sub_orion_ring.max(), sub_orion_ring.mean())
np.median(sub_orion_ring)
import numpy as np
# npad is a tuple of (n_before, n_after) for each dimension
npad_cube = np.pad(ring_cube_8bit, pad_width=((0, 152), (0,56), (0,2)), mode='constant', constant_values=0)
npad_cube.transpose().tofile('orion_rings.raw')
npad_cube.mean()

ipyvolume.quickvolshow(npad_cube, 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)


    

    
Out[69]:





3.8276898860931396






    




0.0 3.82768988609 0.0339702323395






    
Out[69]:





0.0






    
Out[69]:





0.90210723876953125






    






 
 

In [ ]: