2

Basically all I want is the functionally of hp.cartview, but I don't want my machine to waste memory plotting the actual map every single time I call the cartview function. How can I obtain a cartesian projection in healpy in the form of a 2d array without having to plot the projection every time?

Jsn
  • 107
  • 6

1 Answers1

2

First, let me point out that reproject might be a better tool for the job.

You can build a WCS object or a FITS header and then reproject your HEALPix map onto that, and subsequently plot it with wcsaxes, which gives you full support for real-world coordinate pixels (instead of pixel coordinates only).

If you really want to use healpy for these cartview cutouts instead, you can use the underlying healpy.projector.CartesianProj class:

from functools import partial

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

# Build a map
nside = 64
npix = hp.nside2npix(nside)
hpxmap = np.arange(npix)

# Get the cutout via a cartesian projection
lonra = [30, 40]
latra = [-10, 10]

proj = hp.projector.CartesianProj(
    lonra=lonra, latra=latra,
    coord='G',
    xsize=n_pixels, ysize=n_pixels)
reproj_im = proj.projmap(hpxmap, vec2pix_func=partial(hp.vec2pix, nside))

# Plot the cutout
plt.imshow(reproj_im, origin='lower', interpolation='nearest')

Good luck, let me know if you have any follow-up questions!

Daniel Lenz
  • 3,334
  • 17
  • 36
  • This is great. Thank you! – Jsn Aug 23 '19 at 19:57
  • @Daniel Lenz: How do you handle cases where parts of the patch would cross over the pole and fall on the other side? If it's a point I guess one would have - new lattitude = 180 - old lattitude - new longitude = 180 + old longitude But how does one specify new lonra and latra in that case? – Minh N May 11 '21 at 17:04