hpmoc.utils module

Utility functions used across healpix_skymap classes.

exception hpmoc.utils.EmptyStream

Bases: OSError

Raised when a file stream returns no further content.

class hpmoc.utils.TmpGunzipFits(infile: Union[IO, str])

Bases: object

A context manager that unzips a Gzip file to a temporary .fits file, returning a NamedTemporaryFile pointing to the tempfile, and deletes the tempfile when you’re done with it. Can only do binary reads, so mode cannot be specified. Pass it a file name or a gzip file object that’s already been opened in binary read (‘rb’) mode.

filename = None

infile: Union[IO, str]

hpmoc.utils.check_valid_nside(nside)

Checks whether nside is a valid HEALPix NSIDE value. Raises a ValueError if it is not.

Parameters: nside (int or array) – An integer or array of integers representing HEALPix NSIDE values.
Raises: ValueError – If any of the values are not valid HEALPix NSIDE values.

Examples

Does nothing if you provide a valid NSIDE value (floats are accepted as long as their exact values are valid):

>>> check_valid_nside([16, 4])
>>> check_valid_nside(4)
>>> check_valid_nside(1024)
>>> check_valid_nside([1024., 4096])

A descriptive value error is raised if invalid values are provided:

>>> try:
...     check_valid_nside([4.1, 4])
... except ValueError as err:
...     print("Caught exception:", err)
Caught exception: Not a valid NSIDE value: [4.1, 4]
>>> try:
...     check_valid_nside(17)
... except ValueError as err:
...     print("Caught exception:", err)
Caught exception: Not a valid NSIDE value: [17]

hpmoc.utils.check_valid_nuniq(indices)

Checks that indices are valid NUNIQ indices.

Raises: ValueError – If indices are not valid NUNIQ indices, i.e. they are not integers greater than 3.

hpmoc.utils.dangle_rad(ra, dec, mapra, mapdec)

Get an array of angular distances in radians between the point specified in ra, dec and the points specified in mapra, mapdec. All angles, including the return value, are specified in radians.

sin(DEC1)*sin(DEC2) + cos(DEC1)*cos(DEC2)*cos(RA2-RA1)

So we get the angle by taking the arccos of this function.

Parameters

ra (float) – Right ascension of the single point [radians]
dec (float) – Declination of the single point [radians]
mapra (array) – Right ascensions whose angular distances to the single point will be calculated; must be same length as mapdec. [radians]
mapdec (array) – Declinations whose angular distances to the single point will be calculated; must be same length as mapra. [radians]

Returns

Δθ⃗ – Angular distances between the single point provided and the arrays of points [radians].

Return type

array

Examples

Simple examples of distances to north pole:

>>> from math import pi
>>> import numpy as np
>>> dangle_rad(0, pi/2, np.array([0, 0, 0, 2, 3]),
...            np.array([pi/2, -pi/2, 0, 0, 0]))/pi
array([0. , 1. , 0.5, 0.5, 0.5])

hpmoc.utils.density_from_table(table, indices, nside, degrees=False)

Read probability density from table at the specified indices. Will try to read PROBDENSITY, and if it fails, will try to convert PROB to a probability density by dividing by area-per-pixel. Probability can be provided in a density format (usually for MOC skymaps) or a probability-per-pixel format (fixed-res skymaps).

Parameters

table (astropy.table.Table) – The table to read from. Must have either a PROBDENSITY or a PROB column.
indices (array-like) – Indices of the pixels (i.e. table rows) to read.
nside (int or array-like) – The HEALPix NSIDE parameter for the pixels. If each pixel has a different value, you can specify nside as an array of NSIDE values corresponding to the indices in indices.
degrees (bool, optional) – If True, return densities in inverse square degrees. Otherwise, return densities in inverse steradians.

Returns

values – A view onto the probability density values corresponding to indices.

Return type

np.array

Raises

ValueError – If table does not contain PROBDENSITY or PROB as columns or if any nside values are not valid.

hpmoc.utils.fill(u⃗, x⃗, ns, pad=None)

Rasterize a HEALPix multi-order skymap into a fixed-res full-sky HEALPix nested skymap, filling in missing values with a pad values.

Parameters

u⃗ (array-like) – NUNIQ pixel indices of the input skymap
x⃗ (array-like) – Pixel values of the input skymap
nˢ (int) – NSIDE of the output skymap
pad (int or float, optional) – Value to pad missing indices in the output skymap with. If not provided, use healpy.UNSEEN, which will render as blank space in healpy.visufunc plots.

See also

reraster: Similar, but for outputs that are also UNIQ-indexed.
render: Similar, but for outputs with repreating UNIQ indices.
nest_reres: For changing the resolution of NEST indices.

Returns: x⃗ⁿᵉˢᵗ – Fixed-res full-sky version of the input skymap in NEST ordering with missing values filled by pad.
Return type: np.ndarray

hpmoc.utils.group_slices(*u⃗, f=<function <lambda>>, fi=<function <lambda>>, include_empty=False, return_index=False, return_inverse=False, dtype=None): Group elements of u⃗ inputs using some sort of monotonic step function f: u⃗.dtype -> int codomain and a pseudo-inverse fⁱ mapping to the smallest element of the input domain giving that output value (both identity by default) and return a variety of views and slices into these groups. See nside_slices for documentation and an implementation that groups by HEALPix NSIDE order; this function is the same, but with o⃗ replaced by the result of f on elements of u⃗. You can use this function with the default grouping functions to group integers by value, e.g. for working with δo⃗ arrays from uniq_intersection.

See also

nside_slices

hpmoc.utils.handle_compressed_infile(func): Wrap read_partial_skymap so that it opens zipped fits files by temporarily decompressing them.

hpmoc.utils.interp_wcs(wcs: astropy.wcs.WCS, data: nptyping.types._ndarray.NDArray[typing.Any, Any], interp: Optional[Union[str, Tuple[int, Callable[[nptyping.types._ndarray.NDArray[typing.Any, float], nptyping.types._ndarray.NDArray[typing.Any, float], nptyping.types._ndarray.NDArray[typing.Any, Any]], nptyping.types._ndarray.NDArray[typing.Any, Any]]]]] = 'nearest') → Tuple[nptyping.types._ndarray.NDArray[typing.Any, int], nptyping.types._ndarray.NDArray[typing.Any, float]]

