Basic HealSparse Interface¶
Getting Started¶
To conserve memory, HealSparse uses a dual-map approach, where a low-resolution full-sky “coverage map” is combined with a high resolution map containing the pixel data where it is available.
The resolution of the coverage map is controlled by the nside_coverage
parameter, and the resolution of the high-resolution map is controlled by the nside_sparse
parameter.
Behind the scenes, HealSparse uses clever indexing to allow the user to treat these as contiguous maps with minimal overhead.
All HealSparse maps use HEALPix nest indexing behind the scenes, should be treated as nest-indexed maps.
There are 3 basic ways to make a HealSparseMap
.
First, one can read in an existing HEALPix map; second, one can read in an existing HealSparseMap
; and third, one can create a new map.
import numpy as np
import healsparse
# To read a HEALPix map, the nside_coverage must be specified
map1 = healsparse.HealSparseMap.read('healpix_map.fits', nside_coverage=32)
# To read a healsparse map, no additional keywords are necessary
map2 = healsparse.HealSparseMap.read('healsparse_map.hs')
# To read part of a healsparse map, you can specify the coverage pixels to read
map2_partial = healsparse.HealSparseMap.read('healsparse_map.hs', pixels=[100, 101])
# To create a new map, the resolutions and datatype must be specified
nside_coverage = 32
nside_sparse = 4096
map3 = healsparse.HealSparseMap.make_empty(nside_coverage, nside_sparse, np.float64)
To set values in the map, you can use simple indexing or the explicit API:
map3[0: 1000] = np.arange(1000, dtype=np.float64)
map3.update_values_pix(np.arange(1000, 2000), np.arange(1000, dtype=np.float64), nest=True)
To retrieve values from the map, you can use simple indexing or the explicit API via pixels or positions:
print(map3[0: 1000])
>>> [0. ... 999.]
print(map3.get_values_pix(np.arange(1000, 2000), nest=True))
>>> [0. ... 999.]
print(map3.get_values_pos(45.0, 0.1, lonlat=True))
>>> 51.0
A HealSparseMap
has the concept of “valid pixels”, the pixels over which the map is defined (as opposed to hpgeom.UNSEEN
in the case of floating point maps).
You can retrieve the array of valid pixels or the associated positions of the valid pixels easily:
print(map3.valid_pixels)
>>> [ 0 1 2 ... 1997 1998 1999]
ra, dec = map3.valid_pixels_pos(lonlat=True)
print(ra)
>>> [45. ... 45.3515625 ]
print(dec)
>>> [0.00932548 ... 0.81134431]
You can convert a HealSparseMap
to a healpy
map (numpy
array) either by using a full slice ([:]
) or with the generate_healpix_map()
method.
Do watch out, at high resolution this can blow away your memory!
In these cases, generate_healpix_map()
can degrade the map before conversion, using a reduction function (over valid pixels) of your choosing, including mean
, median
, std
, max
, min
, and
, or
, sum
, prod
(product), and wmean
(weighted mean).
hpmap4096 = map3[:]
hpmap128 = map3.generate_healpix_map(nside=128, reduction='mean')
Integer Maps¶
In addition to floating-point maps, which are natively supported by healpy
, HealSparseMap
supports integer maps.
The “sentinel” value of these maps (equivalent to hpgeom.UNSEEN
) is either -MAXINT
or 0
, depending on the desired use of the map (e.g., integer values or positive bitmasks).
Note that these maps cannot be trivially converted to healpy
maps because HEALPix has no concept of sentinel values that are not hpgeom.UNSEEN
, which is a very large negative floating-point value.
import numpy as np
import healsparse
map_int = healsparse.HealSparseMap.make_empty(32, 4096, np.int32)
print(map_int)
>>> HealSparseMap: nside_coverage = 32, nside_sparse = 4096, int32
map_int[0: 1000] = np.arange(1000, dtype=np.int32)
print(map_int[500])
>>> 500
Recarray Maps¶
HealSparseMap
also supports maps made up of numpy
record arrays.
These recarray maps will have one field that is the “primary” field which is used to test if a pixel has a valid value or not.
Therefore, these recarray maps should be used to describe associated values that share the exact same valid footprint.
Each field in the recarray can be treated as its own HealSparseMap
.
For example,
import numpy as np
import healsparse
dtype = [('a', np.float32), ('b', np.int32)]
map_rec = healsparse.HealSparseMap.make_empty(32, 4096, dtype, primary='a')
map_rec[0: 10000] = np.zeros(10000, dtype=dtype)
print(map_rec.valid_pixels)
>>> [ 0 1 2 ... 9997 9998 9999]
map_rec['a'][0: 5000] = np.arange(5000, dtype=np.float32)
map_rec['b'][5000: 10000] = np.arange(5000, dtype=np.int32)
print(map_rec[map_rec.valid_pixels])
>>> [(0., 0) (1., 0) (2., 0) ... (0., 4997) (0., 4998) (0., 4999)]
Note that the call map_rec['a'][0: 5000] = values
will work, but
map_rec[0: 5000]['a'] = values
will not. Also note that using the
fields of the recarray cannot be used to set new pixels, this construction
can only be used to change pixel values.
Wide Masks¶
HealSparse has support for “wide” bit masks with an arbitrary number of bits that are referred to by bit position rather than value.
This is useful, for example, when constructing a coadd coverage map where every pixel can uniquely identify the set of input exposures that contributed at the location of that pixel.
In the case of >64 input exposures you can no longer use a simple 64-bit integer bit mask.
Wide mask bits are always specified by giving a list of integer positions rather than values (e.g., use 10
to set the 10th bit instead of 1024 = 2**10
).
import numpy as np
import healsparse
map_wide = healsparse.HealSparseMap.make_empty(32, 4096, healsparse.WIDE_MASK, wide_mask_maxbits=128)
pixels = np.arange(10000)
map_wide.set_bits_pix(pixels, [4, 100])
print(map_wide.check_bits_pix(pixels, [2]))
>>> [False False False ... False False False]
print(map_wide.check_bits_pix(pixels, [4]))
>>> [ True True True ... True True True]
print(map_wide.check_bits_pix(pixels, [100]))
>>> [ True True True ... True True True]
print(map_wide.check_bits_pix(pixels, [101]))
>>> [False False False ... False False False]
# Check if any of the bits are set
print(map_wide.check_bits_pos([45.2], [0.2], [100, 101], lonlat=True))
>>> [ True]
Writing Maps¶
Writing a HealSparseMap
is easy. To write a map in the default FITS format:
map3.write('output_file.hs', clobber=False)
And to write a map in the Parquet format with pyarrow
:
map3.write('output_file.hsparquet', clobber=False, format='parquet')
Metadata¶
You can also set key/value metadata to a map that will be stored in the fits header of the file and read back in.
The keys must confirm to FITS header key standards (strings, upper case).
The metadata will be stored as a Python dictionary, and can be accessed with the metadata
property.
metadata = {'KEY1': 5, 'KEY2': 10.0}
map3.metadata = metadata
print(map3.metadata['KEY2'])
>>> 10.0
Coverage Masks¶
A HealSparseMap
contains a coverage map that defines the coarse coverage over the sky.
You can retrieve a boolean array describing which pixels are covered in the map with the coverage_mask
property:
import hpgeom as hpg
import matplotlib.pyplot as plt
cov_mask = map3.coverage_mask
cov_pixels, = np.where(cov_mask)
ra, dec = hpg.pixel_to_angle(map3.nside_coverage, cov_pixels)
plt.plot(ra, dec, 'r.')
plt.show()
It is also possible to read the coverage map of a HealSparseMap
on its own:
cov_map = healsparse.HealSparseCoverage.read('output_file.hs')
cov_mask = cov_map.coverage_mask
In some cases, you may be building a map and you already know the coverage when it will be finished.
In this case, it can be faster to initialize the memory at the beginning.
In this case, you can add cov_pixels
to the make_empty
call.
Be aware this may make the map larger than your actual coverage.
import healsparse
nside_coverage = 32
nside_sparse = 4096
map4 = healsparse.HealSparseMap.make_empty(nside_coverage, nside_sparse, np.float32,
cov_pixels=[5, 10, 20, 21])
Fractional Detection Maps¶
One can compute the fractional detection map of a HealSparseMap
with the fracdet_map()
method.
This method will compute the fractional area covered by the sparse map at an arbitrary resolution (not higher than the native resolution, and not lower than the coverage map nside_coverage
).
This is a count of the fraction of “valid” sub-pixels (those that are not equal to the sentinel value) in the original map.
These maps can be useful in conjunction with a degraded map to easily determine the coverage fraction of each degraded pixel.
In order to translate a fracdet_map
to lower resolution, the degrade()
method should be used with the default “mean” reduction operation.
If one tries to compute the fracdet_map
of an existing fracdet_map
then you will not get the expected output, because this is the fractional coverage of the fracdet_map
itself, not of the original sparse map.
Basic Visualization¶
healsparse
does not provide any built-in visualization tools.
However, it is possible to perform quick visualizations of a HealSparseMap
using the matplotlib
package.
For example, we can take render our map as a collection of hexagonal cells using matplotlib.pyplot.hexbin
:
import healsparse
import matplotlib.pyplot as plt
nside_coverage = 32
nside_sparse = 4096
# Generation of the map
hsp_map = healsparse.HealSparseMap.make_empty(nside_coverage, nside_sparse, np.float32)
idx = np.arange(2000, 6000)
hsp_map[idx] = np.random.uniform(size=idx.size).astype(np.float32)
# Visualization of the map
vpix, ra, dec = hsp_map.valid_pixels_pos(return_pixels=True)
plt.hexbin(ra, dec, C=hsp_map[vpix])
plt.colorbar()
plt.show()