swarp coaddition seems to be working. Helps to use the right sematics…

… for running swarp.

pipeline/coaddition.py

@@ -550,7 +550,7 @@ def _coadd_swarp( self,
                 wtdata = ds.image.weight * 10 ** ( ( ds.zp.zp - targds.zp.zp ) / 1.25 )
                 # Make sure weight is 0 for all bad pixels
                 # (This is what swarp expects.)
-                wtdata[ ds.image.flags ] = 0.
+                wtdata[ ds.image.flags != 0 ] = 0.
 
                 tmpim = tmpdir / f'in{imgdex:03d}_image.fits'
                 tmpwt = tmpdir / f'in{imgdex:03d}_weight.fits'
@@ -567,8 +567,29 @@ def _coadd_swarp( self,
             if try_to_reduce_memory_usage:
                 targds.free()
 
-            # Use swarp to coadd all the source images, aligning with the target image
-            # TODO : think about what swarp does with gain and saturate!
+            # Use swarp to coadd all the source images, aligning with the target image.
+            #
+            # We don't want to mess with gain multiplication, so make
+            #   sure that swarp doesn't try to do anything funny by
+            #   setting an unlikely keyword.  Because our weights
+            #   already include noise from objects, not just from the
+            #   sky background, we want gain=0 for swarp to do the
+            #   right thing.  Likewise, saturated pixels should already
+            #   be marked as such from our preprocessing (I really
+            #   hope), so try to avoid letting swarp doing things there
+            #   too.  Also make sure swarp doesn't try to do fscaling,
+            #   since we already scaled the images to zeorpoints.
+            #   (Never know what's in the header!)
+            #
+            # The swarp manual is incomplete.  I don't know what the
+            #   CLIP_SIGMA et al. parameters do.  Are they only used
+            #   with COMBINE_TYPE=CLIPPED?  We really want to do a
+            #   weighted combination, but also we want to do some
+            #   clipping to reject CRs and the like.
+            #   ...Looking at the swarp source code, it looks like CLIPPED
+            #   is doing a weighted mean of the things it doesn't
+            #   throw out, so CLIPPED should be good to just use.
+            #   https://github.com/astromatic/swarp/blob/3d8ddf1e61718a2ba402473990c6483862671806/src/coadd.c#L1418
 
             command = [ 'swarp' ]
             command.extend( sumimgs )
@@ -581,18 +602,28 @@ def _coadd_swarp( self,
                               '-VMEM_MAX', '1024',
                               '-MEM_MAX', '1024',
                               '-WRITE_XML', 'N',
+                              '-INTERPOLATE', 'Y',
+                              '-FSCALE_KEYWORD', 'THIS_KEYWORD_WILL_NEVER_EXIST',
+                              '-FSCALE_DEFAULT', '1.0',
+                              '-GAIN_KEYWORD', 'THIS_KEYWORD_WILL_NEVER_EXIT',
+                              '-GAIN_DEFAULT', '0.0',
+                              '-SATLEV_KEYWORD', 'THIS_KEYWORD_WILL_NEVER_EXIST',
+                              '-SATLEV_DEFAULT', '1e10',
                               '-COMBINE', 'Y',
                               '-COMBINE_TYPE', 'CLIPPED',
                               '-WEIGHT_TYPE', 'MAP_WEIGHT',
-                              '-WEIGHT_IMAGE' ] )
-            command.extend( sumwts )
+                              '-WEIGHT_IMAGE', ','.join([ str(s) for s in sumwts ])
+                             ] )
+
 
             # SCLogger.debug( f"Running swarp to coadd {len(sumimgs)} images" )
             SCLogger.debug( f"Running swarp to coadd {len(sumimgs)} images; swarp command is {command}" )
             t0 = time.perf_counter()
             res = subprocess.run( command, capture_output=True, timeout=self.pars.swarp_timeout )
             t1 = time.perf_counter()
             SCLogger.debug( f"Swarp to sum {len(sumimgs)} images took {t1-t0:.2f} seconds" )
+            SCLogger.debug( f"Swarp stdout:\n{res.stdout}" )
+            SCLogger.debug( f"Swarp stderr:\n{res.stderr}" )
             if res.returncode != 0:
                 raise SubprocessFailure( res )
 

tests/pipeline/test_coaddition.py

@@ -547,6 +547,25 @@ def test_coadd_partial_overlap_swarp( decam_four_offset_refs, decam_four_refs_al
     img = coadder.run( data_store_list=decam_four_offset_refs,
                        alignment_target_datastore=decam_four_refs_alignment_target )
 
+    assert img.data.shape == ( 4096, 2048 )
+    assert img.flags.shape == img.data.shape
+    assert img.weight.shape == img.data.shape
 
-    import pdb; pdb.set_trace()
-    pass
+    # What else to check?
+
+    # Spot check a few points on the image.
+    # (I manually looked at the image and picked out a few spots)
+
+    # Check that the weight is higher in a region where two images actually overlapped
+    assert img.weight[ 550:640, 975:1140 ].mean() == pytest.approx( 0.0199, abs=0.0005 )
+    assert img.weight[ 690:770, 930:1050 ].mean() == pytest.approx( 0.0131, abs=0.0005 )
+
+    # Look at a spot with a star, and a nearby sky, in a place where there was only
+    #   one image in the coadd
+    assert img.data[ 3217:3231, 479:491 ].sum() == pytest.approx( 82530., abs=25. )
+    assert img.data[ 3217:3231, 509:521 ].sum() == pytest.approx( 231., abs=25. )
+
+    # Look at a spot with a galaxy and a nearby sky, in a place where there were
+    #   two images in the sum
+    assert img.data[ 237:266, 978:988 ].sum() == pytest.approx( 7918., abs=10. )
+    assert img.data[ 237:266, 1008:1018 ].sum() == pytest.approx( 44., abs=10. )

tests/seechange_config_test.yaml

@@ -6,8 +6,8 @@ augments:
   - local_augments.yaml
 
 path:
-  data_root: 'tests/temp_data'
-  data_temp: 'tests/temp_data'
+  data_root: '/seechange/tests/temp_data'
+  data_temp: '/seechange/tests/temp_data'
 
 db:
   host: postgres