I have a healpy full-sky map, and need to chose some sky patches centered around given pixels. hp.query_disc finds the pixels within a given disc, which is nice. But I would like to be able to add/subtract on different sky patches, and so the fact that different discs cannot be the same size causes a problem.
I have also used
hp.gnomview(return_projected_map=True)
to look at my sky patches, and found that each returned 2D array has the same size. However since I would probably do stuff on many submaps, I do not want the maps to be drawn every time. Is it possible to obtain only the 2D array? And if so, is it fine to use plt.imshow after adding/subtracting those arrays?
I feel like there should be more elegant ways to achieve what I want -- Are there any existing routines that I overlooked? Thanks!