Map functions

Functions like pix2angNest and ang2pixNest fully define the Healpix tessellation scheme. They are however extremely impractical in a number of situations. It happens often that a large fraction of pixels in a map need to be processed together. Healpix.jl introduces the HealpixMap{T, O <: Order} type, which acts as a collection of all the pixels on the sphere. A HealpixMap type holds the value of all the pixels in its pixels field, and it keeps track of the ordering (either RING or NESTED). Here is an example that shows how to create a map and initialize it:

nside = 32
m = HealpixMap{Float64, RingOrder}(nside)
m.pixels[:] = 1.0  # Set all pixels to 1

Healpix.jl defines the basic operations on maps (sum, subtraction, multiplication, division). These operations can either combine two maps or a map and a scalar value:

mollweide(m * 2.0)
mollweide(m * m)

The HealpixMap{T, O <: Order} is derived from the abstract type AbstractHealpixMap{T}, which does not encode the ordering. It is useful for functions that can either work on ring/nested-ordered maps but cannot be executed on plain generic arrays:

# Return the number of pixels in the map, regardless of its ordering
maplength(m::AbstractHealpixMap{T}) where T = length(m)

# This returns 12
maplength(HealpixMap{Float64, RingOrder}(1))

# This too returns 12
maplength(HealpixMap{Float64, NestedOrder}(1))

# This fails
maplength(zeros(Float64, 12))

Healpix.jl implements the PolarizedHealpixMap{T, O <: Order} type as well, which derives from AbstractPolarizedHealpixMap{T}. This encodes three maps containing the I/Q/U signal: the intensity (I), and the Q and U Stokes parameters. The three maps must have the same resolution.

Healpix.AbstractHealpixMap — Type

AbstractHealpixMap{T} <: AbstractArray{T, 1}

An abstract type representing an Healpix map without a specified ordering. This can be used to implement multiple dispatch when you don't care about the ordering of a map.

Healpix.HealpixMap — Type

struct HealpixMap{T, O <: Order, AA <: AbstractArray{T, 1}} <: AbstractHealpixMap{T}

A Healpix map. The type T is used for the value of the pixels in a map, and it can be anything (even a string!). The type O is used to specify the ordering of the pixels, and it can either be RingOrder or NestedOrder. The type AA is used to store the array of pixels; typical types are Vector, CUArray, SharedArray, etc.

A HealpixMap type contains the following fields:

pixels: array of pixels
resolution: instance of a Resolution object

You can construct a map using one of the following forms:

HealpixMap{T, O, AA}(arr) and HealpixMap{T, O, AA}(nside::Number) will use AA as basetype
HealpixMap{T, O}(arr) and HealpixMap{T, O}(nside::Number) will use Array{T, 1} as basetype

Examples

The following example creates a map with NSIDE=32 in RING order, containing integer values starting from 1:

mymap = Healpix.HealpixMap{Int64, Healpix.RingOrder}(1:Healpix.nside2npix(32))

The call to collect is required to convert the range in an array.

This example creates a map in NESTED order, with NSIDE=64, filled with zeroes:

mymap = Healpix.HealpixMap{Float64, Healpix.NestedOrder}(64)

Finally, the following examples show how to use SharedArray:

using SharedArrays

# Create a map with all pixels set to zero
mymap = Healpix.HealpixMap{Float64, Healpix.NestedOrder, SharedArray{Float64, 1}}(64)

# Create a map and initialize pixel values with a SharedArray
pixels = SharedArray{Int64, 1}(1:12 |> collect)
mymap = Healpix.HealpixMap{Int64, Healpix.RingOrder, SharedArray{Int64, 1}}(m)

Healpix.AbstractPolarizedHealpixMap — Type

AbstractPolarizedHealpixMap{T}

An abstract type representing an Healpix polarized map without a specified ordering. This can be used to implement multiple dispatch when you don't care about the ordering of a map.

Healpix.PolarizedHealpixMap — Type

mutable struct PolarizedHealpixMap{T, O <: Healpix.Order, AA <: AbstractArray{T, 1}}

A polarized I/Q/U map. It contains three Healpix maps with the same NSIDE:

i
q
u

You can create an instance of this type using the function PolarizedHealpixMap{T,O}, which comes in three flavours:

PolarizedHealpixMap(i::HealpixMap{T,O,AA}, q::HealpixMap{T,O,AA}, u::HealpixMap{T,O,AA})
PolarizedHealpixMap{T,O}(i::AbstractVector{T}, q::AbstractVector{T}, u::AbstractVector{T})
PolarizedHealpixMap{T,O}(nside::Number)

