ntessore / healpix

Python and C package for HEALPix discretisation of the sphere
BSD 3-Clause "New" or "Revised" License
4 stars 1 forks source link

Create single-resolution maps from MOCs #37

Open ntessore opened 6 months ago

ntessore commented 6 months ago

Add a function to create single-resolution maps from MOCs, along the following lines:

def moc2map(nside, moc, uniq=False): ...

where uniq is a flag for using NUNIQ or RANGE indexing.

ntessore commented 2 weeks ago

Here's a proof-of-concept implementation using RANGE indexing:

def moc2map(nside, moc):
    order = nside.bit_length() - 1
    max_order = 29
    shift = 2 * (max_order - order)
    if shift < 0:
        raise ValueError(f"maximum NSIDE is 2^{max_order}")
    fact = 1 << shift
    ipix, jpix = np.divmod(moc, fact)
    m = np.zeros(12 * nside**2)
    for (i1, i2), (j1, j2) in zip(ipix.reshape(-1, 2), jpix.reshape(-1, 2)):
        if i1 == i2:
            m[i1] += (j2 - j1) / fact
        else:
            m[i1 : i1 + 1] += 1.0 - j1 / fact
            m[i1 + 1 : i2 - 1] += 1.0
            m[i2 : i2 + 1] += j2 / fact
    return m

The inner loop could be done in C on an array allocated at the Python layer.