2014-08-28 2 views
1

Я пытаюсь создать луч на карте healpix, используя healpy. Для начала я хотел бы иметь возможность создавать 2D-гауссово в проекции mollweide, но я действительно не знаю, с чего начать.Построение массива numpy в healpy

я могу определить 2D гауссово:

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)) 

таким образом, что я могу построить 3D 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, поэтому * map * представляет собой 1D-массив, где индексы соответствуют пикселям. Если вам просто нужна проекция Mollweide, вы можете просто использовать «matplotlib», см. Http://matplotlib.org/examples/pylab_examples/geo_demo.html –

ответ

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) вы s et - значение для этой точки - значение из вашей гауссовой функции. то есть:

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