Interpolate data with coordinates specified by wcs FITS world coordinate system into a HEALPix NUNIQ skymap.

Parameters

wcs (astropy.wcs.WCS) – The world coordinate system defining pixel locations. If loading a FITS file as an HDU called hdu, you can get this argument as astropy.wcs.WCS(hdu.header). Note that you will need to manually include units for dimensionful quantities.
data (array-like) – The data corresponding to WCS. Available from an HDU as hdu.data.
interp (str or (int, func), optional) –
The interpolation strategy to use. Can be a string specifying one of the following pre-defined strategies:
- ”nearest” for nearest-neighbor
- ”bilinear” for bicubic
or else a tuple whose first element is the number of orders by which the pixels covering the WCS should have their resolution increased (“nearest” uses a value of 2, “bilinear” a value of 1; heuristically, a more sophisticated interpolation scheme can probably get away with 1), while the second element is a function taking the x, y coordinates of the pixels followed by the pixel values in data and returning the interpolated pixel values (which will form the return value s of this function).

Returns

u (array) – The corresponding NUNIQ HEALPix indices of the input skymap.
s (array-like) – The pixel-values of the input skymap interpolated at the locations of the pixels in u.

See also

interp_wcs_nn, hpmoc.partial.PartialUniqSkymap, astropy.wcs.WCS

hpmoc.utils.interp_wcs_nn(wcs: astropy.wcs.WCS, data: nptyping.types._ndarray.NDArray[typing.Any, Any]) → Tuple[nptyping.types._ndarray.NDArray[typing.Any, int], nptyping.types._ndarray.NDArray[typing.Any, float]]

Do a nearest-neighbor interpolation of data with coordinates specified by wcs FITS world coordinate system.

Parameters

wcs (astropy.wcs.WCS) – The world coordinate system defining pixel locations. If loading a FITS file as an HDU called hdu, you can get this argument as astropy.wcs.WCS(hdu.header). Note that you will need to manually include units for dimensionful quantities.
data (array-like) – The data corresponding to WCS. Available from an HDU as hdu.data.

Returns

u (array) – The corresponding NUNIQ HEALPix indices of the input skymap.
s (array-like) – The pixel-values of the input skymap interpolated at the locations of the pixels in u.

See also

interp_wcs, hpmoc.partial.PartialUniqSkymap, astropy.wcs.WCS

hpmoc.utils.is_gz(infile: Union[IO, str]): Check whether infile is a GZip file, in which case it should be unzipped.

hpmoc.utils.max_uint_type(largest)

hpmoc.utils.min_int_dtype(vmin, vmax): Find the smallest integer datatype that can represent a given range of values.

hpmoc.utils.monochrome_opacity_colormap(name, color): Get a monochrome matplotlib.colors.LinearSegmentedColormap with color defined by rgba (values between zero and one). Opacity will range from full transparency for the minimum value to the alpha value set in rgba.

hpmoc.utils.nest2ang(n⃗, N⃗s)

Get the angles corresponding to these nested indices n⃗ and NSIDE values N⃗ˢ.

Parameters

n⃗ (array-like) – HEALPix NESTED indices of pixels of interest.
N⃗ˢ (int or array-like) – The NSIDE values corresponding to each pixel. Can be a scalar if all pixels are at the same NSIDE.

Returns

ra_dec – 2D array whose first row is the right-ascension and second row is the declination (in degrees) of each pixel. You can get each of these individually with ra, dec = nest2ang(...).

Return type

astropy.units.Quantity

hpmoc.utils.nest2dangle(n⃗, ns, ra, dec, degrees=True, in_place=False)

Get the angular distance between the pixels defined in n⃗, nˢ and the point located at ra, dec.

Parameters

n⃗ (array-like) – HEALPix NEST indices
nˢ (int) – HEALPix NSIDE value
ra (float or astropy.units.Quantity) – Right ascension of the point; assumed degrees if no unit given.
dec (float or astropy.units.Quantity) – Declination of the point; assumed degrees if no unit given.
degrees (bool, optional) – If True, and assumed degrees if ra and/or dec are not astropy.units.Quantity instances with angular unit defined. If False, assume radians. Ignored if a unit is already specified.
in_place (bool, optional) – If True, store the result in n⃗ to reduce memory usage. Requires n⃗.dtype == np.float64.

Returns

Δθ⃗ – Angular distance in radians between each pixel in n⃗ and the point at ra, dec.

Return type

astropy.units.Quantity

Examples

The 12 base healpix pixels’ distances from north pole should all be equal to 90 minus their declinations:

>>> import numpy as np
>>> from hpmoc.healpy import healpy as hp
>>> Δθ⃗ = nest2dangle(range(12), 1, 32, 90).to('deg')
>>> Δθ⃗
<Quantity [ 48.1896851,  48.1896851,  48.1896851,  48.1896851,  90.       ,
            90.       ,  90.       ,  90.       , 131.8103149, 131.8103149,
           131.8103149, 131.8103149] deg>
>>> np.all(abs(Δθ⃗.value - (90 - hp.pix2ang(1, np.arange(12), nest=True,
...                                        lonlat=True)[1])) < 1e-13)
True

You can run the same check for larger skymaps, too (though note that precision drops for very nearby pixels due to the O(x²) behavior of cos() for small angles and the fact that a dot-product and arccos are used to compute the result):

>>> nside = 2**10
>>> nest = np.arange(12*nside**2)
>>> np.around(
...     (
...         nest2dangle(nest, nside, 32, 90).to('deg').value
...         - (90 - hp.pix2ang(nside, nest, nest=True, lonlat=True)[1])
...     ),
...     11
... ).ptp()
0.0

hpmoc.utils.nest2uniq(indices, nside, in_place=False)

Return the NUNIQ pixel indices for nested indices with NSIDE=``nside``.

Parameters

indices (array) – Indices of HEALPix pixels in NEST ordering.
nside (int) – A valid NSIDE value not greater than any of the NSIDE values of the indices.
in_place (bool, optional) – If True, perform the conversion in-place and return the converted indices object.

Returns

uniq – An array of HEALPix pixels in NEST ordering corresponding to the input indices and nside.

Return type

array

Raises

ValueError – If nside is an invalid NSIDE value, i.e. not a power of 2.

Examples

>>> import numpy as np
>>> nest2uniq(np.array([284, 286, 287,  10,   8,   2, 279, 278, 285]), 8)
array([540, 542, 543, 266, 264, 258, 535, 534, 541])