Unseen pixels and nothingness

You can use the constant UNSEEN to mark unseen pixels, i.e., pixels that lack data associated with them, in a way that is compatible with other Healpix libraries.

m = HealpixMap{Float32, RingOrder}(32)

# Mark all the pixels in the map as «unseen» (missing)
m[:] .= UNSEEN

12288-element view(::HealpixMap{Float32, RingOrder, Vector{Float32}}, :) with eltype Float32:
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
  ⋮
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30
 -1.6375f30

However, Julia provides a nicer way to denote missing pixels through the use of Nothing (type) and nothing (value). Whenever you pass a Union{Nothing, T} type to a Healpix map, the map will be initialized to nothing, and you can test if a pixel has been observer or not using Julia's function isnothing:

m = HealpixMap{Union{Int32, Nothing}, RingOrder}(1)

m[:] = 1:12
m[5] = nothing
@assert isnothing(m[5])

Note that, unlike UNSEEN, this mechanism permits to signal «missing» pixels even for maps that do not use floating-point numbers.

Type stability and Nothing

Using Union{Nothing, T} as the base type of a Healpix map can lead to elegant code, but it is not likely to be efficient!

For instance, consider two implementations of the same code, which sums all the pixels in a map that are not marked as nothing or UNSEEN:

function sumpixels(m::HealpixMap{Union{Nothing, T}, O}) where {T <: Real, O <: Order}
   cumsum = zero(Float64)
   @inbounds for i in eachindex(m)
       (!isnothing(m[i])) && (cumsum += m[i])
   end
   cumsum
end

function sumpixels(m::HealpixMap{T, O}) where {T <: Real, O <: Order}
   cumsum = zero(Float64)
   @inbounds for i in eachindex(m)
       (m[i] != UNSEEN) && (cumsum += m[i])
   end
   cumsum
end

Let's now create two maps with random values and 50% of their pixels marked either as nothing or UNSEEN:

import Random
mnothing = HealpixMap{Union{Float64, Nothing}, RingOrder}(1024)
mnothing[:] = rand(length(mnothing))
@. mnothing[mnothing < 0.5] = nothing

# Create a new map identical to `mnothing`, but use UNSEEN instead of nothing
m = HealpixMap{Float64, RingOrder}(mnothing.resolution.nside)
m[:] = [isnothing(x) ? UNSEEN : x for x in mnothing]

Running sumpixels over the two maps shows that the version with UNSEEN is three times faster.

julia> @benchmark sumpixels(mnothing)
@benchmark sumpixels(m)
BenchmarkTools.Trial: 56 samples with 1 evaluation.
 Range (min … max):  88.233 ms …  91.043 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     89.187 ms               ┊ GC (median):    0.00%
 Time  (mean ± σ):   89.315 ms ± 502.629 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

                  █ █  ▂ █▅▅                                    
  ▅▁▁▁▁▁▁▁▁█▁▁▁▁▁▁█▁█▅██▁███▅█▅██▅█▅▁▅▁▅▁▁▁▁▁▅█▁█▅▁▁▁▁▅▅▁▁▁▁▁▅ ▁
  88.2 ms         Histogram: frequency by time         90.5 ms <

 Memory estimate: 16 bytes, allocs estimate: 1.

julia> @benchmark sumpixels(m)
BenchmarkTools.Trial: 201 samples with 1 evaluation.
 Range (min … max):  22.579 ms … 29.701 ms  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     24.213 ms              ┊ GC (median):    0.00%
 Time  (mean ± σ):   24.906 ms ±  1.939 ms  ┊ GC (mean ± σ):  0.00% ± 0.00%

   ▁ ▃█▁▁▅▁▃▁       ▂                                  ▁       
  ▅█▆█████████▄▆█▁█▅█▆█▃▄█▄▄▃▁▃▆▁▃▃▁▄▃▃▄▁▄▇█▅▁▁▁▁▆▃▃▄▅▄█▅▃▃▁▃ ▄
  22.6 ms         Histogram: frequency by time        29.2 ms <

 Memory estimate: 16 bytes, allocs estimate: 1.

This happens because the implementation of sumpixels that uses nothing values is not type-stable. You should decide if the elegance of using nothing in your code is worth this degradation in performance or not.

You can use saveToFITS or readMapFromFITS on maps whose base type is Union{Nothing, T} only if T is a floating-point number, because for the sake of compatibility with other Healpix libraries the FITS file will use UNSEEN to mark missing values.

Healpix.UNSEEN — Constant

A constant commonly used by Healpix libraries to mark «missing» pixels.

This constant is useful if you need compatibility with other Healpix libraries.

Encoding the order

Healpix.jl distinguishes between RING and NEST orderings using Julia's typesystem. The abstract type Order has two descendeants, RingOrder and NestedOrder, which are used to instantiate objects of type HealpixMap. Applying the functions nest2ring and ring2nest to maps converts those maps to the appropriate orders. In-place nest2ring! and ring2nest! versions are also available.

Healpix.Order — Type

Abstract type representing the ordering of pixels in a Healpix map. See also RingOrder and NestedOrder.

Healpix.RingOrder — Type

The RingOrder type should be used when creating HealpixMap types in order to specify that the pixels in the map are sorted in ring ordering. (See also NestedOrder.)

Healpix.NestedOrder — Type

The NestedOrder type should be used when creating HealpixMap types in order to specify that the pixels in the map are sorted in ring ordering. (See also RingOrder.)

Healpix.nest2ring — Method

nest2ring(m_nest::HealpixMap{T, NestedOrder, AA}) where {T, AA}

Convert a map from nested to ring order. This version allocates a new array of the same array type as the input.

Arguments:

m_nest::HealpixMap{T, NestedOrder, AA}: map of type NestedOrder

Returns:

HealpixMap{T, RingOrder, AA}: the input map converted to RingOrder

Examples

julia> m_nest = HealpixMap{Float64,NestedOrder}(rand(nside2npix(64)));

julia> nest2ring(m_nest)
49152-element HealpixMap{Float64, RingOrder, Vector{Float64}}:
 0.4703834205807309
 ⋮
 0.3945848051663148

Healpix.ring2nest — Method

ring2nest(m_ring::HealpixMap{T, RingOrder, AA}) where {T, AA}

Convert a map from ring to nested order. This version allocates a new array of the same array type as the input.

Arguments:

m_ring::HealpixMap{T, RingOrder, AA}: map of type RingOrder

Returns:

HealpixMap{T, NestedOrder, AA}: the input map converted to NestedOrder

Examples

julia> m_ring = HealpixMap{Float64,RingOrder}(rand(nside2npix(64)));

julia> ring2nest(m_ring)
49152-element HealpixMap{Float64, NestedOrder, Vector{Float64}}:
 0.0673134062168923
 ⋮
 0.703460503535335

Healpix.nest2ring! — Method

nest2ring!(m_ring_dst::HealpixMap{T, RingOrder, AAR},
           m_nest_src::HealpixMap{T, NestedOrder, AAN}) where {T, AAN, AAR}

Convert a map from nested to ring order. This version takes a nested map in the second argument and writes it to the nested map provided in the first argument, following the standard Julia func!(dst, src) convention.

Arguments:

m_ring_dst::HealpixMap{T, NestedOrder, AA}: map of type NestedOrder
m_nest_src::HealpixMap{T, NestedOrder, AAN}: map of type RingOrder

Returns:

HealpixMap{T, RingOrder, AA}: the input map converted to RingOrder

Examples

julia> m_nest = HealpixMap{Float64,NestedOrder}(rand(nside2npix(64)));

julia> m_ring = HealpixMap{Float64,RingOrder}(64);

julia> nest2ring!(m_ring, m_nest)
49152-element HealpixMap{Float64, RingOrder, Vector{Float64}}:
 0.33681791815569895
 ⋮
 0.9092457003948482

Healpix.ring2nest! — Method

ring2nest!(m_nest_dst::HealpixMap{T, NestedOrder, AAN},
           m_ring_src::HealpixMap{T, RingOrder, AAR}) where {T, AAR, AAN}

Convert a map from ring to nested order. This version takes a nested map in the second argument and writes it to the nested map provided in the first argument, following the standard Julia func!(dst, src) convention.

Arguments:

m_nest_dst::HealpixMap{T, NestedOrder, AAN}: map of type RingOrder
m_ring_src::HealpixMap{T, RingOrder, AA}: map of type RingOrder

Returns:

HealpixMap{T, NestedOrder, AA}: the input map converted to NestedOrder

Examples

julia> m_ring = HealpixMap{Float64,RingOrder}(rand(nside2npix(64)));

julia> m_nest = HealpixMap{Float64,RingOrder}(64);

julia> ring2nest!(m_nest, m_ring)
49152-element HealpixMap{Float64, NestedOrder, Vector{Float64}}:
 0.0673134062168923
 ⋮
 0.703460503535335

Pixel functions

When working with maps, it is not needed to pick between ang2pixNest and ang2pixRing because a HealpixMap type already encodes the ordering. Functions pix2ang and ang2pix always choose the correct ordering, but they require a HealpixMap instead of a Resolution as their first argument.

Healpix.pix2ang — Function

pix2ang{T, O <: Order}(map::HealpixMap{T, O}, ipix) -> (Float64, Float64)
pix2ang{T, O <: Order}(map::PolarizedHealpixMap{T, O}, ipix) -> (Float64, Float64)

Return the pair (theta, phi), where theta is the colatitude and phi the longitude of the direction of the pixel center with index ipix.

Healpix.ang2pix — Function

ang2pix{T, O, AA}(map::HealpixMap{T, O}, theta, phi)
ang2pix{T, O, AA}(map::PolarizedHealpixMap{T, O}, theta, phi)

Convert the direction specified by the colatitude theta (∈ [0, π]) and the longitude phi (∈ [0, 2π]) into the index of the pixel in the Healpix map map.

Loading and saving maps

Healpix.jl implements a number of functions to save maps in FITS files.

Healpix.saveToFITS — Function

saveToFITS(map::HealpixMap{T, O}, filename::AbstractString, typechar="D", unit="", extname="MAP") where {T <: Number, O <: Order}
saveToFITS(map::PolarizedHealpixMap{T, O}, filename::AbstractString, typechar="D", unit="", extname="MAP") where {T <: Number, O <: Order}

Save a map into a FITS file. The name of the file is specified in filename; if it begins with !, existing files will be overwritten without warning. The parameter typechar specifies the data type to be used in the FITS file: the default (D) will save 64-bit floating-point values. See the CCFITSIO documentation for other values. The keyword unit specifies the measure unit used for the pixels in the map. The keyword extname specifies the name of the HDU where the map pixels will be written.

Function savePixelsToFITS is a low-level function. It knows nothing about the ordering schema used for the pixels, so the caller should manually write the ORDERING keyword in the HDU header by itself.

Healpix.savePixelsToFITS — Function

savePixelsToFITS(map::HealpixMap{T}, f::CFITSIO.FITSFile, column) where {T <: Number}

Save the pixels of map into the column with index/name column in the FITS file, which must have been already opened.

To load a map from a FITS file, you can either use readMapFromFITS or readPolarizedMapFromFITS.

Healpix.readMapFromFITS — Function

readMapFromFITS{T}(f::CFITSIO.FITSFILE, column, t::Type{T})
readMapFromFITS{T}(fileName::String, column, t::Type{T})

Read a Healpix map from the specified (1-base indexed) column in a FITS file. The values will be read as numbers of type T. If the code fails, CFITSIO will raise an exception. (Refer to the CFITSIO library for more information.)

Healpix.readPolarizedMapFromFITS — Function

readPolarizedMapFromFITS{T}(fileName::AbstractString, column, t::Type{T})

Read a polarized map (I/Q/U) from a FITS file and return a PolarizedHealpixMap object.

The parameter column can be either a number or a 3-element tuple. In the first case, three consecutive columns will be read starting from column (1-based index), otherwise the three column indices will be used.

Testing for conformability

It often happens that two Healpix maps need to be combined together: for instance, pixels on a sky map might need to be masked using a sky mask, or one map might need to be subtracted from another one. «Conformability» means that the operation between the two maps can be done directly on the pixels, without oordering or resolution conversions. The function conformables checks this.


julia> m1 = HealpixMap{Float64, RingOrder}(1)12-element HealpixMap{Float64, RingOrder, Vector{Float64}}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> m2 = HealpixMap{Float64, RingOrder}(1)12-element HealpixMap{Float64, RingOrder, Vector{Float64}}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> m3 = HealpixMap{Float64, NestedOrder}(1)12-element HealpixMap{Float64, NestedOrder, Vector{Float64}}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> m4 = HealpixMap{Float64, NestedOrder}(2)48-element HealpixMap{Float64, NestedOrder, Vector{Float64}}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 ⋮
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0
julia> conformables(m1, m2)true
julia> conformables(m1, m3)false
julia> conformables(m1, m4)false

Healpix.conformables — Function

conformables{T, S, O1, O2}(map1::HealpixMap{T, O1, AA1},
                           map2::HealpixMap{S, O2, AA2}) -> Bool
conformables{T, S, O1, O2}(map1::PolarizedHealpixMap{T, O1, AA1},
                           map2::PolarizedHealpixMap{S, O2, AA2}) -> Bool

Determine if two Healpix maps are "conformables", i.e., if their shape and ordering are the same. The array types AA1 and AA2 are not considered in testing conformability.

