pix2sky breaks at RA = +/- 180 due to reliance on floating-point arithmetic

[zatkins2]

• pixell version: 0.15.1
• Python version: 3.8.3
• Operating System: Springdale 8.6

Description

enmap.pix2sky and enmap.corners call utils.rewind to handle wrapping around the full-sky. The way utils.rewind is implemented is such that if a pixel RA center is within double precision of np.pi, then the pixel returns a RA coordinate of -np.pi (this line). In other words, the bottom-left corner of such a map, enmap.pix2sky(shape, wcs, [[0], [0]]), will have a RA coordinate of -np.pi, not np.pi as is the pixell convention (see comment in fullsky_geometry here that the first column of pixels should correspond to RA=180).

This issue is propagates to other functions that rely on pix2sky, e.g. to enmap.posmap, which will give an RA array spanning (-np.pi, -3*np.pi) instead of (np.pi, -np.pi) as one would like, since enmap.unwind "starts counting" coordinates at the leftmost pixel.

When does this happen?

This issue arises (or not) in interesting cases. Oddly, the dr6v3 maps on-disk have a wcs (likely some floating-point arithmetic stuff with their cdelt etc.) such that this issue just barely doesn't happen:

shape
>>>(10320, 43200)
wcs
>>>car:{cdelt:[-0.008333,0.008333],crval:[0,0],crpix:[21601.00,7561.00]}
wcs.wcs.cdelt[0] # for better printing precision
>>>-0.0083333333333333
enmap.pix2sky(shape, wcs, [[0], [0]])[1, 0] * (180/np.pi)
>>>179.9999999999993

However, if I make a downgraded geometry (e.g. by a factor of 4), whose pixel centers exactly align with pixel centers in the input geometry, this problem occurs:

shape4
>>>(2580, 10800)
wcs4
>>>car:{cdelt:[-0.03333,0.03333],crval:[0.01667,0],crpix:[5400.50,1891.00]}
wcs4.wcs.cdelt[0] # for better printing precision
>>>-0.03333333333333333
enmap.pix2sky(shape4, wcs4, [[0], [0]])[1, 0] * (180/np.pi)
>>>-180.0

If Python could do fractions exactly, we know that these two results should be the same (as a user, I would like them to both be +180.0).

Why does this happen?

Here's what I mean by "if a pixel RA center is within double precision of np.pi":

pix = [[0], [0]]
pix = np.asarray(pix).astype(float)
pflat = pix.reshape(pix.shape[0], -1)

coords = np.asarray(wcsutils.nobcheck(wcs).wcs_pix2world(*(tuple(pflat)[::-1]+(0,)))[::-1])*enmap.get_unit(wcs)
coords4 = np.asarray(wcsutils.nobcheck(wcs4).wcs_pix2world(*(tuple(pflat)[::-1]+(0,)))[::-1])*enmap.get_unit(wcs4)

where I have exactly copied the implementation of pix2sky before the call to rewind. These "true" coordinates give:

coords[1, 0]*180/np.pi
>>>179.9999999999993
coords4[1, 0]*180/np.pi
>>>179.99999999999997

Note the second value is just a smidgen closer to 180 than the dr6v3 geometry value, but it's close enough to cause this issue:

(179.9999999999993 + 180)%360
>>>359.9999999999993
(179.99999999999997 + 180)%360
>>>0.0

where this little arithmetic is what happens in utils.rewind.

How to fix?

I think there are a bunch of ways, but perhaps the most explicit is to add a recentering function (not wcsutils.fix_wcs) which moves the map in steps of 2pi such that pix2sky and corners always give sensible outputs. This could be a new argument like centered=True so that someone could get the old behavior manually if they wanted.

[zatkins2]

