In [15]:
import sys, os
from os.path import join
#from mayavi import mlab
import healpy as hp
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline

#sys.path
ld_lib=os.getenv('LD_LIBRARY_PATH')
scratch='/global/cscratch1/sd/fantaye' #os.getenv('SCRATCH')
home=os.getenv('HOME')
print('home:',home)
print('scratch',scratch)
romepy=home+"/romepy/"
sys.path.append(romepy)


('home:', '/global/homes/f/fantaye')
('scratch', '/global/cscratch1/sd/fantaye')

In [17]:
#my routines
import dotunf as unf

srcdir=join(scratch,'seedlets/fits4096/SurfDensMap/')
rootdir='../'
fdir=rootdir+'figures/'
if not os.path.exists(fdir):
    os.makedirs(fdir)
    

font = {'family' : 'fantasy',
        #'weight' : 'bold',
        'size'   : 16}

matplotlib.rc('font', **font)

#print os.path.abspath(pview.__file__)

nmax=128
lmax=2000
nshell=2*nmax
nside=1024
npix=hp.nside2npix(nside)

In [32]:
#os.listdir(home+'/Radial3Dneedlet/output/seedlets')

In [ ]:
mapj=[]
for j in range(1):
    fname=join(srcdir,'a2beta/mapout/map.unf_j'+str(j)+'_glnpow2.unf')
    m=unf.read_unf(fname).reshape((npix,nshell),order='F')
    mapj.append(m)
    
mapj[0].shape

In [8]:
fig, (ax1, ax2)=plt.subplots(ncols=2,figsize=(15,6))
hp.mollview(mapj[0][:,0],ax=ax1)
hp.mollview(mapj[0][:,1],ax=ax2)


:)

In [34]:
z='''
The following code is taken from 
http://zonca.github.io/2013/03/interactive-3d-plot-of-sky-map.html


fig, (ax1, ax2)=plt.subplots(ncols=2,figsize=(15,6))
# load a WMAP map
npix=hp.nside2npix(512)
m = np.random.randn(npix)
#m = hp.read_map("data/wmap_band_iqumap_r9_7yr_W_v4.fits", 0) * 1e3 # muK
hp.mollview(m,ax=ax1)

nside = hp.npix2nside(len(m))

#value limit
vmin = -1e3; vmax = 1e3
# size of the grid
xsize = ysize = 1000

theta = np.linspace(np.pi, 0, ysize)
phi   = np.linspace(-np.pi, np.pi, xsize)
longitude = np.radians(np.linspace(-180, 180, xsize))
latitude = np.radians(np.linspace(-90, 90, ysize))

# project the map to a rectangular matrix xsize x ysize
PHI, THETA = np.meshgrid(phi, theta)
grid_pix = hp.ang2pix(nside, THETA, PHI)
grid_map = m[grid_pix]

# Create a sphere
r = 0.3
x = r*np.sin(THETA)*np.cos(PHI)
y = r*np.sin(THETA)*np.sin(PHI)
z = r*np.cos(THETA)

mlab.figure(1, bgcolor=(1, 1, 1), fgcolor=(0, 0, 0), size=(400, 300))
mlab.clf()

mlab.mesh(x, y, z, scalars=grid_map, colormap="jet", vmin=vmin, vmax=vmax) 
'''

In [ ]: