How to find all pixel indices within specified coordinates?

Question

There are a number of Healpix IDL routines designed to find the pixel indices belonging to some geometrical area (e.g. spherical triangle, spherical polygon) defined by its vertices. There are the query_* routines (e.g. query_triangle). See the documentation here:

http://healpix.jpl.nasa.gov/html/idlnode45.htm

I wish to use these pixel indices in my healpy program. Either

(A) I could save the output list of pixel indices from the IDL query_* routines in a data_file.save format. You could then import this .sav file of pixel indices into Python using the various modules, e.g. http://www.astropython.org/packages/idlsave94/

(B) It would be far more convenient to somehow not use IDL at all! Healpy has several pixel-related functions, but there seems to be no way to "convert" the IDL query_* routines using healpy alone.

Is there a way to do query_polygon using healpy? Is it possible to do this?

score 1 · Accepted Answer · answered Jun 07 '16 at 07:51

It's a bit unclear what you want precisely, but you can use query_polygon from healpy. For example:

nside = 512
vertices = numpy.array([[....]])
healpy.query_polygon(nside, vertices)

will return the pixels that are inside the polygon. vertices is an (N, 3) array of the polygon's vertices.

From the built-in help:

query_polygon(nside, vertices, inclusive=False, fact=4, nest=False, ndarray buff=None
Returns the pixels whose centers lie within the convex polygon
defined by the *vertices* array (if *inclusive* is False), or which
overlap with this polygon (if *inclusive* is True).

Parameters
----------
nside : int
  The nside of the Healpix map.
vertices : float, array-like
  Vertex array containing the vertices of the polygon, shape (N, 3).
inclusive : bool, optional
  If False, return the exact set of pixels whose pixel centers lie
  within the polygon; if True, return all pixels that overlap with the
  polygon, and maybe a few more. Default: False.
fact : int, optional
  Only used when inclusive=True. The overlapping test will be done at
  the resolution fact*nside. For NESTED ordering, fact must be a power of 2, less than 2**30,
  else it can be any positive integer. Default: 4.
nest: bool, optional
  if True, assume NESTED pixel ordering, otherwise, RING pixel ordering
buff: int array, optional
  if provided, this numpy array is used to contain the return values and must be
  at least long enough to do so

Returns
-------
ipix : int, array
  The pixels which lie within the given polygon.

Note
----
This method is more efficient in the RING scheme.
For inclusive=True, the algorithm may return some pixels which don't overlap
with the disk at all. The higher fact is chosen, the fewer false positives
are returned, at the cost of increased run time.

How to find all pixel indices within specified coordinates?

1 Answers1