Image difference with different filters (HSC-R and R2)

Support
, , ,
#1

Hi all,

We are having an issue to perform difference imaging with HSC R and R2 bands.
We are trying to look at transients in the 2016-2017 period, and we ideally want to form a template with good seeing pre-dating these data.
We thus generated a template using a coadd of the visits of pointing 01172 (2015-03-18) for filter R.

The issue is that the R filter was switched to R2 in-between.
When I try to perform image differencing (using diapipe’s imageDifferenceDriver.py), I get errors like

39112 WARN 2020-08-30T15:38:22.965+0200 imageDifferenceDriver.imageDifference.getTemplate ({'visit': 91552, 'pointing': 1793, 'filter': 'HSC-R2', 'ccd': 57, 'field': 'SSP_UDEEP_COSMOS', 'dateObs': '2016-11-28', 'taiObs': '2016-11-28', 'expTime': 30.0})(getTemplate.py:270)- deepCoadd_sub, tract=9813, patch=40 does not exist

And indeed, if I try in a notebook:

import lsst.ip.diffim as I
mytask=I.GetCoaddAsTemplateTask()
calexprepo = '/sps/lsst/HSC/bracine/Pipeline_runs/SSP_UDEEP_COSMOS_20161101_20170501/DATA/rerun/coaddDriver_R_01172/'
b = Butler(calexprepo)
skymap = b.get('deepCoadd_skyMap')
exposure = b.get('calexp',dataId={'visit': 91552, 'pointing': 1793, 'filter': 'HSC-R2', 'ccd': 1})
tractInfo, patchList, skyCorners = mytask.getOverlapPatchList(exposure, skymap)
# returns an instance of lsst.daf.persistence.butlerSubset.ButlerSubset
availableCoaddRefs = dict()
for patchInfo in patchList:
    patchNumber = tractInfo.getSequentialPatchIndex(patchInfo)
    patchArgDict = dict(
        datasetType=mytask.getCoaddDatasetName() + "_sub",
        bbox=patchInfo.getOuterBBox(),
        tract=tractInfo.getId(),
        patch="%s,%s" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),
        numSubfilters=mytask.config.numSubfilters,
    )
sensorRef = b.dataRef('calexp', dataId={'visit': 91552, 'pointing': 1793, 'filter': 'HSC-R2', 'ccd': 1})
sensorRef.datasetExists(**patchArgDict)

I get False, whereas for the same test on Z band, I get True.

I assume it is because the butler is looking for a coadd in the R2 band, which doesn’t exist.

Is there a way to go around this?
I saw that if I add
patchArgDict[‘filter’]=‘HSC-R’
then
sensorRef.datasetExists(**patchArgDict)
becomes True

Also is there a recommendation not to perform image differencing between R and R2 or I and I2?

Best,

Ben

#2

I suspect there’s no way to get around this without making at least some modifications to the task itself, though there may be lower-level interfaces you could use that let you provide the template to use more directly. I don’t know much about ImageDifferenceTask personally, so I hope someone who does can chime in on that possibility.

As for your last question:

Also is there a recommendation not to perform image differencing between R and R2 or I and I2?

There are significant differences in transmission as a function of wavelength between the versions of the two filters, so there’s at least some danger of getting some “spurious” detections that are really just static objects with SEDs that are sensitive to those differences.

But The R2 and I2 filters exist because the transmission curves of the R and I filters depend strongly on radius, and those internal differences are about as big as the average different between R and R2 or I and I2. So if you really want to avoid problems, you’d need to stick to just R2 and I2, or only work with R and I data with very small dithers (as is true the HSC SSP UltraDeep fields).

I don’t know of any quantitative analysis of the impact of the HSC transmission curve differences on transient detection, though; this is all just making guesses based on what the transmission curves look like.

#3

This is probably a highly unsatisfactory answer, but the gist is that ImageDifferenceTask currently doesn’t support using a template and a science image which have different filters. There is probably a way to hack getTemplateTask in ip_diffim to accept a mismatch filter, but I got around this myself by just building an R2 and an I2 template instead of trying to force it to accept R and I. Half because it was easier and half because of the actual physical filter differences Jim mentions above.

If you do decide to hack around and come up with something that accepts different filers, please do let us know - I can certainly see this being a reasonable use case for various situations, but it’s not at the top of anyone’s priority list for now.

#4

Hi @racine, I’m just following up on Support category posts without a marked solution, and was wondering – did @jbosch or @mrawls’s suggestions of modifying the task parameters or building an R2 filter resolve your issue? Or perhaps did you implement a different solution?

