From e403a2bd39b993bb11ce4d777092400b0a76e48b Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sat, 16 Aug 2025 00:34:55 -0700 Subject: [PATCH 1/8] add analysis and task --- python/lsst/analysis/tools/atools/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/python/lsst/analysis/tools/atools/__init__.py b/python/lsst/analysis/tools/atools/__init__.py index 4eadc6beb..2d00bdac8 100644 --- a/python/lsst/analysis/tools/atools/__init__.py +++ b/python/lsst/analysis/tools/atools/__init__.py @@ -39,4 +39,5 @@ from .skySource import * from .sources import * from .stellarLocus import * +from .skyBrightnessPrecision import * from .wholeSkyPlotTool import * From 84f601af61fed111bfd9336116a00e57579a788f Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sat, 16 Aug 2025 00:35:17 -0700 Subject: [PATCH 2/8] add analysis and task --- .../tools/atools/skyBrightnessPrecision.py | 48 ++ .../tasks/skyBrightnessPrecisionAnalysis.py | 411 ++++++++++++++++++ 2 files changed, 459 insertions(+) create mode 100644 python/lsst/analysis/tools/atools/skyBrightnessPrecision.py create mode 100644 python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py diff --git a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py new file mode 100644 index 000000000..dc293139a --- /dev/null +++ b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py @@ -0,0 +1,48 @@ +# This file is part of analysis_tools. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +__all__ = ("SkyBrightnessPrecisionHistPlot",) + +from ..actions.plot.histPlot import HistPanel, HistPlot +from ..actions.vector import LoadVector +from ..interfaces import AnalysisTool + + +class SkyBrightnessPrecisionHistPlot(AnalysisTool): + parameterizedBand: bool = False + + def setDefaults(self): + super().setDefaults() + self.process.buildActions.hist_sky_brightness_precision = LoadVector( + vectorKey="sky_brightness_precision" + ) + + self.produce.plot = HistPlot() + + self.produce.plot.panels["panel_flux"] = HistPanel() + self.produce.plot.panels["panel_flux"].label = "Sky Brightness Ratio" + self.produce.plot.panels["panel_flux"].hists = dict( + hist_sky_brightness_precision=r"Sky Brightness Ratio", + ) + # self.produce.plot.panels["panel_flux"].rangeType = "sigmaMad" + # self.produce.plot.panels["panel_flux"].lowerRange = 3.5 + # self.produce.plot.panels["panel_flux"].upperRange = 3.5 + # self.produce.plot.panels["panel_flux"].validate() \ No newline at end of file diff --git a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py new file mode 100644 index 000000000..17a5242da --- /dev/null +++ b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py @@ -0,0 +1,411 @@ +# This file is part of analysis_tools. +# +# Developed for the LSST Data Management System. +# This product includes software developed by the LSST Project +# (https://www.lsst.org). +# See the COPYRIGHT file at the top-level directory of this distribution +# for details of code ownership. +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . + +__all__ = ( + "SkyBrightnessPrecisionTask", + "SkyBrightnessPrecisionConfig", + "SkyBrightnessPrecisionConnections", + "SkyBrightnessPrecisionAnalysisConnections", + "SkyBrightnessPrecisionAnalysisConfig", + "SkyBrightnessPrecisionAnalysisTask", +) + +import logging + +import lsst.afw.math as afwMath +import numpy as np +from astropy.table import Table +from lsst.pex.config import Field, ListField +from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct +from lsst.pipe.base.connectionTypes import Input, Output +from lsst.skymap import BaseSkyMap +from lsst.afw.geom import SpanSet + + +from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask + +log = logging.getLogger(__name__) + + +class SkyBrightnessPrecisionConnections( + PipelineTaskConnections, + dimensions=(), + defaultTamplates={ + "detectionTableName": "object_all", + "calexpName": "deep_coadd.calexp", + "backgroundName": "deep_coadd.calexpBackground", + "photoCalibName": "deep_coadd.photoCalib", + "wcsName": "deep_coadd.wcs", + }, +): + """Connections for the SkyBrightnessPrecisionTask.""" + + objTable = Input( + doc="Visit- or coadd-level object table", + name="{detectionTableName}", + storageClass="ArrowAstropy", + multiple=True, + dimensions=(), + deferLoad=True, + ) + + calexp = Input( + name="{calexpName}", + storageClass="ExposureF", + doc="Calibrated exposure (used for photoCalib and image array).", + multiple=True, + dimensions=(), + deferLoad=True, + ) + + calexpBackground = Input( + name="{backgroundName}", + storageClass="Background", + doc="Background model for calexp.", + multiple=True, + dimensions=(), + deferLoad=True, + ) + + skymap = Input( + doc="The skymap that covers the originating data's tract.", + name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, + storageClass="SkyMap", + dimensions=("skymap",), + ) + + photoCalib = Input( + name="{photoCalibName}", + storageClass="PhotoCalib", + doc="Photometric calibration associated with the originating data image.", + multiple=True, + dimensions=(), + deferLoad=True, + ) + + wcs = Input( + name="{wcsName}", + storageClass="Wcs", + doc="WCS header associated with the originating data image.", + multiple=True, + dimensions=(), + deferLoad=True, + ) + + sky_brightness_precision_table = Output( + name="sky_brightness_precision_table", + storageClass="ArrowAstropy", + doc="A table containing two columns: the detector or patch IDs and the values of limiting surface " + "brightness derived for those detectors or patches.", # TODO: change doc + dimensions=(), + ) + + def __init__(self, *, config=None): + super().__init__(config=config) + # Update output table name for configurable dimensions + dimen = "_visit" if "visit" in config.inputTableDimensions else "_tract" + self.objTable = Input( + doc=self.data.doc, + name=self.data.name, + storageClass=self.data.storageClass, + deferLoad=self.data.deferLoad, + dimensions=frozenset(sorted(config.inputTableDimensions)), + multiple=self.data.multiple, + ) + self.calexp = Input( + doc=self.calexp.doc, + name=self.calexp.name, + storageClass=self.calexp.storageClass, + deferLoad=self.calexp.deferLoad, + dimensions=frozenset(sorted(config.inputCalibDimensions)), + multiple=self.calexp.multiple, + ) + self.background = Input( + doc=self.background.doc, + name=self.background.name, + storageClass=self.background.storageClass, + deferLoad=self.background.deferLoad, + dimensions=frozenset(sorted(config.inputCalibDimensions)), + multiple=self.background.multiple, + ) + self.photoCalib = Input( + doc=self.photoCalib.doc, + name=self.photoCalib.name, + storageClass=self.photoCalib.storageClass, + deferLoad=self.photoCalib.deferLoad, + dimensions=frozenset(sorted(config.inputCalibDimensions)), + multiple=self.photoCalib.multiple, + ) + self.wcs = Input( + doc=self.wcs.doc, + name=self.wcs.name, + storageClass=self.wcs.storageClass, + deferLoad=self.wcs.deferLoad, + dimensions=frozenset(sorted(config.inputCalibDimensions)), + multiple=self.wcs.multiple, + ) + self.sky_brightness_precision_table = Output( + doc=self.sky_brightness_precision_table.doc, + name=self.sky_brightness_precision_table.name + dimen, + storageClass=self.sky_brightness_precision_table.storageClass, + dimensions=frozenset(sorted(config.outputDataDimensions)), + ) + + assert config is not None, "Missing required config object." + + if "tract" not in config.inputTableDimensions: + del self.skymap + + +class SkyBrightnessPrecisionConfig( + PipelineTaskConfig, + pipelineConnections=SkyBrightnessPrecisionConnections, +): + """CConfig class for SkyBrightnessPrecisionTask""" + + inputTableDimensions = ListField[str]( + doc="Dimensions of the input object table data.", + default=("skymap", "tract", "patch", "band"), + optional=False, + ) + inputCalibDimensions = ListField[str]( + doc="Dimensions of the input calibration data.", + default=("skymap", "tract", "patch", "band"), + optional=False, + ) + outputDataDimensions = ListField[str]( + doc="Dimensions of the output table data.", + default=("skymap", "tract", "patch", "band"), + optional=False, + ) + apertureSize = Field[int]( + doc="The size of the sky objects photometry aperture.", + default=9, + ) + tolerance = Field[float]( + doc="Fractional tolerance for the precision check (1% = 0.01).", + default=0.01, + ) + + +class SkyBrightnessPrecisionTask(PipelineTask): + """A task for measuring the error in the precision of the sky brightness specified by OSS-REQ-0387-V-05""" + + ConfigClass = SkyBrightnessPrecisionConfig + _DefaultName = "skyBrightnessPrecision" + + def __init__(self, initInputs=None, *args, **kwargs): + super().__init__(*args, **kwargs) + + @staticmethod + def _disk_mask(radius_px: int): + r = int(radius_px) + yy, xx = np.ogrid[-r : r + 1, -r : r + 1] + return (xx * xx + yy * yy) <= r * r + + @staticmethod + def _mean_in_aperture( + image_array: np.ndarray, x: float, y: float, disk_mask: np.ndarray, nAperPixels: int, aper: int + ) -> float: + """Mean over an *effective* area for circular aperture centered at (x,y).""" + r = disk_mask.shape[0] // 2 + x0 = int(np.round(x)) - r + y0 = int(np.round(y)) - r + x1 = x0 + disk_mask.shape[1] + y1 = y0 + disk_mask.shape[0] + + H, W = image_array.shape + xs0, ys0 = max(0, x0), max(0, y0) + xs1, ys1 = min(W, x1), min(H, y1) + if xs0 >= xs1 or ys0 >= ys1: + return 0.0 + + im_cut = image_array[ys0:ys1, xs0:xs1] + mk_cut = disk_mask[ + (ys0 - y0) : (ys1 - y1 + disk_mask.shape[0]), (xs0 - x0) : (xs1 - x1 + disk_mask.shape[1]) + ] + mk_cut = mk_cut[: im_cut.shape[0], : im_cut.shape[1]] + + nPix = int(mk_cut.sum()) + if nPix == 0: + return 0.0 + + ideal_area = np.pi * (aper**2) + eff_area = ideal_area * (nPix / nAperPixels) + return float(im_cut[mk_cut].sum() / eff_area) + + def runQuantum(self, butlerQC, inputRefs, outputRefs): + inputs = butlerQC.get(inputRefs) + sky_brightness_precision_struct = self.run(**{k: v for k, v in inputs.items() if k != "skymap"}) + butlerQC.put(sky_brightness_precision_struct, outputRefs) + + def run(self, objTable, calexp, background, photoCalib, wcs): + """Output a table of per-detector or per-patch sky brightness precision measurement + for the input visit or tract coadd + + Parameters + ---------- + data : `list` + List of dicts with information respecting the extracted source + detection catalogues + photoCalib : `list` + List of dicts with information respecting the extracted image + photometric calibrations + wcs : `list` + List of dicts with information respecting the extracted image + world coordinate systems + + Returns + ------- + `pipe.base.Struct` containing `astropy.table.Table` + Table containing per-detector or per-patch sky brightness precision values + """ + + # lookup_photoCalib = {x.dataId: x for x in photoCalib} + # lookup_wcs = {x.dataId: x for x in wcs} + lookup_calexp = {x.dataId: x for x in calexp} + lookup_bg = {x.dataId: x for x in background} + + source_catalog = objTable[0].get() + + if "band" in self.config.inputCalibDimensions and len(photoCalib) > 0: + band = photoCalib[0].dataId["band"] + "_" + else: + band = "" + + if "visit" in self.config.inputTableDimensions: + skyKeyName = "sky_source" + else: + skyKeyName = "sky_object" + + out = Table( + names=["imageID", "fracWithin1pct", "maxAbsErr", "medianRatio", "nSky"], + dtype=[np.int64, np.float64, np.float64, np.float64, np.int64], + ) + + aper = int(self.config.apertureSize) + ap_col = f"{band}ap{aper:02d}Flux" + disk = self._disk_mask(aper) + nAperPixels = int(SpanSet.fromShape(aper).asArray().sum()) + nPixIdeal = np.pi * aper**2 + + # Iterate over images (detectors or patches) + for dataId in lookup_calexp.keys(): + idKey = "detector" if "detector" in dataId else "patch" + + # Slice catalogue rows for this image and sky sources + isImage = source_catalog[idKey] == dataId[idKey] + if skyKeyName not in source_catalog.colnames: + raise ValueError(f"No sky flag column {skyKeyName} found for {dataId}.") + + isSky = source_catalog[skyKey] > 0 + rows = source_catalog[isImage & isSky] + + if ap_col not in rows.colnames: + raise ValueError(f"No aperture flux column {ap_col} found for {dataId}.") + if len(rows) == 0: + out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, 0]) + continue + + cal = lookup_calexp[dataId].get() + bg = lookup_bg[dataId].get() + + nano = cal.getPhotoCalib().instFluxToNanojansky(1.0) + calimg = cal.image.array.astype("f8", copy=False) * nano + bgimg = bg.getImage().array.astype("f8", copy=False) * nano + + x = np.asarray(rows["x"], dtype="f8") + y = np.asarray(rows["y"], dtype="f8") + mean_flux_sky = np.asarray(rows[ap_col], dtype="f8") / nPixIdeal + + mean_img = np.empty_like(mean_flux_sky) + mean_bg = np.empty_like(mean_flux_sky) + for i in range(mean_img.size): + mean_img[i] = self._mean_in_aperture(calimg, x[i], y[i], disk, nAperPixels, aper) + mean_bg[i] = self._mean_in_aperture(bgimg, x[i], y[i], disk, nAperPixels, aper) + + good = mean_bg != 0.0 + if not np.any(good): + out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, int(mean_bg.size)]) + continue + + # SBRatio per the notebook: (background + skyFlux) / background + sb_ratio = (mean_bg[good] + mean_flux_sky[good]) / mean_bg[good] + + tol = float(self.config.tolerance) + fracWithin = ( + np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size + ) + maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) + medianRatio = np.median(sb_ratio) + + out.add_row( + [ + int(dataId[idKey]), + float(fracWithin), + float(maxAbsErr), + float(medianRatio), + int(sb_ratio.size), + ] + ) + + return Struct(sky_brightness_precision_table=out) + + +class SkyBrightnessPrecisionAnalysisConnections( + AnalysisBaseConnections, + defaultTemplates={"outputName": "skyBrightnessPrecision"}, +): + data = Input( + name="sky_brightness_precision_table", + storageClass="ArrowAstropy", + doc="A table containing sky brightness precision metrics.", + deferLoad=True, + dimensions=(), + ) + + def __init__(self, *, config=None): + super().__init__(config=config) + + dimen = "_visit" if "visit" in config.inputDataDimensions else "_tract" + + self.data = Input( + name=self.data.name + dimen, + storageClass=self.data.storageClass, + doc=self.data.doc, + deferLoad=self.data.deferLoad, + dimensions=frozenset(sorted(config.inputDataDimensions)), + ) + + +class SkyBrightnessPrecisionAnalysisConfig( + AnalysisBaseConfig, pipelineConnections=SkyBrightnessPrecisionAnalysisConnections +): + inputDataDimensions = ListField[str]( + doc="Dimensions of the input table data.", + default=(), + optional=False, + ) + + +class SkyBrightnessPrecisionAnalysisTask(AnalysisPipelineTask): + ConfigClass = SkyBrightnessPrecisionAnalysisConfig + _DefaultName = "skyBrightnessPrecisionAnalysis" From 66d96f15b3b56fad80f932e6fc69380f1d358435 Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sat, 16 Aug 2025 00:47:47 -0700 Subject: [PATCH 3/8] add metric information --- metricInfo/metricInformation.yaml | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/metricInfo/metricInformation.yaml b/metricInfo/metricInformation.yaml index 5a9cc345c..82ed5a6ac 100644 --- a/metricInfo/metricInformation.yaml +++ b/metricInfo/metricInformation.yaml @@ -354,3 +354,13 @@ diaSourcesGoodVsBadRatio: description: | The ratio of the number of diaSources with high reliability (>0.9) over the number of diaSources with low reliability (<0.1) + +skyBrightnessPrecision: + lowThreshold: 0.95 + highThreshold: 1.00 + divergent: false + plot: + debugGroup:surfaceBrightness + atoolsFile: skyBrightnessPrecision.py + description: | + Fraction of sky brightness that is within 1% \ No newline at end of file From 71256d41c769776bf04707b62756198a30b5127f Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sat, 16 Aug 2025 22:10:59 -0700 Subject: [PATCH 4/8] running for se and coadd --- metricInfo/metricInformation.yaml | 2 +- python/lsst/analysis/tools/tasks/__init__.py | 1 + .../tasks/skyBrightnessPrecisionAnalysis.py | 249 +++++++++--------- 3 files changed, 129 insertions(+), 123 deletions(-) diff --git a/metricInfo/metricInformation.yaml b/metricInfo/metricInformation.yaml index 82ed5a6ac..f255aa5f0 100644 --- a/metricInfo/metricInformation.yaml +++ b/metricInfo/metricInformation.yaml @@ -363,4 +363,4 @@ skyBrightnessPrecision: debugGroup:surfaceBrightness atoolsFile: skyBrightnessPrecision.py description: | - Fraction of sky brightness that is within 1% \ No newline at end of file + Fraction of sky brightness precision that is within 1% \ No newline at end of file diff --git a/python/lsst/analysis/tools/tasks/__init__.py b/python/lsst/analysis/tools/tasks/__init__.py index bfa9b9d66..9c81062d0 100644 --- a/python/lsst/analysis/tools/tasks/__init__.py +++ b/python/lsst/analysis/tools/tasks/__init__.py @@ -29,6 +29,7 @@ from .refCatObjectPhotometricAnalysis import * from .refCatSourceAnalysis import * from .refCatSourcePhotometricAnalysis import * +from .skyBrightnessPrecisionAnalysis import * from .sourceObjectTableAnalysis import * from .sourceTableVisitAnalysis import * from .trailedDiaSrcDetectorVisitAnalysis import * diff --git a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py index 17a5242da..e370ba0b5 100644 --- a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py +++ b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py @@ -47,14 +47,14 @@ class SkyBrightnessPrecisionConnections( PipelineTaskConnections, - dimensions=(), - defaultTamplates={ - "detectionTableName": "object_all", - "calexpName": "deep_coadd.calexp", - "backgroundName": "deep_coadd.calexpBackground", - "photoCalibName": "deep_coadd.photoCalib", - "wcsName": "deep_coadd.wcs", + defaultTemplates={ + "detectionTableName": "sourceTable", + "calexpName": "calexp", + "backgroundName": "calexpBackground", + "photoCalibName": "calexp.photoCalib", + "wcsName": "calexp.wcs", }, + dimensions=(), ): """Connections for the SkyBrightnessPrecisionTask.""" @@ -76,7 +76,7 @@ class SkyBrightnessPrecisionConnections( deferLoad=True, ) - calexpBackground = Input( + background = Input( name="{backgroundName}", storageClass="Background", doc="Background model for calexp.", @@ -120,15 +120,17 @@ class SkyBrightnessPrecisionConnections( def __init__(self, *, config=None): super().__init__(config=config) + # Update output table name for configurable dimensions dimen = "_visit" if "visit" in config.inputTableDimensions else "_tract" + self.objTable = Input( - doc=self.data.doc, - name=self.data.name, - storageClass=self.data.storageClass, - deferLoad=self.data.deferLoad, + doc=self.objTable.doc, + name=self.objTable.name, + storageClass=self.objTable.storageClass, + deferLoad=self.objTable.deferLoad, dimensions=frozenset(sorted(config.inputTableDimensions)), - multiple=self.data.multiple, + multiple=self.objTable.multiple, ) self.calexp = Input( doc=self.calexp.doc, @@ -174,6 +176,8 @@ def __init__(self, *, config=None): if "tract" not in config.inputTableDimensions: del self.skymap + self.dimensions.update(frozenset(sorted(config.outputDataDimensions))) + class SkyBrightnessPrecisionConfig( PipelineTaskConfig, @@ -257,117 +261,118 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs): sky_brightness_precision_struct = self.run(**{k: v for k, v in inputs.items() if k != "skymap"}) butlerQC.put(sky_brightness_precision_struct, outputRefs) - def run(self, objTable, calexp, background, photoCalib, wcs): - """Output a table of per-detector or per-patch sky brightness precision measurement - for the input visit or tract coadd - - Parameters - ---------- - data : `list` - List of dicts with information respecting the extracted source - detection catalogues - photoCalib : `list` - List of dicts with information respecting the extracted image - photometric calibrations - wcs : `list` - List of dicts with information respecting the extracted image - world coordinate systems - - Returns - ------- - `pipe.base.Struct` containing `astropy.table.Table` - Table containing per-detector or per-patch sky brightness precision values - """ - - # lookup_photoCalib = {x.dataId: x for x in photoCalib} - # lookup_wcs = {x.dataId: x for x in wcs} - lookup_calexp = {x.dataId: x for x in calexp} - lookup_bg = {x.dataId: x for x in background} - - source_catalog = objTable[0].get() - - if "band" in self.config.inputCalibDimensions and len(photoCalib) > 0: - band = photoCalib[0].dataId["band"] + "_" - else: - band = "" - - if "visit" in self.config.inputTableDimensions: - skyKeyName = "sky_source" - else: - skyKeyName = "sky_object" - - out = Table( - names=["imageID", "fracWithin1pct", "maxAbsErr", "medianRatio", "nSky"], - dtype=[np.int64, np.float64, np.float64, np.float64, np.int64], - ) + def run(self, objTable, calexp, background, photoCalib, wcs): + """Output a table of per-detector or per-patch sky brightness precision measurement + for the input visit or tract coadd + + Parameters + ---------- + data : `list` + List of dicts with information respecting the extracted source + detection catalogues + photoCalib : `list` + List of dicts with information respecting the extracted image + photometric calibrations + wcs : `list` + List of dicts with information respecting the extracted image + world coordinate systems + + Returns + ------- + `pipe.base.Struct` containing `astropy.table.Table` + Table containing per-detector or per-patch sky brightness precision values + """ + + # lookup_photoCalib = {x.dataId: x for x in photoCalib} + # lookup_wcs = {x.dataId: x for x in wcs} + lookup_calexp = {x.dataId: x for x in calexp} + lookup_bg = {x.dataId: x for x in background} + + source_catalog = objTable[0].get() + + if "band" in self.config.inputCalibDimensions and len(photoCalib) > 0: + band = photoCalib[0].dataId["band"] + "_" + else: + band = "" + + + out = Table( + names=["imageID", "fracWithin1pct", "maxAbsErr", "medianRatio", "nSky"], + dtype=[np.int64, np.float64, np.float64, np.float64, np.int64], + ) - aper = int(self.config.apertureSize) + aper = int(self.config.apertureSize) + disk = self._disk_mask(aper) + nAperPixels = int(SpanSet.fromShape(aper).asArray().sum()) + nPixIdeal = np.pi * aper**2 + + if "visit" in self.config.inputTableDimensions: + skyKeyName = "sky_source" + ap_col = f"ap{aper:02d}Flux" + else: + skyKeyName = "sky_object" ap_col = f"{band}ap{aper:02d}Flux" - disk = self._disk_mask(aper) - nAperPixels = int(SpanSet.fromShape(aper).asArray().sum()) - nPixIdeal = np.pi * aper**2 - - # Iterate over images (detectors or patches) - for dataId in lookup_calexp.keys(): - idKey = "detector" if "detector" in dataId else "patch" - - # Slice catalogue rows for this image and sky sources - isImage = source_catalog[idKey] == dataId[idKey] - if skyKeyName not in source_catalog.colnames: - raise ValueError(f"No sky flag column {skyKeyName} found for {dataId}.") - - isSky = source_catalog[skyKey] > 0 - rows = source_catalog[isImage & isSky] - - if ap_col not in rows.colnames: - raise ValueError(f"No aperture flux column {ap_col} found for {dataId}.") - if len(rows) == 0: - out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, 0]) - continue - - cal = lookup_calexp[dataId].get() - bg = lookup_bg[dataId].get() - - nano = cal.getPhotoCalib().instFluxToNanojansky(1.0) - calimg = cal.image.array.astype("f8", copy=False) * nano - bgimg = bg.getImage().array.astype("f8", copy=False) * nano - - x = np.asarray(rows["x"], dtype="f8") - y = np.asarray(rows["y"], dtype="f8") - mean_flux_sky = np.asarray(rows[ap_col], dtype="f8") / nPixIdeal - - mean_img = np.empty_like(mean_flux_sky) - mean_bg = np.empty_like(mean_flux_sky) - for i in range(mean_img.size): - mean_img[i] = self._mean_in_aperture(calimg, x[i], y[i], disk, nAperPixels, aper) - mean_bg[i] = self._mean_in_aperture(bgimg, x[i], y[i], disk, nAperPixels, aper) - - good = mean_bg != 0.0 - if not np.any(good): - out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, int(mean_bg.size)]) - continue - - # SBRatio per the notebook: (background + skyFlux) / background - sb_ratio = (mean_bg[good] + mean_flux_sky[good]) / mean_bg[good] - - tol = float(self.config.tolerance) - fracWithin = ( - np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size - ) - maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) - medianRatio = np.median(sb_ratio) - - out.add_row( - [ - int(dataId[idKey]), - float(fracWithin), - float(maxAbsErr), - float(medianRatio), - int(sb_ratio.size), - ] - ) - - return Struct(sky_brightness_precision_table=out) + + # Iterate over images (detectors or patches) + for dataId in lookup_calexp.keys(): + idKey = "detector" if "detector" in dataId else "patch" + + # Slice catalogue rows for this image and sky sources + isImage = (source_catalog[idKey] == dataId[idKey]) + assert np.any(isImage), f"No sources found for {dataId}." + if skyKeyName not in source_catalog.colnames: + raise ValueError(f"No sky flag column {skyKeyName} found for {dataId}.") + + isSky = source_catalog[skyKeyName] > 0 + rows = source_catalog[isImage & isSky] + + if ap_col not in rows.colnames: + raise ValueError(f"No aperture flux column {ap_col} found for {dataId}.") + if len(rows) == 0: + out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, 0]) + continue + + cal = lookup_calexp[dataId].get() + bg = lookup_bg[dataId].get() + + nano = cal.getPhotoCalib().instFluxToNanojansky(1.0) + calimg = cal.image.array.astype("f8", copy=False) * nano + bgimg = bg.getImage().array.astype("f8", copy=False) * nano + + x = np.asarray(rows["x"], dtype="f8") + y = np.asarray(rows["y"], dtype="f8") + mean_flux_sky = np.asarray(rows[ap_col], dtype="f8") / nPixIdeal + + mean_img = np.empty_like(mean_flux_sky) + mean_bg = np.empty_like(mean_flux_sky) + for i in range(mean_img.size): + mean_img[i] = self._mean_in_aperture(calimg, x[i], y[i], disk, nAperPixels, aper) + mean_bg[i] = self._mean_in_aperture(bgimg, x[i], y[i], disk, nAperPixels, aper) + + good = mean_bg != 0.0 + if not np.any(good): + out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, int(mean_bg.size)]) + continue + + # SBRatio per the notebook: (background + skyFlux) / background + sb_ratio = (mean_bg[good] + mean_flux_sky[good]) / mean_bg[good] + + tol = float(self.config.tolerance) + fracWithin = np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size + maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) + medianRatio = np.median(sb_ratio) + + out.add_row( + [ + int(dataId[idKey]), + float(fracWithin), + float(maxAbsErr), + float(medianRatio), + int(sb_ratio.size), + ] + ) + + return Struct(sky_brightness_precision_table=out) class SkyBrightnessPrecisionAnalysisConnections( From 510929d61b2bd430fe7b95e9f7d53e3e047d7d89 Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sun, 17 Aug 2025 20:28:59 -0700 Subject: [PATCH 5/8] running se --- .../tools/atools/skyBrightnessPrecision.py | 2 +- .../tasks/skyBrightnessPrecisionAnalysis.py | 76 ++++++++----------- 2 files changed, 34 insertions(+), 44 deletions(-) diff --git a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py index dc293139a..3b355eeda 100644 --- a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py +++ b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py @@ -32,7 +32,7 @@ class SkyBrightnessPrecisionHistPlot(AnalysisTool): def setDefaults(self): super().setDefaults() self.process.buildActions.hist_sky_brightness_precision = LoadVector( - vectorKey="sky_brightness_precision" + vectorKey="fracWithin1pct" ) self.produce.plot = HistPlot() diff --git a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py index e370ba0b5..c9eb2d714 100644 --- a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py +++ b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py @@ -220,41 +220,31 @@ def __init__(self, initInputs=None, *args, **kwargs): super().__init__(*args, **kwargs) @staticmethod - def _disk_mask(radius_px: int): - r = int(radius_px) - yy, xx = np.ogrid[-r : r + 1, -r : r + 1] - return (xx * xx + yy * yy) <= r * r + def _mean_in_aperture(image, x, y, aper, npix_unmasked): + # Ensure x and y are integers and is a numpy array + x, y = int(x), int(y) + img = np.asarray(image) + rows, cols = img.shape - @staticmethod - def _mean_in_aperture( - image_array: np.ndarray, x: float, y: float, disk_mask: np.ndarray, nAperPixels: int, aper: int - ) -> float: - """Mean over an *effective* area for circular aperture centered at (x,y).""" - r = disk_mask.shape[0] // 2 - x0 = int(np.round(x)) - r - y0 = int(np.round(y)) - r - x1 = x0 + disk_mask.shape[1] - y1 = y0 + disk_mask.shape[0] - - H, W = image_array.shape - xs0, ys0 = max(0, x0), max(0, y0) - xs1, ys1 = min(W, x1), min(H, y1) - if xs0 >= xs1 or ys0 >= ys1: - return 0.0 - - im_cut = image_array[ys0:ys1, xs0:xs1] - mk_cut = disk_mask[ - (ys0 - y0) : (ys1 - y1 + disk_mask.shape[0]), (xs0 - x0) : (xs1 - x1 + disk_mask.shape[1]) - ] - mk_cut = mk_cut[: im_cut.shape[0], : im_cut.shape[1]] - - nPix = int(mk_cut.sum()) + # Create grid arrays for the entire image using ogrid (memory efficient) + X, Y = np.ogrid[:rows, :cols] + + # Create the mask for points within the circle (distance from (x, y) <= aper) + mask = (X - x) ** 2 + (Y - y) ** 2 <= aper**2 + + # Compute the mean using the actual number of pixels within the mask + nPix = np.sum(mask) + + # 9 pixel aperture has 253 pixels if unmasked, so we scale the area accordingly + area = nPix / npix_unmasked * np.pi * aper**2 + + # There are some sky sources that are out of bound, return 0 in this case if nPix == 0: - return 0.0 + meanFluxInApe = 0 + else: + meanFluxInApe = img[mask].sum() / area - ideal_area = np.pi * (aper**2) - eff_area = ideal_area * (nPix / nAperPixels) - return float(im_cut[mk_cut].sum() / eff_area) + return meanFluxInApe def runQuantum(self, butlerQC, inputRefs, outputRefs): inputs = butlerQC.get(inputRefs) @@ -295,17 +285,15 @@ def run(self, objTable, calexp, background, photoCalib, wcs): else: band = "" - out = Table( names=["imageID", "fracWithin1pct", "maxAbsErr", "medianRatio", "nSky"], dtype=[np.int64, np.float64, np.float64, np.float64, np.int64], ) aper = int(self.config.apertureSize) - disk = self._disk_mask(aper) nAperPixels = int(SpanSet.fromShape(aper).asArray().sum()) nPixIdeal = np.pi * aper**2 - + if "visit" in self.config.inputTableDimensions: skyKeyName = "sky_source" ap_col = f"ap{aper:02d}Flux" @@ -318,7 +306,7 @@ def run(self, objTable, calexp, background, photoCalib, wcs): idKey = "detector" if "detector" in dataId else "patch" # Slice catalogue rows for this image and sky sources - isImage = (source_catalog[idKey] == dataId[idKey]) + isImage = source_catalog[idKey] == dataId[idKey] assert np.any(isImage), f"No sources found for {dataId}." if skyKeyName not in source_catalog.colnames: raise ValueError(f"No sky flag column {skyKeyName} found for {dataId}.") @@ -329,14 +317,15 @@ def run(self, objTable, calexp, background, photoCalib, wcs): if ap_col not in rows.colnames: raise ValueError(f"No aperture flux column {ap_col} found for {dataId}.") if len(rows) == 0: - out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, 0]) - continue + raise ValueError(f"No valid rows found for {dataId}.") cal = lookup_calexp[dataId].get() bg = lookup_bg[dataId].get() nano = cal.getPhotoCalib().instFluxToNanojansky(1.0) calimg = cal.image.array.astype("f8", copy=False) * nano + x0, y0 = cal.getXY0() + bgimg = bg.getImage().array.astype("f8", copy=False) * nano x = np.asarray(rows["x"], dtype="f8") @@ -346,15 +335,14 @@ def run(self, objTable, calexp, background, photoCalib, wcs): mean_img = np.empty_like(mean_flux_sky) mean_bg = np.empty_like(mean_flux_sky) for i in range(mean_img.size): - mean_img[i] = self._mean_in_aperture(calimg, x[i], y[i], disk, nAperPixels, aper) - mean_bg[i] = self._mean_in_aperture(bgimg, x[i], y[i], disk, nAperPixels, aper) + mean_img[i] = self._mean_in_aperture(calimg, x[i] - x0, y[i] - y0, aper, nAperPixels) + mean_bg[i] = self._mean_in_aperture(bgimg, x[i] - x0, y[i] - y0, aper, nAperPixels) good = mean_bg != 0.0 if not np.any(good): - out.add_row([int(dataId[idKey]), np.nan, np.nan, np.nan, int(mean_bg.size)]) - continue + raise ValueError(f"No valid background flux found for {dataId}.") - # SBRatio per the notebook: (background + skyFlux) / background + # SBRatio: (background + skyFlux) / background sb_ratio = (mean_bg[good] + mean_flux_sky[good]) / mean_bg[good] tol = float(self.config.tolerance) @@ -399,6 +387,8 @@ def __init__(self, *, config=None): deferLoad=self.data.deferLoad, dimensions=frozenset(sorted(config.inputDataDimensions)), ) + + self.dimensions.update(frozenset(sorted(config.inputDataDimensions))) class SkyBrightnessPrecisionAnalysisConfig( From 8483780a0a10521ff55117c3d8713372bcf39d81 Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Sun, 17 Aug 2025 20:57:53 -0700 Subject: [PATCH 6/8] running se --- .../tools/atools/skyBrightnessPrecision.py | 2 +- .../tasks/skyBrightnessPrecisionAnalysis.py | 44 ++++++++++++------- 2 files changed, 28 insertions(+), 18 deletions(-) diff --git a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py index 3b355eeda..280f0a908 100644 --- a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py +++ b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py @@ -32,7 +32,7 @@ class SkyBrightnessPrecisionHistPlot(AnalysisTool): def setDefaults(self): super().setDefaults() self.process.buildActions.hist_sky_brightness_precision = LoadVector( - vectorKey="fracWithin1pct" + vectorKey="sb_ratio" ) self.produce.plot = HistPlot() diff --git a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py index c9eb2d714..2da4ac271 100644 --- a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py +++ b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py @@ -114,7 +114,7 @@ class SkyBrightnessPrecisionConnections( name="sky_brightness_precision_table", storageClass="ArrowAstropy", doc="A table containing two columns: the detector or patch IDs and the values of limiting surface " - "brightness derived for those detectors or patches.", # TODO: change doc + "brightness derived for those detectors or patches.", dimensions=(), ) @@ -300,6 +300,9 @@ def run(self, objTable, calexp, background, photoCalib, wcs): else: skyKeyName = "sky_object" ap_col = f"{band}ap{aper:02d}Flux" + + ids_all: list[int] = [] + ratios_all: list[float] = [] # Iterate over images (detectors or patches) for dataId in lookup_calexp.keys(): @@ -344,21 +347,28 @@ def run(self, objTable, calexp, background, photoCalib, wcs): # SBRatio: (background + skyFlux) / background sb_ratio = (mean_bg[good] + mean_flux_sky[good]) / mean_bg[good] - - tol = float(self.config.tolerance) - fracWithin = np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size - maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) - medianRatio = np.median(sb_ratio) - - out.add_row( - [ - int(dataId[idKey]), - float(fracWithin), - float(maxAbsErr), - float(medianRatio), - int(sb_ratio.size), - ] - ) + img_id = int(dataId[idKey]) + ids_all.extend([img_id] * sb_ratio.size) + ratios_all.extend(sb_ratio.astype("f8").tolist()) + + out["imageID"] = np.array(ids_all, dtype=np.int64) + out["sb_ratio"] = np.array(ratios_all, dtype=np.float64) + + + # tol = float(self.config.tolerance) + # fracWithin = np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size + # maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) + # medianRatio = np.median(sb_ratio) + + # out.add_row( + # [ + # int(dataId[idKey]), + # float(fracWithin), + # float(maxAbsErr), + # float(medianRatio), + # int(sb_ratio.size), + # ] + # ) return Struct(sky_brightness_precision_table=out) @@ -387,7 +397,7 @@ def __init__(self, *, config=None): deferLoad=self.data.deferLoad, dimensions=frozenset(sorted(config.inputDataDimensions)), ) - + self.dimensions.update(frozenset(sorted(config.inputDataDimensions))) From 03a7e434f85fc7794cb1eb2f45c67845bbb921fd Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Mon, 18 Aug 2025 11:20:08 -0700 Subject: [PATCH 7/8] finalize changes --- .../tools/atools/skyBrightnessPrecision.py | 42 ++++-- .../tasks/skyBrightnessPrecisionAnalysis.py | 141 +++++++++++------- 2 files changed, 121 insertions(+), 62 deletions(-) diff --git a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py index 280f0a908..0817ca898 100644 --- a/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py +++ b/python/lsst/analysis/tools/atools/skyBrightnessPrecision.py @@ -21,7 +21,9 @@ __all__ = ("SkyBrightnessPrecisionHistPlot",) -from ..actions.plot.histPlot import HistPanel, HistPlot +from ..actions.plot.histPlot import HistPanel, HistPlot, HistStatsPanel +from ..actions.scalar import FracInRange, MedianAction, MaxAction +from ..actions.keyedData import KeyedScalars from ..actions.vector import LoadVector from ..interfaces import AnalysisTool @@ -31,18 +33,34 @@ class SkyBrightnessPrecisionHistPlot(AnalysisTool): def setDefaults(self): super().setDefaults() - self.process.buildActions.hist_sky_brightness_precision = LoadVector( - vectorKey="sb_ratio" + + self.process.buildActions.hist_all = LoadVector(vectorKey="SBRatio") + + self.process.calculateActions.frac_1pct = FracInRange( + vectorKey="SBRatio", + minimum=0.99, + maximum=1.01, + percent=False, ) + self.process.calculateActions.median = MedianAction(vectorKey="SBRatio") + self.process.calculateActions.maximum = MaxAction(vectorKey="SBRatio") self.produce.plot = HistPlot() - self.produce.plot.panels["panel_flux"] = HistPanel() - self.produce.plot.panels["panel_flux"].label = "Sky Brightness Ratio" - self.produce.plot.panels["panel_flux"].hists = dict( - hist_sky_brightness_precision=r"Sky Brightness Ratio", - ) - # self.produce.plot.panels["panel_flux"].rangeType = "sigmaMad" - # self.produce.plot.panels["panel_flux"].lowerRange = 3.5 - # self.produce.plot.panels["panel_flux"].upperRange = 3.5 - # self.produce.plot.panels["panel_flux"].validate() \ No newline at end of file + p = HistPanel() + p.label = "All SBRatio" + p.hists = dict(hist_all="All SBRatio") + p.rangeType = "sigmaMad" + p.lowerRange = 3.5 + p.upperRange = 3.5 + p.referenceValue = 1.0 + + p.statsPanel = HistStatsPanel() + p.statsPanel.statsLabels = ["Frac in [0.99,1.01]", "Median", "Max"] + p.statsPanel.stat1 = ["frac_1pct"] + p.statsPanel.stat2 = ["median"] + p.statsPanel.stat3 = ["maximum"] + + + p.validate() + self.produce.plot.panels["panel_all"] = p diff --git a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py index 2da4ac271..b9d9a2c4c 100644 --- a/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py +++ b/python/lsst/analysis/tools/tasks/skyBrightnessPrecisionAnalysis.py @@ -23,6 +23,9 @@ "SkyBrightnessPrecisionTask", "SkyBrightnessPrecisionConfig", "SkyBrightnessPrecisionConnections", + "SkyBrightnessPrecisionAggregateConnections", + "SkyBrightnessPrecisionAggregateConfig", + "SkyBrightnessPrecisionAggregateTask", "SkyBrightnessPrecisionAnalysisConnections", "SkyBrightnessPrecisionAnalysisConfig", "SkyBrightnessPrecisionAnalysisTask", @@ -35,9 +38,10 @@ from astropy.table import Table from lsst.pex.config import Field, ListField from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections, Struct -from lsst.pipe.base.connectionTypes import Input, Output +from lsst.pipe.base.connectionTypes import Input, Output, PrerequisiteInput from lsst.skymap import BaseSkyMap from lsst.afw.geom import SpanSet +from astropy.table import vstack, Table from ..interfaces import AnalysisBaseConfig, AnalysisBaseConnections, AnalysisPipelineTask @@ -273,8 +277,7 @@ def run(self, objTable, calexp, background, photoCalib, wcs): Table containing per-detector or per-patch sky brightness precision values """ - # lookup_photoCalib = {x.dataId: x for x in photoCalib} - # lookup_wcs = {x.dataId: x for x in wcs} + lookup_calexp = {x.dataId: x for x in calexp} lookup_bg = {x.dataId: x for x in background} @@ -285,11 +288,6 @@ def run(self, objTable, calexp, background, photoCalib, wcs): else: band = "" - out = Table( - names=["imageID", "fracWithin1pct", "maxAbsErr", "medianRatio", "nSky"], - dtype=[np.int64, np.float64, np.float64, np.float64, np.int64], - ) - aper = int(self.config.apertureSize) nAperPixels = int(SpanSet.fromShape(aper).asArray().sum()) nPixIdeal = np.pi * aper**2 @@ -300,9 +298,11 @@ def run(self, objTable, calexp, background, photoCalib, wcs): else: skyKeyName = "sky_object" ap_col = f"{band}ap{aper:02d}Flux" - + ids_all: list[int] = [] ratios_all: list[float] = [] + if "detector" in self.config.inputTableDimensions: + visits_all: list[float] = [] # Iterate over images (detectors or patches) for dataId in lookup_calexp.keys(): @@ -351,66 +351,107 @@ def run(self, objTable, calexp, background, photoCalib, wcs): ids_all.extend([img_id] * sb_ratio.size) ratios_all.extend(sb_ratio.astype("f8").tolist()) - out["imageID"] = np.array(ids_all, dtype=np.int64) - out["sb_ratio"] = np.array(ratios_all, dtype=np.float64) - - - # tol = float(self.config.tolerance) - # fracWithin = np.count_nonzero((sb_ratio >= (1 - tol)) & (sb_ratio <= (1 + tol))) / sb_ratio.size - # maxAbsErr = np.max(np.abs(sb_ratio - 1.0)) - # medianRatio = np.median(sb_ratio) - - # out.add_row( - # [ - # int(dataId[idKey]), - # float(fracWithin), - # float(maxAbsErr), - # float(medianRatio), - # int(sb_ratio.size), - # ] - # ) + if "detector" in dataId: + visit = dataId["visit"] + visits_all.extend([visit] * sb_ratio.size) + + ids = np.array(ids_all, dtype=np.int64) + ratios = np.array(ratios_all, dtype=np.float64) + if visits_all: + visits = np.array(visits_all, dtype=np.int64) + out = Table({"imageID": ids, "visit": visits, "SBRatio": ratios}) + else: + out = Table({"imageID": ids, "SBRatio": ratios}) return Struct(sky_brightness_precision_table=out) -class SkyBrightnessPrecisionAnalysisConnections( - AnalysisBaseConnections, - defaultTemplates={"outputName": "skyBrightnessPrecision"}, +class SkyBrightnessPrecisionAggregateConnections( + PipelineTaskConnections, + dimensions=(), ): - data = Input( - name="sky_brightness_precision_table", + + inputTables = Input( + name="sky_brightness_precision_table_visit", storageClass="ArrowAstropy", - doc="A table containing sky brightness precision metrics.", + doc="Per-visit SB ratio tables to aggregate.", + dimensions=("instrument", "visit", "detector", "band"), + multiple=True, deferLoad=True, + ) + + # Dimensionless outputs for plotting + stacked = Output( + name="sky_brightness_precision_table_all_visit", + storageClass="ArrowAstropy", + doc="All SB ratios stacked (imageID, visit, SBRatio).", + dimensions=(), + ) + visitStats = Output( + name="sky_brightness_precision_visit_stats", + storageClass="ArrowAstropy", + doc="Per-visit summary (visit, median, max).", dimensions=(), ) - def __init__(self, *, config=None): - super().__init__(config=config) - dimen = "_visit" if "visit" in config.inputDataDimensions else "_tract" +class SkyBrightnessPrecisionAggregateConfig( + PipelineTaskConfig, pipelineConnections=SkyBrightnessPrecisionAggregateConnections +): + pass + + +class SkyBrightnessPrecisionAggregateTask(PipelineTask): + ConfigClass = SkyBrightnessPrecisionAggregateConfig + _DefaultName = "skyBrightnessPrecisionAggregate" + + def runQuantum(self, butlerQC, inputRefs, outputRefs): + inputs = butlerQC.get(inputRefs) + result = self.run(**inputs) + butlerQC.put(result, outputRefs) + + def run(self, inputTables): + tables = [ref.get() for ref in inputTables] + if not tables: + return Struct() + + big = vstack(tables, join_type="exact") + + if "visit" in big.colnames: + grouped = big.group_by("visit") + visits = grouped.groups.keys["visit"].astype("i8") + med = grouped.groups.aggregate(np.median)["SBRatio"].astype("f8") + mx = grouped.groups.aggregate(np.max)["SBRatio"].astype("f8") + stats = Table({"visit": visits, "median": med, "max": mx}) + else: + stats = Table(names=("visit", "median", "max"), dtype=("i8", "f8", "f8")) + + return Struct(stacked=big, visitStats=stats) - self.data = Input( - name=self.data.name + dimen, - storageClass=self.data.storageClass, - doc=self.data.doc, - deferLoad=self.data.deferLoad, - dimensions=frozenset(sorted(config.inputDataDimensions)), - ) - self.dimensions.update(frozenset(sorted(config.inputDataDimensions))) +class SkyBrightnessPrecisionAnalysisConnections( + AnalysisBaseConnections, defaultTemplates={"outputName": "skyBrightnessPrecision"} +): + data = Input( + name="sky_brightness_precision_table_all_visit", + storageClass="ArrowAstropy", + deferLoad=True, + dimensions=(), # dimensionless + ) + stats = Input( + name="sky_brightness_precision_visit_stats", + storageClass="ArrowAstropy", + deferLoad=True, + dimensions=(), + ) class SkyBrightnessPrecisionAnalysisConfig( AnalysisBaseConfig, pipelineConnections=SkyBrightnessPrecisionAnalysisConnections ): - inputDataDimensions = ListField[str]( - doc="Dimensions of the input table data.", - default=(), - optional=False, - ) + inputDataDimensions = ListField[str](doc="Dims for analysis", default=(), optional=False) class SkyBrightnessPrecisionAnalysisTask(AnalysisPipelineTask): ConfigClass = SkyBrightnessPrecisionAnalysisConfig - _DefaultName = "skyBrightnessPrecisionAnalysis" + _DefaultName = "skyBrightnessPrecisionAnalysis" \ No newline at end of file From e6bd1136a267cc6bb5db2aa91ba564a8fb78f626 Mon Sep 17 00:00:00 2001 From: Conghao Zhou Date: Mon, 18 Aug 2025 11:20:50 -0700 Subject: [PATCH 8/8] remove metrics --- metricInfo/metricInformation.yaml | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/metricInfo/metricInformation.yaml b/metricInfo/metricInformation.yaml index f255aa5f0..1f705dd92 100644 --- a/metricInfo/metricInformation.yaml +++ b/metricInfo/metricInformation.yaml @@ -353,14 +353,4 @@ diaSourcesGoodVsBadRatio: atoolsFile: diaSourceMetrics.py description: | The ratio of the number of diaSources with high reliability (>0.9) - over the number of diaSources with low reliability (<0.1) - -skyBrightnessPrecision: - lowThreshold: 0.95 - highThreshold: 1.00 - divergent: false - plot: - debugGroup:surfaceBrightness - atoolsFile: skyBrightnessPrecision.py - description: | - Fraction of sky brightness precision that is within 1% \ No newline at end of file + over the number of diaSources with low reliability (<0.1) \ No newline at end of file