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

Question

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:

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

Answer 1

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

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

Answer 2

your mask was apparently defined in Galactic coordinates, with the A part centered on South Galactic pole, and the B part around the NGP. It should appear if you run mollview with coord=['C','G']. Are you sure the mask is stored in Ecliptic coordinates ? It looks like Celestial (or Equatorial) coordinates instead.

To split the mask, get the celestial 3D coordinates of each valid pixel with pix2vec, than rotate them to Galactic using healpy.rotate tools and put those with glat <0 in file A, and those with glat > 0 in file B.