Open zatkins2 opened 1 year ago
(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??)
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.
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]])
This should be fixed with commit #204
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. ]]
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.
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. ]]
Description
enmap.pix2sky
andenmap.corners
callutils.rewind
to handle wrapping around the full-sky. The wayutils.rewind
is implemented is such that if a pixel RA center is within double precision ofnp.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
, notnp.pi
as is thepixell
convention (see comment infullsky_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. toenmap.posmap
, which will give an RA array spanning(-np.pi, -3*np.pi)
instead of(np.pi, -np.pi)
as one would like, sinceenmap.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 theircdelt
etc.) such that this issue just barely doesn't happen: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:
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
":where I have exactly copied the implementation of
pix2sky
before the call torewind
. These "true" coordinates give: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:
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 (notwcsutils.fix_wcs
) which moves the map in steps of2pi
such thatpix2sky
andcorners
always give sensible outputs. This could be a new argument likecentered=True
so that someone could get the old behavior manually if they wanted.