Healpy coordinate error after interpolation: appearance of bisector

Question

I have a coarse skymap made up of 128 points, of which I would like to make a smooth healpix map (see attached Figure, LHS). Figures referenced in the text:

I load my data, then make new longitude and latitude arrays of the appropriate pixel length for the final map (with e.g. nside=32).

My input data are:

lats = pi/2 + ths   # theta from 0, pi, size 8
lons = phs          # phi from 0, 2pi, size 16
data = sky_data[0]  # shape (8,16)

New lon/lat array size based on number of pixels from nside:

nside = 32 
pixIdx = hp.nside2npix(nside) # number of pixels I can get from this nside 
pixIdx = np.arange(pixIdx) # pixel index numbers

I then find the new data values for those pixels by interpolation, and then convert back from angles to pixels.

# new lon/lat
new_lats = hp.pix2ang(nside, pixIdx)[0] # thetas I need to populate with interpolated theta values
new_lons = hp.pix2ang(nside, pixIdx)[1] # phis, same

# interpolation
lut = RectSphereBivariateSpline(lats, lons, data, pole_values=4e-14)
data_interp = lut.ev(new_lats.ravel(), new_lons.ravel()) #interpolate the data
pix = hp.ang2pix(nside, new_lats, new_lons) # convert latitudes and longitudes back to pixels

Then, I construct a healpy map with the interpolated values:

healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double) # create empty map
healpix_map[pix] = data_interp # assign pixels to new interpolated values
testmap = hp.mollview(healpix_map)

The result of the map is the upper RHS of the attached Figure.

(Forgive the use of jet -- viridis doesn't have a "white" zero, so using that colormap adds a blue background.)

The map doesn't look right: you can see from the coarse map in the Figure that there should be a "hotspot" on the lower RHS, but here it appears in the upper left.

As a sanity check, I used matplotlib to make a scatter plot of the interpolated points in a mollview projection, Figure 2, where I removed the edges of the markers to make it look like a map ;)

ax = plt.subplot(111, projection='astro mollweide')
ax.grid()
colors = data_interp
sky=plt.scatter(new_lons, new_lats-pi/2, c = colors, edgecolors='none', cmap ='jet')
plt.colorbar(sky, orientation = 'horizontal')

You can see that this map, lower RHS of attached Figure, produces exactly what I expect! So the coordinates are ok, and I am completely confused.

Has anyone encountered this before? What can I do? I'd like to use the healpy functions on this and future maps, so just using matplotlib isn't an option.

Thanks!

Answer 1

I'm no expert in healpix (actually I've never used it before - I'm a particle physicist), but as far as I can tell it's just a matter of conventions: in a Mollweide projection, healpy places the north pole (positive latitude) at the bottom of the map, for some reason. I'm not sure why it would do that, or whether this is intentional behavior, but it seems pretty clear that's what is happening. If I mask out everything below the equator, i.e. keep only the positive-latitude points

mask = new_lats - pi/2 > 0
pix = hp.ang2pix(nside, new_lats[mask], new_lons[mask])
healpix_map = np.zeros(hp.nside2npix(nside), dtype=np.double)
healpix_map[pix] = data_interp[mask]
testmap = hp.mollview(healpix_map)

it comes up with a plot with no data above the center line:

At least it's easy enough to fix. mollview admits a rot parameter that will effectively rotate the sphere around the viewing axis before projecting it, and a flip parameter which can be set to 'astro' (default) or 'geo' to set whether east is shown at the left or right. A little experimentation shows that you get the coordinate system you want with

hp.mollview(healpix_map, rot=(180, 0, 180), flip='geo')

In the tuple, the first two elements are longitude and latitude of the point to set in the center of the plot, and the third element is the rotation. All are in degrees. With no mask it gives this:

which I believe is just what you're looking for.