Depth1 wrap fix

sotodlib/coords/__init__.py

@@ -1,7 +1,7 @@
 from .pmat import P
 from .helpers import (
     get_radec, get_horiz, get_footprint, get_wcs_kernel, get_supergeom,
-    Timer, DEG, ScalarLastQuat
+    Timer, DEG, ScalarLastQuat, normalize_geometry
 )
 from . import planets
 from . import local

sotodlib/coords/helpers.py

@@ -355,6 +355,7 @@ def get_footprint(tod, wcs_kernel, dets=None, timestamps=None, boresight=None,
         # Works whether rot is a quat or a vector of them.
         asm.Q = rot * asm.Q
     proj.get_planar(asm, output=planar)
+
     # planar is now [ndet,nsamp,{ira,idec}] in intermediate
     # coordinates in radians. These will be 0 at the
     # reference point, and span at most [-pi,pi] in ira
@@ -404,13 +405,8 @@ def get_footprint(tod, wcs_kernel, dets=None, timestamps=None, boresight=None,
     # or add a special case. NB! This case *will* trigger regularly
     # due to wrapping if recentering is turned off.
     shape = (p2-p1+1)[::-1]
-
-    if recenter:
-        # Second part of recentering. Modify crval[0]
-        x_mid  = shape[-1]//2+1
-        ra_mid = w.wcs.crval[0] + (x_mid-w.wcs.crpix[0])*w.wcs.cdelt[0]
-        w.wcs.crpix[0] = x_mid
-        w.wcs.crval[0] = ra_mid
+    # Make sure wcs crval follows so3g pointing matrix assumptions
+    shape, w = normalize_geometry(shape, w)
 
     return (shape, w)
 
@@ -796,3 +792,17 @@ def get_deflected_sightline(aman, wobble_meta, site='so', weather='typical'):
     sight.Q = sight.Q * ~deflq
     return sight
 
+def normalize_geometry(shape, wcs):
+    # Can't freely change wcs for non-separable geometries
+    # (so non-cylindrical ones), as this would change the geometry
+    # in an incompatible way
+    if not wcsutils.is_separable(wcs): return shape, wcs
+    # The pointing matrix assumes that crval ra is in [-180,180],
+    # and that all points are within 180° of this
+    wcs    = wcs.deepcopy()
+    x_mid  = shape[-1]/2+1
+    ra_mid = wcs.wcs.crval[0] + (x_mid-wcs.wcs.crpix[0])*wcs.wcs.cdelt[0]
+    ra_mid = utils.rewind(ra_mid, ref=0, period=360)
+    wcs.wcs.crpix[0] = x_mid
+    wcs.wcs.crval[0] = ra_mid
+    return shape, wcs

sotodlib/site_pipeline/make_depth1_map.py

@@ -167,15 +167,17 @@ def find_footprint(context, tods, ref_wcs, comm=mpi.COMM_WORLD, return_pixboxes=
     # Could be done more generally, but would be much more involved, and this should be
     # good enough
     nphi     = utils.nint(np.abs(360/ref_wcs.wcs.cdelt[0]))
-    widths   = pixboxes[:,1,0]-pixboxes[:,0,0]
-    pixboxes[:,0,0] = utils.rewind(pixboxes[:,0,0], ref=pixboxes[0,0,0], period=nphi)
-    pixboxes[:,1,0] = pixboxes[:,0,0] + widths
+    widths   = pixboxes[:,1,1]-pixboxes[:,0,1]
+    pixboxes[:,0,1] = utils.rewind(pixboxes[:,0,1], ref=pixboxes[0,0,1], period=nphi)
+    pixboxes[:,1,1] = pixboxes[:,0,1] + widths
     # It's now safe to find the total pixel bounding box
     union_pixbox = np.array([np.min(pixboxes[:,0],0)-pad,np.max(pixboxes[:,1],0)+pad])
     # Use this to construct the output geometry
     shape = union_pixbox[1]-union_pixbox[0]
     wcs   = ref_wcs.deepcopy()
     wcs.wcs.crpix -= union_pixbox[0,::-1]
+    # Make sure wcs crval follows so3g pointing matrix assumptions
+    shape, wcs = coords.normalize_geometry(shape, wcs)
     if return_pixboxes: return shape, wcs, pixboxes
     else: return shape, wcs
 
@@ -240,6 +242,7 @@ def calibrate_obs(obs, band, site='so', dtype_tod=np.float32, nocal=True, unit='
 def make_depth1_map(context, obslist, shape, wcs, noise_model, L, preproc, comps="TQU", t0=0, dtype_tod=np.float32, dtype_map=np.float64, comm=mpi.COMM_WORLD, tag="", niter=100, site='so', tiled=0, verbose=0, downsample=1, interpol='nearest', srcsamp_mask=None, unit='K', min_dets=50):
     pre = "" if tag is None else tag + " "
     if comm.rank == 0: L.info(pre + "Initializing equation system")
+
     # Set up our mapmaking equation
     signal_cut = mapmaking.SignalCut(comm, dtype=dtype_tod)
     signal_map = mapmaking.SignalMap(shape, wcs, comm, comps=comps, dtype=dtype_map, tiled=tiled>0, interpol=interpol, ofmt="")
@@ -361,7 +364,7 @@ def main(config_file=None, defaults=defaults, **args):
     wcs = wcsutils.WCS(wcs.to_header())
     # Set shape to None to allow map to fit these TODs exactly.
     #shape = None
-    
+
     comps = args['comps']
     ncomp = len(comps)
     dtype_tod = np.float32
@@ -393,7 +396,7 @@ def main(config_file=None, defaults=defaults, **args):
     else: raise ValueError("Unrecognized noise model '%s'" % args['nmat'])
 
     obslists, obskeys, periods, obs_infos = mapmaking.build_obslists(context, args['query'], mode='depth_1', nset=args['nset'], ntod=args['ntod'], tods=args['tods'], freq=args['freq'],per_tube=True)
-    
+
     for oi in range(comm_inter.rank, len(obskeys), comm_inter.size):
         pid, detset, band = obskeys[oi]
         obslist = obslists[obskeys[oi]]
@@ -457,4 +460,5 @@ def main(config_file=None, defaults=defaults, **args):
     return True
 
 if __name__ == '__main__':
+    from sotodlib.site_pipeline import util
     util.main_launcher(main, get_parser)