hpmoc.utils.nest_reres(nest, nside_in, nside_out)

Change the NSIDE of nest indices. If decreasing resolution, partially-filled pixels will be included but marked in full.

Parameters

nest (array-like) – HEALPix NEST indices.
nside_in (int) – The NSIDE of the provided HEALPix indices.
nside_out (int) – The NSIDE of the output indices.

Returns

reres_nest (array-like) – The smallest set of HEALPix NEST indices at NSIDE = nside_out fully covering nest.
full (bool or array) – If resolution is not decreased, equals True. If resolution is decreased, a boolean array that is True for all indices in reres_nest whose subpixels are all included in nest.

See also

fill: For converting skymap values (not indices) from UNIQ (not NEST) skymaps.

Examples

Double the resolution:

>>> n, full = nest_reres([0, 1], 2, 4)
>>> full
True
>>> print(n)
[0 1 2 3 4 5 6 7]

No effect:

>>> n, full = nest_reres([0, 1], 2, 2)
>>> full
True
>>> print(n)
[0 1]

Halve the resolution: >>> n, full = nest_reres([0, 1], 2, 1) >>> print(full) [False] >>> print(n) [0]

hpmoc.utils.nside2pixarea(nside, degrees=False)

Get the area per-pixel at nside. nside can also be a HEALPix array here.

Parameters

nside (int or array-like) – The NSIDE value or values you would like areas for.
degrees (bool, optional) – If True, return areas in inverse square degrees. Otherwise, return areas in inverse steradians.

Returns

pixarea – The area-per-pixel in the specified units. If nside was a scalar, pixarea is given as a scalar; if it was an array, pixarea is returned as an array whose values correspond to those in nside.

Return type

float or array-like

Raises

ValueError – If any nside values are not valid.

Examples

At NSIDE = 1, we should get 1/12 of the sky, or about 1/steradian:

>>> from math import pi
>>> allsky = 4*pi
>>> nside2pixarea(1) == allsky / 12
True

This should work for a list of NSIDES as well.

>>> import numpy as np
>>> nsides = np.array([1, 1, 2, 4])
>>> areas = np.array([allsky/12, allsky/12, allsky/48, allsky/192])
>>> np.all(nside2pixarea(nsides) == areas)
True

hpmoc.utils.nside_quantile_indices(nside, skymap, quantiles)

Get the indices and cumulative values of pixels falling between the given quantiles (expressed as a fraction of unity) for the skymap value skymap. Also known as confidence/credible region (CR) or percentiles (when using percentages rather than fractions of unity).

Pixels will be sorted by value (density), but quantiles will be taken by total integrated value (i.e. density times area, i.e. in proportion to density over the square of NSIDE). This puts the highest density regions in the upper quantiles while still normalizing for area in a multi-resolution skymap.

Note that this will work perfectly well for partial skymaps, though in that case (as one would expect) the quantiles will be taken with respect to the remaining pixels.

Parameters

nside (int or array) – NSIDE value (or values for a multi-order skymap) for the skymap.
skymap (array) – Skymap values. If nside is an array, values will correspond to those NSIDEs.
quantiles (array) – Quantiles from which to select pixels. Must be in ascending order with values in the interval [0, 1]. These will form endpoints for partitions of the skymap. For example, [0.1, 0.9] will omit the lowest and highest value pixels, giving the intermediate pixels accounting for 80% of the integrated skymap. Note that quantiles returned by this function are non-intersecting and half-open on the right (as with python indices), with the exception of 1 for the last pixel; for example, [0, 1] will include all pixels, [0.5, 1] will include the highest density pixels accounting for 50% of the integrated skymap value, [0, 0.5, 1] will partition the skymap into non-intersecting sets of pixels accounting for the high- and low-density partitions of the skymap by integrated value, etc.

Returns

indices (generator) – An iterator of arrays of indices suitable for selecting the values from skymap corresponing to the selected quantiles, sorted in order of increasing values of skymap. Use these to select the values from skymap corresponding to each of the partitions defined in quantiles. For example, quantiles=[0, 0.5, 1] will return a generator yielding two arrays of indices for accessing the values in the [0, 0.5] and [0.5, 1] quantiles, respectively. For applications where you want to just mask skymap and preserve sort order, you will want to sort this returned quantity before using it.
levels (astropy.units.Quantity or array) – Values of the skymap at each quantile. Useful for e.g. contour plots (though PartialUniqSkymap.plot will handle this automatically).
norm (array) – The total integral of skymap. The integrated region in a partition defined by quantiles [0.1, 0.2], for example, will be (0.2-0.1)*norm.

Raises

ValueError – If quantiles has length less than 2; if its values are not in order and contained in the interval [0, 1]; if nside and skymap cannot be broadcast together; if any values in skymap are negative; or if the total integrated skymap equals zero, in which case quantiles are undefined.

Examples

Get the pixel indices for pixels between the 0.1 (10%) and 0.9 (90%) quantiles on a fixed-resolution full skymap:

>>> import numpy as np
>>> skymap = np.array([ 9, 10, 11,  0,  1,  2,  3,  4,  5,  6,  7,  8])
>>> i, levels, norm = nside_quantile_indices(1, skymap, [0.1, 0.9])
>>> [ii.astype(int) for ii in i]
[array([ 7,  8,  9, 10, 11,  0,  1])]

The levels will be the values of the array at the 0.1 and 0.9 quantiles

>>> print(levels)
[ 4 11]

The norm in this case is just the average pixel value times the total solid angle of the sky:

>>> np.abs(4*np.pi*skymap.mean() - norm) < 1e-13
True

Get the pixel indices for the same interval on a mixed resolution partial skymap (where the last six pixels contribute far less area despite having high density):

>>> nside = np.array(6*[1]+6*[64])
>>> [*nside_quantile_indices(nside, skymap, [0.1, 0.9])[0]
...  ][0].astype(int)
array([0, 1])

Equal lower and upper bounds give empty quantiles:

>>> [*nside_quantile_indices(nside, skymap, [0.5, 0.5])[0]
...  ][0].astype(np.int64)
array([], dtype=int64)

Recover all indices (sorted by density):

>>> [*nside_quantile_indices(nside, skymap, [0, 1])[0]][0].astype(int)
array([ 3,  4,  5,  6,  7,  8,  9, 10, 11,  0,  1,  2])

Pick the 90% CR:

