amaurea / enlib

5 stars 12 forks source link

aberration crashing #41

Open ajvanengelen opened 6 years ago

ajvanengelen commented 6 years ago

I am getting the following error within coordinates.py, using a fresh git pull of enlib. I am getting this when using the aberrate routines -- these have not changed in a couple of years, but coordinates.py has been updated within the last 26 days, evidently.

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
/global/u2/e/engelen/actsims/bin/testAberrate.py in <module>()
     31 aberratedMap = aberration.aberrate(inMap,
     32                                    aberration.dir_equ,
---> 33                                    aberration.beta)#, modulation = True)
     34 

/global/u2/e/engelen/enlib/aberration.py in aberrate(imap, dir, beta, mode, order, recenter)
     40         # It is a bit confusing to have different conventions in enmap
     41         # and coordinates.
---> 42         pos = remap(pos[::-1], dir, beta, pol=pol, recenter=recenter)
     43         pos[:2] = pos[1::-1]
     44         pix = imap.sky2pix(pos[:2], corner=True) # interpol needs corners

/global/u2/e/engelen/enlib/aberration.py in remap(pos, dir, beta, pol, modulation, recenter)
     13         will have three columns, with the third column being
     14     the aberration-induced rotation of the polarization angle."""
---> 15         pos = coordinates.transform("equ",["equ",dir],pos,pol=pol)
     16         if recenter: before = np.mean(pos[1,::10])
     17         # Use -beta here to get the original position from the deflection,

/global/u2/e/engelen/enlib/coordinates/coordinates.pyc in transform(from_sys, to_sys, coords, time, site, pol, mag, bore)
     40                 if len(coords) > 2: fields.append("ang")
     41                 if len(coords) > 3: fields.append("mag")
---> 42         meta = transform_meta(transfunc, coords[:2], fields=fields)
     43 
     44         # Fix the polarization convention. We use healpix

/global/u2/e/engelen/enlib/coordinates/coordinates.pyc in transform_meta(transfun, coords, fields, offset)
     81         for i in range(ntrans):
     82                 # Transpose to get broadcasting right
---> 83                 a = transfun((coords.T + offsets[i].T).T)
     84                 if ocoords is None:
     85                         ocoords = np.zeros((ntrans,)+a.shape, a.dtype)

/global/u2/e/engelen/enlib/coordinates/coordinates.pyc in transfunc(coords)
     33         # polarization rotation and apparent magnification
     34         def transfunc(coords):
---> 35                 return transform_raw(from_info, to_info, coords, time=time, site=site, bore=bore)
     36         fields = []
     37         if pol: fields.append("ang")

/global/u2/e/engelen/enlib/coordinates/coordinates.pyc in transform_raw(from_sys, to_sys, coords, time, site, bore)
    175                         coords[:] = transform_astropy(from_sys, to_sys_astropy, coords)
    176                         from_sys = to_sys_astropy
--> 177         if to_ref is not None: coords[:] = recenter(coords, to_ref[0], restore=to_ref[1])
    178         return coords.reshape(oshape)
    179 

/global/u2/e/engelen/enlib/coordinates/coordinates.pyc in recenter(angs, center, restore)
    275         # Now supports specifying where to recenter by specifying center as
    276         # lon_from,lat_from,lon_to,lat_to
--> 277         if len(center) == 4: ra0, dec0, ra1, dec1 = center
    278         elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
    279         if restore: ra1 += ra0

TypeError: object of type 'numpy.float64' has no len()
> /global/u2/e/engelen/enlib/coordinates/coordinates.py(277)recenter()
    275         # Now supports specifying where to recenter by specifying center as
    276         # lon_from,lat_from,lon_to,lat_to
--> 277         if len(center) == 4: ra0, dec0, ra1, dec1 = center
    278         elif len(center) == 2: ra0, dec0, ra1, dec1 = center[0], center[1], center[0]*0, center[1]*0+np.pi/2
    279         if restore: ra1 += ra0
ajvanengelen commented 6 years ago

Thanks -- this is working now for a stack of three maps (I, Q, U), which is my nominal usage case.

I am finding, though, that it is not working for my initial test run, where I sent it a single map. In that case there is a shape mismatch at some point:

/global/u2/e/engelen/actsims/bin/testAberrate.py in <module>()
     31 aberratedMap = aberration.aberrate(inMap,
     32                                    aberration.dir_equ,
---> 33                                    aberration.beta)#, modulation = True)
     34 

/global/u2/e/engelen/enlib/aberration.py in aberrate(imap, dir, beta, mode, order, recenter)
     48                 omap[1] = c*omap[1] + s*omap[2]
     49                 omap[2] =-s*omap[1] + c*omap[2]
---> 50         omap *= pos[2+pol,None]
     51         return omap
     52 

ValueError: non-broadcastable output operand with shape (10801,21600) doesn't match the broadcast shape (1,10801,21600)

Not sure if you'd call this a bug or not, technically.

Alex