Enlib write_map bug



In [1]:

    
%load_ext autoreload
%autoreload 2



In [2]:

    
from enlib import enmap,wcs as mwcs
import numpy as np
import sys,os

We make a full-sky arcminute resolution geometry. I've only been able to reproduce this bug for res=1.0.



In [3]:

    
res = 1.0
shape, wcs = enmap.fullsky_geometry(res=res*np.pi/180./60., proj="car")
shape = (3,)+shape

We do a pix2sky that is needed by map2alm and make sure it gives a sensible result.



In [4]:

    
ny, nx = shape[-2:]
vy,vx = enmap.pix2sky(shape, wcs, [np.arange(ny),np.zeros(ny)])
hy,hx = enmap.pix2sky(shape, wcs, [np.zeros(nx),np.arange(nx)])
print(vy,vx,hy,hx)









    



(array([-1.57079633, -1.57050544, -1.57021455, ...,  1.57021455,
        1.57050544,  1.57079633]), array([-3.14159265, -3.14159265, -3.14159265, ..., -3.14159265,
       -3.14159265, -3.14159265]), array([-1.57079633, -1.57079633, -1.57079633, ..., -1.57079633,
       -1.57079633, -1.57079633]), array([-3.14159265, -3.14188354, -3.14217443, ..., -9.4239053 ,
       -9.42419618, -9.42448707]))

It makes sense. We now save a map that has this geometry and load it back.



In [5]:

    
root = os.environ['WORK']+"/"
enmap.write_map(root+"temp.fits",enmap.zeros(shape,wcs,dtype=np.uint8))
lshape,lwcs = enmap.read_map_geometry(root+"temp.fits")

print(shape,wcs)
print(lshape,lwcs)
print(mwcs.equal(wcs,lwcs))









    



((3, 10801, 21600), car:{cdelt:[-0.01667,0.01667],crval:[0,0],crpix:[1.08e+04,5401]})
((3, 10801, 21600), car:{cdelt:[-0.01667,0.01667],crval:[0,0],crpix:[1.08e+04,5401]})
True

The shapes and wcs of the geometry we originally made and of the saved map seem to agree. So we proceed to do the same pix2sky operation on the loaded geometry.



In [6]:

    
ny, nx = lshape[-2:]
vy2,vx2 = enmap.pix2sky(lshape, lwcs, [np.arange(ny),np.zeros(ny)])
hy2,hx2 = enmap.pix2sky(lshape, lwcs, [np.zeros(nx),np.arange(nx)])
print(vy2,vx2,hy2,hx2)









    



(array([-1.57079633, -1.57050544, -1.57021455, ...,  1.57021455,
        1.57050544,  1.57079633]), array([-3.14159265, -3.14159265, -3.14159265, ..., -3.14159265,
       -3.14159265, -3.14159265]), array([-1.57079633, -1.57079633, -1.57079633, ..., -1.57079633,
       -1.57079633, -1.57079633]), array([-3.14159265, -3.14188354, -3.14217443, ..., -9.4239053 ,
       -9.42419618, -9.42448707]))

The results are all nans.



In [8]:

    
print(np.all(np.isclose(vy,vy2)))
print(np.all(np.isclose(vx,vx2)))
print(np.all(np.isclose(hy,hy2)))
print(np.all(np.isclose(hx,hx2)))









    



True
True
True
True