plotmaskhealpy

In healpy how to split mask/map into 2 parts?

I have a mask/map as a fits file, downloadable here, and I can plot this mask with healpy using this code:

import healpy as hp
from healpy.newvisufunc import projview, newprojplot

wmap_map_I = hp.read_map('mask.fits')
projview(
    wmap_map_I,
    graticule=True, graticule_labels=True, xlabel="longitude", ylabel="latitude",
    projection_type="mollweide",
    coord=["E"],
    title="Histogram equalized Ecliptic",
    unit="mK",
    norm="hist",
    min=-1,
    max=1,
);

It produces this figure:

enter image description here

What I want to do is split the mask into 2 new masks - one should contain only part A and the other should contain only part B. I don't know how to do this, I thought of making some conditional expression about all points for which longitude + latitude < xyz and then the opposite, but the fits file doesn't contain longitude and latitude.

So how can I split the mask into 2? Thx

Solution

  • I followed EHivon's answer and digged into healpy. I came up with a solution following his/her advice, but first I post 2 solutions that I discovered on the way on my own

    ### 1. SOLUTION: We convert to galactic coordinates and use the order of the pixels to split

# read in the map
map = hp.read_map('map.fits')

# first we define a operator that will rotate from celestial coordinates to galactic
rot_gal2eq = hp.Rotator(coord=['G', 'C'], inv=True)

# rotate the whole map
m_rotated = rot_gal2eq.rotate_map_alms(map)

# there are now some artefacts in the map so I round it
m_rotated = np.rint(m_rotated)

# set the second half of the rotated mask to 0, since we want to split it in half
# For info, the indexing of the elements in a map follows a ring order going down from the most norther point
# to the most souther point. So the first half of the map element are all in the northern sky, the second half
# are all in the souther sky
m_rot_N = m_rotated.copy()
m_rot_N[len(m_rot_N)//2:] =0

# do the same for the south sky
m_rot_S = m_rotated.copy()
m_rot_S[:len(m_rot_N)//2] =0

# rotate back both maps
rot_inv = hp.Rotator(coord=['C', 'G'], inv=True)
m_rot_N_inv = rot_inv.rotate_map_alms(m_rot_N)
m_rot_S_inv = rot_inv.rotate_map_alms(m_rot_S)

# again, round the elements of the masks because of artefacrts
m_rot_N_inv = np.rint(m_rot_N_inv)
m_rot_S_inv = np.rint(m_rot_S_inv)

# eventually, saved the new masks
hp.fitsfunc.write_map('map_N.fits', m_rot_N_inv)
hp.fitsfunc.write_map('map_S.fits', m_rot_S_inv)

    The above we can plot like this:

    hp.mollview(m_rot_N, max=1)
hp.graticule();
hp.mollview(m_rot_S, max=1)
hp.graticule();

    enter image description here

    Here is another way to do the split, substitute the last lines above with these:

    ### 2. SOLUTION: Draw around a point a disc that splits the map as you like
# define a vector pointing at a point on the sky that you want to center on
# they we will draw a half-sky circle around that point and so we can split the map in 2
vec = hp.ang2vec(0, 0) # (0,0) are the angles on the sky pointing at the north pole
ipix_disc = hp.query_disc(nside=hp.pixelfunc.get_nside(m_rotated), vec=vec, radius=np.radians(90))

# set all the values inside the circle to 0, giving you a condition to effectively split the map
m_rot_N = m_rotated[ipix_disc] = 0

    And the solution suggest by EHivon is below, I actually made a slightly more complicated structure with for demonstration purposes:

    ### 3. SOLUTION: turn the pixels to angles and make conditions for the angles
# get the angles of all the pixels on the map
angles = hp.pixelfunc.pix2ang(nside=hp.pixelfunc.get_nside(m_rotated), ipix=np.arange(len(m_rotated)), lonlat=True)
# set all the values 0 that fullfill a condition
m_rotated[(angles[0]<60)|(angles[0]>315)] = 0.5 #angles[0] are the longitude angles, angles[1] are the latitudes of all the pixels

    enter image description here