Enmap modifies wcs when read/writing from file

[damonge]

I believe this is mostly a floating-point error, and it only applies for full-sky maps in CAR.

The issue is that the WCS object created by enmap for a full-sky map differs slightly from the WCS that would be read from a fits file after the same map has been written to file. The difference is in cdelt, and I think it's mostly numerical, but it does prevent you from, for example, finding the coordinates of pixels around the poles.

This was first noted by @kwolz

See example below.

from pixell import enmap
import numpy as np

# Create map
res = 10. * np.pi/180/60
shape, wcs = enmap.fullsky_geometry(res=res, proj='car', variant='CC')
mask = enmap.ones((*shape,), wcs)
map = enmap.ones((3, *shape,), wcs)

# Write map and retrieve wcs
enmap.write_map("dummy.fits", mask)
mask_ones = enmap.read_map("dummy.fits")
wcs_b = mask_ones.wcs

# Compare WCSs
print("Difference in cdelt")
print(wcs.wcs.cdelt-wcs_b.wcs.cdelt)
print(" ")

# Now try to compute the coordinates of one corner of this map
print("Pole coordinates according to each WCS:")
print(wcs.wcs_pix2world([[0, 0]], 0))
print(wcs_b.wcs_pix2world([[0, 0]], 0))

The output of this is:

Difference in cdelt
[ 3.33066907e-15 -3.33066907e-15]
 
Pole coordinates according to each WCS:
[[180. -90.]]
[[nan nan]]

[zatkins2]

(Related: https://github.com/simonsobs/pixell/issues/202)

Just out of curiosity, @damonge what happens if instead of calling the wcs_pix2world method of wcs, you use enmap.pix2sky instead?

[damonge]

Thanks @zatkins2 .

OK, interesting. pix2sky does not crash, although it does put this point outside of the sphere by a little bit. What I mean by that, is that it assigns a value of the latitude that is smaller than -pi/2:

print(enmap.pix2sky(mask_ones.shape, wcs, [[0], [0]]).flatten()-np.array([-np.pi/2, -np.pi]))
print(enmap.pix2sky(mask_ones.shape, wcs_b, [[0], [0]]).flatten()-np.array([-np.pi/2, -np.pi]))

yields

[0. 0.]
[-3.13082893e-14  6.30606678e-14]

What's different here wrt #202 is that the issue is with latitude (which you can't wrap around), not longitude (which you can),

[damonge]

This is similar to the comment by @jadejmvo