2014-08-28 98 views
1

我试图在healpix地图上使用healpy来产生一个波束。对于初学者,我希望能够在mollweide投影中产生2D高斯,但我真的不知道从哪里开始。在healpy中绘制一个numpy数组

我可以定义一个二维高斯:

import numpy as np 
def gaussian_2D(x,y,mu_x=0.,mu_y=0.,sig_x=1.,sig_y=1.): 
    return np.exp(-0.5*(((x-mu_x)/sig_x)**2 + ((y-mu_y)/sig_y)**2)) 

,这样我可以建立类似的三维X,Y,Z空间:

delta = 0.025 
x = np.arange(-4, 4, delta) 
y = np.arange(-4, 4, delta) 
X, Y = np.meshgrid(x,y) 
Z = gaussian_2D(X,Y) 

,但在这里我几乎失去了,并且无法追踪有关如何和/或投影什么的许多有用的文档。任何有关攻击方向的建议都将非常感谢!

+0

'healpy'使用HEALPix像素化,所以a * map *是一维数组,其中索引对应于像素。如果你只需要一个Mollweide投影,你可以使用'matplotlib',参见http://matplotlib.org/examples/pylab_examples/geo_demo.html – 2014-08-28 22:28:56

回答

0

这里是我如何做到这一点:

使用一个小窍门。我在所需的高斯中心插入一个点,然后使用“涂抹”创建一些西格玛高斯。

下面是一些例子:

#!/usr/bin/env python 
import numpy as np 
import healpy as hp 
import pylab as pl 

NSIDE=512 #the map garannularity 

m_sm=np.arange(hp.nside2npix(NSIDE)) # creates the map 
m_sm=m_sm*0. # sets all values to zero 

theta=np.radians(80.) # coordinates for the gaussian 
phi=np.radians(20.) 

indx=hp.pixelfunc.ang2pix(NSIDE,theta,phi) # getting the index of the point corresponding to the coordinates 
m_sm[indx]=1. # setting that point value to 1. 

gmap=hp.smoothing(m_sm, sigma=np.radians(20.),verbose=False,lmax=1024) # creating a new map, smmeared version of m_sm 

hp.mollview(gmap, title="Gaussian Map") #draw it 
pl.show() 

现在如果你想这样做手工,你可以使用一个功能对于高斯

1)你给它的一些坐标

2 )您检索对应于该坐标的索引使用:

indx=hp.pixelfunc.ang2pix(NSIDE,theta,phi) 

3)您将该点的值与高斯函数的值相比较。即:

my_healpy_map[indx]=my_gauss(theta, phy, mean_theta, mean_phy, sigma_theta, sigma_phy)