From 9ff2e6d20c7ac8d8940083121c3a116eeb82f3d2 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 20:31:37 -0700 Subject: [PATCH 1/7] Ignore E203 --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.cfg b/setup.cfg index 52616bd1..7bc9eb18 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,7 +1,7 @@ [flake8] max-line-length = 110 max-doc-length = 79 -ignore = E133, E226, E228, N802, N803, N806, N812, N813, N815, N816, W503 +ignore = E133, E203, E226, E228, N802, N803, N806, N812, N813, N815, N816, W503 exclude = bin, doc/conf.py, From 1e3942e01981a043baa79ae4141d6e80db11588e Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 12:50:28 -0700 Subject: [PATCH 2/7] Add a pyproject.toml file --- pyproject.toml | 9 +++++++++ 1 file changed, 9 insertions(+) create mode 100644 pyproject.toml diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 00000000..aa57a70b --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,9 @@ +[tool.isort] +profile = "black" +line_length = 110 + +[tool.black] +line_length = 110 + +[tool.lsst_versions] +write_to = "python/lsst/drp_tasks/version.py" From 120ea6448d8019d7d4f5d23ff76160f1c296befd Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 17:03:48 -0700 Subject: [PATCH 3/7] Add a pre-commit configuration --- .pre-commit-config.yaml | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 .pre-commit-config.yaml diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 00000000..06f77c67 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,17 @@ +repos: + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.5.0 + hooks: + - id: check-yaml + - id: end-of-file-fixer + - id: trailing-whitespace + - repo: https://github.com/psf/black + rev: 23.10.1 + hooks: + - id: black + language_version: python3.11 + - repo: https://github.com/pycqa/isort + rev: 5.12.0 + hooks: + - id: isort + name: isort (python) From 929500d35832c1d3925e8a91aaca136b59652824 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 16:57:49 -0700 Subject: [PATCH 4/7] Apply black to everything --- python/lsst/drp/tasks/assemble_chi2_coadd.py | 37 +- python/lsst/drp/tasks/assemble_coadd.py | 537 ++++++++----- python/lsst/drp/tasks/dcr_assemble_coadd.py | 522 +++++++----- python/lsst/drp/tasks/forcedPhotCoadd.py | 114 +-- python/lsst/drp/tasks/gbdesAstrometricFit.py | 760 ++++++++++-------- python/lsst/drp/tasks/update_visit_summary.py | 48 +- tests/assemble_coadd_test_utils.py | 134 +-- tests/test_assemble_coadd.py | 76 +- tests/test_dcr_assemble_coadd.py | 5 +- tests/test_gbdesAstrometricFit.py | 290 ++++--- tests/test_update_visit_summary.py | 10 +- 11 files changed, 1450 insertions(+), 1083 deletions(-) diff --git a/python/lsst/drp/tasks/assemble_chi2_coadd.py b/python/lsst/drp/tasks/assemble_chi2_coadd.py index 9d01589e..214ac63a 100644 --- a/python/lsst/drp/tasks/assemble_chi2_coadd.py +++ b/python/lsst/drp/tasks/assemble_chi2_coadd.py @@ -55,7 +55,7 @@ def calculateKernelSize(sigma: float, nSigmaForKernel: float = 7) -> int: size: Size of the smoothing kernel. """ - return (int(sigma * nSigmaForKernel + 0.5)//2)*2 + 1 # make sure it is odd + return (int(sigma * nSigmaForKernel + 0.5) // 2) * 2 + 1 # make sure it is odd def convolveImage(image: afwImage.Image, psf: Psf) -> afwImage.Image: @@ -97,16 +97,17 @@ def convolveImage(image: afwImage.Image, psf: Psf) -> afwImage.Image: return convolvedImage.Factory(convolvedImage, bbox, afwImage.PARENT, False) -class AssembleChi2CoaddConnections(pipeBase.PipelineTaskConnections, - dimensions=("tract", "patch", "skymap"), - defaultTemplates={"inputCoaddName": "deep", - "outputCoaddName": "deepChi2"}): +class AssembleChi2CoaddConnections( + pipeBase.PipelineTaskConnections, + dimensions=("tract", "patch", "skymap"), + defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deepChi2"}, +): inputCoadds = cT.Input( doc="Exposure on which to run deblending", name="{inputCoaddName}Coadd_calexp", storageClass="ExposureF", multiple=True, - dimensions=("tract", "patch", "band", "skymap") + dimensions=("tract", "patch", "band", "skymap"), ) chi2Coadd = cT.Output( doc="Chi^2 exposure, produced by merging multiband coadds", @@ -116,21 +117,20 @@ class AssembleChi2CoaddConnections(pipeBase.PipelineTaskConnections, ) -class AssembleChi2CoaddConfig(pipeBase.PipelineTaskConfig, - pipelineConnections=AssembleChi2CoaddConnections): +class AssembleChi2CoaddConfig(pipeBase.PipelineTaskConfig, pipelineConnections=AssembleChi2CoaddConnections): outputPixelatedVariance = pexConfig.Field( dtype=bool, default=False, doc="Whether to output a pixelated variance map for the generated " - "chi^2 coadd, or to have a flat variance map defined by combining " - "the inverse variance maps of the coadds that were combined." + "chi^2 coadd, or to have a flat variance map defined by combining " + "the inverse variance maps of the coadds that were combined.", ) useUnionForMask = pexConfig.Field( dtype=bool, default=True, doc="Whether to calculate the union of the mask plane in each band, " - "or the intersection of the mask plane in each band." + "or the intersection of the mask plane in each band.", ) @@ -151,6 +151,7 @@ class AssembleChi2CoaddTask(pipeBase.PipelineTask): .. [4] https://project.lsst.org/meetings/law/sites/lsst.org.meetings.law/files/Building%20and%20using%20coadds.pdf # noqa: E501, W505 """ + ConfigClass = AssembleChi2CoaddConfig _DefaultName = "assembleChi2Coadd" @@ -224,10 +225,10 @@ def run(self, inputCoadds: list[afwImage.Exposure]) -> pipeBase.Struct: variance = refExp.variance.Factory(bbox) if self.config.outputPixelatedVariance: # Write the per pixel variance to the output coadd - variance.array[:] = np.sum([1/coadd.variance for coadd in inputCoadds], axis=0) + variance.array[:] = np.sum([1 / coadd.variance for coadd in inputCoadds], axis=0) else: # Use a flat variance in each band - variance.array[:] = np.sum(1/np.array(variance_list)) + variance.array[:] = np.sum(1 / np.array(variance_list)) # Combine the masks planes to calculate the mask plae of the new coadd mask = self.combinedMasks([coadd.mask for coadd in inputCoadds]) # Create the exposure @@ -240,10 +241,7 @@ def run(self, inputCoadds: list[afwImage.Exposure]) -> pipeBase.Struct: class DetectChi2SourcesConnections( pipeBase.PipelineTaskConnections, dimensions=("tract", "patch", "skymap"), - defaultTemplates={ - "inputCoaddName": "deepChi2", - "outputCoaddName": "deepChi2" - } + defaultTemplates={"inputCoaddName": "deepChi2", "outputCoaddName": "deepChi2"}, ): detectionSchema = cT.InitOutput( doc="Schema of the detection catalog", @@ -265,10 +263,7 @@ class DetectChi2SourcesConnections( class DetectChi2SourcesConfig(pipeBase.PipelineTaskConfig, pipelineConnections=DetectChi2SourcesConnections): - detection = pexConfig.ConfigurableField( - target=SourceDetectionTask, - doc="Detect sources in chi2 coadd" - ) + detection = pexConfig.ConfigurableField(target=SourceDetectionTask, doc="Detect sources in chi2 coadd") idGenerator = SkyMapIdGeneratorConfig.make_field() diff --git a/python/lsst/drp/tasks/assemble_coadd.py b/python/lsst/drp/tasks/assemble_coadd.py index 66bf2402..56623668 100644 --- a/python/lsst/drp/tasks/assemble_coadd.py +++ b/python/lsst/drp/tasks/assemble_coadd.py @@ -19,8 +19,13 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -__all__ = ["AssembleCoaddTask", "AssembleCoaddConnections", "AssembleCoaddConfig", - "CompareWarpAssembleCoaddTask", "CompareWarpAssembleCoaddConfig"] +__all__ = [ + "AssembleCoaddTask", + "AssembleCoaddConnections", + "AssembleCoaddConfig", + "CompareWarpAssembleCoaddTask", + "CompareWarpAssembleCoaddConfig", +] import copy import numpy @@ -51,40 +56,45 @@ log = logging.getLogger(__name__) -class AssembleCoaddConnections(pipeBase.PipelineTaskConnections, - dimensions=("tract", "patch", "band", "skymap"), - defaultTemplates={"inputCoaddName": "deep", - "outputCoaddName": "deep", - "warpType": "direct", - "warpTypeSuffix": ""}): - +class AssembleCoaddConnections( + pipeBase.PipelineTaskConnections, + dimensions=("tract", "patch", "band", "skymap"), + defaultTemplates={ + "inputCoaddName": "deep", + "outputCoaddName": "deep", + "warpType": "direct", + "warpTypeSuffix": "", + }, +): inputWarps = pipeBase.connectionTypes.Input( - doc=("Input list of warps to be assemebled i.e. stacked." - "WarpType (e.g. direct, psfMatched) is controlled by the " - "warpType config parameter"), + doc=( + "Input list of warps to be assemebled i.e. stacked." + "WarpType (e.g. direct, psfMatched) is controlled by the " + "warpType config parameter" + ), name="{inputCoaddName}Coadd_{warpType}Warp", storageClass="ExposureF", dimensions=("tract", "patch", "skymap", "visit", "instrument"), deferLoad=True, - multiple=True + multiple=True, ) skyMap = pipeBase.connectionTypes.Input( - doc="Input definition of geometry/bbox and projection/wcs for coadded " - "exposures", + doc="Input definition of geometry/bbox and projection/wcs for coadded " "exposures", name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, storageClass="SkyMap", - dimensions=("skymap", ), + dimensions=("skymap",), ) selectedVisits = pipeBase.connectionTypes.Input( doc="Selected visits to be coadded.", name="{outputCoaddName}Visits", storageClass="StructuredDataDict", - dimensions=("instrument", "tract", "patch", "skymap", "band") + dimensions=("instrument", "tract", "patch", "skymap", "band"), ) brightObjectMask = pipeBase.connectionTypes.PrerequisiteInput( - doc=("Input Bright Object Mask mask produced with external catalogs " - "to be applied to the mask plane BRIGHT_OBJECT." - ), + doc=( + "Input Bright Object Mask mask produced with external catalogs " + "to be applied to the mask plane BRIGHT_OBJECT." + ), name="brightObjectMask", storageClass="ObjectMaskCatalog", dimensions=("tract", "patch", "skymap", "band"), @@ -125,8 +135,9 @@ def __init__(self, *, config=None): self.outputs.remove("inputMap") -class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig, - pipelineConnections=AssembleCoaddConnections): +class AssembleCoaddConfig( + CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig, pipelineConnections=AssembleCoaddConnections +): warpType = pexConfig.Field( doc="Warp name: one of 'direct' or 'psfMatched'", dtype=str, @@ -147,7 +158,7 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig ) doOnlineForMean = pexConfig.Field( dtype=bool, - doc="Perform online coaddition when statistic=\"MEAN\" to save memory?", + doc='Perform online coaddition when statistic="MEAN" to save memory?', default=False, ) doSigmaClip = pexConfig.Field( @@ -158,31 +169,27 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig ) sigmaClip = pexConfig.Field( dtype=float, - doc="Sigma for outlier rejection; ignored if non-clipping statistic " - "selected.", + doc="Sigma for outlier rejection; ignored if non-clipping statistic " "selected.", default=3.0, ) clipIter = pexConfig.Field( dtype=int, - doc="Number of iterations of outlier rejection; ignored if " - "non-clipping statistic selected.", + doc="Number of iterations of outlier rejection; ignored if " "non-clipping statistic selected.", default=2, ) calcErrorFromInputVariance = pexConfig.Field( dtype=bool, doc="Calculate coadd variance from input variance by stacking " - "statistic. Passed to " - "StatisticsControl.setCalcErrorFromInputVariance()", + "statistic. Passed to " + "StatisticsControl.setCalcErrorFromInputVariance()", default=True, ) scaleZeroPoint = pexConfig.ConfigurableField( target=ScaleZeroPointTask, - doc="Task to adjust the photometric zero point of the coadd temp " - "exposures", + doc="Task to adjust the photometric zero point of the coadd temp " "exposures", ) doInterp = pexConfig.Field( - doc="Interpolate over NaN pixels? Also extrapolate, if necessary, but " - "the results are ugly.", + doc="Interpolate over NaN pixels? Also extrapolate, if necessary, but " "the results are ugly.", dtype=bool, default=True, ) @@ -202,17 +209,19 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig ) doUsePsfMatchedPolygons = pexConfig.Field( doc="Use ValidPolygons from shrunk Psf-Matched Calexps? Should be set " - "to True by CompareWarp only.", + "to True by CompareWarp only.", dtype=bool, default=False, ) maskPropagationThresholds = pexConfig.DictField( keytype=str, itemtype=float, - doc=("Threshold (in fractional weight) of rejection at which we " - "propagate a mask plane to the coadd; that is, we set the mask " - "bit on the coadd if the fraction the rejected frames " - "would have contributed exceeds this value."), + doc=( + "Threshold (in fractional weight) of rejection at which we " + "propagate a mask plane to the coadd; that is, we set the mask " + "bit on the coadd if the fraction the rejected frames " + "would have contributed exceeds this value." + ), default={"SAT": 0.1}, ) removeMaskPlanes = pexConfig.ListField( @@ -238,13 +247,15 @@ class AssembleCoaddConfig(CoaddBaseTask.ConfigClass, pipeBase.PipelineTaskConfig dtype=bool, default=False, optional=False, - doc=("Attach a piecewise TransmissionCurve for the coadd? " - "(requires all input Exposures to have TransmissionCurves).") + doc=( + "Attach a piecewise TransmissionCurve for the coadd? " + "(requires all input Exposures to have TransmissionCurves)." + ), ) hasFakes = pexConfig.Field( dtype=bool, default=False, - doc="Should be set to True if fake sources have been inserted into the input data." + doc="Should be set to True if fake sources have been inserted into the input data.", ) doSelectVisits = pexConfig.Field( doc="Coadd only visits selected by a SelectVisitsTask", @@ -271,20 +282,24 @@ def validate(self): # Backwards compatibility. # Configs do not have loggers log.warning("Config doPsfMatch deprecated. Setting warpType='psfMatched'") - self.warpType = 'psfMatched' + self.warpType = "psfMatched" if self.doSigmaClip and self.statistic != "MEANCLIP": log.warning('doSigmaClip deprecated. To replicate behavior, setting statistic to "MEANCLIP"') self.statistic = "MEANCLIP" - if self.doInterp and self.statistic not in ['MEAN', 'MEDIAN', 'MEANCLIP', 'VARIANCE', 'VARIANCECLIP']: - raise ValueError("Must set doInterp=False for statistic=%s, which does not " - "compute and set a non-zero coadd variance estimate." % (self.statistic)) + if self.doInterp and self.statistic not in ["MEAN", "MEDIAN", "MEANCLIP", "VARIANCE", "VARIANCECLIP"]: + raise ValueError( + "Must set doInterp=False for statistic=%s, which does not " + "compute and set a non-zero coadd variance estimate." % (self.statistic) + ) - unstackableStats = ['NOTHING', 'ERROR', 'ORMASK'] + unstackableStats = ["NOTHING", "ERROR", "ORMASK"] if not hasattr(afwMath.Property, self.statistic) or self.statistic in unstackableStats: - stackableStats = [str(k) for k in afwMath.Property.__members__.keys() - if str(k) not in unstackableStats] - raise ValueError("statistic %s is not allowed. Please choose one of %s." - % (self.statistic, stackableStats)) + stackableStats = [ + str(k) for k in afwMath.Property.__members__.keys() if str(k) not in unstackableStats + ] + raise ValueError( + "statistic %s is not allowed. Please choose one of %s." % (self.statistic, stackableStats) + ) class AssembleCoaddTask(CoaddBaseTask, pipeBase.PipelineTask): @@ -332,9 +347,12 @@ def __init__(self, *args, **kwargs): if args: argNames = ["config", "name", "parentTask", "log"] kwargs.update({k: v for k, v in zip(argNames, args)}) - warnings.warn("AssembleCoadd received positional args, and casting them as kwargs: %s. " - "PipelineTask will not take positional args" % argNames, FutureWarning, - stacklevel=2) + warnings.warn( + "AssembleCoadd received positional args, and casting them as kwargs: %s. " + "PipelineTask will not take positional args" % argNames, + FutureWarning, + stacklevel=2, + ) super().__init__(**kwargs) self.makeSubtask("interpImage") @@ -345,8 +363,10 @@ def __init__(self, *args, **kwargs): try: self.brightObjectBitmask = 1 << mask.addMaskPlane(self.config.brightObjectMaskName) except pexExceptions.LsstCppException: - raise RuntimeError("Unable to define mask plane for bright objects; planes used are %s" % - mask.getMaskPlaneDict().keys()) + raise RuntimeError( + "Unable to define mask plane for bright objects; planes used are %s" + % mask.getMaskPlaneDict().keys() + ) del mask if self.config.doInputMap: @@ -363,29 +383,33 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): skyMap = inputData["skyMap"] outputDataId = butlerQC.quantum.dataId - inputData['skyInfo'] = makeSkyInfo(skyMap, - tractId=outputDataId['tract'], - patchId=outputDataId['patch']) + inputData["skyInfo"] = makeSkyInfo( + skyMap, tractId=outputDataId["tract"], patchId=outputDataId["patch"] + ) if self.config.doSelectVisits: - warpRefList = self.filterWarps(inputData['inputWarps'], inputData['selectedVisits']) + warpRefList = self.filterWarps(inputData["inputWarps"], inputData["selectedVisits"]) else: - warpRefList = inputData['inputWarps'] + warpRefList = inputData["inputWarps"] inputs = self.prepareInputs(warpRefList) - self.log.info("Found %d %s", len(inputs.tempExpRefList), - self.getTempExpDatasetName(self.warpType)) + self.log.info("Found %d %s", len(inputs.tempExpRefList), self.getTempExpDatasetName(self.warpType)) if len(inputs.tempExpRefList) == 0: raise pipeBase.NoWorkFound("No coadd temporary exposures found") supplementaryData = self._makeSupplementaryData(butlerQC, inputRefs, outputRefs) - retStruct = self.run(inputData['skyInfo'], inputs.tempExpRefList, inputs.imageScalerList, - inputs.weightList, supplementaryData=supplementaryData) + retStruct = self.run( + inputData["skyInfo"], + inputs.tempExpRefList, + inputs.imageScalerList, + inputs.weightList, + supplementaryData=supplementaryData, + ) - inputData.setdefault('brightObjectMask', None) + inputData.setdefault("brightObjectMask", None) if self.config.doMaskBrightObjects and inputData["brightObjectMask"] is None: log.warning("doMaskBrightObjects is set to True, but brightObjectMask not loaded") - self.processResults(retStruct.coaddExposure, inputData['brightObjectMask'], outputDataId) + self.processResults(retStruct.coaddExposure, inputData["brightObjectMask"], outputDataId) if self.config.doWrite: butlerQC.put(retStruct, outputRefs) @@ -441,7 +465,7 @@ def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs): @deprecated( reason="makeSupplementaryDataGen3 is deprecated in favor of _makeSupplementaryData", version="v25.0", - category=FutureWarning + category=FutureWarning, ) def makeSupplementaryDataGen3(self, butlerQC, inputRefs, outputRefs): return self._makeSupplementaryData(butlerQC, inputRefs, outputRefs) @@ -500,8 +524,9 @@ def prepareInputs(self, refList): except Exception as e: self.log.warning("Scaling failed for %s (skipping it): %s", tempExpRef.dataId, e) continue - statObj = afwMath.makeStatistics(maskedImage.getVariance(), maskedImage.getMask(), - afwMath.MEANCLIP, statsCtrl) + statObj = afwMath.makeStatistics( + maskedImage.getVariance(), maskedImage.getMask(), afwMath.MEANCLIP, statsCtrl + ) meanVar, meanVarErr = statObj.getResult(afwMath.MEANCLIP) weight = 1.0 / float(meanVar) if not numpy.isfinite(weight): @@ -516,8 +541,9 @@ def prepareInputs(self, refList): weightList.append(weight) imageScalerList.append(imageScaler) - return pipeBase.Struct(tempExpRefList=tempExpRefList, weightList=weightList, - imageScalerList=imageScalerList) + return pipeBase.Struct( + tempExpRefList=tempExpRefList, weightList=weightList, imageScalerList=imageScalerList + ) def prepareStats(self, mask=None): """Prepare the statistics for coadding images. @@ -554,8 +580,16 @@ def prepareStats(self, mask=None): return pipeBase.Struct(ctrl=statsCtrl, flags=statsFlags) @timeMethod - def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, - altMaskList=None, mask=None, supplementaryData=None): + def run( + self, + skyInfo, + tempExpRefList, + imageScalerList, + weightList, + altMaskList=None, + mask=None, + supplementaryData=None, + ): """Assemble a coadd from input warps. Assemble the coadd using the provided list of coaddTempExps. Since @@ -618,7 +652,7 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, stats = self.prepareStats(mask=mask) if altMaskList is None: - altMaskList = [None]*len(tempExpRefList) + altMaskList = [None] * len(tempExpRefList) coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs) coaddExposure.setPhotoCalib(self.scaleZeroPoint.getPhotoCalib()) @@ -636,24 +670,38 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, # If inputMap is requested, create the initial version that can be # masked in assembleSubregion. if self.config.doInputMap: - self.inputMapper.build_ccd_input_map(skyInfo.bbox, - skyInfo.wcs, - coaddExposure.getInfo().getCoaddInputs().ccds) + self.inputMapper.build_ccd_input_map( + skyInfo.bbox, skyInfo.wcs, coaddExposure.getInfo().getCoaddInputs().ccds + ) if self.config.doOnlineForMean and self.config.statistic == "MEAN": try: - self.assembleOnlineMeanCoadd(coaddExposure, tempExpRefList, imageScalerList, - weightList, altMaskList, stats.ctrl, - nImage=nImage) + self.assembleOnlineMeanCoadd( + coaddExposure, + tempExpRefList, + imageScalerList, + weightList, + altMaskList, + stats.ctrl, + nImage=nImage, + ) except Exception as e: self.log.exception("Cannot compute online coadd %s", e) raise else: for subBBox in subBBoxIter(skyInfo.bbox, subregionSize): try: - self.assembleSubregion(coaddExposure, subBBox, tempExpRefList, imageScalerList, - weightList, altMaskList, stats.flags, stats.ctrl, - nImage=nImage) + self.assembleSubregion( + coaddExposure, + subBBox, + tempExpRefList, + imageScalerList, + weightList, + altMaskList, + stats.flags, + stats.ctrl, + nImage=nImage, + ) except Exception as e: self.log.exception("Cannot compute coadd %s: %s", subBBox, e) raise @@ -671,9 +719,14 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, # pixels: it identifies pixels that didn't receive any unmasked inputs # (as occurs around the edge of the field). coaddUtils.setCoaddEdgeBits(coaddMaskedImage.getMask(), coaddMaskedImage.getVariance()) - return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage, - warpRefList=tempExpRefList, imageScalerList=imageScalerList, - weightList=weightList, inputMap=inputMap) + return pipeBase.Struct( + coaddExposure=coaddExposure, + nImage=nImage, + warpRefList=tempExpRefList, + imageScalerList=imageScalerList, + weightList=weightList, + inputMap=inputMap, + ) def assembleMetadata(self, coaddExposure, tempExpRefList, weightList): """Set the metadata for the coadd. @@ -702,7 +755,7 @@ def assembleMetadata(self, coaddExposure, tempExpRefList, weightList): # (see #2777). bbox = geom.Box2I(coaddExposure.getBBox().getMin(), geom.Extent2I(1, 1)) - tempExpList = [tempExpRef.get(parameters={'bbox': bbox}) for tempExpRef in tempExpRefList] + tempExpList = [tempExpRef.get(parameters={"bbox": bbox}) for tempExpRef in tempExpRefList] numCcds = sum(len(tempExp.getInfo().getCoaddInputs().ccds) for tempExp in tempExpList) @@ -728,22 +781,35 @@ def assembleMetadata(self, coaddExposure, tempExpRefList, weightList): # Likewise, set the PSF of a PSF-Matched Coadd to the modelPsf # having the maximum width (sufficient because square) modelPsfList = [tempExp.getPsf() for tempExp in tempExpList] - modelPsfWidthList = [modelPsf.computeBBox(modelPsf.getAveragePosition()).getWidth() - for modelPsf in modelPsfList] + modelPsfWidthList = [ + modelPsf.computeBBox(modelPsf.getAveragePosition()).getWidth() for modelPsf in modelPsfList + ] psf = modelPsfList[modelPsfWidthList.index(max(modelPsfWidthList))] else: - psf = measAlg.CoaddPsf(coaddInputs.ccds, coaddExposure.getWcs(), - self.config.coaddPsf.makeControl()) + psf = measAlg.CoaddPsf( + coaddInputs.ccds, coaddExposure.getWcs(), self.config.coaddPsf.makeControl() + ) coaddExposure.setPsf(psf) - apCorrMap = measAlg.makeCoaddApCorrMap(coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT), - coaddExposure.getWcs()) + apCorrMap = measAlg.makeCoaddApCorrMap( + coaddInputs.ccds, coaddExposure.getBBox(afwImage.PARENT), coaddExposure.getWcs() + ) coaddExposure.getInfo().setApCorrMap(apCorrMap) if self.config.doAttachTransmissionCurve: transmissionCurve = measAlg.makeCoaddTransmissionCurve(coaddExposure.getWcs(), coaddInputs.ccds) coaddExposure.getInfo().setTransmissionCurve(transmissionCurve) - def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList, weightList, - altMaskList, statsFlags, statsCtrl, nImage=None): + def assembleSubregion( + self, + coaddExposure, + bbox, + tempExpRefList, + imageScalerList, + weightList, + altMaskList, + statsFlags, + statsCtrl, + nImage=None, + ): """Assemble the coadd for a sub-region. For each coaddTempExp, check for (and swap in) an alternative mask @@ -789,8 +855,7 @@ def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList if nImage is not None: subNImage = afwImage.ImageU(bbox.getWidth(), bbox.getHeight()) for tempExpRef, imageScaler, altMask in zip(tempExpRefList, imageScalerList, altMaskList): - - exposure = tempExpRef.get(parameters={'bbox': bbox}) + exposure = tempExpRef.get(parameters={"bbox": bbox}) maskedImage = exposure.getMaskedImage() mask = maskedImage.getMask() @@ -812,15 +877,21 @@ def assembleSubregion(self, coaddExposure, bbox, tempExpRefList, imageScalerList self.inputMapper.mask_warp_bbox(bbox, visit, mask, statsCtrl.getAndMask()) with self.timer("stack"): - coaddSubregion = afwMath.statisticsStack(maskedImageList, statsFlags, statsCtrl, weightList, - clipped, # also set output to CLIPPED if sigma-clipped - maskMap) + coaddSubregion = afwMath.statisticsStack( + maskedImageList, + statsFlags, + statsCtrl, + weightList, + clipped, # also set output to CLIPPED if sigma-clipped + maskMap, + ) coaddExposure.maskedImage.assign(coaddSubregion, bbox) if nImage is not None: nImage.assign(subNImage, bbox) - def assembleOnlineMeanCoadd(self, coaddExposure, tempExpRefList, imageScalerList, weightList, - altMaskList, statsCtrl, nImage=None): + def assembleOnlineMeanCoadd( + self, coaddExposure, tempExpRefList, imageScalerList, weightList, altMaskList, statsCtrl, nImage=None + ): """Assemble the coadd using the "online" method. This method takes a running sum of images and weights to save memory. @@ -862,13 +933,12 @@ def assembleOnlineMeanCoadd(self, coaddExposure, tempExpRefList, imageScalerList mask_map=maskMap, no_good_pixels_mask=statsCtrl.getNoGoodPixelsMask(), calc_error_from_input_variance=self.config.calcErrorFromInputVariance, - compute_n_image=(nImage is not None) + compute_n_image=(nImage is not None), ) - for tempExpRef, imageScaler, altMask, weight in zip(tempExpRefList, - imageScalerList, - altMaskList, - weightList): + for tempExpRef, imageScaler, altMask, weight in zip( + tempExpRefList, imageScalerList, altMaskList, weightList + ): exposure = tempExpRef.get() maskedImage = exposure.getMaskedImage() mask = maskedImage.getMask() @@ -907,8 +977,9 @@ def removeMaskPlanes(self, maskedImage): try: mask &= ~mask.getPlaneBitMask(maskPlane) except pexExceptions.InvalidParameterError: - self.log.debug("Unable to remove mask plane %s: no mask plane with that name was found.", - maskPlane) + self.log.debug( + "Unable to remove mask plane %s: no mask plane with that name was found.", maskPlane + ) @staticmethod def setRejectedMaskMapping(statsCtrl): @@ -934,9 +1005,11 @@ def setRejectedMaskMapping(statsCtrl): noData = afwImage.Mask.getPlaneBitMask("NO_DATA") clipped = afwImage.Mask.getPlaneBitMask("CLIPPED") toReject = statsCtrl.getAndMask() & (~noData) & (~edge) & (~clipped) - maskMap = [(toReject, afwImage.Mask.getPlaneBitMask("REJECTED")), - (edge, afwImage.Mask.getPlaneBitMask("SENSOR_EDGE")), - (clipped, clipped)] + maskMap = [ + (toReject, afwImage.Mask.getPlaneBitMask("REJECTED")), + (edge, afwImage.Mask.getPlaneBitMask("SENSOR_EDGE")), + (clipped, clipped), + ] return maskMap def applyAltMaskPlanes(self, mask, altMaskSpans): @@ -964,7 +1037,7 @@ def applyAltMaskPlanes(self, mask, altMaskSpans): # self.config.doUsePsfMatchedPolygons should be True only in # CompareWarpAssemble. This mask-clearing step must only occur # *before* applying the new masks below. - for spanSet in altMaskSpans['NO_DATA']: + for spanSet in altMaskSpans["NO_DATA"]: spanSet.clippedTo(mask.getBBox()).clearMask(mask, self.getBadPixelMask()) for plane, spanSetList in altMaskSpans.items(): @@ -987,7 +1060,7 @@ def shrinkValidPolygons(self, coaddInputs): for ccd in coaddInputs.ccds: polyOrig = ccd.getValidPolygon() validPolyBBox = polyOrig.getBBox() if polyOrig else ccd.getBBox() - validPolyBBox.grow(-self.config.matchingKernelSize//2) + validPolyBBox.grow(-self.config.matchingKernelSize // 2) if polyOrig: validPolygon = polyOrig.intersectionSingle(validPolyBBox) else: @@ -1017,18 +1090,20 @@ def setBrightObjectMasks(self, exposure, brightObjectMasks, dataId=None): for rec in brightObjectMasks: center = geom.PointI(wcs.skyToPixel(rec.getCoord())) if rec["type"] == "box": - assert rec["angle"] == 0.0, ("Angle != 0 for mask object %s" % rec["id"]) - width = rec["width"].asArcseconds()/plateScale # convert to pixels - height = rec["height"].asArcseconds()/plateScale # convert to pixels + assert rec["angle"] == 0.0, "Angle != 0 for mask object %s" % rec["id"] + width = rec["width"].asArcseconds() / plateScale # convert to pixels + height = rec["height"].asArcseconds() / plateScale # convert to pixels - halfSize = geom.ExtentI(0.5*width, 0.5*height) + halfSize = geom.ExtentI(0.5 * width, 0.5 * height) bbox = geom.Box2I(center - halfSize, center + halfSize) - bbox = geom.BoxI(geom.PointI(int(center[0] - 0.5*width), int(center[1] - 0.5*height)), - geom.PointI(int(center[0] + 0.5*width), int(center[1] + 0.5*height))) + bbox = geom.BoxI( + geom.PointI(int(center[0] - 0.5 * width), int(center[1] - 0.5 * height)), + geom.PointI(int(center[0] + 0.5 * width), int(center[1] + 0.5 * height)), + ) spans = afwGeom.SpanSet(bbox) elif rec["type"] == "circle": - radius = int(rec["radius"].asArcseconds()/plateScale) # convert to pixels + radius = int(rec["radius"].asArcseconds() / plateScale) # convert to pixels spans = afwGeom.SpanSet.fromShape(radius, offset=center) else: self.log.warning("Unexpected region type %s at %s", rec["type"], center) @@ -1070,12 +1145,11 @@ def filterWarps(self, inputs, goodVisits): Returns ------- - filteredInputs : `list` \ - [`~lsst.pipe.base.connections.DeferredDatasetRef`] + filterInputs : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`] Filtered and sorted list of inputRefs with visitId in goodVisits ordered by goodVisit. """ - inputWarpDict = {inputRef.ref.dataId['visit']: inputRef for inputRef in inputs} + inputWarpDict = {inputRef.ref.dataId["visit"]: inputRef for inputRef in inputs} filteredInputs = [] for visit in goodVisits.keys(): if visit in inputWarpDict: @@ -1112,24 +1186,29 @@ def countMaskFromFootprint(mask, footprint, bitmask, ignoreMask): fp = afwImage.Mask(bbox) subMask = mask.Factory(mask, bbox, afwImage.PARENT) footprint.spans.setMask(fp, bitmask) - return numpy.logical_and((subMask.getArray() & fp.getArray()) > 0, - (subMask.getArray() & ignoreMask) == 0).sum() + return numpy.logical_and( + (subMask.getArray() & fp.getArray()) > 0, (subMask.getArray() & ignoreMask) == 0 + ).sum() class CompareWarpAssembleCoaddConnections(AssembleCoaddConnections): psfMatchedWarps = pipeBase.connectionTypes.Input( - doc=("PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. " - "Only PSF-Matched Warps make sense for image subtraction. " - "Therefore, they must be an additional declared input."), + doc=( + "PSF-Matched Warps are required by CompareWarp regardless of the coadd type requested. " + "Only PSF-Matched Warps make sense for image subtraction. " + "Therefore, they must be an additional declared input." + ), name="{inputCoaddName}Coadd_psfMatchedWarp", storageClass="ExposureF", dimensions=("tract", "patch", "skymap", "visit"), deferLoad=True, - multiple=True + multiple=True, ) templateCoadd = pipeBase.connectionTypes.Output( - doc=("Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, " - "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True"), + doc=( + "Model of the static sky, used to find temporal artifacts. Typically a PSF-Matched, " + "sigma-clipped coadd. Written if and only if assembleStaticSkyModel.doWrite=True" + ), name="{outputCoaddName}CoaddPsfMatched", storageClass="ExposureF", dimensions=("tract", "patch", "skymap", "band"), @@ -1142,66 +1221,67 @@ def __init__(self, *, config=None): config.validate() -class CompareWarpAssembleCoaddConfig(AssembleCoaddConfig, - pipelineConnections=CompareWarpAssembleCoaddConnections): +class CompareWarpAssembleCoaddConfig( + AssembleCoaddConfig, pipelineConnections=CompareWarpAssembleCoaddConnections +): assembleStaticSkyModel = pexConfig.ConfigurableField( target=AssembleCoaddTask, doc="Task to assemble an artifact-free, PSF-matched Coadd to serve as " - "a naive/first-iteration model of the static sky.", + "a naive/first-iteration model of the static sky.", ) detect = pexConfig.ConfigurableField( target=SourceDetectionTask, - doc="Detect outlier sources on difference between each psfMatched warp and static sky model" + doc="Detect outlier sources on difference between each psfMatched warp and static sky model", ) detectTemplate = pexConfig.ConfigurableField( target=SourceDetectionTask, - doc="Detect sources on static sky model. Only used if doPreserveContainedBySource is True" + doc="Detect sources on static sky model. Only used if doPreserveContainedBySource is True", ) maskStreaks = pexConfig.ConfigurableField( target=MaskStreaksTask, doc="Detect streaks on difference between each psfMatched warp and static sky model. Only used if " - "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by" - "streakMaskName" - ) - streakMaskName = pexConfig.Field( - dtype=str, - default="STREAK", - doc="Name of mask bit used for streaks" + "doFilterMorphological is True. Adds a mask plane to an exposure, with the mask plane name set by" + "streakMaskName", ) + streakMaskName = pexConfig.Field(dtype=str, default="STREAK", doc="Name of mask bit used for streaks") maxNumEpochs = pexConfig.Field( doc="Charactistic maximum local number of epochs/visits in which an artifact candidate can appear " - "and still be masked. The effective maxNumEpochs is a broken linear function of local " - "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). " - "For each footprint detected on the image difference between the psfMatched warp and static sky " - "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more " - "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather " - "than transient and not masked.", + "and still be masked. The effective maxNumEpochs is a broken linear function of local " + "number of epochs (N): min(maxFractionEpochsLow*N, maxNumEpochs + maxFractionEpochsHigh*N). " + "For each footprint detected on the image difference between the psfMatched warp and static sky " + "model, if a significant fraction of pixels (defined by spatialThreshold) are residuals in more " + "than the computed effective maxNumEpochs, the artifact candidate is deemed persistant rather " + "than transient and not masked.", dtype=int, - default=2 + default=2, ) maxFractionEpochsLow = pexConfig.RangeField( doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for low N. " - "Effective maxNumEpochs = " - "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)", + "Effective maxNumEpochs = " + "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)", dtype=float, default=0.4, - min=0., max=1., + min=0.0, + max=1.0, ) maxFractionEpochsHigh = pexConfig.RangeField( doc="Fraction of local number of epochs (N) to use as effective maxNumEpochs for high N. " - "Effective maxNumEpochs = " - "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)", + "Effective maxNumEpochs = " + "min(maxFractionEpochsLow * N, maxNumEpochs + maxFractionEpochsHigh * N)", dtype=float, default=0.03, - min=0., max=1., + min=0.0, + max=1.0, ) spatialThreshold = pexConfig.RangeField( doc="Unitless fraction of pixels defining how much of the outlier region has to meet the " - "temporal criteria. If 0, clip all. If 1, clip none.", + "temporal criteria. If 0, clip all. If 1, clip none.", dtype=float, default=0.5, - min=0., max=1., - inclusiveMin=True, inclusiveMax=True + min=0.0, + max=1.0, + inclusiveMin=True, + inclusiveMax=True, ) doScaleWarpVariance = pexConfig.Field( doc="Rescale Warp variance plane using empirical noise?", @@ -1214,55 +1294,55 @@ class CompareWarpAssembleCoaddConfig(AssembleCoaddConfig, ) doPreserveContainedBySource = pexConfig.Field( doc="Rescue artifacts from clipping that completely lie within a footprint detected" - "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.", + "on the PsfMatched Template Coadd. Replicates a behavior of SafeClip.", dtype=bool, default=True, ) doPrefilterArtifacts = pexConfig.Field( doc="Ignore artifact candidates that are mostly covered by the bad pixel mask, " - "because they will be excluded anyway. This prevents them from contributing " - "to the outlier epoch count image and potentially being labeled as persistant." - "'Mostly' is defined by the config 'prefilterArtifactsRatio'.", + "because they will be excluded anyway. This prevents them from contributing " + "to the outlier epoch count image and potentially being labeled as persistant." + "'Mostly' is defined by the config 'prefilterArtifactsRatio'.", dtype=bool, - default=True + default=True, ) prefilterArtifactsMaskPlanes = pexConfig.ListField( doc="Prefilter artifact candidates that are mostly covered by these bad mask planes.", dtype=str, - default=('NO_DATA', 'BAD', 'SAT', 'SUSPECT'), + default=("NO_DATA", "BAD", "SAT", "SUSPECT"), ) prefilterArtifactsRatio = pexConfig.Field( doc="Prefilter artifact candidates with less than this fraction overlapping good pixels", dtype=float, - default=0.05 + default=0.05, ) doFilterMorphological = pexConfig.Field( doc="Filter artifact candidates based on morphological criteria, i.g. those that appear to " - "be streaks.", + "be streaks.", dtype=bool, - default=False + default=False, ) growStreakFp = pexConfig.Field( - doc="Grow streak footprints by this number multiplied by the PSF width", - dtype=float, - default=5 + doc="Grow streak footprints by this number multiplied by the PSF width", dtype=float, default=5 ) def setDefaults(self): AssembleCoaddConfig.setDefaults(self) - self.statistic = 'MEAN' + self.statistic = "MEAN" self.doUsePsfMatchedPolygons = True # Real EDGE removed by psfMatched NO_DATA border half the width of the # matching kernel. CompareWarp applies psfMatched EDGE pixels to # directWarps before assembling. if "EDGE" in self.badMaskPlanes: - self.badMaskPlanes.remove('EDGE') - self.removeMaskPlanes.append('EDGE') - self.assembleStaticSkyModel.badMaskPlanes = ["NO_DATA", ] - self.assembleStaticSkyModel.warpType = 'psfMatched' - self.assembleStaticSkyModel.connections.warpType = 'psfMatched' - self.assembleStaticSkyModel.statistic = 'MEANCLIP' + self.badMaskPlanes.remove("EDGE") + self.removeMaskPlanes.append("EDGE") + self.assembleStaticSkyModel.badMaskPlanes = [ + "NO_DATA", + ] + self.assembleStaticSkyModel.warpType = "psfMatched" + self.assembleStaticSkyModel.connections.warpType = "psfMatched" + self.assembleStaticSkyModel.statistic = "MEANCLIP" self.assembleStaticSkyModel.sigmaClip = 2.5 self.assembleStaticSkyModel.clipIter = 3 self.assembleStaticSkyModel.calcErrorFromInputVariance = False @@ -1286,14 +1366,18 @@ def setDefaults(self): def validate(self): super().validate() if self.assembleStaticSkyModel.doNImage: - raise ValueError("No dataset type exists for a PSF-Matched Template N Image." - "Please set assembleStaticSkyModel.doNImage=False") + raise ValueError( + "No dataset type exists for a PSF-Matched Template N Image." + "Please set assembleStaticSkyModel.doNImage=False" + ) if self.assembleStaticSkyModel.doWrite and (self.warpType == self.assembleStaticSkyModel.warpType): - raise ValueError("warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for " - "the same dataset name. Please set assembleStaticSkyModel.doWrite to False " - "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be " - "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType)) + raise ValueError( + "warpType (%s) == assembleStaticSkyModel.warpType (%s) and will compete for " + "the same dataset name. Please set assembleStaticSkyModel.doWrite to False " + "or warpType to 'direct'. assembleStaticSkyModel.warpType should ways be " + "'PsfMatched'" % (self.warpType, self.assembleStaticSkyModel.warpType) + ) class CompareWarpAssembleCoaddTask(AssembleCoaddTask): @@ -1393,22 +1477,25 @@ def _makeSupplementaryData(self, butlerQC, inputRefs, outputRefs): del staticSkyModelOutputRefs.templateCoadd # A PSF-Matched nImage does not exist as a dataset type - if 'nImage' in staticSkyModelOutputRefs.keys(): + if "nImage" in staticSkyModelOutputRefs.keys(): del staticSkyModelOutputRefs.nImage - templateCoadd = self.assembleStaticSkyModel.runQuantum(butlerQC, staticSkyModelInputRefs, - staticSkyModelOutputRefs) + templateCoadd = self.assembleStaticSkyModel.runQuantum( + butlerQC, staticSkyModelInputRefs, staticSkyModelOutputRefs + ) if templateCoadd is None: raise RuntimeError(self._noTemplateMessage(self.assembleStaticSkyModel.warpType)) - return pipeBase.Struct(templateCoadd=templateCoadd.coaddExposure, - nImage=templateCoadd.nImage, - warpRefList=templateCoadd.warpRefList, - imageScalerList=templateCoadd.imageScalerList, - weightList=templateCoadd.weightList) + return pipeBase.Struct( + templateCoadd=templateCoadd.coaddExposure, + nImage=templateCoadd.nImage, + warpRefList=templateCoadd.warpRefList, + imageScalerList=templateCoadd.imageScalerList, + weightList=templateCoadd.weightList, + ) def _noTemplateMessage(self, warpType): - warpName = (warpType[0].upper() + warpType[1:]) + warpName = warpType[0].upper() + warpType[1:] message = """No %(warpName)s warps were found to build the template coadd which is required to run CompareWarpAssembleCoaddTask. To continue assembling this type of coadd, first either rerun makeCoaddTempExp with config.make%(warpName)s=True or @@ -1419,13 +1506,14 @@ def _noTemplateMessage(self, warpType): from lsst.pipe.tasks.assembleCoadd import SafeClipAssembleCoaddTask config.assemble.retarget(SafeClipAssembleCoaddTask) - """ % {"warpName": warpName} + """ % { + "warpName": warpName + } return message @utils.inheritDoc(AssembleCoaddTask) @timeMethod - def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, - supplementaryData): + def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, supplementaryData): """Notes ----- Assemble the coadd. @@ -1443,23 +1531,26 @@ def run(self, skyInfo, tempExpRefList, imageScalerList, weightList, if dataIds != psfMatchedDataIds: self.log.info("Reordering and or/padding PSF-matched visit input list") - supplementaryData.warpRefList = reorderAndPadList(supplementaryData.warpRefList, - psfMatchedDataIds, dataIds) - supplementaryData.imageScalerList = reorderAndPadList(supplementaryData.imageScalerList, - psfMatchedDataIds, dataIds) + supplementaryData.warpRefList = reorderAndPadList( + supplementaryData.warpRefList, psfMatchedDataIds, dataIds + ) + supplementaryData.imageScalerList = reorderAndPadList( + supplementaryData.imageScalerList, psfMatchedDataIds, dataIds + ) # Use PSF-Matched Warps (and corresponding scalers) and coadd to find # artifacts. - spanSetMaskList = self.findArtifacts(supplementaryData.templateCoadd, - supplementaryData.warpRefList, - supplementaryData.imageScalerList) + spanSetMaskList = self.findArtifacts( + supplementaryData.templateCoadd, supplementaryData.warpRefList, supplementaryData.imageScalerList + ) badMaskPlanes = self.config.badMaskPlanes[:] badMaskPlanes.append("CLIPPED") badPixelMask = afwImage.Mask.getPlaneBitMask(badMaskPlanes) - result = AssembleCoaddTask.run(self, skyInfo, tempExpRefList, imageScalerList, weightList, - spanSetMaskList, mask=badPixelMask) + result = AssembleCoaddTask.run( + self, skyInfo, tempExpRefList, imageScalerList, weightList, spanSetMaskList, mask=badPixelMask + ) # Propagate PSF-matched EDGE pixels to coadd SENSOR_EDGE and # INEXACT_PSF. Psf-Matching moves the real edge inwards. @@ -1482,7 +1573,7 @@ def applyAltEdgeMask(self, mask, altMaskList): maskValue = mask.getPlaneBitMask(["SENSOR_EDGE", "INEXACT_PSF"]) for visitMask in altMaskList: if "EDGE" in visitMask: - for spanSet in visitMask['EDGE']: + for spanSet in visitMask["EDGE"]: spanSet.clippedTo(mask.getBBox()).setMask(mask, maskValue) def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList): @@ -1534,8 +1625,9 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList): if warpDiffExp is not None: # This nImage only approximates the final nImage because it # uses the PSF-matched mask. - nImage.array += (numpy.isfinite(warpDiffExp.image.array) - * ((warpDiffExp.mask.array & badPixelMask) == 0)).astype(numpy.uint16) + nImage.array += ( + numpy.isfinite(warpDiffExp.image.array) * ((warpDiffExp.mask.array & badPixelMask) == 0) + ).astype(numpy.uint16) fpSet = self.detect.detectFootprints(warpDiffExp, doSmooth=False, clearMask=True) fpSet.positive.merge(fpSet.negative) footprints = fpSet.positive @@ -1558,8 +1650,9 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList): maskName = self.config.streakMaskName _ = self.maskStreaks.run(warpDiffExp) streakMask = warpDiffExp.mask - spanSetStreak = afwGeom.SpanSet.fromMask(streakMask, - streakMask.getPlaneBitMask(maskName)).split() + spanSetStreak = afwGeom.SpanSet.fromMask( + streakMask, streakMask.getPlaneBitMask(maskName) + ).split() # Pad the streaks to account for low-surface brightness # wings. psf = warpDiffExp.getPsf() @@ -1579,8 +1672,7 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList): nansMask = afwImage.makeMaskFromArray(nans.astype(afwImage.MaskPixel)) nansMask.setXY0(warpDiffExp.getXY0()) edgeMask = warpDiffExp.mask - spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask, - edgeMask.getPlaneBitMask("EDGE")).split() + spanSetEdgeMask = afwGeom.SpanSet.fromMask(edgeMask, edgeMask.getPlaneBitMask("EDGE")).split() else: # If the directWarp has <1% coverage, the psfMatchedWarp can # have 0% and not exist. In this case, mask the whole epoch. @@ -1603,17 +1695,16 @@ def findArtifacts(self, templateCoadd, tempExpRefList, imageScalerList): for i, spanSetList in enumerate(spanSetArtifactList): if spanSetList: - filteredSpanSetList = self.filterArtifacts(spanSetList, epochCountImage, nImage, - templateFootprints) + filteredSpanSetList = self.filterArtifacts( + spanSetList, epochCountImage, nImage, templateFootprints + ) spanSetArtifactList[i] = filteredSpanSetList if self.config.doFilterMorphological: spanSetArtifactList[i] += spanSetBadMorphoList[i] altMasks = [] for artifacts, noData, edge in zip(spanSetArtifactList, spanSetNoDataMaskList, spanSetEdgeList): - altMasks.append({'CLIPPED': artifacts, - 'NO_DATA': noData, - 'EDGE': edge}) + altMasks.append({"CLIPPED": artifacts, "NO_DATA": noData, "EDGE": edge}) return altMasks def prefilterArtifacts(self, spanSetList, exp): @@ -1643,7 +1734,7 @@ def prefilterArtifacts(self, spanSetList, exp): y, x = span.clippedTo(bbox).indices() yIndexLocal = numpy.array(y) - y0 xIndexLocal = numpy.array(x) - x0 - goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal])/span.getArea() + goodRatio = numpy.count_nonzero(goodArr[yIndexLocal, xIndexLocal]) / span.getArea() if goodRatio > self.config.prefilterArtifactsRatio: returnSpanSetList.append(span) return returnSpanSetList @@ -1676,12 +1767,12 @@ def filterArtifacts(self, spanSetList, epochCountImage, nImage, footprintsToExcl # effectiveMaxNumEpochs is broken line (fraction of N) with # characteristic config.maxNumEpochs. - effMaxNumEpochsHighN = (self.config.maxNumEpochs - + self.config.maxFractionEpochsHigh*numpy.mean(totalN)) + effMaxNumEpochsHighN = self.config.maxNumEpochs + self.config.maxFractionEpochsHigh * numpy.mean( + totalN + ) effMaxNumEpochsLowN = self.config.maxFractionEpochsLow * numpy.mean(totalN) effectiveMaxNumEpochs = int(min(effMaxNumEpochsLowN, effMaxNumEpochsHighN)) - nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0) - & (outlierN <= effectiveMaxNumEpochs)) + nPixelsBelowThreshold = numpy.count_nonzero((outlierN > 0) & (outlierN <= effectiveMaxNumEpochs)) percentBelowThreshold = nPixelsBelowThreshold / len(outlierN) if percentBelowThreshold > self.config.spatialThreshold: maskSpanSetList.append(span) diff --git a/python/lsst/drp/tasks/dcr_assemble_coadd.py b/python/lsst/drp/tasks/dcr_assemble_coadd.py index 17ae0c9e..63e4803f 100644 --- a/python/lsst/drp/tasks/dcr_assemble_coadd.py +++ b/python/lsst/drp/tasks/dcr_assemble_coadd.py @@ -44,23 +44,29 @@ from lsst.pipe.tasks.measurePsf import MeasurePsfTask -class DcrAssembleCoaddConnections(AssembleCoaddConnections, - dimensions=("tract", "patch", "band", "skymap"), - defaultTemplates={"inputWarpName": "deep", - "inputCoaddName": "deep", - "outputCoaddName": "dcr", - "warpType": "direct", - "warpTypeSuffix": "", - "fakesType": ""}): +class DcrAssembleCoaddConnections( + AssembleCoaddConnections, + dimensions=("tract", "patch", "band", "skymap"), + defaultTemplates={ + "inputWarpName": "deep", + "inputCoaddName": "deep", + "outputCoaddName": "dcr", + "warpType": "direct", + "warpTypeSuffix": "", + "fakesType": "", + }, +): inputWarps = pipeBase.connectionTypes.Input( - doc=("Input list of warps to be assembled i.e. stacked." - "Note that this will often be different than the inputCoaddName." - "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter"), + doc=( + "Input list of warps to be assembled i.e. stacked." + "Note that this will often be different than the inputCoaddName." + "WarpType (e.g. direct, psfMatched) is controlled by the warpType config parameter" + ), name="{inputWarpName}Coadd_{warpType}Warp", storageClass="ExposureF", dimensions=("tract", "patch", "skymap", "visit", "instrument"), deferLoad=True, - multiple=True + multiple=True, ) templateExposure = pipeBase.connectionTypes.Input( doc="Input coadded exposure, produced by previous call to AssembleCoadd", @@ -95,8 +101,7 @@ def __init__(self, *, config=None): self.outputs.remove("nImage") -class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, - pipelineConnections=DcrAssembleCoaddConnections): +class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, pipelineConnections=DcrAssembleCoaddConnections): dcrNumSubfilters = pexConfig.Field( dtype=int, doc="Number of sub-filters to forward model chromatic effects to fit the supplied exposures.", @@ -122,18 +127,18 @@ class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, useConvergence = pexConfig.Field( dtype=bool, doc="Use convergence test as a forward modeling end condition?" - "If not set, skips calculating convergence and runs for ``maxNumIter`` iterations", + "If not set, skips calculating convergence and runs for ``maxNumIter`` iterations", default=True, ) baseGain = pexConfig.Field( dtype=float, optional=True, doc="Relative weight to give the new solution vs. the last solution when updating the model." - "A value of 1.0 gives equal weight to both solutions." - "Small values imply slower convergence of the solution, but can " - "help prevent overshooting and failures in the fit." - "If ``baseGain`` is None, a conservative gain " - "will be calculated from the number of subfilters. ", + "A value of 1.0 gives equal weight to both solutions." + "Small values imply slower convergence of the solution, but can " + "help prevent overshooting and failures in the fit." + "If ``baseGain`` is None, a conservative gain " + "will be calculated from the number of subfilters. ", default=None, ) useProgressiveGain = pexConfig.Field( @@ -160,37 +165,30 @@ class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, ) splitSubfilters = pexConfig.Field( dtype=bool, - doc="Calculate DCR for two evenly-spaced wavelengths in each subfilter." - "Instead of at the midpoint", + doc="Calculate DCR for two evenly-spaced wavelengths in each subfilter. Instead of at the midpoint", default=True, ) splitThreshold = pexConfig.Field( dtype=float, doc="Minimum DCR difference within a subfilter to use ``splitSubfilters``, in pixels." - "Set to 0 to always split the subfilters.", + "Set to 0 to always split the subfilters.", default=0.1, ) regularizeModelIterations = pexConfig.Field( dtype=float, - doc="Maximum relative change of the model allowed between iterations." - "Set to zero to disable.", - default=2., + doc="Maximum relative change of the model allowed between iterations. Set to zero to disable.", + default=2.0, ) regularizeModelFrequency = pexConfig.Field( dtype=float, - doc="Maximum relative change of the model allowed between subfilters." - "Set to zero to disable.", - default=4., + doc="Maximum relative change of the model allowed between subfilters. Set to zero to disable.", + default=4.0, ) convergenceMaskPlanes = pexConfig.ListField( - dtype=str, - default=["DETECTED"], - doc="Mask planes to use to calculate convergence." + dtype=str, default=["DETECTED"], doc="Mask planes to use to calculate convergence." ) regularizationWidth = pexConfig.Field( - dtype=int, - default=2, - doc="Minimum radius of a region to include in regularization, in pixels." + dtype=int, default=2, doc="Minimum radius of a region to include in regularization, in pixels." ) imageInterpOrder = pexConfig.Field( dtype=int, @@ -214,7 +212,7 @@ class DcrAssembleCoaddConfig(CompareWarpAssembleCoaddConfig, ) measurePsfSources = pexConfig.ConfigurableField( target=SingleFrameMeasurementTask, - doc="Task to measure sources for PSF measurement, if ``doCalculatePsf`` is set." + doc="Task to measure sources for PSF measurement, if ``doCalculatePsf`` is set.", ) measurePsf = pexConfig.ConfigurableField( target=MeasurePsfTask, @@ -252,8 +250,11 @@ def setDefaults(self): # Use the variance plane to calculate signal to noise self.detectPsfSources.thresholdType = "pixel_stdev" # Ensure psf candidate size is as large as piff psf size. - if (self.doCalculatePsf and self.measurePsf.psfDeterminer.name == "piff" - and self.psfDeterminer["piff"].kernelSize > self.makePsfCandidates.kernelSize): + if ( + self.doCalculatePsf + and self.measurePsf.psfDeterminer.name == "piff" + and self.psfDeterminer["piff"].kernelSize > self.makePsfCandidates.kernelSize + ): self.makePsfCandidates.kernelSize = self.psfDeterminer["piff"].kernelSize @@ -340,30 +341,34 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): skyMap = inputData["skyMap"] outputDataId = butlerQC.quantum.dataId - inputData['skyInfo'] = makeSkyInfo(skyMap, - tractId=outputDataId['tract'], - patchId=outputDataId['patch']) + inputData["skyInfo"] = makeSkyInfo( + skyMap, tractId=outputDataId["tract"], patchId=outputDataId["patch"] + ) # Construct list of input Deferred Datasets - warpRefList = inputData['inputWarps'] + warpRefList = inputData["inputWarps"] inputs = self.prepareInputs(warpRefList) - self.log.info("Found %d %s", len(inputs.tempExpRefList), - self.getTempExpDatasetName(self.warpType)) + self.log.info("Found %d %s", len(inputs.tempExpRefList), self.getTempExpDatasetName(self.warpType)) if len(inputs.tempExpRefList) == 0: self.log.warning("No coadd temporary exposures found") return supplementaryData = self._makeSupplementaryData(butlerQC, inputRefs, outputRefs) - retStruct = self.run(inputData['skyInfo'], inputs.tempExpRefList, inputs.imageScalerList, - inputs.weightList, supplementaryData=supplementaryData) + retStruct = self.run( + inputData["skyInfo"], + inputs.tempExpRefList, + inputs.imageScalerList, + inputs.weightList, + supplementaryData=supplementaryData, + ) - inputData.setdefault('brightObjectMask', None) + inputData.setdefault("brightObjectMask", None) for subfilter in range(self.config.dcrNumSubfilters): # Use the PSF of the stacked dcrModel, and do not recalculate the # PSF for each subfilter retStruct.dcrCoadds[subfilter].setPsf(retStruct.coaddExposure.getPsf()) - self.processResults(retStruct.dcrCoadds[subfilter], inputData['brightObjectMask'], outputDataId) + self.processResults(retStruct.dcrCoadds[subfilter], inputData["brightObjectMask"], outputDataId) if self.config.doWrite: butlerQC.put(retStruct, outputRefs) @@ -396,10 +401,7 @@ def measureCoaddPsf(self, coaddExposure): table = afwTable.SourceTable.make(self.schema) detResults = self.detectPsfSources.run(table, coaddExposure, clearMask=False) coaddSources = detResults.sources - self.measurePsfSources.run( - measCat=coaddSources, - exposure=coaddExposure - ) + self.measurePsfSources.run(measCat=coaddSources, exposure=coaddExposure) # Measure the PSF on the stacked subfilter coadds if possible. # We should already have a decent estimate of the coadd PSF, however, # so in case of any errors simply log them as a warning and use the @@ -438,7 +440,7 @@ def prepareDcrInputs(self, templateCoadd, warpRefList, weightList): If ``lambdaMin`` is missing from the Mapper class of the obs package being used. """ - sigma2fwhm = 2.*np.sqrt(2.*np.log(2.)) + sigma2fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0)) filterLabel = templateCoadd.getFilter() dcrShifts = [] airmassDict = {} @@ -450,7 +452,7 @@ def prepareDcrInputs(self, templateCoadd, warpRefList, weightList): visit = warpExpRef.dataId["visit"] # Just need a rough estimate; average positions are fine psfAvgPos = psf.getAveragePosition() - psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius()*sigma2fwhm + psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius() * sigma2fwhm airmass = visitInfo.getBoresightAirmass() parallacticAngle = visitInfo.getBoresightParAngle().asDegrees() airmassDict[visit] = airmass @@ -458,10 +460,19 @@ def prepareDcrInputs(self, templateCoadd, warpRefList, weightList): psfSizeDict[visit] = psfSize if self.config.doAirmassWeight: weightList[visitNum] *= airmass - dcrShifts.append(np.max(np.abs(calculateDcr(visitInfo, templateCoadd.getWcs(), - self.config.effectiveWavelength, - self.config.bandwidth, - self.config.dcrNumSubfilters)))) + dcrShifts.append( + np.max( + np.abs( + calculateDcr( + visitInfo, + templateCoadd.getWcs(), + self.config.effectiveWavelength, + self.config.bandwidth, + self.config.dcrNumSubfilters, + ) + ) + ) + ) self.log.info("Selected airmasses:\n%s", airmassDict) self.log.info("Selected parallactic angles:\n%s", angleDict) self.log.info("Selected PSF sizes:\n%s", psfSizeDict) @@ -472,18 +483,19 @@ def prepareDcrInputs(self, templateCoadd, warpRefList, weightList): self.log.warning("Unable to calculate restricted PSF, using default coadd PSF: %s", e) else: psf = templateCoadd.getPsf() - dcrModels = DcrModel.fromImage(templateCoadd.maskedImage, - self.config.dcrNumSubfilters, - effectiveWavelength=self.config.effectiveWavelength, - bandwidth=self.config.bandwidth, - wcs=templateCoadd.getWcs(), - filterLabel=filterLabel, - psf=psf) + dcrModels = DcrModel.fromImage( + templateCoadd.maskedImage, + self.config.dcrNumSubfilters, + effectiveWavelength=self.config.effectiveWavelength, + bandwidth=self.config.bandwidth, + wcs=templateCoadd.getWcs(), + filterLabel=filterLabel, + psf=psf, + ) return dcrModels @timeMethod - def run(self, skyInfo, warpRefList, imageScalerList, weightList, - supplementaryData=None): + def run(self, skyInfo, warpRefList, imageScalerList, weightList, supplementaryData=None): r"""Assemble the coadd. Requires additional inputs Struct ``supplementaryData`` to contain a @@ -540,7 +552,7 @@ def run(self, skyInfo, warpRefList, imageScalerList, weightList, `list` of exposure count images for each subfilter. """ minNumIter = self.config.minNumIter or self.config.dcrNumSubfilters - maxNumIter = self.config.maxNumIter or self.config.dcrNumSubfilters*3 + maxNumIter = self.config.maxNumIter or self.config.dcrNumSubfilters * 3 templateCoadd = supplementaryData.templateCoadd baseMask = templateCoadd.mask.clone() # The variance plane is for each subfilter @@ -562,8 +574,9 @@ def run(self, skyInfo, warpRefList, imageScalerList, weightList, stats = self.prepareStats(mask=badPixelMask) dcrModels = self.prepareDcrInputs(templateCoadd, warpRefList, weightList) if self.config.doNImage: - dcrNImages, dcrWeights = self.calculateNImage(dcrModels, skyInfo.bbox, warpRefList, - spanSetMaskList, stats.ctrl) + dcrNImages, dcrWeights = self.calculateNImage( + dcrModels, skyInfo.bbox, warpRefList, spanSetMaskList, stats.ctrl + ) nImage = afwImage.ImageU(skyInfo.bbox) # Note that this nImage will be a factor of dcrNumSubfilters higher # than the nImage returned by assembleCoadd for most pixels. This @@ -575,82 +588,134 @@ def run(self, skyInfo, warpRefList, imageScalerList, weightList, dcrNImages = None subregionSize = geom.Extent2I(*self.config.subregionSize) - nSubregions = (ceil(skyInfo.bbox.getHeight()/subregionSize[1]) - * ceil(skyInfo.bbox.getWidth()/subregionSize[0])) + nSubregions = ceil(skyInfo.bbox.getHeight() / subregionSize[1]) * ceil( + skyInfo.bbox.getWidth() / subregionSize[0] + ) subIter = 0 for subBBox in subBBoxIter(skyInfo.bbox, subregionSize): modelIter = 0 subIter += 1 - self.log.info("Computing coadd over patch %s subregion %s of %s: %s", - skyInfo.patchInfo.getIndex(), subIter, nSubregions, subBBox) + self.log.info( + "Computing coadd over patch %s subregion %s of %s: %s", + skyInfo.patchInfo.getIndex(), + subIter, + nSubregions, + subBBox, + ) dcrBBox = geom.Box2I(subBBox) dcrBBox.grow(self.bufferSize) dcrBBox.clip(dcrModels.bbox) modelWeights = self.calculateModelWeights(dcrModels, dcrBBox) - subExposures = self.loadSubExposures(dcrBBox, stats.ctrl, warpRefList, - imageScalerList, spanSetMaskList) - convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox, - warpRefList, weightList, stats.ctrl) + subExposures = self.loadSubExposures( + dcrBBox, stats.ctrl, warpRefList, imageScalerList, spanSetMaskList + ) + convergenceMetric = self.calculateConvergence( + dcrModels, subExposures, subBBox, warpRefList, weightList, stats.ctrl + ) self.log.info("Initial convergence : %s", convergenceMetric) convergenceList = [convergenceMetric] gainList = [] - convergenceCheck = 1. + convergenceCheck = 1.0 refImage = templateCoadd.image - while (convergenceCheck > self.config.convergenceThreshold or modelIter <= minNumIter): + while convergenceCheck > self.config.convergenceThreshold or modelIter <= minNumIter: gain = self.calculateGain(convergenceList, gainList) - self.dcrAssembleSubregion(dcrModels, subExposures, subBBox, dcrBBox, warpRefList, - stats.ctrl, convergenceMetric, gain, - modelWeights, refImage, dcrWeights) + self.dcrAssembleSubregion( + dcrModels, + subExposures, + subBBox, + dcrBBox, + warpRefList, + stats.ctrl, + convergenceMetric, + gain, + modelWeights, + refImage, + dcrWeights, + ) if self.config.useConvergence: - convergenceMetric = self.calculateConvergence(dcrModels, subExposures, subBBox, - warpRefList, weightList, stats.ctrl) + convergenceMetric = self.calculateConvergence( + dcrModels, subExposures, subBBox, warpRefList, weightList, stats.ctrl + ) if convergenceMetric == 0: - self.log.warning("Coadd patch %s subregion %s had convergence metric of 0.0 which is " - "most likely due to there being no valid data in the region.", - skyInfo.patchInfo.getIndex(), subIter) + self.log.warning( + "Coadd patch %s subregion %s had convergence metric of 0.0 which is " + "most likely due to there being no valid data in the region.", + skyInfo.patchInfo.getIndex(), + subIter, + ) break - convergenceCheck = (convergenceList[-1] - convergenceMetric)/convergenceMetric + convergenceCheck = (convergenceList[-1] - convergenceMetric) / convergenceMetric if (convergenceCheck < 0) & (modelIter > minNumIter): - self.log.warning("Coadd patch %s subregion %s diverged before reaching maximum " - "iterations or desired convergence improvement of %s." - " Divergence: %s", - skyInfo.patchInfo.getIndex(), subIter, - self.config.convergenceThreshold, convergenceCheck) + self.log.warning( + "Coadd patch %s subregion %s diverged before reaching maximum " + "iterations or desired convergence improvement of %s." + " Divergence: %s", + skyInfo.patchInfo.getIndex(), + subIter, + self.config.convergenceThreshold, + convergenceCheck, + ) break convergenceList.append(convergenceMetric) if modelIter > maxNumIter: if self.config.useConvergence: - self.log.warning("Coadd patch %s subregion %s reached maximum iterations " - "before reaching desired convergence improvement of %s." - " Final convergence improvement: %s", - skyInfo.patchInfo.getIndex(), subIter, - self.config.convergenceThreshold, convergenceCheck) + self.log.warning( + "Coadd patch %s subregion %s reached maximum iterations " + "before reaching desired convergence improvement of %s." + " Final convergence improvement: %s", + skyInfo.patchInfo.getIndex(), + subIter, + self.config.convergenceThreshold, + convergenceCheck, + ) break if self.config.useConvergence: - self.log.info("Iteration %s with convergence metric %s, %.4f%% improvement (gain: %.2f)", - modelIter, convergenceMetric, 100.*convergenceCheck, gain) + self.log.info( + "Iteration %s with convergence metric %s, %.4f%% improvement (gain: %.2f)", + modelIter, + convergenceMetric, + 100.0 * convergenceCheck, + gain, + ) modelIter += 1 else: if self.config.useConvergence: - self.log.info("Coadd patch %s subregion %s finished with " - "convergence metric %s after %s iterations", - skyInfo.patchInfo.getIndex(), subIter, convergenceMetric, modelIter) + self.log.info( + "Coadd patch %s subregion %s finished with " + "convergence metric %s after %s iterations", + skyInfo.patchInfo.getIndex(), + subIter, + convergenceMetric, + modelIter, + ) else: - self.log.info("Coadd patch %s subregion %s finished after %s iterations", - skyInfo.patchInfo.getIndex(), subIter, modelIter) + self.log.info( + "Coadd patch %s subregion %s finished after %s iterations", + skyInfo.patchInfo.getIndex(), + subIter, + modelIter, + ) if self.config.useConvergence and convergenceMetric > 0: - self.log.info("Final convergence improvement was %.4f%% overall", - 100*(convergenceList[0] - convergenceMetric)/convergenceMetric) - - dcrCoadds = self.fillCoadd(dcrModels, skyInfo, warpRefList, weightList, - calibration=self.scaleZeroPoint.getPhotoCalib(), - coaddInputs=templateCoadd.getInfo().getCoaddInputs(), - mask=baseMask, - variance=baseVariance) + self.log.info( + "Final convergence improvement was %.4f%% overall", + 100 * (convergenceList[0] - convergenceMetric) / convergenceMetric, + ) + + dcrCoadds = self.fillCoadd( + dcrModels, + skyInfo, + warpRefList, + weightList, + calibration=self.scaleZeroPoint.getPhotoCalib(), + coaddInputs=templateCoadd.getInfo().getCoaddInputs(), + mask=baseMask, + variance=baseVariance, + ) coaddExposure = self.stackCoadd(dcrCoadds) - return pipeBase.Struct(coaddExposure=coaddExposure, nImage=nImage, - dcrCoadds=dcrCoadds, dcrNImages=dcrNImages) + return pipeBase.Struct( + coaddExposure=coaddExposure, nImage=nImage, dcrCoadds=dcrCoadds, dcrNImages=dcrNImages + ) def calculateNImage(self, dcrModels, bbox, warpRefList, spanSetMaskList, statsCtrl): """Calculate the number of exposures contributing to each subfilter. @@ -681,37 +746,49 @@ def calculateNImage(self, dcrModels, bbox, warpRefList, spanSetMaskList, statsCt dcrNImages = [afwImage.ImageU(bbox) for subfilter in range(self.config.dcrNumSubfilters)] dcrWeights = [afwImage.ImageF(bbox) for subfilter in range(self.config.dcrNumSubfilters)] for warpExpRef, altMaskSpans in zip(warpRefList, spanSetMaskList): - exposure = warpExpRef.get(parameters={'bbox': bbox}) + exposure = warpExpRef.get(parameters={"bbox": bbox}) visitInfo = exposure.getInfo().getVisitInfo() wcs = exposure.getInfo().getWcs() mask = exposure.mask if altMaskSpans is not None: self.applyAltMaskPlanes(mask, altMaskSpans) weightImage = np.zeros_like(exposure.image.array) - weightImage[(mask.array & statsCtrl.getAndMask()) == 0] = 1. + weightImage[(mask.array & statsCtrl.getAndMask()) == 0] = 1.0 # The weights must be shifted in exactly the same way as the # residuals, because they will be used as the denominator in the # weighted average of residuals. - weightsGenerator = self.dcrResiduals(weightImage, visitInfo, wcs, - dcrModels.effectiveWavelength, dcrModels.bandwidth) + weightsGenerator = self.dcrResiduals( + weightImage, visitInfo, wcs, dcrModels.effectiveWavelength, dcrModels.bandwidth + ) for shiftedWeights, dcrNImage, dcrWeight in zip(weightsGenerator, dcrNImages, dcrWeights): dcrNImage.array += np.rint(shiftedWeights).astype(dcrNImage.array.dtype) dcrWeight.array += shiftedWeights # Exclude any pixels that don't have at least one exposure contributing # in all subfilters - weightsThreshold = 1. + weightsThreshold = 1.0 goodPix = dcrWeights[0].array > weightsThreshold for weights in dcrWeights[1:]: goodPix = (weights.array > weightsThreshold) & goodPix for subfilter in range(self.config.dcrNumSubfilters): - dcrWeights[subfilter].array[goodPix] = 1./dcrWeights[subfilter].array[goodPix] - dcrWeights[subfilter].array[~goodPix] = 0. + dcrWeights[subfilter].array[goodPix] = 1.0 / dcrWeights[subfilter].array[goodPix] + dcrWeights[subfilter].array[~goodPix] = 0.0 dcrNImages[subfilter].array[~goodPix] = 0 return (dcrNImages, dcrWeights) - def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefList, - statsCtrl, convergenceMetric, - gain, modelWeights, refImage, dcrWeights): + def dcrAssembleSubregion( + self, + dcrModels, + subExposures, + bbox, + dcrBBox, + warpRefList, + statsCtrl, + convergenceMetric, + gain, + modelWeights, + refImage, + dcrWeights, + ): """Assemble the DCR coadd for a sub-region. Build a DCR-matched template for each input exposure, then shift the @@ -762,12 +839,14 @@ def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefLi exposure = subExposures[visit] visitInfo = exposure.getInfo().getVisitInfo() wcs = exposure.getInfo().getWcs() - templateImage = dcrModels.buildMatchedTemplate(exposure=exposure, - bbox=exposure.getBBox(), - order=self.config.imageInterpOrder, - splitSubfilters=self.config.splitSubfilters, - splitThreshold=self.config.splitThreshold, - amplifyModel=self.config.accelerateModel) + templateImage = dcrModels.buildMatchedTemplate( + exposure=exposure, + bbox=exposure.getBBox(), + order=self.config.imageInterpOrder, + splitSubfilters=self.config.splitSubfilters, + splitThreshold=self.config.splitThreshold, + amplifyModel=self.config.accelerateModel, + ) residual = exposure.image.array - templateImage.array # Note that the variance plane here is used to store weights, not # the actual variance @@ -776,15 +855,22 @@ def dcrAssembleSubregion(self, dcrModels, subExposures, bbox, dcrBBox, warpRefLi # This allows the residual for a given subfilter and exposure to be # created on the fly, instead of needing to store them all in # memory. - residualGeneratorList.append(self.dcrResiduals(residual, visitInfo, wcs, - dcrModels.effectiveWavelength, - dcrModels.bandwidth)) - - dcrSubModelOut = self.newModelFromResidual(dcrModels, residualGeneratorList, dcrBBox, statsCtrl, - gain=gain, - modelWeights=modelWeights, - refImage=refImage, - dcrWeights=dcrWeights) + residualGeneratorList.append( + self.dcrResiduals( + residual, visitInfo, wcs, dcrModels.effectiveWavelength, dcrModels.bandwidth + ) + ) + + dcrSubModelOut = self.newModelFromResidual( + dcrModels, + residualGeneratorList, + dcrBBox, + statsCtrl, + gain=gain, + modelWeights=modelWeights, + refImage=refImage, + dcrWeights=dcrWeights, + ) dcrModels.assign(dcrSubModelOut, bbox) def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth): @@ -816,14 +902,27 @@ def dcrResiduals(self, residual, visitInfo, wcs, effectiveWavelength, bandwidth) # Note that `splitSubfilters` is always turned off in the reverse # direction. This option introduces additional blurring if applied to # the residuals. - dcrShift = calculateDcr(visitInfo, wcs, effectiveWavelength, bandwidth, self.config.dcrNumSubfilters, - splitSubfilters=False) + dcrShift = calculateDcr( + visitInfo, + wcs, + effectiveWavelength, + bandwidth, + self.config.dcrNumSubfilters, + splitSubfilters=False, + ) for dcr in dcrShift: - yield applyDcr(filteredResidual, dcr, useInverse=True, splitSubfilters=False, - doPrefilter=False, order=self.config.imageInterpOrder) - - def newModelFromResidual(self, dcrModels, residualGeneratorList, dcrBBox, statsCtrl, - gain, modelWeights, refImage, dcrWeights): + yield applyDcr( + filteredResidual, + dcr, + useInverse=True, + splitSubfilters=False, + doPrefilter=False, + order=self.config.imageInterpOrder, + ) + + def newModelFromResidual( + self, dcrModels, residualGeneratorList, dcrBBox, statsCtrl, gain, modelWeights, refImage, dcrWeights + ): """Calculate a new DcrModel from a set of image residuals. Parameters @@ -866,19 +965,33 @@ def newModelFromResidual(self, dcrModels, residualGeneratorList, dcrBBox, statsC badPixels = ~np.isfinite(newModel.array) newModel.array[badPixels] = model[dcrBBox].array[badPixels] if self.config.regularizeModelIterations > 0: - dcrModels.regularizeModelIter(subfilter, newModel, dcrBBox, - self.config.regularizeModelIterations, - self.config.regularizationWidth) + dcrModels.regularizeModelIter( + subfilter, + newModel, + dcrBBox, + self.config.regularizeModelIterations, + self.config.regularizationWidth, + ) newModelImages.append(newModel) if self.config.regularizeModelFrequency > 0: - dcrModels.regularizeModelFreq(newModelImages, dcrBBox, statsCtrl, - self.config.regularizeModelFrequency, - self.config.regularizationWidth) + dcrModels.regularizeModelFreq( + newModelImages, + dcrBBox, + statsCtrl, + self.config.regularizeModelFrequency, + self.config.regularizationWidth, + ) dcrModels.conditionDcrModel(newModelImages, dcrBBox, gain=gain) self.applyModelWeights(newModelImages, refImage[dcrBBox], modelWeights) - return DcrModel(newModelImages, dcrModels.filter, dcrModels.effectiveWavelength, - dcrModels.bandwidth, dcrModels.psf, - dcrModels.mask, dcrModels.variance) + return DcrModel( + newModelImages, + dcrModels.filter, + dcrModels.effectiveWavelength, + dcrModels.bandwidth, + dcrModels.psf, + dcrModels.mask, + dcrModels.variance, + ) def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weightList, statsCtrl): """Calculate a quality of fit metric for the matched templates. @@ -905,13 +1018,14 @@ def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weigh sub-region. """ significanceImage = np.abs(dcrModels.getReferenceImage(bbox)) - nSigma = 3. - significanceImage += nSigma*dcrModels.calculateNoiseCutoff(dcrModels[1], statsCtrl, - bufferSize=self.bufferSize) + nSigma = 3.0 + significanceImage += nSigma * dcrModels.calculateNoiseCutoff( + dcrModels[1], statsCtrl, bufferSize=self.bufferSize + ) if np.max(significanceImage) == 0: - significanceImage += 1. + significanceImage += 1.0 weight = 0 - metric = 0. + metric = 0.0 metricList = {} for warpExpRef, expWeight in zip(warpRefList, weightList): visit = warpExpRef.dataId["visit"] @@ -919,9 +1033,9 @@ def calculateConvergence(self, dcrModels, subExposures, bbox, warpRefList, weigh singleMetric = self.calculateSingleConvergence(dcrModels, exposure, significanceImage, statsCtrl) metric += singleMetric metricList[visit] = singleMetric - weight += 1. + weight += 1.0 self.log.info("Individual metrics:\n%s", metricList) - return 1.0 if weight == 0.0 else metric/weight + return 1.0 if weight == 0.0 else metric / weight def calculateSingleConvergence(self, dcrModels, exposure, significanceImage, statsCtrl): """Calculate a quality of fit metric for a single matched template. @@ -944,25 +1058,27 @@ def calculateSingleConvergence(self, dcrModels, exposure, significanceImage, sta Quality of fit metric for one exposure, within the sub-region. """ convergeMask = exposure.mask.getPlaneBitMask(self.config.convergenceMaskPlanes) - templateImage = dcrModels.buildMatchedTemplate(exposure=exposure, - bbox=exposure.getBBox(), - order=self.config.imageInterpOrder, - splitSubfilters=self.config.splitSubfilters, - splitThreshold=self.config.splitThreshold, - amplifyModel=self.config.accelerateModel) - diffVals = np.abs(exposure.image.array - templateImage.array)*significanceImage - refVals = np.abs(exposure.image.array + templateImage.array)*significanceImage/2. + templateImage = dcrModels.buildMatchedTemplate( + exposure=exposure, + bbox=exposure.getBBox(), + order=self.config.imageInterpOrder, + splitSubfilters=self.config.splitSubfilters, + splitThreshold=self.config.splitThreshold, + amplifyModel=self.config.accelerateModel, + ) + diffVals = np.abs(exposure.image.array - templateImage.array) * significanceImage + refVals = np.abs(exposure.image.array + templateImage.array) * significanceImage / 2.0 finitePixels = np.isfinite(diffVals) goodMaskPixels = (exposure.mask.array & statsCtrl.getAndMask()) == 0 convergeMaskPixels = exposure.mask.array & convergeMask > 0 usePixels = finitePixels & goodMaskPixels & convergeMaskPixels if np.sum(usePixels) == 0: - metric = 0. + metric = 0.0 else: diffUse = diffVals[usePixels] refUse = refVals[usePixels] - metric = np.sum(diffUse/np.median(diffUse))/np.sum(refUse/np.median(diffUse)) + metric = np.sum(diffUse / np.median(diffUse)) / np.sum(refUse / np.median(diffUse)) return metric def stackCoadd(self, dcrCoadds): @@ -984,8 +1100,17 @@ def stackCoadd(self, dcrCoadds): coaddExposure.maskedImage += coadd.maskedImage return coaddExposure - def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=None, coaddInputs=None, - mask=None, variance=None): + def fillCoadd( + self, + dcrModels, + skyInfo, + warpRefList, + weightList, + calibration=None, + coaddInputs=None, + mask=None, + variance=None, + ): """Create a list of coadd exposures from a list of masked images. Parameters @@ -1017,7 +1142,7 @@ def fillCoadd(self, dcrModels, skyInfo, warpRefList, weightList, calibration=Non refModel = dcrModels.getReferenceImage() for model in dcrModels: if self.config.accelerateModel > 1: - model.array = (model.array - refModel)*self.config.accelerateModel + refModel + model.array = (model.array - refModel) * self.config.accelerateModel + refModel coaddExposure = afwImage.ExposureF(skyInfo.bbox, skyInfo.wcs) if calibration is not None: coaddExposure.setPhotoCalib(calibration) @@ -1074,14 +1199,16 @@ def calculateGain(self, convergenceList, gainList): """ nIter = len(convergenceList) if nIter != len(gainList) + 1: - raise ValueError("convergenceList (%d) must be one element longer than gainList (%d)." - % (len(convergenceList), len(gainList))) + raise ValueError( + "convergenceList (%d) must be one element longer than gainList (%d)." + % (len(convergenceList), len(gainList)) + ) if self.config.baseGain is None: # If ``baseGain`` is not set, calculate it from the number of DCR # subfilters. The more subfilters being modeled, the lower the gain # should be. - baseGain = 1./(self.config.dcrNumSubfilters - 1) + baseGain = 1.0 / (self.config.dcrNumSubfilters - 1) else: baseGain = self.config.baseGain @@ -1093,15 +1220,17 @@ def calculateGain(self, convergenceList, gainList): # should asymptotically approach a final value. We can estimate # that value from the measured changes in convergence weighted by # the gains used in each previous iteration. - estFinalConv = [((1 + gainList[i])*convergenceList[i + 1] - convergenceList[i])/gainList[i] - for i in range(nIter - 1)] + estFinalConv = [ + ((1 + gainList[i]) * convergenceList[i + 1] - convergenceList[i]) / gainList[i] + for i in range(nIter - 1) + ] # The convergence metric is strictly positive, so if the estimated # final convergence is less than zero, force it to zero. estFinalConv = np.array(estFinalConv) estFinalConv[estFinalConv < 0] = 0 # Because the estimate may slowly change over time, only use the # most recent measurements. - estFinalConv = np.median(estFinalConv[max(nIter - 5, 0):]) + estFinalConv = np.median(estFinalConv[max(nIter - 5, 0) :]) lastGain = gainList[-1] lastConv = convergenceList[-2] newConv = convergenceList[-1] @@ -1110,17 +1239,17 @@ def calculateGain(self, convergenceList, gainList): # that the updated model that is actually used is the gain-weighted # average of the new and old model, so the convergence would be # similarly weighted. - predictedConv = (estFinalConv*lastGain + lastConv)/(1. + lastGain) + predictedConv = (estFinalConv * lastGain + lastConv) / (1.0 + lastGain) # If the measured and predicted convergence are very close, that # indicates that our forward model is accurate and we can use a # more aggressive gain. If the measured convergence is # significantly worse (or better!) than predicted, that indicates # that the model is not converging as expected and we should use a # more conservative gain. - delta = (predictedConv - newConv)/((lastConv - estFinalConv)/(1 + lastGain)) + delta = (predictedConv - newConv) / ((lastConv - estFinalConv) / (1 + lastGain)) newGain = 1 - abs(delta) # Average the gains to prevent oscillating solutions. - newGain = (newGain + lastGain)/2. + newGain = (newGain + lastGain) / 2.0 gain = max(baseGain, newGain) else: gain = baseGain @@ -1157,7 +1286,7 @@ def calculateModelWeights(self, dcrModels, dcrBBox): convergeMask = dcrModels.mask.getPlaneBitMask(self.config.convergenceMaskPlanes) convergeMaskPixels = dcrModels.mask[dcrBBox].array & convergeMask > 0 weights = np.zeros_like(dcrModels[0][dcrBBox].array) - weights[convergeMaskPixels] = 1. + weights[convergeMaskPixels] = 1.0 weights = ndimage.gaussian_filter(weights, self.config.modelWeightsWidth) weights /= np.max(weights) return weights @@ -1181,7 +1310,7 @@ def applyModelWeights(self, modelImages, refImage, modelWeights): if self.config.useModelWeights: for model in modelImages: model.array *= modelWeights - model.array += refImage.array*(1. - modelWeights)/self.config.dcrNumSubfilters + model.array += refImage.array * (1.0 - modelWeights) / self.config.dcrNumSubfilters def loadSubExposures(self, bbox, statsCtrl, warpRefList, imageScalerList, spanSetMaskList): """Pre-load sub-regions of a list of exposures. @@ -1211,20 +1340,20 @@ def loadSubExposures(self, bbox, statsCtrl, warpRefList, imageScalerList, spanSe zipIterables = zip(warpRefList, imageScalerList, spanSetMaskList) subExposures = {} for warpExpRef, imageScaler, altMaskSpans in zipIterables: - exposure = warpExpRef.get(parameters={'bbox': bbox}) + exposure = warpExpRef.get(parameters={"bbox": bbox}) visit = warpExpRef.dataId["visit"] if altMaskSpans is not None: self.applyAltMaskPlanes(exposure.mask, altMaskSpans) imageScaler.scaleMaskedImage(exposure.maskedImage) # Note that the variance plane here is used to store weights, not # the actual variance - exposure.variance.array[:, :] = 0. + exposure.variance.array[:, :] = 0.0 # Set the weight of unmasked pixels to 1. - exposure.variance.array[(exposure.mask.array & statsCtrl.getAndMask()) == 0] = 1. + exposure.variance.array[(exposure.mask.array & statsCtrl.getAndMask()) == 0] = 1.0 # Set the image value of masked pixels to zero. # This eliminates needing the mask plane when stacking images in # ``newModelFromResidual`` - exposure.image.array[(exposure.mask.array & statsCtrl.getAndMask()) > 0] = 0. + exposure.image.array[(exposure.mask.array & statsCtrl.getAndMask()) > 0] = 0.0 subExposures[visit] = exposure return subExposures @@ -1244,7 +1373,7 @@ def selectCoaddPsf(self, templateCoadd, warpRefList): psf : `lsst.meas.algorithms.CoaddPsf` The average PSF of the input exposures with the best seeing. """ - sigma2fwhm = 2.*np.sqrt(2.*np.log(2.)) + sigma2fwhm = 2.0 * np.sqrt(2.0 * np.log(2.0)) # Note: ``ccds`` is a `lsst.afw.table.ExposureCatalog` with one entry # per ccd and per visit. If there are multiple ccds, it will have that # many times more elements than ``warpExpRef``. @@ -1252,14 +1381,14 @@ def selectCoaddPsf(self, templateCoadd, warpRefList): templatePsf = templateCoadd.getPsf() # Just need a rough estimate; average positions are fine templateAvgPos = templatePsf.getAveragePosition() - psfRefSize = templatePsf.computeShape(templateAvgPos).getDeterminantRadius()*sigma2fwhm + psfRefSize = templatePsf.computeShape(templateAvgPos).getDeterminantRadius() * sigma2fwhm psfSizes = np.zeros(len(ccds)) ccdVisits = np.array(ccds["visit"]) for warpExpRef in warpRefList: psf = warpExpRef.get(component="psf") visit = warpExpRef.dataId["visit"] psfAvgPos = psf.getAveragePosition() - psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius()*sigma2fwhm + psfSize = psf.computeShape(psfAvgPos).getDeterminantRadius() * sigma2fwhm psfSizes[ccdVisits == visit] = psfSize # Note that the input PSFs include DCR, which should be absent from the # DcrCoadd. The selected PSFs are those that have a FWHM less than or @@ -1267,6 +1396,5 @@ def selectCoaddPsf(self, templateCoadd, warpRefList): # exposures. sizeThreshold = min(np.median(psfSizes), psfRefSize) goodPsfs = psfSizes <= sizeThreshold - psf = measAlg.CoaddPsf(ccds[goodPsfs], templateCoadd.getWcs(), - self.config.coaddPsf.makeControl()) + psf = measAlg.CoaddPsf(ccds[goodPsfs], templateCoadd.getWcs(), self.config.coaddPsf.makeControl()) return psf diff --git a/python/lsst/drp/tasks/forcedPhotCoadd.py b/python/lsst/drp/tasks/forcedPhotCoadd.py index dc3b032d..8a39edbc 100644 --- a/python/lsst/drp/tasks/forcedPhotCoadd.py +++ b/python/lsst/drp/tasks/forcedPhotCoadd.py @@ -33,10 +33,11 @@ __all__ = ("ForcedPhotCoaddConfig", "ForcedPhotCoaddTask") -class ForcedPhotCoaddConnections(pipeBase.PipelineTaskConnections, - dimensions=("band", "skymap", "tract", "patch"), - defaultTemplates={"inputCoaddName": "deep", - "outputCoaddName": "deep"}): +class ForcedPhotCoaddConnections( + pipeBase.PipelineTaskConnections, + dimensions=("band", "skymap", "tract", "patch"), + defaultTemplates={"inputCoaddName": "deep", "outputCoaddName": "deep"}, +): inputSchema = pipeBase.connectionTypes.InitInput( doc="Schema for the input measurement catalogs.", name="{inputCoaddName}Coadd_ref_schema", @@ -63,13 +64,13 @@ class ForcedPhotCoaddConnections(pipeBase.PipelineTaskConnections, doc="Catalog of shapes and positions in the band having forced photometry done", name="{inputCoaddName}Coadd_meas", storageClass="SourceCatalog", - dimensions=("band", "skymap", "tract", "patch") + dimensions=("band", "skymap", "tract", "patch"), ) footprintCatInBand = pipeBase.connectionTypes.Input( doc="Catalog of footprints to attach to sources", name="{inputCoaddName}Coadd_deblendedFlux", storageClass="SourceCatalog", - dimensions=("band", "skymap", "tract", "patch") + dimensions=("band", "skymap", "tract", "patch"), ) scarletModels = pipeBase.connectionTypes.Input( doc="Multiband scarlet models produced by the deblender", @@ -98,11 +99,9 @@ def __init__(self, *, config=None): self.inputs.remove("footprintCatInBand") -class ForcedPhotCoaddConfig(pipeBase.PipelineTaskConfig, - pipelineConnections=ForcedPhotCoaddConnections): +class ForcedPhotCoaddConfig(pipeBase.PipelineTaskConfig, pipelineConnections=ForcedPhotCoaddConnections): measurement = lsst.pex.config.ConfigurableField( - target=ForcedMeasurementTask, - doc="subtask to do forced measurement" + target=ForcedMeasurementTask, doc="subtask to do forced measurement" ) coaddName = lsst.pex.config.Field( doc="coadd name: typically one of deep or goodSeeing", @@ -110,42 +109,40 @@ class ForcedPhotCoaddConfig(pipeBase.PipelineTaskConfig, default="deep", ) doApCorr = lsst.pex.config.Field( - dtype=bool, - default=True, - doc="Run subtask to apply aperture corrections" + dtype=bool, default=True, doc="Run subtask to apply aperture corrections" ) applyApCorr = lsst.pex.config.ConfigurableField( - target=ApplyApCorrTask, - doc="Subtask to apply aperture corrections" + target=ApplyApCorrTask, doc="Subtask to apply aperture corrections" ) catalogCalculation = lsst.pex.config.ConfigurableField( - target=CatalogCalculationTask, - doc="Subtask to run catalogCalculation plugins on catalog" + target=CatalogCalculationTask, doc="Subtask to run catalogCalculation plugins on catalog" ) footprintDatasetName = lsst.pex.config.Field( doc="Dataset (without coadd prefix) that should be used to obtain (Heavy)Footprints for sources. " - "Must have IDs that match those of the reference catalog." - "If None, Footprints will be generated by transforming the reference Footprints.", + "Must have IDs that match those of the reference catalog." + "If None, Footprints will be generated by transforming the reference Footprints.", dtype=str, default="ScarletModelData", - optional=True + optional=True, ) doConserveFlux = lsst.pex.config.Field( dtype=bool, default=True, doc="Whether to use the deblender models as templates to re-distribute the flux " - "from the 'exposure' (True), or to perform measurements on the deblender model footprints. " - "If footprintDatasetName != 'ScarletModelData' then this field is ignored.") + "from the 'exposure' (True), or to perform measurements on the deblender model footprints. " + "If footprintDatasetName != 'ScarletModelData' then this field is ignored.", + ) doStripFootprints = lsst.pex.config.Field( dtype=bool, default=True, doc="Whether to strip footprints from the output catalog before " - "saving to disk. " - "This is usually done when using scarlet models to save disk space.") + "saving to disk. " + "This is usually done when using scarlet models to save disk space.", + ) hasFakes = lsst.pex.config.Field( dtype=bool, default=False, - doc="Should be set to True if fake sources have been inserted into the input data." + doc="Should be set to True if fake sources have been inserted into the input data.", ) idGenerator = SkyMapIdGeneratorConfig.make_field() @@ -158,13 +155,21 @@ def setDefaults(self): self.catalogCalculation.plugins.names = [] self.measurement.copyColumns["id"] = "id" self.measurement.copyColumns["parent"] = "parent" - self.measurement.plugins.names |= ['base_InputCount', 'base_Variance'] - self.measurement.plugins['base_PixelFlags'].masksFpAnywhere = ['CLIPPED', 'SENSOR_EDGE', - 'REJECTED', 'INEXACT_PSF', - 'STREAK'] - self.measurement.plugins['base_PixelFlags'].masksFpCenter = ['CLIPPED', 'SENSOR_EDGE', - 'REJECTED', 'INEXACT_PSF', - 'STREAK'] + self.measurement.plugins.names |= ["base_InputCount", "base_Variance"] + self.measurement.plugins["base_PixelFlags"].masksFpAnywhere = [ + "CLIPPED", + "SENSOR_EDGE", + "REJECTED", + "INEXACT_PSF", + "STREAK", + ] + self.measurement.plugins["base_PixelFlags"].masksFpCenter = [ + "CLIPPED", + "SENSOR_EDGE", + "REJECTED", + "INEXACT_PSF", + "STREAK", + ] class ForcedPhotCoaddTask(pipeBase.PipelineTask): @@ -191,7 +196,7 @@ def __init__(self, refSchema=None, initInputs=None, **kwds): super().__init__(**kwds) if initInputs is not None: - refSchema = initInputs['inputSchema'].schema + refSchema = initInputs["inputSchema"].schema if refSchema is None: raise ValueError("No reference schema provided.") @@ -201,25 +206,27 @@ def __init__(self, refSchema=None, initInputs=None, **kwds): # measurement task, but is passed in by an external caller. if self.config.doApCorr: self.makeSubtask("applyApCorr", schema=self.measurement.schema) - self.makeSubtask('catalogCalculation', schema=self.measurement.schema) + self.makeSubtask("catalogCalculation", schema=self.measurement.schema) self.outputSchema = lsst.afw.table.SourceCatalog(self.measurement.schema) def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) - refCatInBand = inputs.pop('refCatInBand') + refCatInBand = inputs.pop("refCatInBand") if self.config.footprintDatasetName == "ScarletModelData": footprintData = inputs.pop("scarletModels") elif self.config.footprintDatasetName == "DeblendedFlux": footprintData = inputs.pop("footprintCatIndBand") else: footprintData = None - inputs['measCat'], inputs['exposureId'] = self.generateMeasCat(inputRefs.exposure.dataId, - inputs['exposure'], - inputs['refCat'], - refCatInBand, - inputs['refWcs'], - footprintData) + inputs["measCat"], inputs["exposureId"] = self.generateMeasCat( + inputRefs.exposure.dataId, + inputs["exposure"], + inputs["refCat"], + refCatInBand, + inputs["refWcs"], + footprintData, + ) outputs = self.run(**inputs) # Strip HeavyFootprints to save space on disk if self.config.footprintDatasetName == "ScarletModelData" and self.config.doStripFootprints: @@ -264,16 +271,14 @@ def generateMeasCat(self, dataId, exposure, refCat, refCatInBand, refWcs, footpr some sort of mismatch in the two input catalogs) """ id_generator = self.config.idGenerator.apply(dataId) - measCat = self.measurement.generateMeasCat(exposure, refCat, refWcs, - idFactory=id_generator.make_table_id_factory()) + measCat = self.measurement.generateMeasCat( + exposure, refCat, refWcs, idFactory=id_generator.make_table_id_factory() + ) # attach footprints here as this can naturally live inside this method if self.config.footprintDatasetName == "ScarletModelData": # Load the scarlet models self._attachScarletFootprints( - catalog=measCat, - modelData=footprintData, - exposure=exposure, - band=dataId["band"] + catalog=measCat, modelData=footprintData, exposure=exposure, band=dataId["band"] ) else: if self.config.footprintDatasetName is None: @@ -283,9 +288,10 @@ def generateMeasCat(self, dataId, exposure, refCat, refCatInBand, refWcs, footpr for srcRecord in measCat: fpRecord = footprintCat.find(srcRecord.getId()) if fpRecord is None: - raise LookupError("Cannot find Footprint for source {}; please check that {} " - "IDs are compatible with reference source IDs" - .format(srcRecord.getId(), footprintCat)) + raise LookupError( + "Cannot find Footprint for source {}; please check that {} " + "IDs are compatible with reference source IDs".format(srcRecord.getId(), footprintCat) + ) srcRecord.setFootprint(fpRecord.getFootprint()) return measCat, id_generator.catalog_id @@ -318,17 +324,13 @@ def run(self, measCat, exposure, refCat, refWcs, exposureId=None): """ self.measurement.run(measCat, exposure, refCat, refWcs, exposureId=exposureId) if self.config.doApCorr: - self.applyApCorr.run( - catalog=measCat, - apCorrMap=exposure.getInfo().getApCorrMap() - ) + self.applyApCorr.run(catalog=measCat, apCorrMap=exposure.getInfo().getApCorrMap()) self.catalogCalculation.run(measCat) return pipeBase.Struct(measCat=measCat) def _attachScarletFootprints(self, catalog, modelData, exposure, band): - """Attach scarlet models as HeavyFootprints - """ + """Attach scarlet models as HeavyFootprints""" if self.config.doConserveFlux: redistributeImage = exposure else: diff --git a/python/lsst/drp/tasks/gbdesAstrometricFit.py b/python/lsst/drp/tasks/gbdesAstrometricFit.py index e7d8a9c2..7b0c8fe8 100644 --- a/python/lsst/drp/tasks/gbdesAstrometricFit.py +++ b/python/lsst/drp/tasks/gbdesAstrometricFit.py @@ -37,11 +37,12 @@ ReferenceSourceSelectorTask) from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry -__all__ = ['GbdesAstrometricFitConnections', 'GbdesAstrometricFitConfig', 'GbdesAstrometricFitTask'] +__all__ = ["GbdesAstrometricFitConnections", "GbdesAstrometricFitConfig", "GbdesAstrometricFitTask"] -def _make_ref_covariance_matrix(refCat, inputUnit=u.radian, outputCoordUnit=u.marcsec, - outputPMUnit=u.marcsec, version=1): +def _make_ref_covariance_matrix( + refCat, inputUnit=u.radian, outputCoordUnit=u.marcsec, outputPMUnit=u.marcsec, version=1 +): """Make a covariance matrix for the reference catalog including proper motion and parallax. @@ -76,16 +77,18 @@ def _make_ref_covariance_matrix(refCat, inputUnit=u.radian, outputCoordUnit=u.ma # the string in Gaia column names for this # the ordering in the Gaia catalog # and the ordering of the tuples is the order we want in our cov matrix - raErr = (refCat['coord_raErr'] * inputUnit).to(outputCoordUnit).to_value() - decErr = (refCat['coord_decErr'] * inputUnit).to(outputCoordUnit).to_value() - raPMErr = (refCat['pm_raErr'] * inputUnit).to(outputPMUnit).to_value() - decPMErr = (refCat['pm_decErr'] * inputUnit).to(outputPMUnit).to_value() - parallaxErr = (refCat['parallaxErr'] * inputUnit).to(outputPMUnit).to_value() - stdOrder = ((raErr, 'ra', 0), - (decErr, 'dec', 1), - (raPMErr, 'pmra', 3), - (decPMErr, 'pmdec', 4), - (parallaxErr, 'parallax', 2)) + raErr = (refCat["coord_raErr"] * inputUnit).to(outputCoordUnit).to_value() + decErr = (refCat["coord_decErr"] * inputUnit).to(outputCoordUnit).to_value() + raPMErr = (refCat["pm_raErr"] * inputUnit).to(outputPMUnit).to_value() + decPMErr = (refCat["pm_decErr"] * inputUnit).to(outputPMUnit).to_value() + parallaxErr = (refCat["parallaxErr"] * inputUnit).to(outputPMUnit).to_value() + stdOrder = ( + (raErr, "ra", 0), + (decErr, "dec", 1), + (raPMErr, "pmra", 3), + (decPMErr, "pmdec", 4), + (parallaxErr, "parallax", 2), + ) k = 0 for i, pr1 in enumerate(stdOrder): @@ -97,20 +100,20 @@ def _make_ref_covariance_matrix(refCat, inputUnit=u.radian, outputCoordUnit=u.ma else: # diagnonal element cov[:, k] = pr1[0] * pr2[0] - k = k+1 + k = k + 1 elif version == 2: - positionParameters = ['coord_ra', 'coord_dec', 'pm_ra', 'pm_dec', 'parallax'] + positionParameters = ["coord_ra", "coord_dec", "pm_ra", "pm_dec", "parallax"] units = [outputCoordUnit, outputCoordUnit, outputPMUnit, outputPMUnit, outputPMUnit] k = 0 for i, pi in enumerate(positionParameters): for j, pj in enumerate(positionParameters): if i == j: - cov[:, k] = (refCat[f'{pi}Err']**2 * inputUnit**2).to_value(units[j] * units[j]) + cov[:, k] = (refCat[f"{pi}Err"] ** 2 * inputUnit**2).to_value(units[j] * units[j]) elif i > j: - cov[:, k] = (refCat[f'{pj}_{pi}_Cov'] * inputUnit**2).to_value(units[i] * units[j]) + cov[:, k] = (refCat[f"{pj}_{pi}_Cov"] * inputUnit**2).to_value(units[i] * units[j]) else: - cov[:, k] = (refCat[f'{pi}_{pj}_Cov'] * inputUnit**2).to_value(units[i] * units[j]) + cov[:, k] = (refCat[f"{pi}_{pj}_Cov"] * inputUnit**2).to_value(units[i] * units[j]) k += 1 return cov @@ -138,14 +141,14 @@ def _convert_to_ast_polymap_coefficients(coefficients): N = len(coefficients) / 2 # Get the degree of the polynomial by applying the quadratic formula to the # formula for calculating the number of coefficients of the polynomial. - degree = int(-1.5 + 0.5 * (1 + 8 * N)**0.5) + degree = int(-1.5 + 0.5 * (1 + 8 * N) ** 0.5) for outVar in [1, 2]: for i in range(degree + 1): for j in range(degree + 1): if (i + j) > degree: continue - vectorIndex = int(((i+j)*(i+j+1))/2+j + N * (outVar - 1)) + vectorIndex = int(((i + j) * (i + j + 1)) / 2 + j + N * (outVar - 1)) polyArray[vectorIndex, 0] = coefficients[vectorIndex] polyArray[vectorIndex, 1] = outVar polyArray[vectorIndex, 2] = i @@ -169,197 +172,198 @@ def _get_wcs_from_sip(butlerWcs): WCS object in TPV format. """ fits_metadata = butlerWcs.getFitsMetadata() - if not ((fits_metadata.get('CTYPE1') == 'RA---TAN-SIP') - and (fits_metadata.get('CTYPE2') == 'DEC--TAN-SIP')): - raise ValueError(f"CTYPES {fits_metadata.get('CTYPE1')} and {fits_metadata.get('CTYPE2')}" - "do not match SIP convention") + if not ( + (fits_metadata.get("CTYPE1") == "RA---TAN-SIP") and (fits_metadata.get("CTYPE2") == "DEC--TAN-SIP") + ): + raise ValueError( + f"CTYPES {fits_metadata.get('CTYPE1')} and {fits_metadata.get('CTYPE2')}" + "do not match SIP convention" + ) # Correct CRPIX values to correspond to source table pixel indexing # convention - crpix1 = fits_metadata.get('CRPIX1') - crpix2 = fits_metadata.get('CRPIX2') - fits_metadata.set('CRPIX1', crpix1 - 1) - fits_metadata.set('CRPIX2', crpix2 - 1) + crpix1 = fits_metadata.get("CRPIX1") + crpix2 = fits_metadata.get("CRPIX2") + fits_metadata.set("CRPIX1", crpix1 - 1) + fits_metadata.set("CRPIX2", crpix2 - 1) floatDict = {k: fits_metadata[k] for k in fits_metadata if isinstance(fits_metadata[k], (int, float))} - wcs = wcsfit.readTPVFromSIP(floatDict, 'SIP') + wcs = wcsfit.readTPVFromSIP(floatDict, "SIP") return wcs -class GbdesAstrometricFitConnections(pipeBase.PipelineTaskConnections, - dimensions=('skymap', 'tract', 'instrument', 'physical_filter')): +class GbdesAstrometricFitConnections( + pipeBase.PipelineTaskConnections, dimensions=("skymap", "tract", "instrument", "physical_filter") +): """Middleware input/output connections for task data.""" + inputCatalogRefs = pipeBase.connectionTypes.Input( doc="Source table in parquet format, per visit.", - name='preSourceTable_visit', - storageClass='DataFrame', - dimensions=('instrument', 'visit'), + name="preSourceTable_visit", + storageClass="DataFrame", + dimensions=("instrument", "visit"), deferLoad=True, multiple=True, ) inputVisitSummaries = pipeBase.connectionTypes.Input( - doc=("Per-visit consolidated exposure metadata built from calexps. " - "These catalogs use detector id for the id and must be sorted for " - "fast lookups of a detector."), - name='visitSummary', - storageClass='ExposureCatalog', - dimensions=('instrument', 'visit'), + doc=( + "Per-visit consolidated exposure metadata built from calexps. " + "These catalogs use detector id for the id and must be sorted for " + "fast lookups of a detector." + ), + name="visitSummary", + storageClass="ExposureCatalog", + dimensions=("instrument", "visit"), multiple=True, ) referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput( doc="The astrometry reference catalog to match to loaded input catalog sources.", - name='gaia_dr3_20230707', - storageClass='SimpleCatalog', - dimensions=('skypix',), + name="gaia_dr3_20230707", + storageClass="SimpleCatalog", + dimensions=("skypix",), deferLoad=True, multiple=True, ) outputWcs = pipeBase.connectionTypes.Output( - doc=("Per-tract, per-visit world coordinate systems derived from the fitted model." - " These catalogs only contain entries for detectors with an output, and use" - " the detector id for the catalog id, sorted on id for fast lookups of a detector."), - name='gbdesAstrometricFitSkyWcsCatalog', - storageClass='ExposureCatalog', - dimensions=('instrument', 'visit', 'skymap', 'tract'), - multiple=True + doc=( + "Per-tract, per-visit world coordinate systems derived from the fitted model." + " These catalogs only contain entries for detectors with an output, and use" + " the detector id for the catalog id, sorted on id for fast lookups of a detector." + ), + name="gbdesAstrometricFitSkyWcsCatalog", + storageClass="ExposureCatalog", + dimensions=("instrument", "visit", "skymap", "tract"), + multiple=True, ) outputCatalog = pipeBase.connectionTypes.Output( - doc=("Source table with stars used in fit, along with residuals in pixel coordinates and tangent " - "plane coordinates and chisq values."), - name='gbdesAstrometricFit_fitStars', - storageClass='ArrowNumpyDict', - dimensions=('instrument', 'skymap', 'tract', 'physical_filter'), + doc=( + "Source table with stars used in fit, along with residuals in pixel coordinates and tangent " + "plane coordinates and chisq values." + ), + name="gbdesAstrometricFit_fitStars", + storageClass="ArrowNumpyDict", + dimensions=("instrument", "skymap", "tract", "physical_filter"), ) starCatalog = pipeBase.connectionTypes.Output( doc="Star catalog.", - name='gbdesAstrometricFit_starCatalog', - storageClass='ArrowNumpyDict', - dimensions=('instrument', 'skymap', 'tract', 'physical_filter') + name="gbdesAstrometricFit_starCatalog", + storageClass="ArrowNumpyDict", + dimensions=("instrument", "skymap", "tract", "physical_filter"), ) def getSpatialBoundsConnections(self): return ("inputVisitSummaries",) -class GbdesAstrometricFitConfig(pipeBase.PipelineTaskConfig, - pipelineConnections=GbdesAstrometricFitConnections): +class GbdesAstrometricFitConfig( + pipeBase.PipelineTaskConfig, pipelineConnections=GbdesAstrometricFitConnections +): """Configuration for GbdesAstrometricFitTask""" + sourceSelector = sourceSelectorRegistry.makeField( - doc="How to select sources for cross-matching.", - default='science' + doc="How to select sources for cross-matching.", default="science" ) referenceSelector = pexConfig.ConfigurableField( target=ReferenceSourceSelectorTask, doc="How to down-select the loaded astrometry reference catalog.", ) matchRadius = pexConfig.Field( - doc="Matching tolerance between associated objects (arcseconds).", - dtype=float, - default=1.0 + doc="Matching tolerance between associated objects (arcseconds).", dtype=float, default=1.0 ) minMatches = pexConfig.Field( - doc="Number of matches required to keep a source object.", - dtype=int, - default=2 + doc="Number of matches required to keep a source object.", dtype=int, default=2 ) allowSelfMatches = pexConfig.Field( doc="Allow multiple sources from the same visit to be associated with the same object.", dtype=bool, - default=False + default=False, ) sourceFluxType = pexConfig.Field( dtype=str, doc="Source flux field to use in source selection and to get fluxes from the catalog.", - default='apFlux_12_0' + default="apFlux_12_0", ) systematicError = pexConfig.Field( dtype=float, - doc=("Systematic error padding added in quadrature for the science catalogs (marcsec). The default" - "value is equivalent to 0.02 pixels for HSC."), - default=0.0034 + doc=( + "Systematic error padding added in quadrature for the science catalogs (marcsec). The default" + "value is equivalent to 0.02 pixels for HSC." + ), + default=0.0034, ) referenceSystematicError = pexConfig.Field( dtype=float, doc="Systematic error padding added in quadrature for the reference catalog (marcsec).", - default=0.0 + default=0.0, ) modelComponents = pexConfig.ListField( dtype=str, - doc=("List of mappings to apply to transform from pixels to sky, in order of their application." - "Supported options are 'INSTRUMENT/DEVICE' and 'EXPOSURE'."), - default=['INSTRUMENT/DEVICE', 'EXPOSURE'] + doc=( + "List of mappings to apply to transform from pixels to sky, in order of their application." + "Supported options are 'INSTRUMENT/DEVICE' and 'EXPOSURE'." + ), + default=["INSTRUMENT/DEVICE", "EXPOSURE"], ) deviceModel = pexConfig.ListField( dtype=str, - doc=("List of mappings to apply to transform from detector pixels to intermediate frame. Map names" - "should match the format 'BAND/DEVICE/'."), - default=['BAND/DEVICE/poly'] + doc=( + "List of mappings to apply to transform from detector pixels to intermediate frame. Map names" + "should match the format 'BAND/DEVICE/'." + ), + default=["BAND/DEVICE/poly"], ) exposureModel = pexConfig.ListField( dtype=str, - doc=("List of mappings to apply to transform from intermediate frame to sky coordinates. Map names" - "should match the format 'EXPOSURE/'."), - default=['EXPOSURE/poly'] - ) - devicePolyOrder = pexConfig.Field( - dtype=int, - doc="Order of device polynomial model.", - default=4 - ) - exposurePolyOrder = pexConfig.Field( - dtype=int, - doc="Order of exposure polynomial model.", - default=6 - ) - fitProperMotion = pexConfig.Field( - dtype=bool, - doc="Fit the proper motions of the objects.", - default=False + doc=( + "List of mappings to apply to transform from intermediate frame to sky coordinates. Map names" + "should match the format 'EXPOSURE/'." + ), + default=["EXPOSURE/poly"], ) + devicePolyOrder = pexConfig.Field(dtype=int, doc="Order of device polynomial model.", default=4) + exposurePolyOrder = pexConfig.Field(dtype=int, doc="Order of exposure polynomial model.", default=6) + fitProperMotion = pexConfig.Field(dtype=bool, doc="Fit the proper motions of the objects.", default=False) excludeNonPMObjects = pexConfig.Field( - dtype=bool, - doc="Exclude reference objects without proper motion/parallax information.", - default=True + dtype=bool, doc="Exclude reference objects without proper motion/parallax information.", default=True ) fitReserveFraction = pexConfig.Field( - dtype=float, - default=0.2, - doc="Fraction of objects to reserve from fit for validation." + dtype=float, default=0.2, doc="Fraction of objects to reserve from fit for validation." ) fitReserveRandomSeed = pexConfig.Field( dtype=int, doc="Set the random seed for selecting data points to reserve from the fit for validation.", - default=1234 + default=1234, ) def setDefaults(self): # Use only stars because aperture fluxes of galaxies are biased and # depend on seeing. - self.sourceSelector['science'].doUnresolved = True - self.sourceSelector['science'].unresolved.name = 'extendedness' + self.sourceSelector["science"].doUnresolved = True + self.sourceSelector["science"].unresolved.name = "extendedness" # Use only isolated sources. - self.sourceSelector['science'].doIsolated = True - self.sourceSelector['science'].isolated.parentName = 'parentSourceId' - self.sourceSelector['science'].isolated.nChildName = 'deblend_nChild' + self.sourceSelector["science"].doIsolated = True + self.sourceSelector["science"].isolated.parentName = "parentSourceId" + self.sourceSelector["science"].isolated.nChildName = "deblend_nChild" # Do not use either flux or centroid measurements with flags, # chosen from the usual QA flags for stars. - self.sourceSelector['science'].doFlags = True - badFlags = ['pixelFlags_edge', - 'pixelFlags_saturated', - 'pixelFlags_interpolatedCenter', - 'pixelFlags_interpolated', - 'pixelFlags_crCenter', - 'pixelFlags_bad', - 'hsmPsfMoments_flag', - f'{self.sourceFluxType}_flag', - ] - self.sourceSelector['science'].flags.bad = badFlags + self.sourceSelector["science"].doFlags = True + badFlags = [ + "pixelFlags_edge", + "pixelFlags_saturated", + "pixelFlags_interpolatedCenter", + "pixelFlags_interpolated", + "pixelFlags_crCenter", + "pixelFlags_bad", + "hsmPsfMoments_flag", + f"{self.sourceFluxType}_flag", + ] + self.sourceSelector["science"].flags.bad = badFlags # Use only primary sources. - self.sourceSelector['science'].doRequirePrimary = True + self.sourceSelector["science"].doRequirePrimary = True def validate(self): super().validate() @@ -367,14 +371,20 @@ def validate(self): # Check if all components of the device and exposure models are # supported. for component in self.deviceModel: - if not (('poly' in component.lower()) or ('identity' in component.lower())): - raise pexConfig.FieldValidationError(GbdesAstrometricFitConfig.deviceModel, self, - f'deviceModel component {component} is not supported.') + if not (("poly" in component.lower()) or ("identity" in component.lower())): + raise pexConfig.FieldValidationError( + GbdesAstrometricFitConfig.deviceModel, + self, + f"deviceModel component {component} is not supported.", + ) for component in self.exposureModel: - if not (('poly' in component.lower()) or ('identity' in component.lower())): - raise pexConfig.FieldValidationError(GbdesAstrometricFitConfig.exposureModel, self, - f'exposureModel component {component} is not supported.') + if not (("poly" in component.lower()) or ("identity" in component.lower())): + raise pexConfig.FieldValidationError( + GbdesAstrometricFitConfig.exposureModel, + self, + f"exposureModel component {component} is not supported.", + ) class GbdesAstrometricFitTask(pipeBase.PipelineTask): @@ -383,52 +393,55 @@ class GbdesAstrometricFitTask(pipeBase.PipelineTask): """ ConfigClass = GbdesAstrometricFitConfig - _DefaultName = 'gbdesAstrometricFit' + _DefaultName = "gbdesAstrometricFit" def __init__(self, **kwargs): super().__init__(**kwargs) - self.makeSubtask('sourceSelector') - self.makeSubtask('referenceSelector') + self.makeSubtask("sourceSelector") + self.makeSubtask("referenceSelector") def runQuantum(self, butlerQC, inputRefs, outputRefs): # We override runQuantum to set up the refObjLoaders inputs = butlerQC.get(inputRefs) - instrumentName = butlerQC.quantum.dataId['instrument'] + instrumentName = butlerQC.quantum.dataId["instrument"] # Ensure the inputs are in a consistent order - inputCatVisits = np.array([inputCat.dataId['visit'] for inputCat in inputs['inputCatalogRefs']]) - inputs['inputCatalogRefs'] = [inputs['inputCatalogRefs'][v] for v in inputCatVisits.argsort()] - inputSumVisits = np.array([inputSum[0]['visit'] for inputSum in inputs['inputVisitSummaries']]) - inputs['inputVisitSummaries'] = [inputs['inputVisitSummaries'][v] for v in inputSumVisits.argsort()] - inputRefHtm7s = np.array([inputRefCat.dataId['htm7'] for inputRefCat in inputRefs.referenceCatalog]) + inputCatVisits = np.array([inputCat.dataId["visit"] for inputCat in inputs["inputCatalogRefs"]]) + inputs["inputCatalogRefs"] = [inputs["inputCatalogRefs"][v] for v in inputCatVisits.argsort()] + inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]]) + inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()] + inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog]) inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()] - inputRefCats = np.array([inputRefCat.dataId['htm7'] for inputRefCat in inputs['referenceCatalog']]) - inputs['referenceCatalog'] = [inputs['referenceCatalog'][v] for v in inputRefCats.argsort()] + inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]]) + inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()] - sampleRefCat = inputs['referenceCatalog'][0].get() - refEpoch = sampleRefCat[0]['epoch'] + sampleRefCat = inputs["referenceCatalog"][0].get() + refEpoch = sampleRefCat[0]["epoch"] refConfig = LoadReferenceObjectsConfig() - refConfig.anyFilterMapsToThis = 'phot_g_mean' + refConfig.anyFilterMapsToThis = "phot_g_mean" refConfig.requireProperMotion = True - refObjectLoader = ReferenceObjectLoader(dataIds=[ref.datasetRef.dataId - for ref in inputRefCatRefs], - refCats=inputs.pop('referenceCatalog'), - config=refConfig, - log=self.log) + refObjectLoader = ReferenceObjectLoader( + dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs], + refCats=inputs.pop("referenceCatalog"), + config=refConfig, + log=self.log, + ) - output = self.run(**inputs, instrumentName=instrumentName, refEpoch=refEpoch, - refObjectLoader=refObjectLoader) + output = self.run( + **inputs, instrumentName=instrumentName, refEpoch=refEpoch, refObjectLoader=refObjectLoader + ) for outputRef in outputRefs.outputWcs: - visit = outputRef.dataId['visit'] + visit = outputRef.dataId["visit"] butlerQC.put(output.outputWCSs[visit], outputRef) butlerQC.put(output.outputCatalog, outputRefs.outputCatalog) butlerQC.put(output.starCatalog, outputRefs.starCatalog) - def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch=None, - refObjectLoader=None): + def run( + self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch=None, refObjectLoader=None + ): """Run the WCS fit for a given set of visits Parameters @@ -462,8 +475,9 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch instrument = wcsfit.Instrument(instrumentName) # Get RA, Dec, MJD, etc., for the input visits - exposureInfo, exposuresHelper, extensionInfo = self._get_exposure_info(inputVisitSummaries, - instrument) + exposureInfo, exposuresHelper, extensionInfo = self._get_exposure_info( + inputVisitSummaries, instrument + ) # Get information about the extent of the input visits fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummaries, exposureInfo.medianEpoch) @@ -471,18 +485,24 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch self.log.info("Load catalogs and associate sources") # Set up class to associate sources into matches using a # friends-of-friends algorithm - associations = wcsfit.FoFClass(fields, [instrument], exposuresHelper, - [fieldRadius.asDegrees()], - (self.config.matchRadius * u.arcsec).to(u.degree).value) + associations = wcsfit.FoFClass( + fields, + [instrument], + exposuresHelper, + [fieldRadius.asDegrees()], + (self.config.matchRadius * u.arcsec).to(u.degree).value, + ) # Add the reference catalog to the associator - medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format='decimalyear').mjd - refObjects, refCovariance = self._load_refcat(associations, refObjectLoader, fieldCenter, fieldRadius, - extensionInfo, epoch=medianEpoch) + medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="decimalyear").mjd + refObjects, refCovariance = self._load_refcat( + associations, refObjectLoader, fieldCenter, fieldRadius, extensionInfo, epoch=medianEpoch + ) # Add the science catalogs and associate new sources as they are added - sourceIndices, usedColumns = self._load_catalogs_and_associate(associations, inputCatalogRefs, - extensionInfo) + sourceIndices, usedColumns = self._load_catalogs_and_associate( + associations, inputCatalogRefs, extensionInfo + ) self.log.info("Fit the WCSs") # Set up a YAML-type string using the config variables and a sample @@ -501,33 +521,42 @@ def run(self, inputCatalogRefs, inputVisitSummaries, instrumentName="", refEpoch verbose = 2 # Set up the WCS-fitting class using the results of the FOF associator - wcsf = wcsfit.WCSFit(fields, [instrument], exposuresHelper, - extensionInfo.visitIndex, extensionInfo.detectorIndex, - inputYAML, extensionInfo.wcs, associations.sequence, associations.extn, - associations.obj, sysErr=self.config.systematicError, - refSysErr=self.config.referenceSystematicError, - usePM=self.config.fitProperMotion, - verbose=verbose) + wcsf = wcsfit.WCSFit( + fields, + [instrument], + exposuresHelper, + extensionInfo.visitIndex, + extensionInfo.detectorIndex, + inputYAML, + extensionInfo.wcs, + associations.sequence, + associations.extn, + associations.obj, + sysErr=self.config.systematicError, + refSysErr=self.config.referenceSystematicError, + usePM=self.config.fitProperMotion, + verbose=verbose, + ) # Add the science and reference sources self._add_objects(wcsf, inputCatalogRefs, sourceIndices, extensionInfo, usedColumns) self._add_ref_objects(wcsf, refObjects, refCovariance, extensionInfo) # Do the WCS fit - wcsf.fit(reserveFraction=self.config.fitReserveFraction, - randomNumberSeed=self.config.fitReserveRandomSeed) + wcsf.fit( + reserveFraction=self.config.fitReserveFraction, randomNumberSeed=self.config.fitReserveRandomSeed + ) self.log.info("WCS fitting done") outputWCSs = self._make_outputs(wcsf, inputVisitSummaries, exposureInfo) outputCatalog = wcsf.getOutputCatalog() starCatalog = wcsf.getStarCatalog() - return pipeBase.Struct(outputWCSs=outputWCSs, - fitModel=wcsf, - outputCatalog=outputCatalog, - starCatalog=starCatalog) + return pipeBase.Struct( + outputWCSs=outputWCSs, fitModel=wcsf, outputCatalog=outputCatalog, starCatalog=starCatalog + ) - def _prep_sky(self, inputVisitSummaries, epoch, fieldName='Field'): + def _prep_sky(self, inputVisitSummaries, epoch, fieldName="Field"): """Get center and radius of the input tract. This assumes that all visits will be put into the same `wcsfit.Field` and fit together. @@ -551,8 +580,10 @@ def _prep_sky(self, inputVisitSummaries, epoch, fieldName='Field'): """ allDetectorCorners = [] for visSum in inputVisitSummaries: - detectorCorners = [lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector() for (ra, dec) - in zip(visSum['raCorners'].ravel(), visSum['decCorners'].ravel())] + detectorCorners = [ + lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector() + for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel()) + ] allDetectorCorners.extend(detectorCorners) boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(allDetectorCorners).getBoundingCircle() center = lsst.geom.SpherePoint(boundingCircle.getCenter()) @@ -566,8 +597,9 @@ def _prep_sky(self, inputVisitSummaries, epoch, fieldName='Field'): return fields, center, radius - def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, instrumentNumber=0, - refEpoch=None): + def _get_exposure_info( + self, inputVisitSummaries, instrument, fieldNumber=0, instrumentNumber=0, refEpoch=None + ): """Get various information about the input visits to feed to the fitting routines. @@ -636,7 +668,7 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins # Get information for all the science visits for v, visitSummary in enumerate(inputVisitSummaries): visitInfo = visitSummary[0].getVisitInfo() - visit = visitSummary[0]['visit'] + visit = visitSummary[0]["visit"] visits.append(visit) exposureNames.append(str(visit)) raDec = visitInfo.getBoresightRaDec() @@ -652,17 +684,18 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins obsLat = visitInfo.observatory.getLatitude().asDegrees() obsElev = visitInfo.observatory.getElevation() earthLocation = astropy.coordinates.EarthLocation.from_geodetic(obsLon, obsLat, obsElev) - observatory_gcrs = earthLocation.get_gcrs(astropy.time.Time(obsMJD, format='mjd')) + observatory_gcrs = earthLocation.get_gcrs(astropy.time.Time(obsMJD, format="mjd")) observatory_icrs = observatory_gcrs.transform_to(astropy.coordinates.ICRS()) # We want the position in AU in Cartesian coordinates observatories.append(observatory_icrs.cartesian.xyz.to(u.AU).value) for row in visitSummary: - detector = row['id'] + detector = row["id"] if detector not in detectors: detectors.append(detector) - detectorBounds = wcsfit.Bounds(row['bbox_min_x'], row['bbox_max_x'], - row['bbox_min_y'], row['bbox_max_y']) + detectorBounds = wcsfit.Bounds( + row["bbox_min_x"], row["bbox_max_x"], row["bbox_min_y"], row["bbox_max_y"] + ) instrument.addDevice(str(detector), detectorBounds) detectorIndex = np.flatnonzero(detector == np.array(detectors))[0] @@ -670,7 +703,7 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins extensionDetectorIndices.append(detectorIndex) extensionVisits.append(visit) extensionDetectors.append(detector) - extensionType.append('SCIENCE') + extensionType.append("SCIENCE") wcs = row.getWcs() wcss.append(_get_wcs_from_sip(wcs)) @@ -681,11 +714,11 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins # Set the reference epoch to be the median of the science visits. # The reference catalog will be shifted to this date. medianMJD = np.median(mjds) - medianEpoch = astropy.time.Time(medianMJD, format='mjd').decimalyear + medianEpoch = astropy.time.Time(medianMJD, format="mjd").decimalyear # Add information for the reference catalog. Most of the values are # not used. - exposureNames.append('REFERENCE') + exposureNames.append("REFERENCE") visits.append(-1) fieldNumbers.append(0) if self.config.fitProperMotion: @@ -700,44 +733,47 @@ def _get_exposure_info(self, inputVisitSummaries, instrument, fieldNumber=0, ins observatories.append(np.array([0, 0, 0])) identity = wcsfit.IdentityMap() icrs = wcsfit.SphericalICRS() - refWcs = wcsfit.Wcs(identity, icrs, 'Identity', np.pi / 180.) + refWcs = wcsfit.Wcs(identity, icrs, "Identity", np.pi / 180.0) wcss.append(refWcs) extensionVisitIndices.append(len(exposureNames) - 1) extensionDetectorIndices.append(-1) # REFERENCE device must be -1 extensionVisits.append(-1) extensionDetectors.append(-1) - extensionType.append('REFERENCE') + extensionType.append("REFERENCE") # Make a table of information to use elsewhere in the class - extensionInfo = pipeBase.Struct(visit=np.array(extensionVisits), - detector=np.array(extensionDetectors), - visitIndex=np.array(extensionVisitIndices), - detectorIndex=np.array(extensionDetectorIndices), - wcs=np.array(wcss), - extensionType=np.array(extensionType)) + extensionInfo = pipeBase.Struct( + visit=np.array(extensionVisits), + detector=np.array(extensionDetectors), + visitIndex=np.array(extensionVisitIndices), + detectorIndex=np.array(extensionDetectorIndices), + wcs=np.array(wcss), + extensionType=np.array(extensionType), + ) # Make the exposureHelper object to use in the fitting routines - exposuresHelper = wcsfit.ExposuresHelper(exposureNames, - fieldNumbers, - instrumentNumbers, - ras, - decs, - airmasses, - exposureTimes, - mjds, - observatories) - - exposureInfo = pipeBase.Struct(visits=visits, - detectors=detectors, - ras=ras, - decs=decs, - medianEpoch=medianEpoch) + exposuresHelper = wcsfit.ExposuresHelper( + exposureNames, + fieldNumbers, + instrumentNumbers, + ras, + decs, + airmasses, + exposureTimes, + mjds, + observatories, + ) + + exposureInfo = pipeBase.Struct( + visits=visits, detectors=detectors, ras=ras, decs=decs, medianEpoch=medianEpoch + ) return exposureInfo, exposuresHelper, extensionInfo - def _load_refcat(self, associations, refObjectLoader, center, radius, extensionInfo, epoch=None, - fieldIndex=0): + def _load_refcat( + self, associations, refObjectLoader, center, radius, extensionInfo, epoch=None, fieldIndex=0 + ): """Load the reference catalog and add reference objects to the `wcsfit.FoFClass` object. @@ -766,7 +802,7 @@ def _load_refcat(self, associations, refObjectLoader, center, radius, extensionI refCovariance : `list` of `float` Flattened output covariance matrix. """ - formattedEpoch = astropy.time.Time(epoch, format='mjd') + formattedEpoch = astropy.time.Time(epoch, format="mjd") refFilter = refObjectLoader.config.anyFilterMapsToThis skyCircle = refObjectLoader.loadSkyCircle(center, radius, refFilter, epoch=formattedEpoch) @@ -779,54 +815,68 @@ def _load_refcat(self, associations, refObjectLoader, center, radius, extensionI refCat = selected.sourceCat # In Gaia DR3, missing values are denoted by NaNs. - finiteInd = np.isfinite(refCat['coord_ra']) & np.isfinite(refCat['coord_dec']) + finiteInd = np.isfinite(refCat["coord_ra"]) & np.isfinite(refCat["coord_dec"]) refCat = refCat[finiteInd] if self.config.excludeNonPMObjects: # Gaia DR2 has zeros for missing data, while Gaia DR3 has NaNs: - hasPM = ((refCat['pm_raErr'] != 0) & np.isfinite(refCat['pm_raErr']) - & np.isfinite(refCat['pm_decErr'])) + hasPM = ( + (refCat["pm_raErr"] != 0) & np.isfinite(refCat["pm_raErr"]) & np.isfinite(refCat["pm_decErr"]) + ) refCat = refCat[hasPM] - ra = (refCat['coord_ra'] * u.radian).to(u.degree).to_value().tolist() - dec = (refCat['coord_dec'] * u.radian).to(u.degree).to_value().tolist() - raCov = ((refCat['coord_raErr'] * u.radian).to(u.degree).to_value()**2).tolist() - decCov = ((refCat['coord_decErr'] * u.radian).to(u.degree).to_value()**2).tolist() + ra = (refCat["coord_ra"] * u.radian).to(u.degree).to_value().tolist() + dec = (refCat["coord_dec"] * u.radian).to(u.degree).to_value().tolist() + raCov = ((refCat["coord_raErr"] * u.radian).to(u.degree).to_value() ** 2).tolist() + decCov = ((refCat["coord_decErr"] * u.radian).to(u.degree).to_value() ** 2).tolist() # Get refcat version from refcat metadata refCatMetadata = refObjectLoader.refCats[0].get().getMetadata() - refCatVersion = refCatMetadata['REFCAT_FORMAT_VERSION'] + refCatVersion = refCatMetadata["REFCAT_FORMAT_VERSION"] if refCatVersion == 2: - raDecCov = (refCat['coord_ra_coord_dec_Cov'] * u.radian**2).to(u.degree**2).to_value().tolist() + raDecCov = ( + (refCat["coord_ra_coord_dec_Cov"] * u.radian**2).to(u.degree**2).to_value().tolist() + ) else: raDecCov = np.zeros(len(ra)) - refObjects = {'ra': ra, 'dec': dec, 'raCov': raCov, 'decCov': decCov, 'raDecCov': raDecCov} + refObjects = {"ra": ra, "dec": dec, "raCov": raCov, "decCov": decCov, "raDecCov": raDecCov} refCovariance = [] if self.config.fitProperMotion: - raPM = (refCat['pm_ra'] * u.radian).to(u.marcsec).to_value().tolist() - decPM = (refCat['pm_dec'] * u.radian).to(u.marcsec).to_value().tolist() - parallax = (refCat['parallax'] * u.radian).to(u.marcsec).to_value().tolist() + raPM = (refCat["pm_ra"] * u.radian).to(u.marcsec).to_value().tolist() + decPM = (refCat["pm_dec"] * u.radian).to(u.marcsec).to_value().tolist() + parallax = (refCat["parallax"] * u.radian).to(u.marcsec).to_value().tolist() cov = _make_ref_covariance_matrix(refCat, version=refCatVersion) - pmDict = {'raPM': raPM, 'decPM': decPM, 'parallax': parallax} + pmDict = {"raPM": raPM, "decPM": decPM, "parallax": parallax} refObjects.update(pmDict) refCovariance = cov - extensionIndex = np.flatnonzero(extensionInfo.extensionType == 'REFERENCE')[0] + extensionIndex = np.flatnonzero(extensionInfo.extensionType == "REFERENCE")[0] visitIndex = extensionInfo.visitIndex[extensionIndex] detectorIndex = extensionInfo.detectorIndex[extensionIndex] instrumentIndex = -1 # -1 indicates the reference catalog refWcs = extensionInfo.wcs[extensionIndex] - associations.addCatalog(refWcs, 'STELLAR', visitIndex, fieldIndex, instrumentIndex, detectorIndex, - extensionIndex, np.ones(len(refCat), dtype=bool), - ra, dec, np.arange(len(ra))) + associations.addCatalog( + refWcs, + "STELLAR", + visitIndex, + fieldIndex, + instrumentIndex, + detectorIndex, + extensionIndex, + np.ones(len(refCat), dtype=bool), + ra, + dec, + np.arange(len(ra)), + ) return refObjects, refCovariance - def _load_catalogs_and_associate(self, associations, inputCatalogRefs, extensionInfo, - fieldIndex=0, instrumentIndex=0): + def _load_catalogs_and_associate( + self, associations, inputCatalogRefs, extensionInfo, fieldIndex=0, instrumentIndex=0 + ): """Load the science catalogs and add the sources to the associator class `wcsfit.FoFClass`, associating them into matches as you go. @@ -853,8 +903,19 @@ class `wcsfit.FoFClass`, associating them into matches as you go. columns : `list` of `str` List of columns needed from source tables. """ - columns = ['detector', 'sourceId', 'x', 'xErr', 'y', 'yErr', 'ixx', 'iyy', 'ixy', - f'{self.config.sourceFluxType}_instFlux', f'{self.config.sourceFluxType}_instFluxErr'] + columns = [ + "detector", + "sourceId", + "x", + "xErr", + "y", + "yErr", + "ixx", + "iyy", + "ixy", + f"{self.config.sourceFluxType}_instFlux", + f"{self.config.sourceFluxType}_instFluxErr", + ] if self.sourceSelector.config.doFlags: columns.extend(self.sourceSelector.config.flags.bad) if self.sourceSelector.config.doUnresolved: @@ -867,25 +928,27 @@ class `wcsfit.FoFClass`, associating them into matches as you go. sourceIndices = [None] * len(extensionInfo.visit) for inputCatalogRef in inputCatalogRefs: - visit = inputCatalogRef.dataId['visit'] - inputCatalog = inputCatalogRef.get(parameters={'columns': columns}) + visit = inputCatalogRef.dataId["visit"] + inputCatalog = inputCatalogRef.get(parameters={"columns": columns}) # Get a sorted array of detector names - detectors = np.unique(inputCatalog['detector']) + detectors = np.unique(inputCatalog["detector"]) for detector in detectors: - detectorSources = inputCatalog[inputCatalog['detector'] == detector] - xCov = detectorSources['xErr']**2 - yCov = detectorSources['yErr']**2 - xyCov = (detectorSources['ixy'] * (xCov + yCov) - / (detectorSources['ixx'] + detectorSources['iyy'])) + detectorSources = inputCatalog[inputCatalog["detector"] == detector] + xCov = detectorSources["xErr"] ** 2 + yCov = detectorSources["yErr"] ** 2 + xyCov = ( + detectorSources["ixy"] * (xCov + yCov) / (detectorSources["ixx"] + detectorSources["iyy"]) + ) # Remove sources with bad shape measurements goodShapes = xyCov**2 <= (xCov * yCov) selected = self.sourceSelector.run(detectorSources) goodInds = selected.selected & goodShapes isStar = np.ones(goodInds.sum()) - extensionIndex = np.flatnonzero((extensionInfo.visit == visit) - & (extensionInfo.detector == detector))[0] + extensionIndex = np.flatnonzero( + (extensionInfo.visit == visit) & (extensionInfo.detector == detector) + )[0] detectorIndex = extensionInfo.detectorIndex[extensionIndex] visitIndex = extensionInfo.visitIndex[extensionIndex] @@ -894,14 +957,23 @@ class `wcsfit.FoFClass`, associating them into matches as you go. wcs = extensionInfo.wcs[extensionIndex] associations.reprojectWCS(wcs, fieldIndex) - associations.addCatalog(wcs, 'STELLAR', visitIndex, fieldIndex, - instrumentIndex, detectorIndex, extensionIndex, isStar, - detectorSources[goodInds]['x'].to_list(), - detectorSources[goodInds]['y'].to_list(), - np.arange(goodInds.sum())) - - associations.sortMatches(fieldIndex, minMatches=self.config.minMatches, - allowSelfMatches=self.config.allowSelfMatches) + associations.addCatalog( + wcs, + "STELLAR", + visitIndex, + fieldIndex, + instrumentIndex, + detectorIndex, + extensionIndex, + isStar, + detectorSources[goodInds]["x"].to_list(), + detectorSources[goodInds]["y"].to_list(), + np.arange(goodInds.sum()), + ) + + associations.sortMatches( + fieldIndex, minMatches=self.config.minMatches, allowSelfMatches=self.config.allowSelfMatches + ) return sourceIndices, columns @@ -922,50 +994,53 @@ def make_yaml(self, inputVisitSummary, inputFile=None): YAML object containing the model description. """ if inputFile is not None: - inputYAML = wcsfit.YAMLCollector(inputFile, 'PixelMapCollection') + inputYAML = wcsfit.YAMLCollector(inputFile, "PixelMapCollection") else: - inputYAML = wcsfit.YAMLCollector('', 'PixelMapCollection') + inputYAML = wcsfit.YAMLCollector("", "PixelMapCollection") inputDict = {} - modelComponents = ['INSTRUMENT/DEVICE', 'EXPOSURE'] - baseMap = {'Type': 'Composite', 'Elements': modelComponents} - inputDict['EXPOSURE/DEVICE/base'] = baseMap + modelComponents = ["INSTRUMENT/DEVICE", "EXPOSURE"] + baseMap = {"Type": "Composite", "Elements": modelComponents} + inputDict["EXPOSURE/DEVICE/base"] = baseMap - xMin = str(inputVisitSummary['bbox_min_x'].min()) - xMax = str(inputVisitSummary['bbox_max_x'].max()) - yMin = str(inputVisitSummary['bbox_min_y'].min()) - yMax = str(inputVisitSummary['bbox_max_y'].max()) + xMin = str(inputVisitSummary["bbox_min_x"].min()) + xMax = str(inputVisitSummary["bbox_max_x"].max()) + yMin = str(inputVisitSummary["bbox_min_y"].min()) + yMax = str(inputVisitSummary["bbox_max_y"].max()) - deviceModel = {'Type': 'Composite', 'Elements': self.config.deviceModel.list()} - inputDict['INSTRUMENT/DEVICE'] = deviceModel + deviceModel = {"Type": "Composite", "Elements": self.config.deviceModel.list()} + inputDict["INSTRUMENT/DEVICE"] = deviceModel for component in self.config.deviceModel: - if 'poly' in component.lower(): - componentDict = {'Type': 'Poly', - 'XPoly': {'OrderX': self.config.devicePolyOrder, - 'SumOrder': True}, - 'YPoly': {'OrderX': self.config.devicePolyOrder, - 'SumOrder': True}, - 'XMin': xMin, 'XMax': xMax, 'YMin': yMin, 'YMax': yMax} - elif 'identity' in component.lower(): - componentDict = {'Type': 'Identity'} + if "poly" in component.lower(): + componentDict = { + "Type": "Poly", + "XPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True}, + "YPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True}, + "XMin": xMin, + "XMax": xMax, + "YMin": yMin, + "YMax": yMax, + } + elif "identity" in component.lower(): + componentDict = {"Type": "Identity"} inputDict[component] = componentDict - exposureModel = {'Type': 'Composite', 'Elements': self.config.exposureModel.list()} - inputDict['EXPOSURE'] = exposureModel + exposureModel = {"Type": "Composite", "Elements": self.config.exposureModel.list()} + inputDict["EXPOSURE"] = exposureModel for component in self.config.exposureModel: - if 'poly' in component.lower(): - componentDict = {'Type': 'Poly', - 'XPoly': {'OrderX': self.config.exposurePolyOrder, - 'SumOrder': 'true'}, - 'YPoly': {'OrderX': self.config.exposurePolyOrder, - 'SumOrder': 'true'}} - elif 'identity' in component.lower(): - componentDict = {'Type': 'Identity'} + if "poly" in component.lower(): + componentDict = { + "Type": "Poly", + "XPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"}, + "YPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"}, + } + elif "identity" in component.lower(): + componentDict = {"Type": "Identity"} inputDict[component] = componentDict inputYAML.addInput(yaml.dump(inputDict)) - inputYAML.addInput('Identity:\n Type: Identity\n') + inputYAML.addInput("Identity:\n Type: Identity\n") return inputYAML @@ -987,27 +1062,32 @@ def _add_objects(self, wcsf, inputCatalogRefs, sourceIndices, extensionInfo, col List of columns needed from source tables. """ for inputCatalogRef in inputCatalogRefs: - visit = inputCatalogRef.dataId['visit'] - inputCatalog = inputCatalogRef.get(parameters={'columns': columns}) - detectors = np.unique(inputCatalog['detector']) + visit = inputCatalogRef.dataId["visit"] + inputCatalog = inputCatalogRef.get(parameters={"columns": columns}) + detectors = np.unique(inputCatalog["detector"]) for detector in detectors: - detectorSources = inputCatalog[inputCatalog['detector'] == detector] + detectorSources = inputCatalog[inputCatalog["detector"] == detector] - extensionIndex = np.flatnonzero((extensionInfo.visit == visit) - & (extensionInfo.detector == detector))[0] + extensionIndex = np.flatnonzero( + (extensionInfo.visit == visit) & (extensionInfo.detector == detector) + )[0] sourceCat = detectorSources[sourceIndices[extensionIndex]] - xCov = sourceCat['xErr']**2 - yCov = sourceCat['yErr']**2 - xyCov = (sourceCat['ixy'] * (xCov + yCov) - / (sourceCat['ixx'] + sourceCat['iyy'])) + xCov = sourceCat["xErr"] ** 2 + yCov = sourceCat["yErr"] ** 2 + xyCov = sourceCat["ixy"] * (xCov + yCov) / (sourceCat["ixx"] + sourceCat["iyy"]) # TODO: add correct xyErr if DM-7101 is ever done. - d = {'x': sourceCat['x'].to_numpy(), 'y': sourceCat['y'].to_numpy(), - 'xCov': xCov.to_numpy(), 'yCov': yCov.to_numpy(), 'xyCov': xyCov.to_numpy()} + d = { + "x": sourceCat["x"].to_numpy(), + "y": sourceCat["y"].to_numpy(), + "xCov": xCov.to_numpy(), + "yCov": yCov.to_numpy(), + "xyCov": xyCov.to_numpy(), + } - wcsf.setObjects(extensionIndex, d, 'x', 'y', ['xCov', 'yCov', 'xyCov']) + wcsf.setObjects(extensionIndex, d, "x", "y", ["xCov", "yCov", "xyCov"]) def _add_ref_objects(self, wcsf, refObjects, refCovariance, extensionInfo): """Add reference sources to the wcsfit.WCSFit object. @@ -1023,14 +1103,23 @@ def _add_ref_objects(self, wcsf, refObjects, refCovariance, extensionInfo): extensionInfo : `lsst.pipe.base.Struct` Struct containing properties for each extension. """ - extensionIndex = np.flatnonzero(extensionInfo.extensionType == 'REFERENCE')[0] + extensionIndex = np.flatnonzero(extensionInfo.extensionType == "REFERENCE")[0] if self.config.fitProperMotion: - wcsf.setObjects(extensionIndex, refObjects, 'ra', 'dec', ['raCov', 'decCov', 'raDecCov'], - pmDecKey='decPM', pmRaKey='raPM', parallaxKey='parallax', pmCovKey='fullCov', - pmCov=refCovariance) + wcsf.setObjects( + extensionIndex, + refObjects, + "ra", + "dec", + ["raCov", "decCov", "raDecCov"], + pmDecKey="decPM", + pmRaKey="raPM", + parallaxKey="parallax", + pmCovKey="fullCov", + pmCov=refCovariance, + ) else: - wcsf.setObjects(extensionIndex, refObjects, 'ra', 'dec', ['raCov', 'decCov', 'raDecCov']) + wcsf.setObjects(extensionIndex, refObjects, "ra", "dec", ["raCov", "decCov", "raDecCov"]) def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, xScale=1, yScale=1): """Make an `lsst.afw.geom.SkyWcs` from a dictionary of mappings. @@ -1058,13 +1147,12 @@ def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, x WCS constructed from the input mappings """ # Set up pixel frames - pixelFrame = astshim.Frame(2, 'Domain=PIXELS') - normedPixelFrame = astshim.Frame(2, 'Domain=NORMEDPIXELS') + pixelFrame = astshim.Frame(2, "Domain=PIXELS") + normedPixelFrame = astshim.Frame(2, "Domain=NORMEDPIXELS") if doNormalizePixels: # Pixels will need to be rescaled before going into the mappings - normCoefficients = [-1.0, 2.0/xScale, 0, - -1.0, 0, 2.0/yScale] + normCoefficients = [-1.0, 2.0 / xScale, 0, -1.0, 0, 2.0 / yScale] normMap = _convert_to_ast_polymap_coefficients(normCoefficients) else: normMap = astshim.UnitMap(2) @@ -1073,34 +1161,34 @@ def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, x tangentPoint = lsst.geom.SpherePoint(centerRA, centerDec) cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0 * lsst.geom.degrees, True) iwcToSkyWcs = afwgeom.makeSkyWcs(lsst.geom.Point2D(0, 0), tangentPoint, cdMatrix) - iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping('PIXELS', 'SKY') - skyFrame = iwcToSkyWcs.getFrameDict().getFrame('SKY') + iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY") + skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY") frameDict = astshim.FrameDict(pixelFrame) - frameDict.addFrame('PIXELS', normMap, normedPixelFrame) + frameDict.addFrame("PIXELS", normMap, normedPixelFrame) - currentFrameName = 'NORMEDPIXELS' + currentFrameName = "NORMEDPIXELS" # Dictionary values are ordered according to the maps' application. for m, mapElement in enumerate(mapDict.values()): - mapType = mapElement['Type'] + mapType = mapElement["Type"] - if mapType == 'Poly': - mapCoefficients = mapElement['Coefficients'] + if mapType == "Poly": + mapCoefficients = mapElement["Coefficients"] astMap = _convert_to_ast_polymap_coefficients(mapCoefficients) - elif mapType == 'Identity': + elif mapType == "Identity": astMap = astshim.UnitMap(2) else: raise ValueError(f"Converting map type {mapType} to WCS is not supported") if m == len(mapDict) - 1: - newFrameName = 'IWC' + newFrameName = "IWC" else: - newFrameName = 'INTERMEDIATE' + str(m) - newFrame = astshim.Frame(2, f'Domain={newFrameName}') + newFrameName = "INTERMEDIATE" + str(m) + newFrame = astshim.Frame(2, f"Domain={newFrameName}") frameDict.addFrame(currentFrameName, astMap, newFrame) currentFrameName = newFrameName - frameDict.addFrame('IWC', iwcToSkyMap, skyFrame) + frameDict.addFrame("IWC", iwcToSkyMap, skyFrame) outWCS = afwgeom.SkyWcs(frameDict) return outWCS @@ -1129,40 +1217,44 @@ def _make_outputs(self, wcsf, visitSummaryTables, exposureInfo): # Set up the schema for the output catalogs schema = lsst.afw.table.ExposureTable.makeMinimalSchema() - schema.addField('visit', type='L', doc='Visit number') + schema.addField("visit", type="L", doc="Visit number") # Pixels will need to be rescaled before going into the mappings sampleDetector = visitSummaryTables[0][0] - xscale = sampleDetector['bbox_max_x'] - sampleDetector['bbox_min_x'] - yscale = sampleDetector['bbox_max_y'] - sampleDetector['bbox_min_y'] + xscale = sampleDetector["bbox_max_x"] - sampleDetector["bbox_min_x"] + yscale = sampleDetector["bbox_max_y"] - sampleDetector["bbox_min_y"] catalogs = {} for v, visitSummary in enumerate(visitSummaryTables): - visit = visitSummary[0]['visit'] + visit = visitSummary[0]["visit"] catalog = lsst.afw.table.ExposureCatalog(schema) catalog.resize(len(exposureInfo.detectors)) - catalog['visit'] = visit + catalog["visit"] = visit - for d, detector in enumerate(visitSummary['id']): - mapName = f'{visit}/{detector}' + for d, detector in enumerate(visitSummary["id"]): + mapName = f"{visit}/{detector}" - mapElements = wcsf.mapCollection.orderAtoms(f'{mapName}/base') + mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base") mapDict = {} for m, mapElement in enumerate(mapElements): mapType = wcsf.mapCollection.getMapType(mapElement) - mapDict[mapElement] = {'Type': mapType} + mapDict[mapElement] = {"Type": mapType} - if mapType == 'Poly': + if mapType == "Poly": mapCoefficients = mapParams[mapElement] - mapDict[mapElement]['Coefficients'] = mapCoefficients + mapDict[mapElement]["Coefficients"] = mapCoefficients # The RA and Dec of the visit are needed for the last step of # the mapping from the visit tangent plane to RA and Dec - outWCS = self._make_afw_wcs(mapDict, exposureInfo.ras[v] * lsst.geom.radians, - exposureInfo.decs[v] * lsst.geom.radians, - doNormalizePixels=True, - xScale=xscale, yScale=yscale) + outWCS = self._make_afw_wcs( + mapDict, + exposureInfo.ras[v] * lsst.geom.radians, + exposureInfo.decs[v] * lsst.geom.radians, + doNormalizePixels=True, + xScale=xscale, + yScale=yscale, + ) catalog[d].setId(detector) catalog[d].setWcs(outWCS) diff --git a/python/lsst/drp/tasks/update_visit_summary.py b/python/lsst/drp/tasks/update_visit_summary.py index 703120e8..29377923 100644 --- a/python/lsst/drp/tasks/update_visit_summary.py +++ b/python/lsst/drp/tasks/update_visit_summary.py @@ -194,9 +194,7 @@ def best_for_detector( if record is None: continue if center is None: - center_for_record = compute_center_for_detector_record( - record, bbox=bbox - ) + center_for_record = compute_center_for_detector_record(record, bbox=bbox) if center_for_record is None: continue else: @@ -357,9 +355,7 @@ def __init__(self, *, config: UpdateVisitSummaryConfig | None = None): case "global": self.inputs.remove("wcs_overrides_tract") case bad: - raise ValueError( - f"Invalid value wcs_provider={bad!r}; config was not validated." - ) + raise ValueError(f"Invalid value wcs_provider={bad!r}; config was not validated.") match self.config.photo_calib_provider: case "input_summary": self.inputs.remove("photo_calib_overrides_tract") @@ -369,9 +365,7 @@ def __init__(self, *, config: UpdateVisitSummaryConfig | None = None): case "global": self.inputs.remove("photo_calib_overrides_tract") case bad: - raise ValueError( - f"Invalid value photo_calib_provider={bad!r}; config was not validated." - ) + raise ValueError(f"Invalid value photo_calib_provider={bad!r}; config was not validated.") match self.config.background_provider: case "input_summary": self.inputs.remove("background_originals") @@ -379,14 +373,10 @@ def __init__(self, *, config: UpdateVisitSummaryConfig | None = None): case "replacement": pass case bad: - raise ValueError( - f"Invalid value background_provider={bad!r}; config was not validated." - ) + raise ValueError(f"Invalid value background_provider={bad!r}; config was not validated.") -class UpdateVisitSummaryConfig( - PipelineTaskConfig, pipelineConnections=UpdateVisitSummaryConnections -): +class UpdateVisitSummaryConfig(PipelineTaskConfig, pipelineConnections=UpdateVisitSummaryConnections): """Configuration for UpdateVisitSummaryTask. Notes @@ -520,9 +510,7 @@ def __init__(self, *, initInputs: dict[str, Any] | None = None, **kwargs: Any): self.schema_mapper.addMinimalSchema(input_summary_schema) self.schema = self.schema_mapper.getOutputSchema() if self.config.wcs_provider == "tract": - self.schema.addField( - "wcsTractId", type="L", doc="ID of the tract that provided the WCS." - ) + self.schema.addField("wcsTractId", type="L", doc="ID of the tract that provided the WCS.") if self.config.photo_calib_provider == "tract": self.schema.addField( "photoCalibTractId", @@ -546,14 +534,10 @@ def runQuantum( # objects). match self.config.wcs_provider: case "tract": - inputs["wcs_overrides"] = PerTractInput.load( - butlerQC, sky_map, inputRefs.wcs_overrides_tract - ) + inputs["wcs_overrides"] = PerTractInput.load(butlerQC, sky_map, inputRefs.wcs_overrides_tract) del inputRefs.wcs_overrides_tract case "global": - inputs["wcs_overrides"] = GlobalInput( - butlerQC.get(inputRefs.wcs_overrides_global) - ) + inputs["wcs_overrides"] = GlobalInput(butlerQC.get(inputRefs.wcs_overrides_global)) del inputRefs.wcs_overrides_global case "input_summary": inputs["wcs_overrides"] = None @@ -582,9 +566,7 @@ def runQuantum( # deferLoad=True connections into mappings keyed by detector ID. for name in deferred_dataset_types: handles_list = inputs[name] - inputs[name] = { - handle.dataId["detector"]: handle for handle in handles_list - } + inputs[name] = {handle.dataId["detector"]: handle for handle in handles_list} for record in inputs["input_summary_catalog"]: detector_id = record.getId() if detector_id not in inputs[name]: @@ -703,9 +685,7 @@ def run( bbox = exposure.getBBox() if wcs_overrides: - wcs_tract, wcs_record = wcs_overrides.best_for_detector( - detector_id, bbox=bbox - ) + wcs_tract, wcs_record = wcs_overrides.best_for_detector(detector_id, bbox=bbox) if wcs_record is not None: wcs = wcs_record.getWcs() else: @@ -755,9 +735,7 @@ def run( if self.config.photo_calib_provider == "tract": output_record["photoCalibTractId"] = photo_calib_tract output_record.setPhotoCalib(photo_calib) - self.compute_summary_stats.update_photo_calib_stats( - summary_stats, photo_calib - ) + self.compute_summary_stats.update_photo_calib_stats(summary_stats, photo_calib) if background_overrides is not None: if (handle := background_overrides.get(detector_id)) is not None: @@ -771,9 +749,7 @@ def run( for layer in new_bkg: full_bkg.append(layer) exposure.image -= new_bkg.getImage() - self.compute_summary_stats.update_background_stats( - summary_stats, full_bkg - ) + self.compute_summary_stats.update_background_stats(summary_stats, full_bkg) self.compute_summary_stats.update_masked_image_stats( summary_stats, exposure.getMaskedImage() ) diff --git a/tests/assemble_coadd_test_utils.py b/tests/assemble_coadd_test_utils.py index 087a7543..fd885adc 100644 --- a/tests/assemble_coadd_test_utils.py +++ b/tests/assemble_coadd_test_utils.py @@ -60,8 +60,10 @@ def makeMockSkyInfo(bbox, wcs, patch): skyInfo : `lsst.pipe.base.Struct` Patch geometry information. """ + def getIndex(): return patch + patchInfo = pipeBase.Struct(getIndex=getIndex) skyInfo = pipeBase.Struct(bbox=bbox, wcs=wcs, patchInfo=patchInfo) return skyInfo @@ -123,7 +125,8 @@ class MockCoaddTestData: `~lsst.meas.algorithms.testUtils.plantSources` lacking the option to specify the pixel origin. """ - rotAngle = 0.*degrees + + rotAngle = 0.0 * degrees """Rotation of the pixel grid on the sky, East from North (`lsst.geom.Angle`). """ @@ -161,12 +164,25 @@ class MockCoaddTestData: detector = None "Properties of the CCD for the exposure (`lsst.afw.cameraGeom.Detector`)." - def __init__(self, shape=geom.Extent2I(201, 301), offset=geom.Point2I(-123, -45), - backgroundLevel=314.592, seed=42, nSrc=37, - fluxRange=2., noiseLevel=5, sourceSigma=200., - minPsfSize=1.5, maxPsfSize=3., - pixelScale=0.2*arcseconds, ra=209.*degrees, dec=-20.25*degrees, - ccd=37, patch=42, tract=0): + def __init__( + self, + shape=geom.Extent2I(201, 301), + offset=geom.Point2I(-123, -45), + backgroundLevel=314.592, + seed=42, + nSrc=37, + fluxRange=2.0, + noiseLevel=5, + sourceSigma=200.0, + minPsfSize=1.5, + maxPsfSize=3.0, + pixelScale=0.2 * arcseconds, + ra=209.0 * degrees, + dec=-20.25 * degrees, + ccd=37, + patch=42, + tract=0, + ): self.ra = ra self.dec = dec self.pixelScale = pixelScale @@ -182,15 +198,15 @@ def __init__(self, shape=geom.Extent2I(201, 301), offset=geom.Point2I(-123, -45) # Set up properties of the simulations nSigmaForKernel = 5 - self.kernelSize = (int(maxPsfSize*nSigmaForKernel + 0.5)//2)*2 + 1 # make sure it is odd + self.kernelSize = (int(maxPsfSize * nSigmaForKernel + 0.5) // 2) * 2 + 1 # make sure it is odd - bufferSize = self.kernelSize//2 + bufferSize = self.kernelSize // 2 x0, y0 = self.bbox.getBegin() xSize, ySize = self.bbox.getDimensions() # Set the pixel coordinates and fluxes of the simulated sources. - self.xLoc = self.rngData.random(nSrc)*(xSize - 2*bufferSize) + bufferSize + x0 - self.yLoc = self.rngData.random(nSrc)*(ySize - 2*bufferSize) + bufferSize + y0 - self.flux = (self.rngData.random(nSrc)*(fluxRange - 1.) + 1.)*sourceSigma*noiseLevel + self.xLoc = self.rngData.random(nSrc) * (xSize - 2 * bufferSize) + bufferSize + x0 + self.yLoc = self.rngData.random(nSrc) * (ySize - 2 * bufferSize) + bufferSize + y0 + self.flux = (self.rngData.random(nSrc) * (fluxRange - 1.0) + 1.0) * sourceSigma * noiseLevel self.backgroundLevel = backgroundLevel self.noiseLevel = noiseLevel @@ -290,42 +306,54 @@ def makeDummyVisitInfo(self, exposureId, randomizeTime=False): visitInfo : `lsst.afw.image.VisitInfo` VisitInfo for the exposure. """ - lsstLat = -30.244639*u.degree - lsstLon = -70.749417*u.degree - lsstAlt = 2663.*u.m - lsstTemperature = 20.*u.Celsius - lsstHumidity = 40. # in percent - lsstPressure = 73892.*u.pascal - loc = EarthLocation(lat=lsstLat, - lon=lsstLon, - height=lsstAlt) + lsstLat = -30.244639 * u.degree + lsstLon = -70.749417 * u.degree + lsstAlt = 2663.0 * u.m + lsstTemperature = 20.0 * u.Celsius + lsstHumidity = 40.0 # in percent + lsstPressure = 73892.0 * u.pascal + loc = EarthLocation(lat=lsstLat, lon=lsstLon, height=lsstAlt) time = Time(2000.0, format="jyear", scale="tt") if randomizeTime: # Pick a random time within a 6 hour window - time += 6*u.hour*(self.rngMods.random() - 0.5) - radec = SkyCoord(dec=self.dec.asDegrees(), ra=self.ra.asDegrees(), - unit='deg', obstime=time, frame='icrs', location=loc) - airmass = float(1.0/np.sin(radec.altaz.alt)) - obsInfo = makeObservationInfo(location=loc, - detector_exposure_id=exposureId, - datetime_begin=time, - datetime_end=time, - boresight_airmass=airmass, - boresight_rotation_angle=Angle(0.*u.degree), - boresight_rotation_coord='sky', - temperature=lsstTemperature, - pressure=lsstPressure, - relative_humidity=lsstHumidity, - tracking_radec=radec, - altaz_begin=radec.altaz, - observation_type='science', - ) + time += 6 * u.hour * (self.rngMods.random() - 0.5) + radec = SkyCoord( + dec=self.dec.asDegrees(), + ra=self.ra.asDegrees(), + unit="deg", + obstime=time, + frame="icrs", + location=loc, + ) + airmass = float(1.0 / np.sin(radec.altaz.alt)) + obsInfo = makeObservationInfo( + location=loc, + detector_exposure_id=exposureId, + datetime_begin=time, + datetime_end=time, + boresight_airmass=airmass, + boresight_rotation_angle=Angle(0.0 * u.degree), + boresight_rotation_coord="sky", + temperature=lsstTemperature, + pressure=lsstPressure, + relative_humidity=lsstHumidity, + tracking_radec=radec, + altaz_begin=radec.altaz, + observation_type="science", + ) visitInfo = MakeRawVisitInfoViaObsInfo.observationInfo2visitInfo(obsInfo) return visitInfo - def makeTestImage(self, expId, noiseLevel=None, psfSize=None, backgroundLevel=None, - detectionSigma=5., badRegionBox=None): + def makeTestImage( + self, + expId, + noiseLevel=None, + psfSize=None, + backgroundLevel=None, + detectionSigma=5.0, + badRegionBox=None, + ): """Make a reproduceable PSF-convolved masked image for testing. Parameters @@ -346,29 +374,31 @@ def makeTestImage(self, expId, noiseLevel=None, psfSize=None, backgroundLevel=No if backgroundLevel is None: backgroundLevel = self.backgroundLevel if noiseLevel is None: - noiseLevel = 5. + noiseLevel = 5.0 visitInfo = self.makeDummyVisitInfo(expId, randomizeTime=True) if psfSize is None: - psfSize = self.rngMods.random()*(self.maxPsfSize - self.minPsfSize) + self.minPsfSize + psfSize = self.rngMods.random() * (self.maxPsfSize - self.minPsfSize) + self.minPsfSize nSrc = len(self.flux) sigmas = [psfSize for src in range(nSrc)] sigmasPsfMatched = [self.maxPsfSize for src in range(nSrc)] coordList = list(zip(self.xLoc, self.yLoc, self.flux, sigmas)) coordListPsfMatched = list(zip(self.xLoc, self.yLoc, self.flux, sigmasPsfMatched)) xSize, ySize = self.bbox.getDimensions() - model = plantSources(self.bbox, self.kernelSize, self.backgroundLevel, - coordList, addPoissonNoise=False) - modelPsfMatched = plantSources(self.bbox, self.kernelSize, self.backgroundLevel, - coordListPsfMatched, addPoissonNoise=False) + model = plantSources( + self.bbox, self.kernelSize, self.backgroundLevel, coordList, addPoissonNoise=False + ) + modelPsfMatched = plantSources( + self.bbox, self.kernelSize, self.backgroundLevel, coordListPsfMatched, addPoissonNoise=False + ) model.variance.array = np.abs(model.image.array) + noiseLevel modelPsfMatched.variance.array = np.abs(modelPsfMatched.image.array) + noiseLevel - noise = self.rngData.random((ySize, xSize))*noiseLevel + noise = self.rngData.random((ySize, xSize)) * noiseLevel noise -= np.median(noise) model.image.array += noise modelPsfMatched.image.array += noise detectedMask = afwImage.Mask.getPlaneBitMask("DETECTED") - detectionThreshold = self.backgroundLevel + detectionSigma*noiseLevel + detectionThreshold = self.backgroundLevel + detectionSigma * noiseLevel model.mask.array[model.image.array > detectionThreshold] += detectedMask if badRegionBox is not None: @@ -406,9 +436,9 @@ def makeDataRefList(exposures, matchedExposures, warpType, tract=0, patch=42, co """ dataRefList = [] for expId in exposures: - if warpType == 'direct': + if warpType == "direct": exposure = exposures[expId] - elif warpType == 'psfMatched': + elif warpType == "psfMatched": exposure = matchedExposures[expId] else: raise ValueError("warpType must be one of 'direct' or 'psfMatched'") @@ -419,7 +449,7 @@ def makeDataRefList(exposures, matchedExposures, warpType, tract=0, patch=42, co tract=tract, patch=patch, visit=expId, - coaddName=coaddName + coaddName=coaddName, ) dataRefList.append(dataRef) return dataRefList diff --git a/tests/test_assemble_coadd.py b/tests/test_assemble_coadd.py index e0304a19..f0ee8d67 100644 --- a/tests/test_assemble_coadd.py +++ b/tests/test_assemble_coadd.py @@ -33,12 +33,15 @@ from lsst.drp.tasks.dcr_assemble_coadd import DcrAssembleCoaddTask, DcrAssembleCoaddConfig from assemble_coadd_test_utils import makeMockSkyInfo, MockCoaddTestData -__all__ = ["MockAssembleCoaddConfig", "MockAssembleCoaddTask", - "MockCompareWarpAssembleCoaddConfig", "MockCompareWarpAssembleCoaddTask"] +__all__ = [ + "MockAssembleCoaddConfig", + "MockAssembleCoaddTask", + "MockCompareWarpAssembleCoaddConfig", + "MockCompareWarpAssembleCoaddTask", +] class MockAssembleCoaddConfig(AssembleCoaddConfig): - def setDefaults(self): super().setDefaults() self.doWrite = False @@ -51,6 +54,7 @@ class MockAssembleCoaddTask(AssembleCoaddTask): up the Task, and instead supply in-memory mock data references to the `run` method so that the coaddition algorithms can be tested without a Butler. """ + ConfigClass = MockAssembleCoaddConfig def __init__(self, **kwargs): @@ -84,13 +88,17 @@ def runQuantum(self, mockSkyInfo, warpRefList, *args): """ inputs = self.prepareInputs(warpRefList) - retStruct = self.run(mockSkyInfo, inputs.tempExpRefList, inputs.imageScalerList, - inputs.weightList, supplementaryData=pipeBase.Struct()) + retStruct = self.run( + mockSkyInfo, + inputs.tempExpRefList, + inputs.imageScalerList, + inputs.weightList, + supplementaryData=pipeBase.Struct(), + ) return retStruct class MockCompareWarpAssembleCoaddConfig(CompareWarpAssembleCoaddConfig): - def setDefaults(self): super().setDefaults() self.assembleStaticSkyModel.retarget(MockAssembleCoaddTask) @@ -106,6 +114,7 @@ class MockCompareWarpAssembleCoaddTask(MockAssembleCoaddTask, CompareWarpAssembl up the Task, and instead supply in-memory mock data references to the `run` method so that the coaddition algorithms can be tested without a Butler. """ + ConfigClass = MockCompareWarpAssembleCoaddConfig _DefaultName = "compareWarpAssembleCoadd" @@ -123,22 +132,27 @@ def runQuantum(self, mockSkyInfo, warpRefList, *args): nImage=templateCoadd.nImage, warpRefList=templateCoadd.warpRefList, imageScalerList=templateCoadd.imageScalerList, - weightList=templateCoadd.weightList) - - retStruct = self.run(mockSkyInfo, inputs.tempExpRefList, inputs.imageScalerList, - inputs.weightList, supplementaryData=supplementaryData) + weightList=templateCoadd.weightList, + ) + + retStruct = self.run( + mockSkyInfo, + inputs.tempExpRefList, + inputs.imageScalerList, + inputs.weightList, + supplementaryData=supplementaryData, + ) return retStruct class MockDcrAssembleCoaddConfig(DcrAssembleCoaddConfig): - def setDefaults(self): super().setDefaults() self.assembleStaticSkyModel.retarget(MockCompareWarpAssembleCoaddTask) self.assembleStaticSkyModel.doWrite = False self.doWrite = False self.effectiveWavelength = 476.31 # Use LSST g band values for the test. - self.bandwidth = 552. - 405. + self.bandwidth = 552.0 - 405.0 class MockDcrAssembleCoaddTask(MockCompareWarpAssembleCoaddTask, DcrAssembleCoaddTask): @@ -149,6 +163,7 @@ class MockDcrAssembleCoaddTask(MockCompareWarpAssembleCoaddTask, DcrAssembleCoad up the Task, and instead supply in-memory mock data references to the `run` method so that the coaddition algorithms can be tested without a Butler. """ + ConfigClass = MockDcrAssembleCoaddConfig _DefaultName = "dcrAssembleCoadd" @@ -157,7 +172,6 @@ def __init__(self, *args, **kwargs): class MockInputMapAssembleCoaddConfig(MockCompareWarpAssembleCoaddConfig): - def setDefaults(self): super().setDefaults() self.doInputMap = True @@ -171,6 +185,7 @@ class MockInputMapAssembleCoaddTask(MockCompareWarpAssembleCoaddTask): up the Task, and instead supply in-memory mock data references to the `run` method so that the coaddition algorithms can be tested without a Butler. """ + ConfigClass = MockInputMapAssembleCoaddConfig _DefaultName = "inputMapAssembleCoadd" @@ -193,10 +208,12 @@ def setUp(self): matchedExposures = {} for expId in range(100, 110): exposures[expId], matchedExposures[expId] = testData.makeTestImage(expId) - self.dataRefList = testData.makeDataRefList(exposures, matchedExposures, - 'direct', patch=patch, tract=tract) - self.dataRefListPsfMatched = testData.makeDataRefList(exposures, matchedExposures, - 'psfMatched', patch=patch, tract=tract) + self.dataRefList = testData.makeDataRefList( + exposures, matchedExposures, "direct", patch=patch, tract=tract + ) + self.dataRefListPsfMatched = testData.makeDataRefList( + exposures, matchedExposures, "psfMatched", patch=patch, tract=tract + ) self.skyInfo = makeMockSkyInfo(testData.bbox, testData.wcs, patch=patch) def checkRun(self, assembleTask, warpType="direct"): @@ -245,10 +262,8 @@ def testOnlineCoadd(self): resultsOnline = assembleTaskOnline.runQuantum(self.skyInfo, dataRefList) coaddOnline = resultsOnline.coaddExposure - self.assertFloatsAlmostEqual(coaddOnline.image.array, - coadd.image.array, rtol=1e-3) - self.assertFloatsAlmostEqual(coaddOnline.variance.array, - coadd.variance.array, rtol=1e-6) + self.assertFloatsAlmostEqual(coaddOnline.image.array, coadd.image.array, rtol=1e-3) + self.assertFloatsAlmostEqual(coaddOnline.variance.array, coadd.variance.array, rtol=1e-6) self.assertMasksEqual(coaddOnline.mask, coadd.mask) def testInputMap(self): @@ -263,15 +278,16 @@ def testInputMap(self): matchedExposures = {} for expId in range(100, 110): if expId == 105: - badBox = lsst.geom.Box2I(lsst.geom.Point2I(testData.bbox.beginX + 10, - testData.bbox.beginY + 10), - lsst.geom.Extent2I(100, 100)) + badBox = lsst.geom.Box2I( + lsst.geom.Point2I(testData.bbox.beginX + 10, testData.bbox.beginY + 10), + lsst.geom.Extent2I(100, 100), + ) else: badBox = None - exposures[expId], matchedExposures[expId] = testData.makeTestImage(expId, - badRegionBox=badBox) - dataRefList = testData.makeDataRefList(exposures, matchedExposures, - 'direct', patch=patch, tract=tract) + exposures[expId], matchedExposures[expId] = testData.makeTestImage(expId, badRegionBox=badBox) + dataRefList = testData.makeDataRefList( + exposures, matchedExposures, "direct", patch=patch, tract=tract + ) results = assembleTask.runQuantum(self.skyInfo, dataRefList) @@ -292,8 +308,8 @@ def testInputMap(self): metadata = inputMap.metadata visitBitDict = {} for bit in range(inputMap.wide_mask_maxbits): - if f'B{bit:04d}VIS' in metadata: - visitBitDict[metadata[f'B{bit:04d}VIS']] = bit + if f"B{bit:04d}VIS" in metadata: + visitBitDict[metadata[f"B{bit:04d}VIS"]] = bit for expId in range(100, 110): if expId == 105: self.assertFalse(np.all(inputMap.check_bits_pix(validPix, [visitBitDict[expId]]))) diff --git a/tests/test_dcr_assemble_coadd.py b/tests/test_dcr_assemble_coadd.py index 31876f4c..916d329c 100644 --- a/tests/test_dcr_assemble_coadd.py +++ b/tests/test_dcr_assemble_coadd.py @@ -28,6 +28,7 @@ class DcrAssembleCoaddCalculateGainTestCase(lsst.utils.tests.TestCase): """Tests of dcrAssembleCoaddTask.calculateGain().""" + def setUp(self): self.baseGain = 0.5 self.gainList = [self.baseGain, self.baseGain] @@ -36,7 +37,7 @@ def setUp(self): # perfectly, so that the improvement is limited only by our # conservative gain. for i in range(2): - self.convergenceList.append(self.convergenceList[i]/(self.baseGain + 1)) + self.convergenceList.append(self.convergenceList[i] / (self.baseGain + 1)) self.nextGain = (1 + self.baseGain) / 2 self.config = DcrAssembleCoaddConfig() @@ -107,7 +108,7 @@ def testProgressiveGainBadFit(self): gainList = [self.baseGain, self.baseGain] convergenceList = [0.2] for i in range(2): - convergenceList.append(convergenceList[i]/(wrongGain + 1)) + convergenceList.append(convergenceList[i] / (wrongGain + 1)) # The below math is a simplified version of the full algorithm, # assuming the predicted convergence is zero. # Note that in this case, nextGain is smaller than wrongGain. diff --git a/tests/test_gbdesAstrometricFit.py b/tests/test_gbdesAstrometricFit.py index 6f8d99f2..ddc0beae 100644 --- a/tests/test_gbdesAstrometricFit.py +++ b/tests/test_gbdesAstrometricFit.py @@ -42,10 +42,8 @@ class TestGbdesAstrometricFit(lsst.utils.tests.TestCase): - @classmethod def setUpClass(cls): - # Set random seed np.random.seed(1234) @@ -58,7 +56,7 @@ def setUpClass(cls): cls.datadir = os.path.join(TESTDIR, "data") cls.fieldNumber = 0 - cls.instrumentName = 'HSC' + cls.instrumentName = "HSC" cls.instrument = wcsfit.Instrument(cls.instrumentName) cls.refEpoch = 57205.5 @@ -67,8 +65,9 @@ def setUpClass(cls): cls.testVisits = [1176, 17900, 17930, 17934] cls.inputVisitSummary = [] for testVisit in cls.testVisits: - visSum = afwTable.ExposureCatalog.readFits(os.path.join(cls.datadir, - f'visitSummary_{testVisit}.fits')) + visSum = afwTable.ExposureCatalog.readFits( + os.path.join(cls.datadir, f"visitSummary_{testVisit}.fits") + ) cls.inputVisitSummary.append(visSum) cls.config = GbdesAstrometricFitConfig() @@ -80,28 +79,32 @@ def setUpClass(cls): cls.task = GbdesAstrometricFitTask(config=cls.config) cls.exposureInfo, cls.exposuresHelper, cls.extensionInfo = cls.task._get_exposure_info( - cls.inputVisitSummary, cls.instrument, refEpoch=cls.refEpoch) + cls.inputVisitSummary, cls.instrument, refEpoch=cls.refEpoch + ) cls.fields, cls.fieldCenter, cls.fieldRadius = cls.task._prep_sky( - cls.inputVisitSummary, cls.exposureInfo.medianEpoch) + cls.inputVisitSummary, cls.exposureInfo.medianEpoch + ) # Bounding box of observations: raMins, raMaxs = [], [] decMins, decMaxs = [], [] for visSum in cls.inputVisitSummary: - raMins.append(visSum['raCorners'].min()) - raMaxs.append(visSum['raCorners'].max()) - decMins.append(visSum['decCorners'].min()) - decMaxs.append(visSum['decCorners'].max()) + raMins.append(visSum["raCorners"].min()) + raMaxs.append(visSum["raCorners"].max()) + decMins.append(visSum["decCorners"].min()) + decMaxs.append(visSum["decCorners"].max()) raMin = min(raMins) raMax = max(raMaxs) decMin = min(decMins) decMax = max(decMaxs) - corners = [lsst.geom.SpherePoint(raMin, decMin, lsst.geom.degrees).getVector(), - lsst.geom.SpherePoint(raMax, decMin, lsst.geom.degrees).getVector(), - lsst.geom.SpherePoint(raMax, decMax, lsst.geom.degrees).getVector(), - lsst.geom.SpherePoint(raMin, decMax, lsst.geom.degrees).getVector()] + corners = [ + lsst.geom.SpherePoint(raMin, decMin, lsst.geom.degrees).getVector(), + lsst.geom.SpherePoint(raMax, decMin, lsst.geom.degrees).getVector(), + lsst.geom.SpherePoint(raMax, decMax, lsst.geom.degrees).getVector(), + lsst.geom.SpherePoint(raMin, decMax, lsst.geom.degrees).getVector(), + ] cls.boundingPolygon = sphgeom.ConvexPolygon(corners) # Make random set of data in a bounding box determined by input visits @@ -115,23 +118,26 @@ def setUpClass(cls): refDataId, deferredRefCat = cls._make_refCat(starIds, starRAs, starDecs, inReferenceFraction) cls.refObjectLoader = ReferenceObjectLoader([refDataId], [deferredRefCat]) cls.refObjectLoader.config.requireProperMotion = False - cls.refObjectLoader.config.anyFilterMapsToThis = 'test_filter' + cls.refObjectLoader.config.anyFilterMapsToThis = "test_filter" cls.task.refObjectLoader = cls.refObjectLoader # Get True WCS for stars: - with open(os.path.join(cls.datadir, 'sample_wcs.yaml'), 'r') as f: + with open(os.path.join(cls.datadir, "sample_wcs.yaml"), "r") as f: cls.trueModel = yaml.load(f, Loader=yaml.Loader) trueWCSs = cls._make_wcs(cls.trueModel, cls.inputVisitSummary) # Make source catalogs: - cls.inputCatalogRefs = cls._make_sourceCat(starIds, starRAs, starDecs, trueWCSs, - inScienceFraction) + cls.inputCatalogRefs = cls._make_sourceCat(starIds, starRAs, starDecs, trueWCSs, inScienceFraction) - cls.outputs = cls.task.run(cls.inputCatalogRefs, cls.inputVisitSummary, - instrumentName=cls.instrumentName, refEpoch=cls.refEpoch, - refObjectLoader=cls.refObjectLoader) + cls.outputs = cls.task.run( + cls.inputCatalogRefs, + cls.inputVisitSummary, + instrumentName=cls.instrumentName, + refEpoch=cls.refEpoch, + refObjectLoader=cls.refObjectLoader, + ) @classmethod def _make_refCat(cls, starIds, starRas, starDecs, inReferenceFraction): @@ -161,7 +167,7 @@ def _make_refCat(cls, starIds, starRas, starDecs, inReferenceFraction): # determined by bounding box used in above simulate. refSchema = afwTable.SimpleTable.makeMinimalSchema() idKey = refSchema.addField("sourceId", type="I") - fluxKey = refSchema.addField("test_filter_flux", units='nJy', type=np.float64) + fluxKey = refSchema.addField("test_filter_flux", units="nJy", type=np.float64) raErrKey = refSchema.addField("coord_raErr", type=np.float64) decErrKey = refSchema.addField("coord_decErr", type=np.float64) pmraErrKey = refSchema.addField("pm_raErr", type=np.float64) @@ -211,10 +217,14 @@ def _make_sourceCat(cls, starIds, starRas, starDecs, trueWCSs, inScienceFraction inputCatalogRefs = [] # Take a subset of the simulated data # Use true wcs objects to put simulated data into ccds - bbox = lsst.geom.BoxD(lsst.geom.Point2D(cls.inputVisitSummary[0][0]['bbox_min_x'], - cls.inputVisitSummary[0][0]['bbox_min_y']), - lsst.geom.Point2D(cls.inputVisitSummary[0][0]['bbox_max_x'], - cls.inputVisitSummary[0][0]['bbox_max_y'])) + bbox = lsst.geom.BoxD( + lsst.geom.Point2D( + cls.inputVisitSummary[0][0]["bbox_min_x"], cls.inputVisitSummary[0][0]["bbox_min_y"] + ), + lsst.geom.Point2D( + cls.inputVisitSummary[0][0]["bbox_max_x"], cls.inputVisitSummary[0][0]["bbox_max_y"] + ), + ) bboxCorners = bbox.getCorners() cls.inputCatalogRefs = [] for v, visit in enumerate(cls.testVisits): @@ -226,51 +236,64 @@ def _make_sourceCat(cls, starIds, starRas, starDecs, trueWCSs, inScienceFraction sourceCats = [] for detector in trueWCSs[visit]: detWcs = detector.getWcs() - detectorId = detector['id'] + detectorId = detector["id"] radecCorners = detWcs.pixelToSky(bboxCorners) detectorFootprint = sphgeom.ConvexPolygon([rd.getVector() for rd in radecCorners]) - detectorIndices = detectorFootprint.contains((visitStarRas*u.degree).to(u.radian), - (visitStarDecs*u.degree).to(u.radian)) + detectorIndices = detectorFootprint.contains( + (visitStarRas * u.degree).to(u.radian), (visitStarDecs * u.degree).to(u.radian) + ) nDetectorStars = detectorIndices.sum() - detectorArray = np.ones(nDetectorStars, dtype=bool) * detector['id'] + detectorArray = np.ones(nDetectorStars, dtype=bool) * detector["id"] ones_like = np.ones(nDetectorStars) zeros_like = np.zeros(nDetectorStars, dtype=bool) - x, y = detWcs.skyToPixelArray(visitStarRas[detectorIndices], visitStarDecs[detectorIndices], - degrees=True) + x, y = detWcs.skyToPixelArray( + visitStarRas[detectorIndices], visitStarDecs[detectorIndices], degrees=True + ) - origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]['id'] == detectorId])[0].getWcs() + origWcs = (cls.inputVisitSummary[v][cls.inputVisitSummary[v]["id"] == detectorId])[0].getWcs() inputRa, inputDec = origWcs.pixelToSkyArray(x, y, degrees=True) sourceDict = {} - sourceDict['detector'] = detectorArray - sourceDict['sourceId'] = visitStarIds[detectorIndices] - sourceDict['x'] = x - sourceDict['y'] = y - sourceDict['xErr'] = 1e-3 * ones_like - sourceDict['yErr'] = 1e-3 * ones_like - sourceDict['inputRA'] = inputRa - sourceDict['inputDec'] = inputDec - sourceDict['trueRA'] = visitStarRas[detectorIndices] - sourceDict['trueDec'] = visitStarDecs[detectorIndices] - for key in ['apFlux_12_0_flux', 'apFlux_12_0_instFlux', 'ixx', 'iyy']: + sourceDict["detector"] = detectorArray + sourceDict["sourceId"] = visitStarIds[detectorIndices] + sourceDict["x"] = x + sourceDict["y"] = y + sourceDict["xErr"] = 1e-3 * ones_like + sourceDict["yErr"] = 1e-3 * ones_like + sourceDict["inputRA"] = inputRa + sourceDict["inputDec"] = inputDec + sourceDict["trueRA"] = visitStarRas[detectorIndices] + sourceDict["trueDec"] = visitStarDecs[detectorIndices] + for key in ["apFlux_12_0_flux", "apFlux_12_0_instFlux", "ixx", "iyy"]: sourceDict[key] = ones_like - for key in ['pixelFlags_edge', 'pixelFlags_saturated', 'pixelFlags_interpolatedCenter', - 'pixelFlags_interpolated', 'pixelFlags_crCenter', 'pixelFlags_bad', - 'hsmPsfMoments_flag', 'apFlux_12_0_flag', 'extendedness', 'parentSourceId', - 'deblend_nChild', 'ixy']: + for key in [ + "pixelFlags_edge", + "pixelFlags_saturated", + "pixelFlags_interpolatedCenter", + "pixelFlags_interpolated", + "pixelFlags_crCenter", + "pixelFlags_bad", + "hsmPsfMoments_flag", + "apFlux_12_0_flag", + "extendedness", + "parentSourceId", + "deblend_nChild", + "ixy", + ]: sourceDict[key] = zeros_like - sourceDict['apFlux_12_0_instFluxErr'] = 1e-3 * ones_like - sourceDict['detect_isPrimary'] = ones_like.astype(bool) + sourceDict["apFlux_12_0_instFluxErr"] = 1e-3 * ones_like + sourceDict["detect_isPrimary"] = ones_like.astype(bool) sourceCat = pd.DataFrame(sourceDict) sourceCats.append(sourceCat) visitSourceTable = pd.concat(sourceCats) - inputCatalogRef = InMemoryDatasetHandle(visitSourceTable, storageClass="DataFrame", - dataId={"visit": visit}) + inputCatalogRef = InMemoryDatasetHandle( + visitSourceTable, storageClass="DataFrame", dataId={"visit": visit} + ) inputCatalogRefs.append(inputCatalogRef) @@ -293,45 +316,50 @@ def _make_wcs(cls, model, inputVisitSummaries): """ # Pixels will need to be rescaled before going into the mappings - xscale = inputVisitSummaries[0][0]['bbox_max_x'] - inputVisitSummaries[0][0]['bbox_min_x'] - yscale = inputVisitSummaries[0][0]['bbox_max_y'] - inputVisitSummaries[0][0]['bbox_min_y'] + xscale = inputVisitSummaries[0][0]["bbox_max_x"] - inputVisitSummaries[0][0]["bbox_min_x"] + yscale = inputVisitSummaries[0][0]["bbox_max_y"] - inputVisitSummaries[0][0]["bbox_min_y"] catalogs = {} schema = lsst.afw.table.ExposureTable.makeMinimalSchema() - schema.addField('visit', type='L', doc='Visit number') + schema.addField("visit", type="L", doc="Visit number") for visitSum in inputVisitSummaries: - visit = visitSum[0]['visit'] - visitMapName = f'{visit}/poly' + visit = visitSum[0]["visit"] + visitMapName = f"{visit}/poly" visitModel = model[visitMapName] catalog = lsst.afw.table.ExposureCatalog(schema) catalog.resize(len(visitSum)) - catalog['visit'] = visit + catalog["visit"] = visit raDec = visitSum[0].getVisitInfo().getBoresightRaDec() - visitMapType = visitModel['Type'] - visitDict = {'Type': visitMapType} - if visitMapType == 'Poly': - mapCoefficients = (visitModel['XPoly']['Coefficients'] - + visitModel['YPoly']['Coefficients']) + visitMapType = visitModel["Type"] + visitDict = {"Type": visitMapType} + if visitMapType == "Poly": + mapCoefficients = visitModel["XPoly"]["Coefficients"] + visitModel["YPoly"]["Coefficients"] visitDict["Coefficients"] = mapCoefficients for d, detector in enumerate(visitSum): - detectorId = detector['id'] - detectorMapName = f'HSC/{detectorId}/poly' + detectorId = detector["id"] + detectorMapName = f"HSC/{detectorId}/poly" detectorModel = model[detectorMapName] - detectorMapType = detectorModel['Type'] - mapDict = {detectorMapName: {'Type': detectorMapType}, - visitMapName: visitDict} - if detectorMapType == 'Poly': - mapCoefficients = (detectorModel['XPoly']['Coefficients'] - + detectorModel['YPoly']['Coefficients']) - mapDict[detectorMapName]['Coefficients'] = mapCoefficients - - outWCS = cls.task._make_afw_wcs(mapDict, raDec.getRa(), raDec.getDec(), - doNormalizePixels=True, xScale=xscale, yScale=yscale) + detectorMapType = detectorModel["Type"] + mapDict = {detectorMapName: {"Type": detectorMapType}, visitMapName: visitDict} + if detectorMapType == "Poly": + mapCoefficients = ( + detectorModel["XPoly"]["Coefficients"] + detectorModel["YPoly"]["Coefficients"] + ) + mapDict[detectorMapName]["Coefficients"] = mapCoefficients + + outWCS = cls.task._make_afw_wcs( + mapDict, + raDec.getRa(), + raDec.getDec(), + doNormalizePixels=True, + xScale=xscale, + yScale=yscale, + ) catalog[d].setId(detectorId) catalog[d].setWcs(outWCS) @@ -359,25 +387,27 @@ def test_get_exposure_info(self): yy = np.linspace(0, 4000, 6) xgrid, ygrid = np.meshgrid(xx, yy) for visSum in self.inputVisitSummary: - visit = visSum[0]['visit'] + visit = visSum[0]["visit"] for detectorInfo in visSum: - detector = detectorInfo['id'] - extensionIndex = np.flatnonzero((self.extensionInfo.visit == visit) - & (self.extensionInfo.detector == detector))[0] + detector = detectorInfo["id"] + extensionIndex = np.flatnonzero( + (self.extensionInfo.visit == visit) & (self.extensionInfo.detector == detector) + )[0] fitWcs = self.extensionInfo.wcs[extensionIndex] calexpWcs = detectorInfo.getWcs() - tanPlaneXY = np.array([fitWcs.toWorld(x, y) for (x, y) in zip(xgrid.ravel(), - ygrid.ravel())]) + tanPlaneXY = np.array([fitWcs.toWorld(x, y) for (x, y) in zip(xgrid.ravel(), ygrid.ravel())]) calexpra, calexpdec = calexpWcs.pixelToSkyArray(xgrid.ravel(), ygrid.ravel(), degrees=True) - tangentPoint = calexpWcs.pixelToSky(calexpWcs.getPixelOrigin().getX(), - calexpWcs.getPixelOrigin().getY()) + tangentPoint = calexpWcs.pixelToSky( + calexpWcs.getPixelOrigin().getX(), calexpWcs.getPixelOrigin().getY() + ) cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0 * lsst.geom.degrees, True) iwcToSkyWcs = afwgeom.makeSkyWcs(lsst.geom.Point2D(0, 0), tangentPoint, cdMatrix) - newRAdeg, newDecdeg = iwcToSkyWcs.pixelToSkyArray(tanPlaneXY[:, 0], tanPlaneXY[:, 1], - degrees=True) + newRAdeg, newDecdeg = iwcToSkyWcs.pixelToSkyArray( + tanPlaneXY[:, 0], tanPlaneXY[:, 1], degrees=True + ) # One WCS is in SIP and the other is TPV. The pixel-to-sky # conversion is not exactly the same but should be close. @@ -385,20 +415,29 @@ def test_get_exposure_info(self): # particularly for detector # >= 100. See if we can improve/if # improving is necessary. Check if matching in corner detectors # is ok. - rtol = (1e-3 if (detector >= 100) else 1e-5) + rtol = 1e-3 if (detector >= 100) else 1e-5 np.testing.assert_allclose(calexpra, newRAdeg, rtol=rtol) np.testing.assert_allclose(calexpdec, newDecdeg, rtol=rtol) def test_refCatLoader(self): - """Test that we can load objects from refCat - """ - - tmpAssociations = wcsfit.FoFClass(self.fields, [self.instrument], self.exposuresHelper, - [self.fieldRadius.asDegrees()], - (self.task.config.matchRadius * u.arcsec).to(u.degree).value) - - self.task._load_refcat(tmpAssociations, self.refObjectLoader, self.fieldCenter, self.fieldRadius, - self.extensionInfo, epoch=2015) + """Test that we can load objects from refCat""" + + tmpAssociations = wcsfit.FoFClass( + self.fields, + [self.instrument], + self.exposuresHelper, + [self.fieldRadius.asDegrees()], + (self.task.config.matchRadius * u.arcsec).to(u.degree).value, + ) + + self.task._load_refcat( + tmpAssociations, + self.refObjectLoader, + self.fieldCenter, + self.fieldRadius, + self.extensionInfo, + epoch=2015, + ) # We have only loaded one catalog, so getting the 'matches' should just # return the same objects we put in, except some random objects that @@ -411,46 +450,48 @@ def test_refCatLoader(self): self.assertGreater(nMatches, self.nStars * 0.9) def test_load_catalogs_and_associate(self): - - tmpAssociations = wcsfit.FoFClass(self.fields, [self.instrument], self.exposuresHelper, - [self.fieldRadius.asDegrees()], - (self.task.config.matchRadius * u.arcsec).to(u.degree).value) + tmpAssociations = wcsfit.FoFClass( + self.fields, + [self.instrument], + self.exposuresHelper, + [self.fieldRadius.asDegrees()], + (self.task.config.matchRadius * u.arcsec).to(u.degree).value, + ) self.task._load_catalogs_and_associate(tmpAssociations, self.inputCatalogRefs, self.extensionInfo) tmpAssociations.sortMatches(self.fieldNumber, minMatches=2) matchIds = [] correctMatches = [] - for (s, e, o) in zip(tmpAssociations.sequence, tmpAssociations.extn, tmpAssociations.obj): + for s, e, o in zip(tmpAssociations.sequence, tmpAssociations.extn, tmpAssociations.obj): objVisitInd = self.extensionInfo.visitIndex[e] objDet = self.extensionInfo.detector[e] - ExtnInds = self.inputCatalogRefs[objVisitInd].get()['detector'] == objDet + ExtnInds = self.inputCatalogRefs[objVisitInd].get()["detector"] == objDet objInfo = self.inputCatalogRefs[objVisitInd].get()[ExtnInds].iloc[o] if s == 0: if len(matchIds) > 0: correctMatches.append(len(set(matchIds)) == 1) matchIds = [] - matchIds.append(objInfo['sourceId']) + matchIds.append(objInfo["sourceId"]) # A few matches may incorrectly associate sources because of the random # positions self.assertGreater(sum(correctMatches), len(correctMatches) * 0.95) def test_make_outputs(self): - """Test that the run method recovers the input model parameters. - """ + """Test that the run method recovers the input model parameters.""" for v, visit in enumerate(self.testVisits): visitSummary = self.inputVisitSummary[v] outputWcsCatalog = self.outputs.outputWCSs[visit] visitSources = self.inputCatalogRefs[v].get() for d, detectorRow in enumerate(visitSummary): - detectorId = detectorRow['id'] + detectorId = detectorRow["id"] fitwcs = outputWcsCatalog[d].getWcs() - detSources = visitSources[visitSources['detector'] == detectorId] - fitRA, fitDec = fitwcs.pixelToSkyArray(detSources['x'], detSources['y'], degrees=True) - dRA = fitRA - detSources['trueRA'] - dDec = fitDec - detSources['trueDec'] + detSources = visitSources[visitSources["detector"] == detectorId] + fitRA, fitDec = fitwcs.pixelToSkyArray(detSources["x"], detSources["y"], degrees=True) + dRA = fitRA - detSources["trueRA"] + dDec = fitDec - detSources["trueDec"] # Check that input coordinates match the output coordinates self.assertAlmostEqual(np.mean(dRA), 0) self.assertAlmostEqual(np.std(dRA), 0) @@ -458,21 +499,20 @@ def test_make_outputs(self): self.assertAlmostEqual(np.std(dDec), 0) def test_run(self): - """Test that run method recovers the input model parameters - """ + """Test that run method recovers the input model parameters""" outputMaps = self.outputs.fitModel.mapCollection.getParamDict() for v, visit in enumerate(self.testVisits): visitSummary = self.inputVisitSummary[v] - visitMapName = f'{visit}/poly' + visitMapName = f"{visit}/poly" origModel = self.trueModel[visitMapName] - if origModel['Type'] != 'Identity': + if origModel["Type"] != "Identity": fitModel = outputMaps[visitMapName] - origXPoly = origModel['XPoly']['Coefficients'] - origYPoly = origModel['YPoly']['Coefficients'] - fitXPoly = fitModel[:len(origXPoly)] - fitYPoly = fitModel[len(origXPoly):] + origXPoly = origModel["XPoly"]["Coefficients"] + origYPoly = origModel["YPoly"]["Coefficients"] + fitXPoly = fitModel[: len(origXPoly)] + fitYPoly = fitModel[len(origXPoly) :] absDiffX = abs(fitXPoly - origXPoly) absDiffY = abs(fitYPoly - origYPoly) @@ -480,15 +520,15 @@ def test_run(self): np.testing.assert_array_less(absDiffX, 1e-6) np.testing.assert_array_less(absDiffY, 1e-6) for d, detectorRow in enumerate(visitSummary): - detectorId = detectorRow['id'] - detectorMapName = f'HSC/{detectorId}/poly' + detectorId = detectorRow["id"] + detectorMapName = f"HSC/{detectorId}/poly" origModel = self.trueModel[detectorMapName] - if (origModel['Type'] != 'Identity') and (v == 0): + if (origModel["Type"] != "Identity") and (v == 0): fitModel = outputMaps[detectorMapName] - origXPoly = origModel['XPoly']['Coefficients'] - origYPoly = origModel['YPoly']['Coefficients'] - fitXPoly = fitModel[:len(origXPoly)] - fitYPoly = fitModel[len(origXPoly):] + origXPoly = origModel["XPoly"]["Coefficients"] + origYPoly = origModel["YPoly"]["Coefficients"] + fitXPoly = fitModel[: len(origXPoly)] + fitYPoly = fitModel[len(origXPoly) :] absDiffX = abs(fitXPoly - origXPoly) absDiffY = abs(fitYPoly - origYPoly) # Check that input detector model matches fit diff --git a/tests/test_update_visit_summary.py b/tests/test_update_visit_summary.py index 5dc34878..20bd88b8 100644 --- a/tests/test_update_visit_summary.py +++ b/tests/test_update_visit_summary.py @@ -31,15 +31,13 @@ class UpdateVisitSummaryTestCase(unittest.TestCase): - def setUp(self) -> None: self.input_schema = ExposureTable.makeMinimalSchema() ExposureSummaryStats.update_schema(self.input_schema) self.init_inputs = {"input_summary_schema": ExposureCatalog(self.input_schema)} def test_wcs_provider(self) -> None: - """Test the wcs_provider config option's effect on connections. - """ + """Test the wcs_provider config's effect on connections.""" config = UpdateVisitSummaryConfig() config.wcs_provider = "input_summary" connections = UpdateVisitSummaryConnections(config=config) @@ -62,8 +60,7 @@ def test_wcs_provider(self) -> None: self.assertEqual(task.schema, self.input_schema) def test_photo_calib_provider(self) -> None: - """Test the photo_calib_provider config option's effect on connections. - """ + """Test the photo_calib_provider config's effect on connections.""" config = UpdateVisitSummaryConfig() config.photo_calib_provider = "input_summary" connections = UpdateVisitSummaryConnections(config=config) @@ -86,8 +83,7 @@ def test_photo_calib_provider(self) -> None: self.assertEqual(task.schema, self.input_schema) def test_background_provider(self) -> None: - """Test the background_provider config option's effect on connections. - """ + """Test the background_provider config's effect on connections.""" config = UpdateVisitSummaryConfig() config.background_provider = "input_summary" connections = UpdateVisitSummaryConnections(config=config) From f30c70860a3b9b1b8a87d83093c95bd3c837e29a Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 13:09:09 -0700 Subject: [PATCH 5/7] Apply isort to all files --- doc/conf.py | 1 - python/lsst/__init__.py | 1 + python/lsst/drp/__init__.py | 1 + python/lsst/drp/tasks/__init__.py | 3 +-- python/lsst/drp/tasks/assemble_chi2_coadd.py | 11 ++++----- python/lsst/drp/tasks/assemble_coadd.py | 25 ++++++++++---------- python/lsst/drp/tasks/dcr_assemble_coadd.py | 25 +++++++++++--------- python/lsst/drp/tasks/forcedPhotCoadd.py | 6 ++--- python/lsst/drp/tasks/gbdesAstrometricFit.py | 20 +++++++++------- tests/assemble_coadd_test_utils.py | 18 +++++++------- tests/test_assemble_coadd.py | 18 +++++++------- tests/test_dcr_assemble_coadd.py | 3 +-- tests/test_gbdesAstrometricFit.py | 20 ++++++++-------- tests/test_update_visit_summary.py | 2 +- 14 files changed, 78 insertions(+), 76 deletions(-) diff --git a/doc/conf.py b/doc/conf.py index bbbeabc3..4c636f17 100644 --- a/doc/conf.py +++ b/doc/conf.py @@ -7,7 +7,6 @@ from documenteer.conf.pipelinespkg import * - project = "drp_tasks" html_theme_options["logotext"] = project html_title = project diff --git a/python/lsst/__init__.py b/python/lsst/__init__.py index bb61062c..f77af49c 100644 --- a/python/lsst/__init__.py +++ b/python/lsst/__init__.py @@ -1,2 +1,3 @@ import pkgutil + __path__ = pkgutil.extend_path(__path__, __name__) diff --git a/python/lsst/drp/__init__.py b/python/lsst/drp/__init__.py index bb61062c..f77af49c 100644 --- a/python/lsst/drp/__init__.py +++ b/python/lsst/drp/__init__.py @@ -1,2 +1,3 @@ import pkgutil + __path__ = pkgutil.extend_path(__path__, __name__) diff --git a/python/lsst/drp/tasks/__init__.py b/python/lsst/drp/tasks/__init__.py index b8d4bb6f..49b66b0c 100644 --- a/python/lsst/drp/tasks/__init__.py +++ b/python/lsst/drp/tasks/__init__.py @@ -19,6 +19,5 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -from .version import * # Generated by sconsUtils - from .forcedPhotCoadd import * +from .version import * # Generated by sconsUtils diff --git a/python/lsst/drp/tasks/assemble_chi2_coadd.py b/python/lsst/drp/tasks/assemble_chi2_coadd.py index 214ac63a..bb8b76ea 100644 --- a/python/lsst/drp/tasks/assemble_chi2_coadd.py +++ b/python/lsst/drp/tasks/assemble_chi2_coadd.py @@ -21,18 +21,17 @@ import logging -import numpy as np -from lsst.afw.detection import Psf -import lsst.afw.math as afwMath import lsst.afw.image as afwImage +import lsst.afw.math as afwMath import lsst.afw.table as afwTable -from lsst.meas.algorithms import SourceDetectionTask -from lsst.meas.base import SkyMapIdGeneratorConfig import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase import lsst.pipe.base.connectionTypes as cT - +import numpy as np +from lsst.afw.detection import Psf +from lsst.meas.algorithms import SourceDetectionTask +from lsst.meas.base import SkyMapIdGeneratorConfig log = logging.getLogger(__name__) diff --git a/python/lsst/drp/tasks/assemble_coadd.py b/python/lsst/drp/tasks/assemble_coadd.py index 56623668..f56da360 100644 --- a/python/lsst/drp/tasks/assemble_coadd.py +++ b/python/lsst/drp/tasks/assemble_coadd.py @@ -28,30 +28,31 @@ ] import copy -import numpy -import warnings import logging -import lsst.pex.config as pexConfig -import lsst.pex.exceptions as pexExceptions -import lsst.geom as geom +import warnings + import lsst.afw.geom as afwGeom import lsst.afw.image as afwImage import lsst.afw.math as afwMath import lsst.afw.table as afwTable import lsst.coadd.utils as coaddUtils -import lsst.pipe.base as pipeBase +import lsst.geom as geom import lsst.meas.algorithms as measAlg -import lsstDebug +import lsst.pex.config as pexConfig +import lsst.pex.exceptions as pexExceptions +import lsst.pipe.base as pipeBase import lsst.utils as utils -from lsst.skymap import BaseSkyMap +import lsstDebug +import numpy +from deprecated.sphinx import deprecated +from lsst.meas.algorithms import AccumulatorMeanStack, ScaleVarianceTask, SourceDetectionTask from lsst.pipe.tasks.coaddBase import CoaddBaseTask, makeSkyInfo, reorderAndPadList, subBBoxIter +from lsst.pipe.tasks.healSparseMapping import HealSparseInputMapTask from lsst.pipe.tasks.interpImage import InterpImageTask -from lsst.pipe.tasks.scaleZeroPoint import ScaleZeroPointTask from lsst.pipe.tasks.maskStreaks import MaskStreaksTask -from lsst.pipe.tasks.healSparseMapping import HealSparseInputMapTask -from lsst.meas.algorithms import SourceDetectionTask, AccumulatorMeanStack, ScaleVarianceTask +from lsst.pipe.tasks.scaleZeroPoint import ScaleZeroPointTask +from lsst.skymap import BaseSkyMap from lsst.utils.timer import timeMethod -from deprecated.sphinx import deprecated log = logging.getLogger(__name__) diff --git a/python/lsst/drp/tasks/dcr_assemble_coadd.py b/python/lsst/drp/tasks/dcr_assemble_coadd.py index 63e4803f..567d9efd 100644 --- a/python/lsst/drp/tasks/dcr_assemble_coadd.py +++ b/python/lsst/drp/tasks/dcr_assemble_coadd.py @@ -22,26 +22,29 @@ __all__ = ["DcrAssembleCoaddConnections", "DcrAssembleCoaddTask", "DcrAssembleCoaddConfig"] from math import ceil -import numpy as np -from scipy import ndimage -import lsst.geom as geom + import lsst.afw.image as afwImage import lsst.afw.table as afwTable import lsst.coadd.utils as coaddUtils -from lsst.ip.diffim.dcrModel import applyDcr, calculateDcr, DcrModel +import lsst.geom as geom import lsst.meas.algorithms as measAlg -from lsst.meas.base import SingleFrameMeasurementTask import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase import lsst.utils as utils -from lsst.utils.timer import timeMethod -from .assemble_coadd import (AssembleCoaddConnections, - AssembleCoaddTask, - CompareWarpAssembleCoaddConfig, - CompareWarpAssembleCoaddTask, - ) +import numpy as np +from lsst.ip.diffim.dcrModel import DcrModel, applyDcr, calculateDcr +from lsst.meas.base import SingleFrameMeasurementTask from lsst.pipe.tasks.coaddBase import makeSkyInfo, subBBoxIter from lsst.pipe.tasks.measurePsf import MeasurePsfTask +from lsst.utils.timer import timeMethod +from scipy import ndimage + +from .assemble_coadd import ( + AssembleCoaddConnections, + AssembleCoaddTask, + CompareWarpAssembleCoaddConfig, + CompareWarpAssembleCoaddTask, +) class DcrAssembleCoaddConnections( diff --git a/python/lsst/drp/tasks/forcedPhotCoadd.py b/python/lsst/drp/tasks/forcedPhotCoadd.py index 8a39edbc..c8229b16 100644 --- a/python/lsst/drp/tasks/forcedPhotCoadd.py +++ b/python/lsst/drp/tasks/forcedPhotCoadd.py @@ -19,15 +19,13 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import lsst.pex.config import lsst.afw.table - +import lsst.pex.config import lsst.pipe.base as pipeBase - from lsst.meas.base._id_generator import SkyMapIdGeneratorConfig -from lsst.meas.base.forcedMeasurement import ForcedMeasurementTask from lsst.meas.base.applyApCorr import ApplyApCorrTask from lsst.meas.base.catalogCalculation import CatalogCalculationTask +from lsst.meas.base.forcedMeasurement import ForcedMeasurementTask from lsst.meas.extensions.scarlet.io import updateCatalogFootprints __all__ = ("ForcedPhotCoaddConfig", "ForcedPhotCoaddTask") diff --git a/python/lsst/drp/tasks/gbdesAstrometricFit.py b/python/lsst/drp/tasks/gbdesAstrometricFit.py index 7b0c8fe8..223420fe 100644 --- a/python/lsst/drp/tasks/gbdesAstrometricFit.py +++ b/python/lsst/drp/tasks/gbdesAstrometricFit.py @@ -19,22 +19,24 @@ # the GNU General Public License along with this program. If not, # see . # -import numpy as np +import astropy.coordinates import astropy.time import astropy.units as u -import astropy.coordinates -import yaml -import wcsfit import astshim - +import lsst.afw.geom as afwgeom +import lsst.afw.table import lsst.geom import lsst.pex.config as pexConfig import lsst.pipe.base as pipeBase import lsst.sphgeom -import lsst.afw.table -import lsst.afw.geom as afwgeom -from lsst.meas.algorithms import (LoadReferenceObjectsConfig, ReferenceObjectLoader, - ReferenceSourceSelectorTask) +import numpy as np +import wcsfit +import yaml +from lsst.meas.algorithms import ( + LoadReferenceObjectsConfig, + ReferenceObjectLoader, + ReferenceSourceSelectorTask, +) from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry __all__ = ["GbdesAstrometricFitConnections", "GbdesAstrometricFitConfig", "GbdesAstrometricFitTask"] diff --git a/tests/assemble_coadd_test_utils.py b/tests/assemble_coadd_test_utils.py index fd885adc..1a55e5c6 100644 --- a/tests/assemble_coadd_test_utils.py +++ b/tests/assemble_coadd_test_utils.py @@ -25,22 +25,20 @@ This is not intended to test accessing data with the Butler and instead uses mock Butler data references to pass in the simulated data. """ -from astropy.time import Time -from astropy import units as u -from astropy.coordinates import SkyCoord, EarthLocation, Angle -import numpy as np - -from lsst.afw.cameraGeom.testUtils import DetectorWrapper import lsst.afw.geom as afwGeom import lsst.afw.image as afwImage import lsst.geom as geom +import lsst.pipe.base as pipeBase +import numpy as np +from astro_metadata_translator import makeObservationInfo +from astropy import units as u +from astropy.coordinates import Angle, EarthLocation, SkyCoord +from astropy.time import Time +from lsst.afw.cameraGeom.testUtils import DetectorWrapper from lsst.geom import arcseconds, degrees from lsst.meas.algorithms.testUtils import plantSources from lsst.obs.base import MakeRawVisitInfoViaObsInfo -import lsst.pipe.base as pipeBase -from lsst.pipe.tasks.coaddInputRecorder import CoaddInputRecorderTask, CoaddInputRecorderConfig - -from astro_metadata_translator import makeObservationInfo +from lsst.pipe.tasks.coaddInputRecorder import CoaddInputRecorderConfig, CoaddInputRecorderTask __all__ = ["makeMockSkyInfo", "MockCoaddTestData"] diff --git a/tests/test_assemble_coadd.py b/tests/test_assemble_coadd.py index f0ee8d67..ee0ccba0 100644 --- a/tests/test_assemble_coadd.py +++ b/tests/test_assemble_coadd.py @@ -22,16 +22,18 @@ """Test AssembleCoaddTask and its variants. """ import unittest -import numpy as np - -import lsst.utils.tests import lsst.pipe.base as pipeBase -from lsst.drp.tasks.assemble_coadd import (AssembleCoaddTask, AssembleCoaddConfig, - CompareWarpAssembleCoaddTask, CompareWarpAssembleCoaddConfig, - ) -from lsst.drp.tasks.dcr_assemble_coadd import DcrAssembleCoaddTask, DcrAssembleCoaddConfig -from assemble_coadd_test_utils import makeMockSkyInfo, MockCoaddTestData +import lsst.utils.tests +import numpy as np +from assemble_coadd_test_utils import MockCoaddTestData, makeMockSkyInfo +from lsst.drp.tasks.assemble_coadd import ( + AssembleCoaddConfig, + AssembleCoaddTask, + CompareWarpAssembleCoaddConfig, + CompareWarpAssembleCoaddTask, +) +from lsst.drp.tasks.dcr_assemble_coadd import DcrAssembleCoaddConfig, DcrAssembleCoaddTask __all__ = [ "MockAssembleCoaddConfig", diff --git a/tests/test_dcr_assemble_coadd.py b/tests/test_dcr_assemble_coadd.py index 916d329c..076c1515 100644 --- a/tests/test_dcr_assemble_coadd.py +++ b/tests/test_dcr_assemble_coadd.py @@ -22,8 +22,7 @@ import unittest import lsst.utils.tests - -from lsst.drp.tasks.dcr_assemble_coadd import DcrAssembleCoaddTask, DcrAssembleCoaddConfig +from lsst.drp.tasks.dcr_assemble_coadd import DcrAssembleCoaddConfig, DcrAssembleCoaddTask class DcrAssembleCoaddCalculateGainTestCase(lsst.utils.tests.TestCase): diff --git a/tests/test_gbdesAstrometricFit.py b/tests/test_gbdesAstrometricFit.py index ddc0beae..04b69e96 100644 --- a/tests/test_gbdesAstrometricFit.py +++ b/tests/test_gbdesAstrometricFit.py @@ -19,24 +19,24 @@ # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import unittest import os.path -import numpy as np -import yaml +import unittest + import astropy.units as u +import lsst.afw.geom as afwgeom +import lsst.afw.table as afwTable +import lsst.geom +import lsst.utils +import numpy as np import pandas as pd import wcsfit - -import lsst.utils -import lsst.afw.table as afwTable -import lsst.afw.geom as afwgeom -from lsst.drp.tasks.gbdesAstrometricFit import GbdesAstrometricFitConfig, GbdesAstrometricFitTask +import yaml +from lsst import sphgeom from lsst.daf.base import PropertyList +from lsst.drp.tasks.gbdesAstrometricFit import GbdesAstrometricFitConfig, GbdesAstrometricFitTask from lsst.meas.algorithms import ReferenceObjectLoader from lsst.meas.algorithms.testUtils import MockRefcatDataId from lsst.pipe.base import InMemoryDatasetHandle -from lsst import sphgeom -import lsst.geom TESTDIR = os.path.abspath(os.path.dirname(__file__)) diff --git a/tests/test_update_visit_summary.py b/tests/test_update_visit_summary.py index 20bd88b8..5969264b 100644 --- a/tests/test_update_visit_summary.py +++ b/tests/test_update_visit_summary.py @@ -21,8 +21,8 @@ import unittest -from lsst.afw.table import ExposureCatalog, ExposureTable from lsst.afw.image import ExposureSummaryStats +from lsst.afw.table import ExposureCatalog, ExposureTable from lsst.drp.tasks.update_visit_summary import ( UpdateVisitSummaryConfig, UpdateVisitSummaryConnections, From bfd69b4b108ce7b39a5ab5383629095c3ae184e9 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 17:04:29 -0700 Subject: [PATCH 6/7] Add the GitHub Actions workflow for formatting --- .github/workflows/formatting.yaml | 11 +++++++++++ 1 file changed, 11 insertions(+) create mode 100644 .github/workflows/formatting.yaml diff --git a/.github/workflows/formatting.yaml b/.github/workflows/formatting.yaml new file mode 100644 index 00000000..deecc490 --- /dev/null +++ b/.github/workflows/formatting.yaml @@ -0,0 +1,11 @@ +name: check Python formatting + +on: + push: + branches: + - main + pull_request: null + +jobs: + call-workflow: + uses: lsst/rubin_workflows/.github/workflows/formatting.yaml@main From 2f0ae690d6f3040e29012b0dd095a050a90d24b9 Mon Sep 17 00:00:00 2001 From: Arun Kannawadi Date: Thu, 26 Oct 2023 20:24:56 -0700 Subject: [PATCH 7/7] Ignore formatting commits in git blame --- .git-blame-ignore-revs | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 .git-blame-ignore-revs diff --git a/.git-blame-ignore-revs b/.git-blame-ignore-revs new file mode 100644 index 00000000..d8f08c58 --- /dev/null +++ b/.git-blame-ignore-revs @@ -0,0 +1,7 @@ +# Run this command to always ignore formatting commits in `git blame` +# git config blame.ignoreRevsFile .git-blame-ignore-revs + +# Changes to be compliant with black +929500d35832c1d3925e8a91aaca136b59652824 +# Changes to be compliant with isort +f30c70860a3b9b1b8a87d83093c95bd3c837e29a