Interpolation

The fact that a Healpix map is, well, pixelized means that there is a sharp boundary between adjacent pixels. This can lead to undesidable effects, and therefore Healpix.jl provides a function, interpolate, that returns the interpolated value of the map along some direction in the sky:

If the direction (θ, ɸ) passes through the center of a pixel, interpolate returns the value of the pixel itself;
Otherwise, the value of the pixel and its neighbours will be interpolated using a linear function to return a weighted value.

Healpix.interpolate — Function

interpolate(m::HealpixMap{T, RingOrder, AA}, θ, ϕ) -> Value
interpolate(m::HealpixMap{T, RingOrder, AA}, θ, ϕ, pixbuf, weightbuf) -> Value

Return an interpolated value of the map along the specified direction.

When provided, the parameters pixbuf and weightbuf must be 4-element arrays of integer and floating-point values, respectively. They can be reused across multiple calls of interpolate, to save heap allocations, and they do not need to be initialized, as they are used internally:

pixbuf = Array{Int}(undef, 4)
weightbuf = Array{Float64}(undef, 4)

m = HealpixMap{Float64, RingOrder}(1)
for (θ, ϕ) in [(0., 0.), (π/2, π/2)]
    println(interpolate(m, θ, ϕ, pixbuf, weightbuf))
end

Passing pixbuf and weightbuf saves some time, as this simple benchmark shows:

julia> @benchmark interpolate(m, rand(), rand(), pixbuf, weightbuf)
BenchmarkTools.Trial:
  memory estimate:  618 bytes
  allocs estimate:  9
  --------------
  minimum time:     283.184 ns (0.00% GC)
  median time:      296.140 ns (0.00% GC)
  mean time:        348.132 ns (9.55% GC)
  maximum time:     10.627 μs (95.88% GC)
  --------------
  samples:          10000
  evals/sample:     282

julia> @benchmark interpolate(m, rand(), rand())
BenchmarkTools.Trial:
  memory estimate:  837 bytes
  allocs estimate:  11
  --------------
  minimum time:     329.825 ns (0.00% GC)
  median time:      345.504 ns (0.00% GC)
  mean time:        417.004 ns (11.04% GC)
  maximum time:     13.733 μs (96.21% GC)
  --------------
  samples:          10000
  evals/sample:     223

Upgrading and Downgrading

Changing resolution is done with udgrade. This is very fast for nested orderings, but slow for ring ordering.

Healpix.udgrade — Function

udgrade(input_map::HealpixMap{T,O,AA}, output_nside; kw...) where {T,O,AA} -> HealpixMap{T,O,AA}

Upgrades or downgrades a map to a target nside. Always makes a copy. This is very fast for nested orderings, but slow for ring because one needs to transform to nested ordering first.

Arguments:

input_map::HealpixMap{T,O,AA}: the map to upgrade/downgrade
output_nside: desired nside

Keywords:

threshold=abs(1e-6UNSEEN): absolute tolerance for identifying a bad pixel vs UNSEEN
pess=false: if false, estimate pixels from remaining good pixels when downgrading. if true, the entire downgraded pixel is set to UNSEEN.

Returns:

HealpixMap{T,O,AA}: upgraded/downgraded map in the same ordering as the input

Examples

julia> A = HealpixMap{Float64, NestedOrder}(ones(nside2npix(4)))
192-element HealpixMap{Float64, RingOrder, Vector{Float64}}:
 1.0
 ⋮
 1.0

julia> Healpix.udgrade(A, 2)
48-element HealpixMap{Float64, NestedOrder, Vector{Float64}}:
 1.0
 ⋮
 1.0

Map-making

Map-making is the process of converting a time series of measurements into a sky map. The most basic form of map-making is the so-called "binning", where samples in the time stream falling within the same sky pixel are averaged. This map-making algorithm is strictly accurate only if the noise in the time stream is white.

Healpix.jl implements two functions to perform binning, tod2map and combinemaps!.

Healpix.tod2map — Function

tod2map{T,O}(pixidx, tod::Array{T}; nside=128) :: (map, hits)

Create a binned map for a TOD and return a tuple containing the map itself and the hit map.

Healpix.combinemaps! — Function

combinemaps{T, O, H}(destmap::HealpixMap{T, O}, desthitmap::HealpixMap{H, O}, othermap::HealpixMap{T, O}, otherhitmap::HealpixMap{H, O})

Sum "othermap" to "destmap", assuming that both maps have been produced by binning TODs. The parameters desthitmap and otherhitmap are the two hit maps. At the end of the call, destmap and desthitmap are updated.