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 aNamedTemporaryFile
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 inmapra
,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 specifiedindices
. 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 aPROB
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 inindices
.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 containPROBDENSITY
orPROB
as columns or if anynside
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 inhealpy.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 functionf: u⃗.dtype -> int
codomain and a pseudo-inversefⁱ
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. Seenside_slices
for documentation and an implementation that groups by HEALPix NSIDE order; this function is the same, but witho⃗
replaced by the result off
on elements ofu⃗
. You can use this function with the default grouping functions to group integers by value, e.g. for working withδo⃗
arrays fromuniq_intersection
.See also
- 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 bywcs
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 asastropy.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 ashdu.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 indata
and returning the interpolated pixel values (which will form the return values
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 bywcs
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 asastropy.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 ashdu.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 byrgba
(values between zero and one). Opacity will range from full transparency for the minimum value to the alpha value set inrgba
.
- hpmoc.utils.nest2ang(n⃗, N⃗s)
Get the angles corresponding to these nested indices
n⃗
and NSIDE valuesN⃗ˢ
.- 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 atra, 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 ifra
and/ordec
are notastropy.units.Quantity
instances with angular unit defined. IfFalse
, assume radians. Ignored if a unit is already specified.in_place (bool, optional) – If
True
, store the result inn⃗
to reduce memory usage. Requiresn⃗.dtype == np.float64
.
- Returns
Δθ⃗ – Angular distance in radians between each pixel in
n⃗
and the point atra, 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 convertedindices
object.
- Returns
uniq – An array of HEALPix pixels in NEST ordering corresponding to the input
indices
andnside
.- 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 coveringnest
.full (bool or array) – If resolution is not decreased, equals
True
. If resolution is decreased, a boolean array that isTrue
for all indices inreres_nest
whose subpixels are all included innest
.
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 innside
.- 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 valueskymap
. 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 theskymap
. 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 of1
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 ofskymap
. Use these to select the values fromskymap
corresponding to each of the partitions defined inquantiles
. 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 maskskymap
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]
; ifnside
andskymap
cannot be broadcast together; if any values inskymap
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 indicesarray-like –
np.array
instances containing NUNIQ HEALPix indicesinclude_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 ifTrue
.return_inverse (bool, optional) – Whether to return
u⃗̇ˢ
. Only returned ifTrue
.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 ordero⃗ (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 inslices
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
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. Ifinfile
is gzip-compressed andmemmap
isTrue
, theninfile
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 ofTmpGunzipFits
when working with zipped files in order to be able to usememmap=True
.
- Returns
partial_skymap – A partial skymap table in
nested
ordering. Has two columns:UNIQ
andPROBDENSITY
. If the resolution of the original HEALPix skymap file is lower than that of the u⃗, then any pixels overlapping with those inu⃗
will be used; this might result in a larger portion of the skymap being used than that specified inu⃗
. The resolution of this skymap will be the resolution of the smallest pixel loaded from the input file (in the case ofring
ornested
ordering, this is just the resolution of the input skymap).- Return type
astropy.table.Table
See also
- 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-coveredu⃗ᵒ
skymap, e.g. for rendering a plot, thanks to a call tonp.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 anastropy.wcs.WCS
world coordinate system, thenwcs2mask_and_uniq
will be used to get the indices. Non-valid pixels (i.e. pixels outside the projection area) will take onnp.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 ahealpy
plot with missing pixels, passpad=healpy.UNSEEN
.valid (array, optional) – If provided, results will be scattered into an array of the same shape as
valid
, filling the indices wherevalid==True
. The number ofTrue
values invalid
must therefore equal the length ofu⃗ᵒ
. This argument only makes sense ifu⃗ᵒ
is an array of NUNIQ indices; if it is aWCS
instance andvalid
is provided, an error is raised. Usevalid
to produce plots or to reuse indices produced bywcs2mask_and_uniq
in severalrender
invocations. See note on howmask_missing
affects the result.mask_missing (bool) – If
mask_missing=True
, return anp.ma.core.MaskedArray
. Missing values are tolerated and are marked asTrue
in themask_missing
. They will be set topad or 0
in thedata
field. Ifvalid
is also provided, then the output will still be anp.ma.core.MaskedArray
, but will be set toTrue
wherevervalid == False
in addition to wherever pixels are missing (and will still take on masked values ofnp.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 anp.ma.core.MaskedArray
set toTrue
at the missing values in thevalid
field with missingdata
field values set topad or None
.- Return type
array-like
- Raises
ValueError – If
u⃗ᵒ
is aWCS
instance andvalid
is notNone
.
See also
reraster
,hpmoc.partial.PartialUniqSkymap.render
,hpmoc.points.PointsTuple.render
,np.ma.core.MaskedArray
- hpmoc.utils.reraster(u⃗, x⃗, u⃗o, pad=None, mask_missing=False, Ii⃗i⃗o=None)
Rasterize skymap pixel values
x⃗
with NUNIQ indicesu⃗
to match pixelsu⃗ᵒ
, discarding sky areas excluded byu⃗ᵒ
and (optionally) padding missing values withpad
.- 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 coveru⃗ᵒ
. Usehealpy.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 anp.ma.core.MaskedArray
. Missing values are tolerated and are marked asTrue
in themask_missing
. They will be set topad or 0
in thedata
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 inu⃗ᵒ
are not fully covered by pixels fromu⃗
. Any parts of the sky defined inu⃗
that are not covered byu⃗ᵒ
are omitted, so this function can also be used to mask a skymap in a single step. Ifmask_missing=True
, is anp.ma.core.MaskedArray
.- Return type
array-like
- Raises
ValueError – If
pad
is not provided butu⃗
does not cover all pixels inu⃗ᵒ
.
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 inu⃗ᵒ
, 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 anp.ma.core.MaskedArray
which isTrue
for values which were missing (the missing/masked values themselves will be set to zero orpad
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
ofu⃗
andu⃗ᵒ
from a previous calculation, you can avoid recomputing it during rasterization by passing it as theIᵢ⃗ⁱ⃗ᵒ
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. IfTrue
, pick a lower resolution (necessary since HEALPix resolutions increment in discrete steps).degrees (bool, optional) – Whether
res
is specified in degrees. IfFalse
, 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 inputu⃗
.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
- 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 asnˢ
).
- Returns
u⃗ᵒ – Indices covering the same sky region as
u⃗
(possibly a larger region if resolution is reduced) at resolutionnˢ
in either NEST order (ifnest
isTrue
) or NUNIQ order (ifnest
isFalse
).- Return type
array
- Raises
ValueError – If
nˢ
is not a valid NSIDE value oru⃗
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 convertedindices
object along with the calculatednside
.
- 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
- 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
minusorders
, unless this value is negative, in which case the output will only consist of base pixels. If pixels inu
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 inu⃗ⁱ
.pad (float or int, optional) – A pad value to use for parts of the sky contained in either
u⃗ⁱ[0]
oru⃗ⁱ[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 inu⃗ⁱ
.coarse (bool, optional) – If
True
, for sky areas whereu⃗ⁱ[0]
andu⃗ⁱ[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. IfFalse
, split coarse pixels into the higher resolution of those specified inu⃗ⁱ[0]
andu⃗ⁱ[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 bothu⃗ⁱ
inputs.y⃗ (np.ndarray) – Pixel values of the result of
Ω
corresponding to indicesu⃗ʸ
.
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 inu⃗2
and return pairs of indices into both of the input index lists showing which pixel inu⃗2
each downselected pixel fromu⃗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 withu⃗2
.u⃗̇2 (array) – Corresponding indices into
u⃗2
that overlap withu⃗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 ofu⃗̇2
has NSIDE 1024 (order 10), then the first entry ofδo⃗
will be (10 - 4) = 6.
- Raises
ValueError – If either
u⃗1
oru⃗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 ofu⃗1
; and correspondence between index 2 ofu⃗1
and index 2 ofu⃗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 valuesx⃗
) and find the shortest equivalent multi-order pixelation by combining pixels. Ifx⃗
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 allx
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 forx
, as well as a boolean mask of pixelsi
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 tocombined
.
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 valuex2
:>>> 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
andskymap
at the specifiedquantiles
assuming that each pixel inskymap
has a uniform PDF within its borders. The difference between successive values in the density-ordered CDF ofskymap
will be used to select subpixels in NEST ordering from the pixel inuniq
with the greatest CDF value lower than the corresponding quantile, allowing for reproducible results. The returned values have the same order as the inputquantiles
. You can usehealpy.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. IfFalse
, 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 likelynp.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 anastropy.wcs.WCS
instance’s pixels as well as theirx
andy
pixel coordinates. All returned pixels and coordinates will be within the image boundaries and the indices will be non-repeating. Iforder_delta
is provided, then the NEST resolution is doubled the number of times specified thereby. Can only pass one ofnside
ororder_delta
or else aValueError
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 anastropy.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