>>> [*nside_quantile_indices(nside, skymap, [0.1, 1])[0]][0].astype(int)
array([0, 1, 2])

Get the four top 20% quantiles:

>>> for q in nside_quantile_indices(nside, skymap,
...                                 np.arange(0.2, 1.1, 0.2))[0]:
...     print(q)
[0]
[]
[1]
[2]

hpmoc.utils.nside_slices(*u⃗, include_empty=False, return_index=False, return_inverse=False, dtype=None)

Sort and slice up a list of NUNIQ pixel index arrays, returning the sorted arrays as well as slice information for chunking them by NSIDE (pixel resolution), accessing the original array data, and the NSIDE orders of each chunk.

This is just a wrapper around group_slices using HEALPix NSIDE order as the grouping function.

Parameters

*u⃗ – np.array instances containing NUNIQ HEALPix indices
array-like – np.array instances containing NUNIQ HEALPix indices
include_empty (bool, optional) – If True, also include NSIDE orders not included in the input indices. Affects all return values.
return_index (bool, optional) – Whether to return u⃗̇. Only returned if True.
return_inverse (bool, optional) – Whether to return u⃗̇ˢ. Only returned if True.
dtype (int or numpy.dtype, optional) – If provided, cast the returned array to this data type. Useful for pre-allocating output arrays that only depend on spatial information.

Returns

u⃗ˢ (List[array]) – Sorted versions of each input array
s⃗ (List[List[slice]]) – Slices into each u⃗ˢ chunked by NSIDE order
o⃗ (array) – An array of HEALPix NSIDE orders included in the input indices
𝓁⃗ (List[array]) – The lengths of each slice in slice_starts
v⃗ (List[List[array]]) – Lists of array views into each u⃗ˢ corresponding to the slices given in slices
u⃗̇ (List[array], optional) – Indices into the original array that give u⃗ˢ
u⃗̇ˢ (List[Array], optional) – Indices into each u⃗ˢ that give the original arrays

See also

group_slices

Examples

>>> import numpy as np
>>> from pprint import pprint
>>> u⃗1 = np.array([1024, 4100, 1027, 263168, 263169, 1026, 44096])
>>> u⃗2 = np.array([4096, 4097, 1025, 16842752, 1026, 11024])
>>> us, s, o, l, v, ius = nside_slices(u⃗1, u⃗2, return_index=True)
>>> pprint([uu.astype(int) for uu in us])
[array([  1024,   1026,   1027,   4100,  44096, 263168, 263169]),
 array([    1025,     1026,     4096,     4097,    11024, 16842752])]
>>> pprint(s)
[[slice(0, 3, None),
  slice(3, 4, None),
  slice(4, 5, None),
  slice(5, 7, None),
  slice(7, 7, None)],
 [slice(0, 2, None),
  slice(2, 5, None),
  slice(5, 5, None),
  slice(5, 5, None),
  slice(5, 6, None)]]
>>> o.astype(int)
array([ 4,  5,  6,  8, 11])
>>> [ll.astype(int) for ll in l]
[array([3, 1, 1, 2, 0]), array([2, 3, 0, 0, 1])]
>>> for vv in v[0]:
...     print(vv)
[1024 1026 1027]
[4100]
[44096]
[263168 263169]
[]
>>> for vv in v[1]:
...     print(vv)
[1025 1026]
[ 4096  4097 11024]
[]
[]
[16842752]
>>> [ii.astype(int) for ii in ius]
[array([0, 5, 2, 1, 6, 3, 4]), array([2, 4, 0, 1, 5, 3])]

hpmoc.utils.outline_effect(): Get a matplotlib.patheffects.withStroke effect that outlines text nicely to improve plot readability.

hpmoc.utils.read_partial_skymap(infile: Union[IO, str], u⃗, memmap=True)

Read in pixels from a FITS skymap (or a gzip-compressed FITS skymap) that lie in a specific sky region. Attempts to minimize memory usage by memory-mapping pixels in the input file and only loading those specified in u⃗.

Parameters

infile (str or file) – A FITS HEALPix skymap file path or file object opened in binary read mode ‘rb’ (optionally compressed; see note under memmap)
u⃗ (array-like) – HEALPix pixel indices in NUNIQ ordering specifying the region of the skymap that should be loaded.
memmap (bool, optional) – If True, use memory mapping when reading in files. This can VASTLY reduce memory required for high-resolution skymaps. If infile is gzip-compressed and memmap is True, then infile will be decompressed to a temporary file and data will be read from it (necessary to constrain memory usage); for high-resolution skymaps, this can require the availability of several gigabytes of tempfile storage. You will need to make use of TmpGunzipFits when working with zipped files in order to be able to use memmap=True.

Returns

partial_skymap – A partial skymap table in nested ordering. Has two columns: UNIQ and PROBDENSITY. If the resolution of the original HEALPix skymap file is lower than that of the u⃗, then any pixels overlapping with those in u⃗ will be used; this might result in a larger portion of the skymap being used than that specified in u⃗. The resolution of this skymap will be the resolution of the smallest pixel loaded from the input file (in the case of ring or nested ordering, this is just the resolution of the input skymap).

Return type

astropy.table.Table

See also

uniq_minimize, TmpGunzipFits

hpmoc.utils.render(u⃗, x⃗, u⃗o, pad=None, valid=None, mask_missing=False, Ii⃗i⃗o=None)

Like reraster, but allows you to map to a partially-covered u⃗ᵒ skymap, e.g. for rendering a plot, thanks to a call to np.unique(u⃗ᵒ, return_inverse=True) wrapping the whole thing (to take care of scattering values to repeated pixels).

Parameters