(While posmap breaking for my special downgrade-by-4 geometry is annoying, I am more worried about something going wrong in map2alm, alm2map etc. Looking at the code, I think it still calls the buggy pix2sky to get phi0 in minfo, but I'm guessing/hoping phi0=np.pi and phi0=-np.pi should give me the same results??)

[amaurea]

Rotating the sky by 360° has no physical effect, so don't worry about map2alm or alm2map. Similarly, having postmap output coordinates that are 360° off from standard values doesn't really do anything either, it's just annoying. But it's worth it to fix that annoyance. I've made a fix, but I need a self-contained test case to test it.

[zatkins2]

As discussed, the following example should produce an offending geometry:

from pixell import enmap
import numpy as np

def downgrade_geometry_cc_quad(shape, wcs, dg):
    """Get downgraded geometry that adheres to Clenshaw-Curtis quadrature.

    Parameters
    ----------
    shape : tuple
        Shape of map geometry to be downgraded.
    wcs : astropy.wcs.WCS
        Wcs of map geometry to be downgraded
    dg : int or float
        Downgrade factor.

    Returns
    -------
    (tuple, astropy.wcs.WCS)
        The shape and wcs of the downgraded geometry. If dg <= 1, the 
        geometry of the input map.

    Notes
    -----
    Works by generating a fullsky geometry at downgraded resolution and slicing
    out coordinates of original map.
    """
    if dg <= 1:
        return shape[-2:], wcs
    else:
        # get the shape, wcs corresponding to the sliced fullsky geometry, with resolution
        # downgraded vs. imap resolution by factor dg, containg sky footprint of imap
        res = np.deg2rad(np.abs(wcs.wcs.cdelt)) * dg
        full_dshape, full_dwcs = enmap.fullsky_geometry(res=res)
        full_dpixbox = enmap.skybox2pixbox(
            full_dshape, full_dwcs, enmap.corners(shape, wcs, corner=False), corner=False
            )

        # we need to round this to an exact integer in order to guarantee that 
        # the sliced geometry is still clenshaw-curtis compatible. we use 
        # np.round here (as opposed to np.floor) since we want pixel centers,
        # see enmap.sky2pix documemtation
        full_dpixbox = np.round(full_dpixbox).astype(int)

        slice_dshape, slice_dwcs = slice_geometry_by_pixbox(
            full_dshape, full_dwcs, full_dpixbox
            )
        return slice_dshape, slice_dwcs

def slice_geometry_by_pixbox(ishape, iwcs, pixbox):
    pb = np.asarray(pixbox)
    return enmap.slice_geometry(
        ishape[-2:], iwcs, (slice(*pb[:, -2]), slice(*pb[:, -1])), nowrap=True
        )

# assume shape, wcs is the geometry from a raw dr6 map
shape4, wcs4 = downgrade_geometry_cc_quad(shape, wcs, 4)

print(enmap.pix2sky(shape4, wcs4, pix=[[0], [0]], corner=False))

gives

array([[-1.09955743],
       [-3.14159265]])

[amaurea]

This should be fixed with commit #204

[zatkins2]

It seems like this may still be happening:

import pixell
print(pixell.__version__)

from pixell import enmap
import numpy as np

print(np.rad2deg(enmap.pix2sky(*enmap.fullsky_geometry(res=np.deg2rad(1/120), variant='CC'), [[0], [0]])))
print(np.rad2deg(enmap.pix2sky(*enmap.fullsky_geometry(res=np.deg2rad(1/120), variant='fejer1'), [[0], [0]])))

gives

0.17.3
[[ -90.]
 [-180.]]
[[ -89.99583333]
 [-180.        ]]

[amaurea]

What's the bug here, @zatkins2? Fejer1 and CC aren't supposed to have the same coordinates for the first pixel. The output looks correct to me.

[zatkins2]

Yes, the point of the code snippet was not that the Dec coordinates of the two geometries differ (of course they do), but rather that the RA coordinate of the first pixel is still -180 deg regardless of the ring scheme.

I've just tested this again for the most-recent pixell version and it's the same:

import pixell
print(pixell.__version__)

from pixell import enmap
import numpy as np

print(np.rad2deg(enmap.pix2sky(*enmap.fullsky_geometry(res=np.deg2rad(1/120), variant='CC'), [[0], [0]])))
print(np.rad2deg(enmap.pix2sky(*enmap.fullsky_geometry(res=np.deg2rad(1/120), variant='fejer1'), [[0], [0]])))

gives

0.23.8
[[ -90.]
 [-180.]]
[[ -89.99583333]
 [-180.        ]]

[jadejmvo]

I encountered a similar problem. A discrepancy between pixell.enmap.pix2sky and astropy.wcs.all_pix2world when converting pixel coordinates to world coordinates at the borders of a map. The RA and Dec values returned by pix2sky appear to be incorrect (Dec out of [-90, 90]), while astropy produces consistent and expected results.