LSST Applications g0da5cf3356+25b44625d0,g17e5ecfddb+50a5ac4092,g1c76d35bf8+585f0f68a2,g295839609d+8ef6456700,g2e2c1a68ba+cc1f6f037e,g38293774b4+62d12e78cb,g3b44f30a73+2891c76795,g48ccf36440+885b902d19,g4b2f1765b6+0c565e8f25,g5320a0a9f6+bd4bf1dc76,g56364267ca+403c24672b,g56b687f8c9+585f0f68a2,g5c4744a4d9+78cd207961,g5ffd174ac0+bd4bf1dc76,g6075d09f38+3075de592a,g667d525e37+cacede5508,g6f3e93b5a3+da81c812ee,g71f27ac40c+cacede5508,g7212e027e3+eb621d73aa,g774830318a+18d2b9fa6c,g7985c39107+62d12e78cb,g79ca90bc5c+fa2cc03294,g881bdbfe6c+cacede5508,g91fc1fa0cf+82a115f028,g961520b1fb+2534687f64,g96f01af41f+f2060f23b6,g9ca82378b8+cacede5508,g9d27549199+78cd207961,gb065e2a02a+ad48cbcda4,gb1df4690d6+585f0f68a2,gb35d6563ee+62d12e78cb,gbc3249ced9+bd4bf1dc76,gbec6a3398f+bd4bf1dc76,gd01420fc67+bd4bf1dc76,gd59336e7c4+c7bb92e648,gf46e8334de+81c9a61069,gfed783d017+bd4bf1dc76,v25.0.1.rc3
LSST Data Management Base Package
|
detectorId = exposure.getInfo().getDetector().getId() if finalizedPsfApCorrCatalog is not None: row = finalizedPsfApCorrCatalog.find(detectorId) if row is None: self.log.warning("Detector id %s not found in finalizedPsfApCorrCatalog; " "Using original psf.", detectorId) else: psf = row.getPsf() apCorrMap = row.getApCorrMap() if psf is None or apCorrMap is None: self.log.warning("Detector id %s has None for psf/apCorrMap in " "finalizedPsfApCorrCatalog; Using original psf.", detectorId) else: exposure.setPsf(psf) exposure.info.setApCorrMap(apCorrMap) return exposure @timeMethod def run(self, exposure=None, selectSources=None, templateExposure=None, templateSources=None, idFactory=None, calexpBackgroundExposure=None, subtractedExposure=None):
subtractRes = None controlSources = None subtractedExposure = None scoreExposure = None diaSources = None kernelSources = None # We'll clone exposure if modified but will still need the original exposureOrig = exposure if self.config.doAddCalexpBackground: mi = exposure.getMaskedImage() mi += calexpBackgroundExposure.getImage() if not exposure.hasPsf(): raise pipeBase.TaskError("Exposure has no psf") sciencePsf = exposure.getPsf() if self.config.doSubtract: if self.config.doScaleTemplateVariance: self.log.info("Rescaling template variance") templateVarFactor = self.scaleVariance.run( templateExposure.getMaskedImage()) self.log.info("Template variance scaling factor: %.2f", templateVarFactor) self.metadata.add("scaleTemplateVarianceFactor", templateVarFactor) self.metadata.add("psfMatchingAlgorithm", self.config.subtract.name) if self.config.subtract.name == 'zogy': subtractRes = self.subtract.run(exposure, templateExposure, doWarping=True) scoreExposure = subtractRes.scoreExp subtractedExposure = subtractRes.diffExp subtractRes.subtractedExposure = subtractedExposure subtractRes.matchedExposure = None elif self.config.subtract.name == 'al': # compute scienceSigmaOrig: sigma of PSF of science image before pre-convolution # Just need a rough estimate; average positions are fine sciAvgPos = sciencePsf.getAveragePosition() scienceSigmaOrig = sciencePsf.computeShape(sciAvgPos).getDeterminantRadius() templatePsf = templateExposure.getPsf() templateAvgPos = templatePsf.getAveragePosition() templateSigma = templatePsf.computeShape(templateAvgPos).getDeterminantRadius() # if requested, convolve the science exposure with its PSF # (properly, this should be a cross-correlation, but our code does not yet support that) # compute scienceSigmaPost: sigma of science exposure with pre-convolution, if done, # else sigma of original science exposure # TODO: DM-22762 This functional block should be moved into its own method preConvPsf = None if self.config.useScoreImageDetection: self.log.warning("AL likelihood image: pre-convolution of PSF is not implemented.") convControl = afwMath.ConvolutionControl() # cannot convolve in place, so need a new image anyway srcMI = exposure.maskedImage exposure = exposure.clone() # New deep copy srcPsf = sciencePsf if self.config.useGaussianForPreConvolution: self.log.info( "AL likelihood image: Using Gaussian (sigma=%.2f) PSF estimation " "for science image pre-convolution", scienceSigmaOrig) # convolve with a simplified PSF model: a double Gaussian kWidth, kHeight = sciencePsf.getLocalKernel( sciencePsf.getAveragePosition() ).getDimensions() preConvPsf = SingleGaussianPsf(kWidth, kHeight, scienceSigmaOrig) else: # convolve with science exposure's PSF model self.log.info( "AL likelihood image: Using the science image PSF for pre-convolution.") preConvPsf = srcPsf afwMath.convolve( exposure.maskedImage, srcMI, preConvPsf.getLocalKernel(preConvPsf.getAveragePosition()), convControl ) scienceSigmaPost = scienceSigmaOrig*math.sqrt(2) else: scienceSigmaPost = scienceSigmaOrig # If requested, find and select sources from the image # else, AL subtraction will do its own source detection # TODO: DM-22762 This functional block should be moved into its own method if self.config.doSelectSources: if selectSources is None: self.log.warning("Src product does not exist; running detection, measurement," " selection") # Run own detection and measurement; necessary in nightly processing selectSources = self.subtract.getSelectSources( exposure, sigma=scienceSigmaPost, doSmooth=not self.config.useScoreImageDetection, idFactory=idFactory, ) if self.config.doAddMetrics: # Number of basis functions nparam = len(makeKernelBasisList(self.subtract.config.kernel.active, referenceFwhmPix=scienceSigmaPost*FwhmPerSigma, targetFwhmPix=templateSigma*FwhmPerSigma)) # Modify the schema of all Sources # DEPRECATED: This is a data dependent (nparam) output product schema # outside the task constructor. # NOTE: The pre-determination of nparam at this point # may be incorrect as the template psf is warped later in # ImagePsfMatchTask.matchExposures() kcQa = KernelCandidateQa(nparam) selectSources = kcQa.addToSchema(selectSources) if self.config.kernelSourcesFromRef: # match exposure sources to reference catalog astromRet = self.astrometer.loadAndMatch(exposure=exposure, sourceCat=selectSources) matches = astromRet.matches elif templateSources: # match exposure sources to template sources mc = afwTable.MatchControl() mc.findOnlyClosest = False matches = afwTable.matchRaDec(templateSources, selectSources, 1.0*geom.arcseconds, mc) else: raise RuntimeError("doSelectSources=True and kernelSourcesFromRef=False," "but template sources not available. Cannot match science " "sources with template sources. Run process* on data from " "which templates are built.") kernelSources = self.sourceSelector.run(selectSources, exposure=exposure, matches=matches).sourceCat random.shuffle(kernelSources, random.random) controlSources = kernelSources[::self.config.controlStepSize] kernelSources = [k for i, k in enumerate(kernelSources) if i % self.config.controlStepSize] if self.config.doSelectDcrCatalog: redSelector = DiaCatalogSourceSelectorTask( DiaCatalogSourceSelectorConfig(grMin=self.sourceSelector.config.grMax, grMax=99.999)) redSources = redSelector.selectStars(exposure, selectSources, matches=matches).starCat controlSources.extend(redSources) blueSelector = DiaCatalogSourceSelectorTask( DiaCatalogSourceSelectorConfig(grMin=-99.999, grMax=self.sourceSelector.config.grMin)) blueSources = blueSelector.selectStars(exposure, selectSources, matches=matches).starCat controlSources.extend(blueSources) if self.config.doSelectVariableCatalog: varSelector = DiaCatalogSourceSelectorTask( DiaCatalogSourceSelectorConfig(includeVariable=True)) varSources = varSelector.selectStars(exposure, selectSources, matches=matches).starCat controlSources.extend(varSources) self.log.info("Selected %d / %d sources for Psf matching (%d for control sample)", len(kernelSources), len(selectSources), len(controlSources)) allresids = {} # TODO: DM-22762 This functional block should be moved into its own method if self.config.doUseRegister: self.log.info("Registering images") if templateSources is None: # Run detection on the template, which is # temporarily background-subtracted # sigma of PSF of template image before warping templateSources = self.subtract.getSelectSources( templateExposure, sigma=templateSigma, doSmooth=True, idFactory=idFactory ) # Third step: we need to fit the relative astrometry. # wcsResults = self.fitAstrometry(templateSources, templateExposure, selectSources) warpedExp = self.register.warpExposure(templateExposure, wcsResults.wcs, exposure.getWcs(), exposure.getBBox()) templateExposure = warpedExp # Create debugging outputs on the astrometric # residuals as a function of position. Persistence # not yet implemented; expected on (I believe) #2636. if self.config.doDebugRegister: # Grab matches to reference catalog srcToMatch = {x.second.getId(): x.first for x in matches} refCoordKey = wcsResults.matches[0].first.getTable().getCoordKey() inCentroidKey = wcsResults.matches[0].second.getTable().getCentroidSlot().getMeasKey() sids = [m.first.getId() for m in wcsResults.matches] positions = [m.first.get(refCoordKey) for m in wcsResults.matches] residuals = [m.first.get(refCoordKey).getOffsetFrom(wcsResults.wcs.pixelToSky( m.second.get(inCentroidKey))) for m in wcsResults.matches] allresids = dict(zip(sids, zip(positions, residuals))) cresiduals = [m.first.get(refCoordKey).getTangentPlaneOffset( wcsResults.wcs.pixelToSky( m.second.get(inCentroidKey))) for m in wcsResults.matches] colors = numpy.array([-2.5*numpy.log10(srcToMatch[x].get("g")) + 2.5*numpy.log10(srcToMatch[x].get("r")) for x in sids if x in srcToMatch.keys()]) dlong = numpy.array([r[0].asArcseconds() for s, r in zip(sids, cresiduals) if s in srcToMatch.keys()]) dlat = numpy.array([r[1].asArcseconds() for s, r in zip(sids, cresiduals) if s in srcToMatch.keys()]) idx1 = numpy.where(colors < self.sourceSelector.config.grMin) idx2 = numpy.where((colors >= self.sourceSelector.config.grMin) & (colors <= self.sourceSelector.config.grMax)) idx3 = numpy.where(colors > self.sourceSelector.config.grMax) rms1Long = IqrToSigma*( (numpy.percentile(dlong[idx1], 75) - numpy.percentile(dlong[idx1], 25))) rms1Lat = IqrToSigma*(numpy.percentile(dlat[idx1], 75) - numpy.percentile(dlat[idx1], 25)) rms2Long = IqrToSigma*( (numpy.percentile(dlong[idx2], 75) - numpy.percentile(dlong[idx2], 25))) rms2Lat = IqrToSigma*(numpy.percentile(dlat[idx2], 75) - numpy.percentile(dlat[idx2], 25)) rms3Long = IqrToSigma*( (numpy.percentile(dlong[idx3], 75) - numpy.percentile(dlong[idx3], 25))) rms3Lat = IqrToSigma*(numpy.percentile(dlat[idx3], 75) - numpy.percentile(dlat[idx3], 25)) self.log.info("Blue star offsets'': %.3f %.3f, %.3f %.3f", numpy.median(dlong[idx1]), rms1Long, numpy.median(dlat[idx1]), rms1Lat) self.log.info("Green star offsets'': %.3f %.3f, %.3f %.3f", numpy.median(dlong[idx2]), rms2Long, numpy.median(dlat[idx2]), rms2Lat) self.log.info("Red star offsets'': %.3f %.3f, %.3f %.3f", numpy.median(dlong[idx3]), rms3Long, numpy.median(dlat[idx3]), rms3Lat) self.metadata.add("RegisterBlueLongOffsetMedian", numpy.median(dlong[idx1])) self.metadata.add("RegisterGreenLongOffsetMedian", numpy.median(dlong[idx2])) self.metadata.add("RegisterRedLongOffsetMedian", numpy.median(dlong[idx3])) self.metadata.add("RegisterBlueLongOffsetStd", rms1Long) self.metadata.add("RegisterGreenLongOffsetStd", rms2Long) self.metadata.add("RegisterRedLongOffsetStd", rms3Long) self.metadata.add("RegisterBlueLatOffsetMedian", numpy.median(dlat[idx1])) self.metadata.add("RegisterGreenLatOffsetMedian", numpy.median(dlat[idx2])) self.metadata.add("RegisterRedLatOffsetMedian", numpy.median(dlat[idx3])) self.metadata.add("RegisterBlueLatOffsetStd", rms1Lat) self.metadata.add("RegisterGreenLatOffsetStd", rms2Lat) self.metadata.add("RegisterRedLatOffsetStd", rms3Lat) # warp template exposure to match exposure, # PSF match template exposure to exposure, # then return the difference # Return warped template... Construct sourceKernelCand list after subtract self.log.info("Subtracting images") subtractRes = self.subtract.subtractExposures( templateExposure=templateExposure, scienceExposure=exposure, candidateList=kernelSources, convolveTemplate=self.config.convolveTemplate, doWarping=not self.config.doUseRegister ) if self.config.useScoreImageDetection: scoreExposure = subtractRes.subtractedExposure else: subtractedExposure = subtractRes.subtractedExposure if self.config.doDetection: self.log.info("Computing diffim PSF") # Get Psf from the appropriate input image if it doesn't exist if subtractedExposure is not None and not subtractedExposure.hasPsf(): if self.config.convolveTemplate: subtractedExposure.setPsf(exposure.getPsf()) else: subtractedExposure.setPsf(templateExposure.getPsf()) # If doSubtract is False, then subtractedExposure was fetched from disk (above), # thus it may have already been decorrelated. Thus, we do not decorrelate if # doSubtract is False. # NOTE: At this point doSubtract == True if self.config.doDecorrelation and self.config.doSubtract: preConvKernel = None if self.config.useGaussianForPreConvolution: preConvKernel = preConvPsf.getLocalKernel(preConvPsf.getAveragePosition()) if self.config.useScoreImageDetection: scoreExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure, scoreExposure, subtractRes.psfMatchingKernel, spatiallyVarying=self.config.doSpatiallyVarying, preConvKernel=preConvKernel, templateMatched=True, preConvMode=True).correctedExposure # Note that the subtracted exposure is always decorrelated, # even if the score image is used for detection subtractedExposure = self.decorrelate.run(exposureOrig, subtractRes.warpedExposure, subtractedExposure, subtractRes.psfMatchingKernel, spatiallyVarying=self.config.doSpatiallyVarying, preConvKernel=None, templateMatched=self.config.convolveTemplate, preConvMode=False).correctedExposure # END (if subtractAlgorithm == 'AL') # END (if self.config.doSubtract) if self.config.doDetection: self.log.info("Running diaSource detection") # subtractedExposure - reserved for task return value # in zogy, it is always the proper difference image # in AL, it may be (yet) pre-convolved and/or decorrelated # # detectionExposure - controls which exposure to use for detection # in-place modifications will appear in task return if self.config.useScoreImageDetection: # zogy with score image detection enabled self.log.info("Detection, diffim rescaling and measurements are " "on AL likelihood or Zogy score image.") detectionExposure = scoreExposure else: # AL or zogy with no score image detection detectionExposure = subtractedExposure # Rescale difference image variance plane if self.config.doScaleDiffimVariance: self.log.info("Rescaling diffim variance") diffimVarFactor = self.scaleVariance.run(detectionExposure.getMaskedImage()) self.log.info("Diffim variance scaling factor: %.2f", diffimVarFactor) self.metadata.add("scaleDiffimVarianceFactor", diffimVarFactor) # Erase existing detection mask planes mask = detectionExposure.getMaskedImage().getMask() mask &= ~(mask.getPlaneBitMask("DETECTED") | mask.getPlaneBitMask("DETECTED_NEGATIVE")) table = afwTable.SourceTable.make(self.schema, idFactory) table.setMetadata(self.algMetadata) results = self.detection.run( table=table, exposure=detectionExposure, doSmooth=not self.config.useScoreImageDetection ) if self.config.doMerge: fpSet = results.fpSets.positive fpSet.merge(results.fpSets.negative, self.config.growFootprint, self.config.growFootprint, False) diaSources = afwTable.SourceCatalog(table) fpSet.makeSources(diaSources) self.log.info("Merging detections into %d sources", len(diaSources)) else: diaSources = results.sources # Inject skySources before measurement. if self.config.doSkySources: skySourceFootprints = self.skySources.run( mask=detectionExposure.mask, seed=detectionExposure.info.id) if skySourceFootprints: for foot in skySourceFootprints: s = diaSources.addNew() s.setFootprint(foot) s.set(self.skySourceKey, True) if self.config.doMeasurement: newDipoleFitting = self.config.doDipoleFitting self.log.info("Running diaSource measurement: newDipoleFitting=%r", newDipoleFitting) if not newDipoleFitting: # Just fit dipole in diffim self.measurement.run(diaSources, detectionExposure) else: # Use (matched) template and science image (if avail.) to constrain dipole fitting if self.config.doSubtract and 'matchedExposure' in subtractRes.getDict(): self.measurement.run(diaSources, detectionExposure, exposure, subtractRes.matchedExposure) else: self.measurement.run(diaSources, detectionExposure, exposure) if self.config.doApCorr: self.applyApCorr.run( catalog=diaSources, apCorrMap=detectionExposure.getInfo().getApCorrMap() ) if self.config.doForcedMeasurement: # Run forced psf photometry on the PVI at the diaSource locations. # Copy the measured flux and error into the diaSource. forcedSources = self.forcedMeasurement.generateMeasCat( exposure, diaSources, detectionExposure.getWcs()) self.forcedMeasurement.run(forcedSources, exposure, diaSources, detectionExposure.getWcs()) mapper = afwTable.SchemaMapper(forcedSources.schema, diaSources.schema) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFlux")[0], "ip_diffim_forced_PsfFlux_instFlux", True) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_instFluxErr")[0], "ip_diffim_forced_PsfFlux_instFluxErr", True) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_area")[0], "ip_diffim_forced_PsfFlux_area", True) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag")[0], "ip_diffim_forced_PsfFlux_flag", True) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_noGoodPixels")[0], "ip_diffim_forced_PsfFlux_flag_noGoodPixels", True) mapper.addMapping(forcedSources.schema.find("base_PsfFlux_flag_edge")[0], "ip_diffim_forced_PsfFlux_flag_edge", True) for diaSource, forcedSource in zip(diaSources, forcedSources): diaSource.assign(forcedSource, mapper) # Match with the calexp sources if possible if self.config.doMatchSources: if selectSources is not None: # Create key,val pair where key=diaSourceId and val=sourceId matchRadAsec = self.config.diaSourceMatchRadius matchRadPixel = matchRadAsec/exposure.getWcs().getPixelScale().asArcseconds() srcMatches = afwTable.matchXy(selectSources, diaSources, matchRadPixel) srcMatchDict = dict([(srcMatch.second.getId(), srcMatch.first.getId()) for srcMatch in srcMatches]) self.log.info("Matched %d / %d diaSources to sources", len(srcMatchDict), len(diaSources)) else: self.log.warning("Src product does not exist; cannot match with diaSources") srcMatchDict = {} # Create key,val pair where key=diaSourceId and val=refId refAstromConfig = AstrometryConfig() refAstromConfig.matcher.maxMatchDistArcSec = matchRadAsec refAstrometer = AstrometryTask(self.refObjLoader, config=refAstromConfig) astromRet = refAstrometer.run(exposure=exposure, sourceCat=diaSources) refMatches = astromRet.matches if refMatches is None: self.log.warning("No diaSource matches with reference catalog") refMatchDict = {} else: self.log.info("Matched %d / %d diaSources to reference catalog", len(refMatches), len(diaSources)) refMatchDict = dict([(refMatch.second.getId(), refMatch.first.getId()) for refMatch in refMatches]) # Assign source Ids for diaSource in diaSources: sid = diaSource.getId() if sid in srcMatchDict: diaSource.set("srcMatchId", srcMatchDict[sid]) if sid in refMatchDict: diaSource.set("refMatchId", refMatchDict[sid]) if self.config.doAddMetrics and self.config.doSelectSources: self.log.info("Evaluating metrics and control sample") kernelCandList = [] for cell in subtractRes.kernelCellSet.getCellList(): for cand in cell.begin(False): # include bad candidates kernelCandList.append(cand) # Get basis list to build control sample kernels basisList = kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelList() nparam = len(kernelCandList[0].getKernel(KernelCandidateF.ORIG).getKernelParameters()) controlCandList = ( diffimTools.sourceTableToCandidateList(controlSources, subtractRes.warpedExposure, exposure, self.config.subtract.kernel.active, self.config.subtract.kernel.active.detectionConfig, self.log, doBuild=True, basisList=basisList)) KernelCandidateQa.apply(kernelCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel, dof=nparam) KernelCandidateQa.apply(controlCandList, subtractRes.psfMatchingKernel, subtractRes.backgroundModel) if self.config.doDetection: KernelCandidateQa.aggregate(selectSources, self.metadata, allresids, diaSources) else: KernelCandidateQa.aggregate(selectSources, self.metadata, allresids) self.runDebug(exposure, subtractRes, selectSources, kernelSources, diaSources) return pipeBase.Struct( subtractedExposure=subtractedExposure, scoreExposure=scoreExposure, warpedExposure=subtractRes.warpedExposure, matchedExposure=subtractRes.matchedExposure, subtractRes=subtractRes, diaSources=diaSources, selectSources=selectSources ) def fitAstrometry(self, templateSources, templateExposure, selectSources):
Definition at line 1312 of file imageDifference.py.