u⃗ (array) – The indices of the skymap.
x⃗ (array) – The values of the skymap.
u⃗ᵒ (array or astropy.wcs.WCS) – If u⃗ᵒ is an astropy.wcs.WCS world coordinate system, then wcs2mask_and_uniq will be used to get the indices. Non-valid pixels (i.e. pixels outside the projection area) will take on np.nan values, while valid pixels will be rendered as usual.
pad (float, optional) – Pad value for missing pixels. If not provided, will raise an error if missing parts of the skymap fall in u⃗ᵒ. To render a healpy plot with missing pixels, pass pad=healpy.UNSEEN.
valid (array, optional) – If provided, results will be scattered into an array of the same shape as valid, filling the indices where valid==True. The number of True values in valid must therefore equal the length of u⃗ᵒ. This argument only makes sense if u⃗ᵒ is an array of NUNIQ indices; if it is a WCS instance and valid is provided, an error is raised. Use valid to produce plots or to reuse indices produced by wcs2mask_and_uniq in several render invocations. See note on how mask_missing affects the result.
mask_missing (bool) – If mask_missing=True, return a np.ma.core.MaskedArray. Missing values are tolerated and are marked as True in the mask_missing. They will be set to pad or 0 in the data field. If valid is also provided, then the output will still be a np.ma.core.MaskedArray, but will be set to True wherever valid == False in addition to wherever pixels are missing (and will still take on masked values of np.nan in the invalid regions).
Iᵢ⃗ⁱ⃗ᵒ (Tuple[array, array, array]) – Return tuple of uniq_intersection. Use this to save time in repeated invocations.

Returns

s⃗ₒ – The pixel values at locations specified by u⃗ᵒ. If mask_missing=True, will be a np.ma.core.MaskedArray set to True at the missing values in the valid field with missing data field values set to pad or None.

Return type

array-like

Raises

ValueError – If u⃗ᵒ is a WCS instance and valid is not None.

hpmoc.utils.reraster(u⃗, x⃗, u⃗o, pad=None, mask_missing=False, Ii⃗i⃗o=None)

Rasterize skymap pixel values x⃗ with NUNIQ indices u⃗ to match pixels u⃗ᵒ, discarding sky areas excluded by u⃗ᵒ and (optionally) padding missing values with pad.

Parameters

u⃗ (array-like) – NUNIQ indices of the skymap.
x⃗ (array-like) – Pixel values. Must be the same length as u⃗.
u⃗ᵒ (array-like) – NUNIQ indices of the output skymap.
pad (float or int, optional) – A pad value to use for pixels missing from the input skymap. Only used if u⃗ does not fully cover u⃗ᵒ. Use healpy.UNSEEN for this value if you want to mark pixels as not-observed for HEALPy plots etc.
mask_missing (bool) – If mask_missing=True, return a np.ma.core.MaskedArray. Missing values are tolerated and are marked as True in the mask_missing. They will be set to pad or 0 in the data field.
Iᵢ⃗ⁱ⃗ᵒ (Tuple[np.ndarray, np.ndarray, np.ndarray], optional) – If you’ve already computed uniq_intersection(u⃗, u⃗ᵒ), you can pass it as this argument to avoid recomputing it. No checks will be made for correctness if provided.

Returns

x⃗ᵒ – Pixel values of the rasterized skymap corresponding to the indices given in u⃗ᵒ. x⃗ᵒ values are pixel-area-weighted averages of the input pixel values, even if some pixels in u⃗ᵒ are not fully covered by pixels from u⃗. Any parts of the sky defined in u⃗ that are not covered by u⃗ᵒ are omitted, so this function can also be used to mask a skymap in a single step. If mask_missing=True, is a np.ma.core.MaskedArray.

Return type

array-like

Raises

ValueError – If pad is not provided but u⃗ does not cover all pixels in u⃗ᵒ.

See also

render, hpmoc.partial.PartialUniqSkymap.reraster, np.ma.core.MaskedArray

Examples

Create a small partial skymap with example pixel values:

>>> import numpy as np
>>> u⃗ = np.array([1024, 4100, 1027, 1026, 44096])
>>> x⃗ = np.array([1.,   2.,   3.,   4.,   5.])

We will rerasterize this skymap to these sky areas:

>>> u⃗ᵒ = np.array([4096, 4097, 1025, 1026, 11024])
>>> reraster(u⃗, x⃗, u⃗ᵒ)
array([1., 1., 2., 4., 5.])

The third pixel in u⃗ is not present in u⃗ᵒ, so we will need to provide a default pad value for it when rasterizing in the other direction. Note that the first pixel of the result is the average of the first and second pixels in the input map, since both of these have equal area and overlap with the first pixel:

>>> reraster(u⃗ᵒ, x⃗, u⃗, pad=0.)
array([1.5, 3. , 0. , 4. , 5. ])

We can also simply mask that value by passing mask_missing=True, in which case the result will be a np.ma.core.MaskedArray which is True for values which were missing (the missing/masked values themselves will be set to zero or pad if provided):

>>> m = reraster(u⃗ᵒ, x⃗, u⃗, mask_missing=True)
>>> print(m.data)
[1.5 3.  0.  4.  5. ]
>>> print(m.mask)
[False False  True False False]

Note that the values are averages of input pixels; in cases where only one input pixel is sampled from, the value remains unchanged. This makes reraster good for working with densities and other intensive spatial values; extensive values should have their pixel areas divided out before being rasterized.

If you’ve already got the uniq_intersection of u⃗ and u⃗ᵒ from a previous calculation, you can avoid recomputing it during rasterization by passing it as the Iᵢ⃗ⁱ⃗ᵒ argument, though beware it will not be checked for correctness:

>>> reraster(u⃗, x⃗, u⃗ᵒ, Iᵢ⃗ⁱ⃗ᵒ=uniq_intersection(u⃗, u⃗ᵒ))
array([1., 1., 2., 4., 5.])

hpmoc.utils.resol2nside(res, coarse=False, degrees=True)

Get the HEALPix NSIDE value corresponding to a given resolution.

Parameters

res (float or array-like) – The required resolution, measured as the angular width of a pixel.
coarse (bool, optional) – If False, pick a higher resolution than the one specified. If True, pick a lower resolution (necessary since HEALPix resolutions increment in discrete steps).
degrees (bool, optional) – Whether res is specified in degrees. If False, radians are assumed.

hpmoc.utils.set_partial_skymap_metadata(meta, mask, caller): Write metadata to a partial skymap.

hpmoc.utils.sky_area(): Get the area of the entire sky as an Astropy.units.Quantity.

hpmoc.utils.sky_area_deg(): Get the area of the entire sky in square degrees.

hpmoc.utils.uniq2dangle(u⃗, ra, dec, degrees=True)

Like nest2dangle, but takes HEALPix NUNIQ indices as input u⃗.

Examples

The 12 base healpix pixels’ distances from north pole should all be equal to 90 minus their declinations:

