diff --git a/sotodlib/coords/planets.py b/sotodlib/coords/planets.py
index 5b2e4ad83..8089084ab 100644
--- a/sotodlib/coords/planets.py
+++ b/sotodlib/coords/planets.py
@@ -151,6 +151,7 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
     # Get reference elevation...
     el = np.median(tod.boresight.el[::10])
     az = np.median(tod.boresight.az[::10])
+    roll = np.median(tod.boresight.roll[::10])
     t = (tod.timestamps[0] + tod.timestamps[-1]) / 2
     if isinstance(planet, (list, tuple)):
         _, ra, dec = planet
@@ -159,15 +160,15 @@ def get_scan_q(tod, planet, boresight_offset=None, refq=None):
     else:
         planet = SlowSource.for_named_source(planet, t)
 
-    def scan_q_model(t, az, el, planet):
+    def scan_q_model(t, az, el, roll, planet):
         csl = so3g.proj.CelestialSightLine.az_el(
-            t, az, el, weather='typical', site='so')
+           t, az, el, roll=roll, weather='typical', site='so')
         ra0, dec0 = planet.pos(t)
         return csl.Q, ~so3g.proj.quat.rotation_lonlat(ra0, dec0) * csl.Q * q_xieta
 
     def distance(p):
         dt, daz = p
-        q, qnet = scan_q_model(t + dt, az + daz, el, planet)
+        q, qnet = scan_q_model(t + dt, az + daz, el, roll, planet)
         lon, lat, phi = so3g.proj.quat.decompose_lonlat(qnet)
         return 90 * coords.DEG - lat
 
@@ -177,7 +178,7 @@ def distance(p):
     if warnflag != 0:
         logger.warning('Source-scan solver failed to converge or otherwise '
                        f'complained!  warnflag={warnflag}')
-    q, qnet = scan_q_model(t+p[0], az+p[1], el, planet)
+    q, qnet = scan_q_model(t+p[0], az+p[1], el, roll, planet)
     psi = so3g.proj.quat.decompose_xieta(qnet)[2][0]
     ra, dec = planet.pos(t+p[0])
     rot = ~so3g.proj.quat.rotation_lonlat(ra, dec, psi)
@@ -591,19 +592,28 @@ def compute_source_flags(tod=None, P=None, mask=None, wrap=None,
       the center of the circle, all in degrees.
 
     """
-    if P is None:
-        logger.info('Getting Projection Matrix ...')
-        P, X = get_scan_P(tod, center_on, res=res, comps='T')
-        shape, wcs = tuple(P.geom)
-        if shape[0] * shape[1] > max_pix:
-            raise ValueError(f'Mask map too large: {shape}')
-
     if isinstance(mask, str):
         # Assume it's a filename, and file is simple columns of (x, y,
         # radius) in deg.  (Deprecated!)
         mask = [{'xyr': list(map(float, line.split()))}
                 for line in open(mask)]
 
+    if P is None:
+        x, y, r = mask['xyr']
+        shape = 3 * np.ceil(r * coords.DEG / res).astype(int)
+        # ensure odd dimensions
+        if shape % 2 == 0:
+            shape += 1
+        if shape**2 > max_pix:
+            raise ValueError(f'Mask map too large: {shape}')
+        logger.info('Getting Projection Matrix ...')
+        P, X = get_scan_P(tod, center_on, res=res, comps='T')
+
+        geom = enmap.geometry(np.array((0, 0)), shape=(shape, shape),
+                              proj='tan', res=(res, -res))
+        P.geom = enmap.Geometry(*geom)
+        shape, wcs = tuple(P.geom)
+
     mask_map = P.zeros()
     _add_to_mask(mask, mask_map)
     a = P.from_map(mask_map)