diff --git a/pipeline/coaddition.py b/pipeline/coaddition.py index 169e8ad9..d1d510e4 100644 --- a/pipeline/coaddition.py +++ b/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,11 +602,19 @@ 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}" ) @@ -593,6 +622,8 @@ def _coadd_swarp( self, 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 ) diff --git a/tests/pipeline/test_coaddition.py b/tests/pipeline/test_coaddition.py index 43d2dcaa..a195829c 100644 --- a/tests/pipeline/test_coaddition.py +++ b/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. ) diff --git a/tests/seechange_config_test.yaml b/tests/seechange_config_test.yaml index bdcfc8f0..4b268c2c 100644 --- a/tests/seechange_config_test.yaml +++ b/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