>>> import numpy as np
>>> from hpmoc.healpy import healpy as hp
>>> Δθ⃗ = uniq2dangle(range(4, 16), 32, 90).to('deg')
>>> Δθ⃗
<Quantity [ 48.1896851,  48.1896851,  48.1896851,  48.1896851,  90.       ,
            90.       ,  90.       ,  90.       , 131.8103149, 131.8103149,
           131.8103149, 131.8103149] deg>
>>> np.all(abs(Δθ⃗.value - (90 - hp.pix2ang(1, np.arange(12), nest=True,
...                                    lonlat=True)[1])) < 1e-13)
True

See also

nest2dangle

hpmoc.utils.uniq2nest(u⃗, ns, nest=True)

Take a set of HEALPix NUNIQ-ordered indices at arbitrary resolution covering an arbitrary portion of the sky and convert them to non-overlapping pixels at a fixed NSIDE (resolution), returning the indices of the resulting skymap in NUNIQ or NEST ordering.

Parameters

u⃗ (array) – Indices of HEALPix pixels in NUNIQ ordering.
nˢ (int) – HEALPix NSIDE value of the output map.
nest (bool, optional) – Whether to return the fixed-resolution indices in NEST ordering. If False, leave them in NUNIQ ordering (though they will still be at the fixed resolution specified as nˢ).

Returns

u⃗ᵒ – Indices covering the same sky region as u⃗ (possibly a larger region if resolution is reduced) at resolution nˢ in either NEST order (if nest is True) or NUNIQ order (if nest is False).

Return type

array

Raises

ValueError – If nˢ is not a valid NSIDE value or u⃗ are not valid NUNIQ indices, or if the requested resolution is too high too represent with int64.

Examples

Let’s convert this NUNIQ sky region to an NSIDE=32 nested indices (this will split the first, third, and fourth pixels, which are coarser–and select the pixel containing the last pixel, which is smaller than–than the target pixel size):

>>> import numpy as np
>>> u⃗ = np.array([1024, 4100, 1027, 1026, 44096])
>>> uniq2nest(u⃗, 32).astype(int)
array([   0,    1,    2,    3,    4,    8,    9,   10,   11,   12,   13,
         14,   15, 6928])

Same pixel indices, but keep them in NUNIQ format:

>>> uniq2nest(u⃗, 32, nest=False).astype(int)
array([ 4096,  4097,  4098,  4099,  4100,  4104,  4105,  4106,  4107,
        4108,  4109,  4110,  4111, 11024])

Coarsen the pixels to NSIDE=16:

>>> uniq2nest(u⃗, 16).astype(int)
array([   0,    1,    2,    3, 1732])
>>> uniq2nest(u⃗, 16, False).astype(int)
array([1024, 1025, 1026, 1027, 2756])

Increase resolution of all pixels to NSIDE=64:

>>> uniq2nest(u⃗, 64).astype(int)
array([    0,     1,     2,     3,     4,     5,     6,     7,     8,
           9,    10,    11,    12,    13,    14,    15,    16,    17,
          18,    19,    32,    33,    34,    35,    36,    37,    38,
          39,    40,    41,    42,    43,    44,    45,    46,    47,
          48,    49,    50,    51,    52,    53,    54,    55,    56,
          57,    58,    59,    60,    61,    62,    63, 27712])

hpmoc.utils.uniq2nest_and_nside(indices, in_place=False)

Parameters

indices (array) – A scalar or numpy array of NUNIQ indices
in_place (bool, optional) – If True, perform the conversion in-place and return the converted indices object along with the calculated nside.

Returns

nest_indices (array) – The indices in nest ordering at their respective resolutions
nside (array) – The resolution expressed as NSIDE of the indices given in nest_indices

Examples

>>> import numpy as np
>>> from pprint import pprint
>>> nuniq = np.array([540, 542, 543, 266, 264, 258, 535, 534, 541])
>>> pprint([u.astype(int) for u in uniq2nest_and_nside(nuniq)])
[array([284, 286, 287,  10,   8,   2, 279, 278, 285]),
 array([8, 8, 8, 8, 8, 8, 8, 8, 8])]

Confirm that this is the inverse of nest2uniq:

>>> all(nest2uniq(*uniq2nest_and_nside(nuniq)) == nuniq)
True

hpmoc.utils.uniq2nside(indices)

Get the NSIDE value of the given NUNIQ-ordered indices.

Raises: ValueError – If indices are not valid NUNIQ indices, i.e. they are not integers greater than 3.

hpmoc.utils.uniq2order(indices)

Get the HEALPix order of the given NUNIQ-ordered indices.

Raises: ValueError – If indices are not valid NUNIQ indices, i.e. they are not integers greater than 3.

hpmoc.utils.uniq2xyf_nside(u⃗)

Examples

>>> uniq2xyf_nside(8)
(0, 0, 4, 1)

See also

xyf_nside2uniq

hpmoc.utils.uniq_coarsen(u, orders)

Coarsen the pixel indices in u to reduce storage and computation requirements.

Parameters

u (array) – UNIQ indices to coarsen.
orders (int) – How many times the resolution of the smallest pixels will be halved.

Returns

uc – Unique coarsened pixel values in ascending order. All pixels will have a HEALPix order capped at the maximum order present in u minus orders, unless this value is negative, in which case the output will only consist of base pixels. If pixels in u overlap, it is possible that there will also be overlapping pixels in the output; no check is made for this.

Return type

array

Raises

ValueError – If orders < 0.

Examples

Pixels 353, 354, and 355 lie within pixels 88, 22, and 5; pixels 80 and 81 lie within pixels 20 and 5; pixel 21 lies within pixel 5.

>>> u = [4, 21, 80, 81, 353, 354, 355]

Coarsening by zero will have no effect:

>>> print(uniq_coarsen(u, 0))
[  4  21  80  81 353 354 355]

Coarsening by one will only combine the very smallest pixels:

>>> print(uniq_coarsen(u, 1))
[ 4 21 80 81 88]

Coarsening by larger numbers will combine so many higher orders:

>>> print(uniq_coarsen(u, 2))
[ 4 20 21 22]
>>> print(uniq_coarsen(u, 3))
[4 5]

Coarsening by a value greater than the largest order will have no further effect, since the 12 base pixels cannot be combined:

>>> print(uniq_coarsen([4, 21, 80, 81, 353, 354, 355], 4))
[4 5]

hpmoc.utils.uniq_diadic(Ω, u⃗i, x⃗i, pad=None, coarse=True)

