CODEWITHSUNDEEP
pythonmapsastronomyhealpy
Steve
Steve

Reputation: 11

Conserving pixel counts with HEALPy

I have a list of galaxies to plot onto a healpix map (which I'm using healpy to do) each galaxy has a set flux and I need to have them plotted in such a way that the flux for each galaxy is conserved on the map.

This is my code:

import numpy as np
import matplotlib.pyplot as plt
import healpy as hp
pi = np.pi

nside = 8
xsize = 100

ra = np.array([pi/4,pi/3])
dec = np.array([pi/4,pi/3])
flux = np.array([10,20])

hpm = np.zeros(hp.nside2npix(nside)) #Blank healpix map
pixindex = hp.ang2pix(nside, dec, ra)

np.add.at(hpm,pixindex,flux) #Add flux onto correct pixels

img=hp.mollview(hpm,coord=['E'],xsize=xsize,return_projected_map=True)
print(np.sum(img[img>0]))

The result I get is 140 and not 30 which is the true sum of the fluxes.

I get what is going on and that the same flux is being spread over multiple pixels (6 for the first galaxy and 4 pixels for the second) and I'm aware I could just do:

newimg = img * (np.sum(flux)/np.sum(img[img>0]))

and this would conserve the total photon count but it wouldn't necessarily conserve the photon count of each galaxy which is what I need. i.e this method ends up with the first galaxy contributing a flux of 12.86 and the second galaxy a flux of 17.14.

Is there a way of working out before how many pixels each coordinate will take up then altering the amount of flux dumped based on this?

Thanks in advance!

Upvotes: 1

Views: 417

Answers (1)

Marco
Marco

Reputation: 21

The xsize parameter for the function hp.mollview should be used for plotting purpose only. If you would like to manipulate the resolution of the map, use hp.pixelfunc.ud_grade

For example, if you want to go from nside=8 to nside=32,

After the line

np.add.at(hpm, pixindex, flux) #Add flux onto correct pixels

use the ud_grade function with power=-2, so the total flux can be conserved:

hpm_nside_32 = hp.pixelfunc.ud_grade(hpm, power=-2, nside_out=32)

the sum 

np.sum(hpm_nside_32)

will be conserved at 30.

I can't offer a solution if you need the mollview to have the flux conserved. The closest I can get is by scaling the img value down based on the ratio of the number of pixels in the mollview image and the number of pixels in the hpm. The first term xsize * xsize / 2. is the number of pixels in the mollview, the second term 2. * np.pi / 8. is the ratio of the area of an ellipse with semi-minor axis half the length of the semi-major axis pi * r * 2r to the area of a rectangle 4r * 2r.

len(hpm) / ((xsize * xsize / 2.) * (2. * np.pi / 8.)) * np.sum(img[img > 0])

When xsize = 100, the sum will become 27.380; when xsize = 1000, the approximation is much better, at 29.981; when xsize = 10000, total flux becomes 29.988.

An alternative way to give you a good approximation is to compute the ratio of the number of non -inf pixels in img to the number of pixel in your map (which is 768 for nside=8):

float(len(hpm)) / float(np.sum(img>=0)) * np.sum(img[img>0])

At xsize = 100, 1000, 10000, the flux will be at 28.295, 30.075, 29.998 respectively.

Upvotes: 1

Related Questions