diff --git a/sotodlib/coords/__init__.py b/sotodlib/coords/__init__.py index d29e0373c..f72b0cebe 100644 --- a/sotodlib/coords/__init__.py +++ b/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 diff --git a/sotodlib/coords/helpers.py b/sotodlib/coords/helpers.py index 9717bf892..e297cc9df 100644 --- a/sotodlib/coords/helpers.py +++ b/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 diff --git a/sotodlib/site_pipeline/make_depth1_map.py b/sotodlib/site_pipeline/make_depth1_map.py index bf7dacca8..59ee73f90 100644 --- a/sotodlib/site_pipeline/make_depth1_map.py +++ b/sotodlib/site_pipeline/make_depth1_map.py @@ -167,15 +167,28 @@ 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] + # Cap xshape to nphi. To see why, consider this example: + # Sky width: 100 + # box 0: 0 30 + # box 1: -40 10 + # box 2: 20 70 + # box 3: 45 110 + # union: -40 110: Wider than the whole sky! + # But since we use union_pixbox[0] as the zero-pixel in + # our output geometry, this overflow just results in + # unhittable pixels for x >= nphi, which we can just chop off here + shape[-1] = min(shape[-1], nphi) 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 +253,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 +375,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 +407,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 +471,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)