Apply a diadic function Ω(x⃗ᶠ1, x⃗ᶠ2) -> y⃗ᶠ that operates on skymap pixel values of the same resolution to skymaps with arbitrary pixelization/resolution schemes and pixel orders, returning the indices and pixel values of the resulting skymap. Useful for binary operations between arbitrary partial skymaps.

Parameters

Ω (Callable[[np.ndarray, np.ndarray], np.ndarray]) – A binary function operating on two sets of skymap pixel values corresponding elementwise to the same sky locations. Must be a skymap-resolution independent operation for the results to make sense.
u⃗ⁱ (Tuple[np.ndarray, np.ndarray]) – NUNIQ indices of the two skymaps to be passed to Ω.
x⃗ⁱ (List[np.ndarray, np.ndarray]) – Pixel values (corresponding to the locations specified by u⃗ⁱ) of the skymaps to be passed to Ω. Must have same lengths as the arrays in u⃗ⁱ.
pad (float or int, optional) – A pad value to use for parts of the sky contained in either u⃗ⁱ[0] or u⃗ⁱ[1] but not in both (since Ω will be undefined in these regions). If not provided, the returned skymap will only contain the intersection of the sky areas defined in u⃗ⁱ.
coarse (bool, optional) – If True, for sky areas where u⃗ⁱ[0] and u⃗ⁱ[1] have different resolutions, pick the lower resolution for the output map (using pixel-area-weighted averages to decimate the higher-res regions). This produces shorter output arrays. If False, split coarse pixels into the higher resolution of those specified in u⃗ⁱ[0] and u⃗ⁱ[1] for a given sky region; use this if you need to maintain resolution for subsequent calculations, but be aware that this may impact performance without improving accuracy, e.g. if you’re planning to integrate the result of this operation. This can also be useful if you need to cover the exact area defined by the input skymaps.

Returns

u⃗ʸ (np.ndarray) – Sorted NUNIQ indices of the result of Ω. In general, will be different from both u⃗ⁱ inputs.
y⃗ (np.ndarray) – Pixel values of the result of Ω corresponding to indices u⃗ʸ.

Examples

Take the product of two heterogeneous multi-order skymaps:

>>> from pprint import pprint
>>> from operator import mul
>>> import numpy as np
>>> u⃗1 = np.array([1024, 4100, 1027, 1026, 44096])
>>> x⃗1 = np.array([1.,   2.,   3.,   4.,   5.])
>>> u⃗2 = np.array([4096, 4097, 1025, 1026, 11024])
>>> x⃗2 = np.array([0.,   10.,  1.,   100., 1000.])
>>> pprint(uniq_diadic(mul, (u⃗1, u⃗2), (x⃗1, x⃗2)), width=60)
(array([ 1024,  1025,  1026, 11024]),
 array([5.e+00, 2.e+00, 4.e+02, 5.e+03]))

Provide a default pad value for indices non-overlapping parts of the input skymaps:

>>> pprint(uniq_diadic(mul, (u⃗1, u⃗2), (x⃗1, x⃗2), pad=0.), width=60)
(array([ 1024,  1025,  1026,  1027, 11024]),
 array([5.e+00, 2.e+00, 4.e+02, 0.e+00, 5.e+03]))

Increase resolution as necessary (do not combine pixels):

>>> pprint(uniq_diadic(mul, (u⃗1, u⃗2), (x⃗1, x⃗2), coarse=False), width=60)
(array([ 1026,  4096,  4097,  4100, 44096]),
 array([4.e+02, 0.e+00, 1.e+01, 2.e+00, 5.e+03]))

hpmoc.utils.uniq_intersection(u⃗1, u⃗2)

Downselect the pixel indices given in u⃗1 to the set that overlaps with pixels in u⃗2 and return pairs of indices into both of the input index lists showing which pixel in u⃗2 each downselected pixel from u⃗1 overlaps with. Use this to get rid of pixels outside of a given region, or alternatively use it as part of a calculation involving two multi-resolution skymaps whose pixel sizes are non-identical. Written to perform efficiently on arbitrary inputs with O(MlogN+NlogM) performance (where M, N are the lengths of the input arrays).

Parameters

u⃗1 (array) – Indices of HEALPix pixels in NUNIQ ordering. Pixels corresponding to these indices MUST NOT OVERLAP.
u⃗2 (array) – Indices of HEALPix pixels in NUNIQ ordering. Pixels corresponding to these indices MUST NOT OVERLAP.

Returns

u⃗̇1 (array) – Indices into u⃗1 that overlap with u⃗2.
u⃗̇2 (array) – Corresponding indices into u⃗2 that overlap with u⃗1.
δo⃗ (array) – Corresponding differences in order between the indices, e.g. if the first entry of u⃗̇1 has NSIDE 16 (order 4) and the corresponding entry of u⃗̇2 has NSIDE 1024 (order 10), then the first entry of δo⃗ will be (10 - 4) = 6.

Raises

ValueError – If either u⃗1 or u⃗2 contain indices referring to overlapping pixels (note that this may happen even if the inputs do not contain repeating values since different pixel sizes can overlap).

Examples

Some pixels with NSIDE of 16, 32, and 64, respectively:

>>> from pprint import pprint
>>> import numpy as np
>>> u⃗1 = np.array([1024, 4100, 1027, 1026, 44096])

Pixels at NSIDE = 32 that overlap with only the first and last pixels of u⃗1:

>>> u⃗2 = np.array([4096, 4097, 1025, 1026, 11024])

We should see correspondence between index 0 of u⃗1 and indices 0, 1 of u⃗1; and correspondence between index 2 of u⃗1 and index 2 of u⃗2:

>>> pprint(tuple(a.astype(int) for a in uniq_intersection(u⃗1, u⃗2)),
...        width=60)
(array([4, 3, 0, 0, 1]),
 array([4, 3, 0, 1, 2]),
 array([-1,  0,  1,  1, -1]))

hpmoc.utils.uniq_minimize(u, *x, test=<built-in function eq>, combine=<function <lambda>>)

Take a set of HEALPix NUNIQ indices u⃗ (and, optionally, pixel values x⃗) and find the shortest equivalent multi-order pixelation by combining pixels. If x⃗ is provided, only combine pixels whose values are equal. This can also be used if a canonical pixelization is needed for a given mask or skymap.

Parameters

