Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
tbdata = hdulist[1].data
pix=hp.ang2pix(nside,(90-tbdata['DEC'])*np.pi/180,tbdata['RA']*np.pi/180)
n=np.bincount(pix,minlength=npix)
e1=np.bincount(pix,minlength=npix,weights=tbdata['E1'])
e2=np.bincount(pix,minlength=npix,weights=tbdata['E2'])
nmap+=n
e1map+=e1
e2map+=e2
ifile+=1
ndens=(np.sum(nmap)+0.0)/(4*np.pi)
mp_e1=e1map/nmap; mp_e1[nmap<=0]=0
mp_e2=e2map/nmap; mp_e2[nmap<=0]=0
mp_d=(nmap+0.0)/np.mean(nmap+0.0)-1
mp_db,mp_E,mp_B=hp.alm2map(hp.map2alm(np.array([mp_d,mp_e1,mp_e2]),pol=True),pol=False,nside=nside)
lt,cls_dd=np.loadtxt('test/outlj_cl_dd.txt',unpack=True);
lt,clt_dl=np.loadtxt('test/outlj_cl_d1l2.txt',unpack=True);
lt,clt_ll=np.loadtxt('test/outlj_cl_ll.txt',unpack=True);
lt,clt_kd=np.loadtxt('test/outlj_cl_dc.txt',unpack=True);
lt,clt_kk=np.loadtxt('test/outlj_cl_cc.txt',unpack=True);
lt,clt_id=np.loadtxt('test/outlj_cl_di.txt',unpack=True);
lt,clt_ii=np.loadtxt('test/outlj_cl_ii.txt',unpack=True);
cln_dd=np.ones_like(lt)/ndens
clt_dd=cls_dd+cln_dd
cld_dd,cld_ee,cld_bb,cld_de,cld_eb,cld_db=hp.anafast(np.array([mp_d,mp_e1,mp_e2]),pol=True);
ld=np.arange(len(cld_dd));
#Analyze kappa
mp_k=hp.read_map("test/out_kappa_z000.fits")
cld_kk=hp.anafast(mp_k); ld=np.arange(len(cld_kk))
# Call pys2let and compute directional wavelet transform.
# Returns an array of MW maps
# The way to access them is described below.
print 'Running analysis_lm2wav'
f_wav, f_scal = analysis_lm2wav(f_lm, B, L, J_min, N, spin, upsample)
print 'Done'
print 'size f_scal f_wav', f_scal.size, f_wav.size, f_wav.size / f_scal.size
# Uses synthesis to reconstruct the input alms.
print 'Running synthesis_wav2lm'
f_lm_rec = synthesis_wav2lm(f_wav, f_scal, B, L, J_min, N, spin, upsample)
print 'Done'
f_rec = hp.alm2map(f_lm_rec, nside=nside, lmax=L-1)
#plot to compare the input/output maps
#hp.mollview(f)
#hp.mollview(f_rec)
# Convert to MW sampling from spherical harmonics
f_mw = alm2map_mw(f_lm, L, spin)
f_lm_rec = map2alm_mw(f_mw, L, spin) # To check that the analysis routine works
f_mw_rec = alm2map_mw(f_lm, L, spin)
# Home made plotting routine! inputs : function f (1D array of MW signal), bandlimit L, plot axis ax, and title
def myplot(f, L, ax, title=''):
thetas, phis = mw_sampling(L) # Compute the theta and phi values of the MW equiangular grid.
ntheta = len(thetas)
nphi = len(phis)
arr = f.reshape((ntheta, nphi)) # Convert the input MW 1D array into 2D array with rows for theta and columns for phi. As simple as that!
import matplotlib.pyplot as plt
nside = 128
L = 128
J_min = 1
B = 3
J = pys2let_j_max(B, L, J_min)
# The filename of some random healpix map
fname = '/Users/bl/Dropbox/Wavelets/s2let/data/somecmbsimu_hpx_128.fits'
# Read healpix map and compute alms.
# f_lm has size L*(L+1)/2
f_ini = hp.read_map(fname) # Initial map
f_lm = hp.map2alm(f_ini, lmax=L-1) # Its alms
f = hp.alm2map(f_lm, nside=nside, lmax=L-1) # Band limited version
hp.mollview(f)
# Call pys2let and compute wavelet transform. Returns the harmonic coefficients of the wavelets.
# f_scal_lm has size L*(L+1)/2
# f_wav_lm has size L*(L+1)/2 by J-J_min+1
f_wav_lm, f_scal_lm = analysis_axisym_lm_wav(f_lm, B, L, J_min)
# Reconstruct healpix maps on the sphere and plot them
f_scal = hp.alm2map(f_scal_lm, nside=nside, lmax=L-1)
hp.mollview(f_scal)
f_wav = np.empty([12*nside*nside, J-J_min+1])
for j in range(J-J_min+1):
flm = f_wav_lm[:,j].ravel()
f_wav[:,j] = hp.alm2map(flm, nside=nside, lmax=L-1)
hp.mollview(f_wav[:,j])
J_min = 0
L_transitions = np.array([64])
L = 256
hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits, J, L_bounds = construct_hybrid_tiling(L, J_min, Bs, L_transitions)
hybrid_wav_l = hybrid_wav_l.T.ravel()
res = verify_tiling(L, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)
print 'Is the tiling OK and giving an invertible wavelet transform?', res
# The filename of some arbitrary healpix map
fname = '/Users/bl/Dropbox/Astrodata/SDSS/Fields/Planck_EBV_256rQ.fits'
#fname = '/Users/bl/Dropbox/Wavelets/s2let/data/somecmbsimu_hpx_128.fits'
f_ini = hp.read_map(fname) # Initial map
f_lm = hp.map2alm(f_ini, lmax=L-1) # Its alms
f = hp.alm2map(f_lm, nside=nside, lmax=L-1) # Band limited version
print 'Running analysis_lm2wav'
f_wav, f_scal = analysis_lm2wav_manualtiling(f_lm, L, N, spin, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)
# Uses synthesis to reconstruct the input alms.
print 'Running synthesis_wav2lm'
f_lm_rec = synthesis_wav2lm_manualtiling(f_wav, f_scal, L, N, spin, hybrid_scal_l, hybrid_wav_l, hybrid_scal_bandlimit, hybrid_wav_bandlimits)
f_rec = hp.alm2map(f_lm_rec, nside=nside, lmax=L-1)
# Convert to MW sampling from spherical harmonics
f_mw = alm2map_mw(lm_hp2lm(f_lm, L), L, spin)
#f_lm_rec = map2alm_mw(f_mw, L, spin) # To check that the analysis routine works
#f_mw_rec = alm2map_mw(f_lm_rec, L, spin)
f_mw_rec = alm2map_mw(lm_hp2lm(f_lm_rec, L), L, spin)
m_rotated : np.ndarray
The rotated HEALPix array.
"""
npix = len(m)
nside = hp.npix2nside(npix)
theta = 0.5 * np.pi - dec
phi = ra
if method == 'fft':
if nest:
m = hp.reorder(m, n2r=True)
alm = hp.map2alm(m)
hp.rotate_alm(alm, -phi, -theta, 0.0)
ret = hp.alm2map(alm, nside, verbose=False)
if nest:
ret = hp.reorder(ret, r2n=True)
elif method == 'direct':
R = hp.Rotator(
rot=np.asarray([0, theta, -phi]),
deg=False, inv=False, eulertype='Y')
theta, phi = hp.pix2ang(nside, np.arange(npix), nest=nest)
ipix = hp.ang2pix(nside, *R(theta, phi), nest=nest)
ret = m[ipix]
else:
raise ValueError('Unrecognized method: {0}'.format(method))
return ret
m_rotated : np.ndarray
The rotated HEALPix array.
"""
npix = len(m)
nside = hp.npix2nside(npix)
theta = 0.5 * np.pi - dec
phi = ra
if method == 'fft':
if nest:
m = hp.reorder(m, n2r=True)
alm = hp.map2alm(m)
hp.rotate_alm(alm, -phi, -theta, 0.0)
ret = hp.alm2map(alm, nside, verbose=False)
if nest:
ret = hp.reorder(ret, r2n=True)
elif method == 'direct':
R = hp.Rotator(
rot=np.asarray([0, theta, -phi]),
deg=False, inv=False, eulertype='Y')
theta, phi = hp.pix2ang(nside, np.arange(npix), nest=nest)
ipix = hp.ang2pix(nside, *R(theta, phi), nest=nest)
ret = m[ipix]
else:
raise ValueError('Unrecognized method: {0}'.format(method))
return ret
Args:
imap: ndmap of shape (Ny,Nx)
lmax: integer specifying maximum multipole of map
nside: integer specifying nside of healpix map
Returns:
retmap: (Npix,) healpix map as array
"""
import healpy as hp
alm = curvedsky.map2alm(imap, lmax=lmax, spin=0)
if alm.ndim > 1:
assert alm.shape[0] == 1
alm = alm[0]
retmap = hp.alm2map(alm.astype(np.complex128), nside, lmax=lmax)
return retmap
def mapFromClm_fast(clm, nside):
"""
Given an array of C_{lm} values, produce a pixel-power-map (non-Nested) for
healpix pixelation with nside
@param clm: Array of C_{lm} values (inc. 0,0 element)
@param nside: Nside of the healpix pixelation
return: Healpix pixels
Use Healpix spherical harmonics for computational efficiency
"""
maxl = int(np.sqrt(len(clm))) - 1
alm = almFromClm(clm)
h = hp.alm2map(alm, nside, maxl, verbose=False)
return h
sigma : float, scalar, optional
The sigma of the Gaussian used to smooth the map (applied on alm)
[in radians]
pol : bool, optional
If True, assumes input alms are TEB. Output will be TQU maps.
(input must be 1 or 3 alms)
If False, apply spin 0 harmonic transform to each alm.
(input can be any number of alms)
If there is only one input alm, it has no effect. Default: True.
Returns
-------
maps : array or list of arrays
A Healpix map in RING scheme at nside or a list of T,Q,U maps (if
polarized input)'''
return healpy.alm2map(self.get_data(), nside,
lmax=self._lmax, mmax=self._mmax,
pixwin=pixwin, fwhm=fwhm, sigma=sigma, pol=pol, verbose=verbose)
def from_map(self, data, iter=3, pol=True, use_weights=False, gal_cut=0):