#5

Hi @MelissaGraham,
Sorry I only saw your question.

so I tried modifying getTemplate from ip_diffim (with version v20.0.0 I am using) and it worked (see below git diff)

diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py
index fce609a..6deec56 100644
--- a/python/lsst/ip/diffim/getTemplate.py
+++ b/python/lsst/ip/diffim/getTemplate.py
@@ -110,10 +110,28 @@ class GetCoaddAsTemplateTask(pipeBase.Task):
                 patch="%s,%s" % (patchInfo.getIndex()[0], patchInfo.getIndex()[1]),
                 numSubfilters=self.config.numSubfilters,
             )
-
             if sensorRef.datasetExists(**patchArgDict):
                 self.log.info("Reading patch %s" % patchArgDict)
                 availableCoaddRefs[patchNumber] = patchArgDict
+            elif exposure.getInfo().getFilter().getAliases()[1]=='HSC-R2':
+                # note that for some reason, the getAliases returns ['r', 'HSC-R2'], so we need to look at [1], see next
+
+                templateIdList={'filter':'HSC-R'}
+                patchArgDict.update(templateIdList)
+                
+                if sensorRef.datasetExists(**patchArgDict):
+                    self.log.info("Reading patch %s" % patchArgDict)
+                    self.log.info("Note that we used the hack that switches from R2 to R for the template")
+
+                    availableCoaddRefs[patchNumber] = patchArgDict
+            elif exposure.getInfo().getFilter().getAliases()[0]=='HSC-I2':
+                # note that for some reason, the getAliases returns ['HSC-I2','i'], so we need to look at [0]
+                templateIdList={'filter':'HSC-I'}
+                patchArgDict.update(templateIdList)
+                if sensorRef.datasetExists(**patchArgDict):
+                    self.log.info("Reading patch %s" % patchArgDict)
+                    self.log.info("Note that we used the hack that switches from I2 to I for the template")
+                    availableCoaddRefs[patchNumber] = patchArgDict
 
         templateExposure = self.run(
             tractInfo, patchList, skyCorners, availableCoaddRefs,

This is of course not the best way to do that, but it did the job.

1 Like
#6

Hi again,

I just wanted to switch (back) to using ap_pipe, mostly for the use of a database, and I now have an issue with the same I/I2 and R/R2 problem.
I guess I can just continue the discussion on the present posting.

Using the v20.0.0 of ap_association, I had to modify filterName from mapApData.py, as is now done in the latset version:
filterName = exposure.getFilter().getCanonicalName()[0]
to avoid the error: sqlalchemy.exc.OperationalError: (sqlite3.OperationalError) table DiaObject has no column named r2PSFluxLinearSlope.

But I now get a sqlite error
FATAL 2020-11-12T22:23:51.594+0100 apPipe ({'visit': 91554, 'ccd': 42, 'filter': 'HSC-R2', 'field': 'SSP_UDEEP_COSMOS', 'dateObs': '2016-11-28', 'pointing': 1793, 'taiObs': '2016-11-28', 'expTime': 300.0})(cmdLineTask.py:397)- Failed on dataId={'visit': 91554, 'ccd': 42, 'filter': 'HSC-R2', 'field': 'SSP_UDEEP_COSMOS', 'dateObs': '2016-11-28', 'pointing': 1793, 'taiObs': '2016-11-28', 'expTime': 300.0}: IntegrityError: (sqlite3.IntegrityError) UNIQUE constraint failed: DiaSource.diaSourceId
which I don’t get if I run that hacked code on the “HSC-Z” band.

I assume this is still an issue due to the fact that the template and visit don’t have the same band name.
Is there an obvious solution I missed here?

#7

If you’re reasonably sure the reason the unique constraint is failing is because it can’t find any matching HSC-R2 data (only HSC-R data exists), there is currently no obvious solution. You should consider building templates in HSC-R2, or possibly working out a way to lie to the Butler that the existing templates are HSC-R2. I don’t have actual guidance for how to do the lying bit, but it would involve some ugly registry and/or file path hacking at the ingestion phase followed by reprocessing.

Since we are in the process of switching our filter system to distinguish the concept of “band” from “physical filter,” I think in that new reality it will be possible to do difference imaging with matching bands (e.g., “r”) but different filters (e.g., “HSC-R” and “HSC-R2”). RFC-730 describes this.

2 Likes
#8

Thanks so much for posting the solution, @racine, that could help any others who find themselves in a similar situation.