u (array-like) – HEALPix NUNIQ indices of the skymap in question.
x (array-like, optional) – Pixel values of the array. If included, sub-pixels will only be combined if their pixel values are equal accordint to test for all x values provided.
test (func, optional) – An equality test for determining whether adjacent pixels can be combined. Defaults to a standard equality check, operator.eq. Override this if you want, e.g., approximately equal floating point values to be combined, small values rounded to zero, etc.
combine (func, optional) – A function for combining pixels. Expects an argument x, of the form of one of the arrays passed in for x, as well as a boolean mask of pixels i which select the first pixel of each four to be combined. By default, simply selects this first pixel; you could alternatively provide a function which, e.g., averages the combined pixels.

Returns

u⃗ᵐ (numpy.ndarray) – The shortest equivalent NUNIQ indexing that can describe u.
*x⃗ᵐ (numpy.ndarray, optional) – Corresponding pixel values in x, combined according to combined.

Examples

Make a set of pixels corresponding to the first four base pixels as well as the first pixel at NSIDE = 2 lying in the fifth base pixel:

>>> import numpy as np
>>> u = np.concatenate([nest2uniq(np.arange(2), 1),
...                     nest2uniq(np.arange(8, 17), 2)])
>>> print(u)
[ 4  5 24 25 26 27 28 29 30 31 32]

UNIQ indices 24-31 cover the same area as 6-7; uniq_minimize will detect this:

>>> um, = uniq_minimize(u)
>>> print(um)
[ 4  5  6  7 32]

This makes no difference against a constant skymap, with, e.g., values of 1 everywhere:

>>> um, xm = uniq_minimize(u, np.full_like(u, 1))
>>> print(um)
[ 4  5  6  7 32]
>>> print(xm)
[1 1 1 1 1]

If, however, the 4 pixels in the range 28-31 do not have equal values, they will not be combined with the default choice of test:

>>> um, xm = uniq_minimize(u, np.array([1, 2, 3, 3, 3, 3, 4, 5, 6, 7, 8]))
>>> print(um)
[ 4  5  6 28 29 30 31 32]
>>> print(xm)
[1 2 3 4 5 6 7 8]

This can be very effective for combining pixels in skymaps made by algorithms like CWB, which often produce adjacent pixels with zero probability, or BAYESTAR, which is natively MOC but is often distributed at fixed resolution. For example, an NSIDE=128 pixel in a BAYESTAR skymap whose smallest pixel has NSIDE=1024 will be split into 64 NSIDE=1024 subpixels to bring the entire skymap to a single resolution without loss of precision; we can represent such a single pixel using its NUNIQ indices u2 and a constant skymap value x2:

>>> u2 = np.arange(4194304, 4194368)
>>> x2 = np.ones_like(u2)
>>> um2, xm2 = uniq_minimize(u2, x2)

The pixel indices will be combined to give the NUNIQ index of the original pixel at NSIDE=128:

>>> print(um2)
[65536]

The pixel values will remain the same:

>>> print(xm2)
[1]

hpmoc.utils.uniq_sample(uniq, skymap, quantiles)

Do a probability-weighted random sample of locations in a skymap based on quantile. Follows ascending density order, just like nside_quantile_indices.

Parameters

uniq (array-like) – NUNIQ indices of the skymap to be sampled.
skymap (array-like) – Probability density at each of the pixels in uniq.
quantiles (array-like) – Values in the interval [0, 1] specifying which quantile to select from the skymap. For example, [0.5, 0.9] will select locations from within the pixels for which the integral of all higher-probability-density regions are 0.5 and 0.1, respectively.

Returns

locations – Max-resolution (order = 30) NUNIQ indices drawn from the PDF defined by uniq and skymap at the specified quantiles assuming that each pixel in skymap has a uniform PDF within its borders. The difference between successive values in the density-ordered CDF of skymap will be used to select subpixels in NEST ordering from the pixel in uniq with the greatest CDF value lower than the corresponding quantile, allowing for reproducible results. The returned values have the same order as the input quantiles. You can use healpy.pix2ang(*uniq2nest_and_nside(locations)[::-1], nest=True) to recover the angles corresponding to the recovered locations.

Return type

array

See also

nside_quantile_indices, uniq2nest_and_nside, healpy.pix2ang

hpmoc.utils.wcs2ang(wcs: astropy.wcs.WCS, lonlat=True)

Convert an astropy.wcs.WCS world coordinate system’s pixels into ICRS coordinate angles.

Parameters

wcs (astropy.wcs.WCS) – The world coordinate system for whose pixels you want coordinate values.
lonlat (bool, optional) – If True, return right ascension/declination. If False, return (phi, theta) angles.

Returns

valid (array) – Boolean mask indicating which values from the WCS are valid. The rest can be padded with a fill value by the user (most likely np.nan.
ra_or_theta (astropy.units.Quantity) – The right-ascension/longitude if lonlat=True, otherwise the zenith/theta angle.
dec_or_phi (astropy.units.Quantity) – The declination/latitude angle if lonlat=True, otherwise the azimuthal/phi angle.

hpmoc.utils.wcs2mask_and_uniq(wcs): Convert an astropy.wcs.WCS world coordinate system’s pixels into NUNIQ indices for HEALPix pixels of approximately the same size.

hpmoc.utils.wcs2nest(wcs, nside=None, order_delta=None)

Get NEST pixels at nside resolution covering an astropy.wcs.WCS instance’s pixels as well as their x and y pixel coordinates. All returned pixels and coordinates will be within the image boundaries and the indices will be non-repeating. If order_delta is provided, then the NEST resolution is doubled the number of times specified thereby. Can only pass one of nside or order_delta or else a ValueError is raised.

Returns

nside (int) – The NSIDE of the output indices.
nest (NDArray[(Any,), int]) – The HEALPix NEST indices.
x (NDArray[(Any,), float]) – The pixel-space x-coordinates of the points in nest.
y (NDArray[(Any,), float]) – The pixel-space y-coordinates of the points in nest.

hpmoc.utils.wcs2resol(wcs): Get the resolution of an astropy.wcs.WCS coordinate system, i.e. the smallest inter-pixel distance, as an astropy.units.Quantity with angular unit.

hpmoc.utils.xyf_nside2uniq(x, y, f, ns)

Examples

>>> import numpy as np
>>> (xyf_nside2uniq(*uniq2xyf_nside(np.arange(4, 10000))) == np.arange(4, 10000)).all()
True

See also

uniq2xyf_nside