From d61885b9ff364f91cf4460b82f1bf7e527964ed6 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Thu, 30 Jul 2026 11:45:13 -0700 Subject: [PATCH 01/30] Add interface to deep_coadds and cell-based coadd PSFs --- imsim/__init__.py | 1 + imsim/coadd_psf.py | 153 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 154 insertions(+) create mode 100644 imsim/coadd_psf.py diff --git a/imsim/__init__.py b/imsim/__init__.py index 69965189..5022d4af 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -41,3 +41,4 @@ from .table_row import * from .photon_pooling import * from .full_focal_plane import * +from .coadd_psf import * diff --git a/imsim/coadd_psf.py b/imsim/coadd_psf.py new file mode 100644 index 00000000..a599d389 --- /dev/null +++ b/imsim/coadd_psf.py @@ -0,0 +1,153 @@ +""" +Interface to Rubin deep_coadd PSFs +""" + +import galsim +from galsim.config import InputLoader, RegisterInputType, RegisterObjectType +import lsst.daf.butler as daf_butler +import lsst.geom + + +class DeepCoadd: + def __init__(self, butler, band, skymap, dstype="deep_coadd"): + """ + Parameters + ---------- + butler : lsst.daf.butler.Butler + Butler for the data repo and collection containing the deep_coadds. + band : str + Band of interest, i.e., in "ugrizy". + skymap : str + Name of the skymap to use, e.g., "lsst_cells_v2". + dstype : str + The dataset type of the cell-based coadds. Default: "deep_coadd". + """ + self.butler = butler + self.band = band + self.skymap_name = skymap + self.skymap = butler.get("skyMap", name=skymap) + self.dstype = dstype + self._psf_cache = {} + self._butler_cache = {} + + def getPSF(self, ra, dec): + """Return the cell coadd PSF, evaluated at the center of + the cell containing this sky position. + """ + # Find the tract, patch for this location. + sky_coords = lsst.geom.SpherePoint( + lsst.geom.Angle(ra*lsst.geom.degrees), + lsst.geom.Angle(dec*lsst.geom.degrees) + ) + tract_info = self.skymap.findTract(sky_coords) + tract = tract_info.getId() + patch = tract_info.findPatch(sky_coords).getSequentialIndex() + + # Cache info from the butler for the corresponding deep_coadd. + butler_key = tract, patch + if butler_key not in self._butler_cache: + dataId = dict(skymap=self.skymap_name, tract=tract, patch=patch, + band=self.band) + deep_coadd = self.butler.get(self.dstype, **dataId) + if hasattr(deep_coadd, "to_legacy"): + # This is a `lsst.images` object. Use the `.to_legacy` + # function to convert to a `lsst.afw.image.Exposure` so + # that we can access the functions to generate a PSF + # image array to pass to galsim. + deep_coadd = deep_coadd.to_legacy() + wcs = deep_coadd.getWcs() + psf_grid = deep_coadd.getPsf() # This is the grid of coadd psfs. + self._butler_cache[butler_key] = wcs, psf_grid + else: + wcs, psf_grid = self._butler_cache[butler_key] + + # Cache cell PSFs based on tract, patch, psf_grid index. + pixel_coords = wcs.skyToPixel(sky_coords) + psf_grid_index = psf_grid.grid.index(lsst.geom.Point2I(pixel_coords)) + psf_key = tract, patch, psf_grid_index + if psf_key not in self._psf_cache: + # Evaluate PSF at cell center. + x0 = ((pixel_coords.x // psf_grid.grid.cell_size.x) + * psf_grid.grid.cell_size.x + psf_grid.grid.shape.x/2) + y0 = ((pixel_coords.y // psf_grid.grid.cell_size.y) + * psf_grid.grid.cell_size.y + psf_grid.grid.shape.y/2) + pixel_cell_center = lsst.geom.Point2D(x0, y0) + sky_cell_center = wcs.pixelToSky(pixel_cell_center) + # The following is based on the lsst/source_injection + # implementation. See https://github.com/lsst/source_injection/blob/w.2026.31/python/lsst/source/injection/inject_engine.py#L705 + mat = wcs.linearizePixelToSky( + sky_cell_center, lsst.geom.arcseconds).getMatrix() + galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], + mat[1, 0], mat[1, 1]) + psf_array = psf_grid.computeKernelImage(pixel_cell_center).array + self._psf_cache[psf_key] = galsim.InterpolatedImage( + galsim.Image(psf_array), wcs=galsim_wcs) + + return self._psf_cache[psf_key] + + +class DeepCoaddLoader(InputLoader): + """ + Load the deep_coadd input object. Here's an example yaml entry: + + input.deep_coadd: + repo: dp2 + collection: LSSTCam/runs/DRP/DP2 + band: $band # This will be obtained from the opsim metadata entry. + """ + def __init__(self): + super().__init__(init_func=DeepCoadd, takes_logger=True, + use_proxy=False) + self.butler = None + + def getKwargs(self, config, base, logger): + logger.debug("Get kwargs for DeepCoadd") + req = { + "repo": str, + "collection": str, + "band": str + } + opt = { + "skymap": str, + "dstype": str + } + params, _ = galsim.config.GetAllParams(config, base, req=req, opt=opt) + if self.butler is None: + self.butler = daf_butler.Butler(params["repo"], + collections=[params["collection"]]) + + kwargs = { + "butler": self.butler, + "band": params["band"], + "skymap": params.get("skymap", "lsst_cells_v2"), + "dstype": params.get("dstype", "deep_coadd"), + } + safe = True + return kwargs, safe + + +def BuildRubinCoaddPSF(config, base, ignore, gsparams, logger): + """ + Build PSFs from Rubin deep_coadds. The deep_coadd image will be + retrieved from the data repository using the butler for the + tract-patch combination corresponding to the object's sky position. + + Assuming the input.deep_coadd object is defined, to use this + PSF, add the following to the config yaml: + + input.atm_psf: "" # disable the atmospheric PSF + psf: + type: RubinCoaddPSF + """ + deep_coadd = galsim.config.GetInputObj('deep_coadd', config, base, + 'RubinCoaddPSF') + image_pos = base['image_pos'] + celestial_coord = base['wcs'].toWorld(image_pos) + ra = celestial_coord.ra / galsim.degrees + dec = celestial_coord.dec / galsim.degrees + safe = False + return deep_coadd.getPSF(ra, dec), safe + + +RegisterInputType("deep_coadd", DeepCoaddLoader()) +RegisterObjectType("RubinCoaddPSF", BuildRubinCoaddPSF, input_type="deep_coadd") From c4cb7ab08254ffecd0de52e64a62be3f4bc4d001 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 8 Sep 2026 15:22:43 -0700 Subject: [PATCH 02/30] rename coadd_psf.py -> deep_coadd.py --- imsim/__init__.py | 2 +- imsim/{coadd_psf.py => deep_coadd.py} | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename imsim/{coadd_psf.py => deep_coadd.py} (100%) diff --git a/imsim/__init__.py b/imsim/__init__.py index 5022d4af..df5d52f8 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -41,4 +41,4 @@ from .table_row import * from .photon_pooling import * from .full_focal_plane import * -from .coadd_psf import * +from .deep_coadd import * diff --git a/imsim/coadd_psf.py b/imsim/deep_coadd.py similarity index 100% rename from imsim/coadd_psf.py rename to imsim/deep_coadd.py From 489317e2c5fe61fbf7ee799e1177c198cb6ca677 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 9 Sep 2026 12:27:21 -0700 Subject: [PATCH 03/30] provide info on deep_coadds more generally --- imsim/deep_coadd.py | 49 ++++++++++++++++++++++++++++++++++++++------- 1 file changed, 42 insertions(+), 7 deletions(-) diff --git a/imsim/deep_coadd.py b/imsim/deep_coadd.py index a599d389..61fc3982 100644 --- a/imsim/deep_coadd.py +++ b/imsim/deep_coadd.py @@ -1,7 +1,7 @@ """ -Interface to Rubin deep_coadd PSFs +Interface to Rubin deep_coadds """ - +import pandas as pd import galsim from galsim.config import InputLoader, RegisterInputType, RegisterObjectType import lsst.daf.butler as daf_butler @@ -9,7 +9,8 @@ class DeepCoadd: - def __init__(self, butler, band, skymap, dstype="deep_coadd"): + def __init__(self, butler, band, skymap_name, data_ids=None, + dstype="deep_coadd"): """ Parameters ---------- @@ -17,19 +18,43 @@ def __init__(self, butler, band, skymap, dstype="deep_coadd"): Butler for the data repo and collection containing the deep_coadds. band : str Band of interest, i.e., in "ugrizy". - skymap : str + skymap_name : str Name of the skymap to use, e.g., "lsst_cells_v2". + data_ids: list + List of data_ids (dicts) for identifying deep_coadds to simulate. dstype : str The dataset type of the cell-based coadds. Default: "deep_coadd". """ self.butler = butler self.band = band - self.skymap_name = skymap - self.skymap = butler.get("skyMap", name=skymap) + self.skymap_name = skymap_name + self.skymap = butler.get("skyMap", skymap=skymap_name) + self.data_ids = data_ids self.dstype = dstype self._psf_cache = {} self._butler_cache = {} + def get(self, index=None, data_id=None): + """Return the deep_coadd using the index of the self.data_ids list, + or if provided, using the data_id. + """ + if (data_id, index) == (None, None): + return None + if data_id is None and index >= 0 and index < len(self.data_ids): + data_id = self.data_ids[index] + elif data_id not in self.data_ids: + return None + data_id['skymap'] = self.skymap_name + deep_coadd = self.butler.get(self.dstype, **data_id) + + # Ensure this is a lsst.images version: + assert hasattr(deep_coadd, 'to_legacy') + return deep_coadd + + def getWcs(self, data_id): + data_id['skymap'] = self.skymap_name + return galsim.AstropyWCS(wcs=self.get(data_id=data_id).fits_wcs) + def getPSF(self, ra, dec): """Return the cell coadd PSF, evaluated at the center of the cell containing this sky position. @@ -99,6 +124,7 @@ def __init__(self): super().__init__(init_func=DeepCoadd, takes_logger=True, use_proxy=False) self.butler = None + self.deep_coadd_list = None def getKwargs(self, config, base, logger): logger.debug("Get kwargs for DeepCoadd") @@ -109,13 +135,22 @@ def getKwargs(self, config, base, logger): } opt = { "skymap": str, - "dstype": str + "dstype": str, + "deep_coadd_list_file": str } params, _ = galsim.config.GetAllParams(config, base, req=req, opt=opt) if self.butler is None: self.butler = daf_butler.Butler(params["repo"], collections=[params["collection"]]) + deep_coadd_list_file = params.get("deep_coadd_list_file", None) + + if self.deep_coadd_list is None and deep_coadd_list_file is not None: + df = pd.read_parquet("deep_coadd_list_file") + columns = ["band", "tract", "patch"] + self.deep_coadd_list = [dict(zip(columns, row)) + for row in zip(*[df[_] for _ in columns])] + kwargs = { "butler": self.butler, "band": params["band"], From e160fd23222184f9b1fc44d8c4653d3a3def7553 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 9 Sep 2026 12:28:24 -0700 Subject: [PATCH 04/30] initial implementation of coadd image builder --- imsim/__init__.py | 1 + imsim/lsst_coadd_image.py | 83 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+) create mode 100644 imsim/lsst_coadd_image.py diff --git a/imsim/__init__.py b/imsim/__init__.py index df5d52f8..be9e2dad 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -33,6 +33,7 @@ from .bandpass import * from .telescope_loader import * from .lsst_image import * +from .lsst_coadd_image import * from .checkpoint import * from .opd import * from .vignetting import * diff --git a/imsim/lsst_coadd_image.py b/imsim/lsst_coadd_image.py new file mode 100644 index 00000000..565bb47f --- /dev/null +++ b/imsim/lsst_coadd_image.py @@ -0,0 +1,83 @@ +import galsim +from galsim.config import GetAllParams, AddNoise, RegisterImageType +from .lsst_image import LSST_ImageBuilderBase + + +__all__ = ['LSST_CoaddImageBuilder'] + + +class LSST_CoaddImageBuilder(LSST_ImageBuilderBase): # noqa: N801 + def setup(self, config, base, image_num, obj_num, ignore, logger): + """ + Do the initialization and setup for a Rubin/LSST deep coadd image. + + Parameters: + config: The configuration dict for the image field. + base: The base configuration dict. + image_num: The current image number. + obj_num: The first object number in the image. + ignore: A list of parameters that are allowed to be in config + that we can ignore here. + logger: A logger object to log progress. + + Returns: + xsize, ysize + """ + self.nobjects = self.getNObj(config, base, image_num, logger=logger) + if 'nobjects' in config: + # User specified nobjects. + # Make sure it's not more than what any input catalog can + # handle (we don't want repeated objects). + input_nobj = galsim.config.ProcessInputNObjects(base) + if input_nobj is not None: + self.nobjects = min(self.nobjects, input_nobj) + logger.info('image %d: nobj = %d', image_num, self.nobjects) + + req = { 'tract': int, 'patch': int } + opt = { 'nbatch': int, 'nsubbatch': int, 'nbatch_fft': int, + 'npbatch_per_checkpoint': int} + extra_ignore = ['image_pos', 'world_pos', 'stamp_size', + 'stamp_xsize', 'stamp_ysize', 'nobjects' ] + params = GetAllParams(config, base, req=req, opt=opt, + ignore=ignore+extra_ignore)[0] + + self.tract = params['tract'] + self.patch = params['patch'] + + self.add_noise = True + + self.nbatch = params.get('nbatch', 10) + self.nbatch_per_checkpoint = params.get('nbatch_per_checkpoint', 1) + self.nsubbatch = params.get('nsubbatch', 50) + self.nbatch_fft = params.get('nbatch_fft', 1) + try: + self.checkpoint = galsim.config.GetInputObj( + 'checkpoint', config, base, 'LSST_CoaddImageBuilder') + except galsim.config.GalSimConfigError: + self.checkpoint = None + + deep_coadd = galsim.config.GetInputObj( + 'deep_coadd', config, base, 'LSSTC_CoaddImageBuilder') + bbox = deep_coadd.skymap[self.tract][self.patch].getOuterBBox() + xsize, ysize = bbox.width, bbox.height + + return xsize, ysize + + def addNoise(self, image, config, base, image_num, obj_num, current_var, + logger): + """Add image noise + + Parameters: + image: The image onto which to add the noise. + config: The configuration dict for the image field. + base: The base configuration dict. + image_num: The current image number. + obj_num: The first object number in the image. + current_var: The current noise variance in each postage stamps. + logger: A logger object to log progress. + """ + if self.add_noise: + AddNoise(base, image, current_var, logger) + + +RegisterImageType('LSST_DeepCoadd', LSST_CoaddImageBuilder()) From fffa0207e296d951791db4fa4b53c747f7305113 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 9 Sep 2026 12:28:46 -0700 Subject: [PATCH 05/30] initial implementation of coadd wcs builder --- imsim/__init__.py | 1 + imsim/deep_coadd_wcs.py | 35 +++++++++++++++++++++++++++++++++++ 2 files changed, 36 insertions(+) create mode 100644 imsim/deep_coadd_wcs.py diff --git a/imsim/__init__.py b/imsim/__init__.py index be9e2dad..4329ad2e 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -43,3 +43,4 @@ from .photon_pooling import * from .full_focal_plane import * from .deep_coadd import * +from .deep_coadd_wcs import * diff --git a/imsim/deep_coadd_wcs.py b/imsim/deep_coadd_wcs.py new file mode 100644 index 00000000..004eafae --- /dev/null +++ b/imsim/deep_coadd_wcs.py @@ -0,0 +1,35 @@ +""" +WCS builder for LSST deep_coadds +""" +from galsim.config import WCSBuilder, GetAllParams, GetInputObj, RegisterWCSType + + +__all__ = ["LSSTDeepCoaddWCSBuilder"] + + +class LSSTDeepCoaddWCSBuilder(WCSBuilder): + + def buildWCS(self, config, base, logger): + """Build an AstropyWCS from a deep_coadd + + Parameters: + config: The configuration dict for the wcs type. + base: The base configuration dict. + logger: If provided, a logger for logging debug statements. + + Returns: + The constructed WCS object (a galsim.AstropyWCS instance). + """ + req = { + "tract": int, + "patch": int, + "band": str, + } + params, _ = GetAllParams(config, base, req=req) + data_id = { key: params[key] for key in req } + deep_coadd = GetInputObj("deep_coadd", config, base, + "LSSTDeepCoaddWCSBuilder") + return deep_coadd.getWcs(data_id) + + +RegisterWCSType("LSSTDeepCoaddWcs", LSSTDeepCoaddWCSBuilder()) From ef0fae0e6a36e5f92e20f8437b6bfb0a0d425f4c Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 9 Sep 2026 17:59:55 -0700 Subject: [PATCH 06/30] add RubinDeepCoaddOutput type; adjust various configs --- imsim/__init__.py | 3 +- imsim/deep_coadd.py | 63 ++++++++++++------- ...sst_coadd_image.py => deep_coadd_image.py} | 13 ++-- imsim/deep_coadd_output.py | 36 +++++++++++ imsim/deep_coadd_wcs.py | 10 +-- imsim/skycat.py | 2 +- imsim/templates.py | 1 + 7 files changed, 91 insertions(+), 37 deletions(-) rename imsim/{lsst_coadd_image.py => deep_coadd_image.py} (89%) create mode 100644 imsim/deep_coadd_output.py diff --git a/imsim/__init__.py b/imsim/__init__.py index 4329ad2e..bfdb02a1 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -33,7 +33,6 @@ from .bandpass import * from .telescope_loader import * from .lsst_image import * -from .lsst_coadd_image import * from .checkpoint import * from .opd import * from .vignetting import * @@ -43,4 +42,6 @@ from .photon_pooling import * from .full_focal_plane import * from .deep_coadd import * +from .deep_coadd_image import * from .deep_coadd_wcs import * +from .deep_coadd_output import * diff --git a/imsim/deep_coadd.py b/imsim/deep_coadd.py index 61fc3982..3ebc5675 100644 --- a/imsim/deep_coadd.py +++ b/imsim/deep_coadd.py @@ -3,21 +3,19 @@ """ import pandas as pd import galsim -from galsim.config import InputLoader, RegisterInputType, RegisterObjectType +from galsim.config import (InputLoader, RegisterInputType, RegisterObjectType, + RegisterValueType, GetAllParams, GetInputObj) import lsst.daf.butler as daf_butler import lsst.geom class DeepCoadd: - def __init__(self, butler, band, skymap_name, data_ids=None, - dstype="deep_coadd"): + def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """ Parameters ---------- butler : lsst.daf.butler.Butler Butler for the data repo and collection containing the deep_coadds. - band : str - Band of interest, i.e., in "ugrizy". skymap_name : str Name of the skymap to use, e.g., "lsst_cells_v2". data_ids: list @@ -26,7 +24,6 @@ def __init__(self, butler, band, skymap_name, data_ids=None, The dataset type of the cell-based coadds. Default: "deep_coadd". """ self.butler = butler - self.band = band self.skymap_name = skymap_name self.skymap = butler.get("skyMap", skymap=skymap_name) self.data_ids = data_ids @@ -52,10 +49,10 @@ def get(self, index=None, data_id=None): return deep_coadd def getWcs(self, data_id): - data_id['skymap'] = self.skymap_name - return galsim.AstropyWCS(wcs=self.get(data_id=data_id).fits_wcs) + deep_coadd = self.get(data_id=data_id) + return galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) - def getPSF(self, ra, dec): + def getPSF(self, ra, dec, band): """Return the cell coadd PSF, evaluated at the center of the cell containing this sky position. """ @@ -72,7 +69,7 @@ def getPSF(self, ra, dec): butler_key = tract, patch if butler_key not in self._butler_cache: dataId = dict(skymap=self.skymap_name, tract=tract, patch=patch, - band=self.band) + band=band) deep_coadd = self.butler.get(self.dstype, **dataId) if hasattr(deep_coadd, "to_legacy"): # This is a `lsst.images` object. Use the `.to_legacy` @@ -117,8 +114,7 @@ class DeepCoaddLoader(InputLoader): input.deep_coadd: repo: dp2 - collection: LSSTCam/runs/DRP/DP2 - band: $band # This will be obtained from the opsim metadata entry. + collection: dp2 """ def __init__(self): super().__init__(init_func=DeepCoadd, takes_logger=True, @@ -131,14 +127,13 @@ def getKwargs(self, config, base, logger): req = { "repo": str, "collection": str, - "band": str } opt = { - "skymap": str, + "skymap_name": str, "dstype": str, - "deep_coadd_list_file": str + "deep_coadd_list_file": str, } - params, _ = galsim.config.GetAllParams(config, base, req=req, opt=opt) + params, _ = GetAllParams(config, base, req=req, opt=opt) if self.butler is None: self.butler = daf_butler.Butler(params["repo"], collections=[params["collection"]]) @@ -146,15 +141,15 @@ def getKwargs(self, config, base, logger): deep_coadd_list_file = params.get("deep_coadd_list_file", None) if self.deep_coadd_list is None and deep_coadd_list_file is not None: - df = pd.read_parquet("deep_coadd_list_file") + df = pd.read_parquet(deep_coadd_list_file) columns = ["band", "tract", "patch"] self.deep_coadd_list = [dict(zip(columns, row)) for row in zip(*[df[_] for _ in columns])] kwargs = { "butler": self.butler, - "band": params["band"], - "skymap": params.get("skymap", "lsst_cells_v2"), + "skymap_name": params.get("skymap", "lsst_cells_v2"), + "data_ids": self.deep_coadd_list, "dstype": params.get("dstype", "deep_coadd"), } safe = True @@ -174,15 +169,35 @@ def BuildRubinCoaddPSF(config, base, ignore, gsparams, logger): psf: type: RubinCoaddPSF """ - deep_coadd = galsim.config.GetInputObj('deep_coadd', config, base, - 'RubinCoaddPSF') + deep_coadd = GetInputObj('deep_coadd', config, base, 'RubinCoaddPSF') + coadd_num = base['coadd_num'] + assert (coadd_num >= 0 and coadd_num < len(deep_coadd.data_ids)) + data_id = deep_coadd.data_ids[coadd_num] image_pos = base['image_pos'] celestial_coord = base['wcs'].toWorld(image_pos) ra = celestial_coord.ra / galsim.degrees dec = celestial_coord.dec / galsim.degrees + band = data_id['band'] safe = False - return deep_coadd.getPSF(ra, dec), safe + return deep_coadd.getPSF(ra, dec, band), safe + + +def DeepCoaddData(config, base, value_type): + deep_coadd = GetInputObj('deep_coadd', config, base, 'DeepCoaddData') + req = { 'field': str } + params, safe = GetAllParams(config, base, req=req) + field = params['field'] + + coadd_num = base['coadd_num'] + assert (coadd_num >= 0 and coadd_num < len(deep_coadd.data_ids)) + data_id = deep_coadd.data_ids[coadd_num] + + val = value_type(data_id.get(field, None)) + + return val, safe -RegisterInputType("deep_coadd", DeepCoaddLoader()) -RegisterObjectType("RubinCoaddPSF", BuildRubinCoaddPSF, input_type="deep_coadd") +RegisterInputType('deep_coadd', DeepCoaddLoader()) +RegisterObjectType('RubinCoaddPSF', BuildRubinCoaddPSF, input_type='deep_coadd') +RegisterValueType('DeepCoaddData', DeepCoaddData, [int, str], + input_type='deep_coadd') diff --git a/imsim/lsst_coadd_image.py b/imsim/deep_coadd_image.py similarity index 89% rename from imsim/lsst_coadd_image.py rename to imsim/deep_coadd_image.py index 565bb47f..2f505a29 100644 --- a/imsim/lsst_coadd_image.py +++ b/imsim/deep_coadd_image.py @@ -3,10 +3,10 @@ from .lsst_image import LSST_ImageBuilderBase -__all__ = ['LSST_CoaddImageBuilder'] +__all__ = ['RubinCoaddImageBuilder'] -class LSST_CoaddImageBuilder(LSST_ImageBuilderBase): # noqa: N801 +class RubinCoaddImageBuilder(LSST_ImageBuilderBase): # noqa: N801 def setup(self, config, base, image_num, obj_num, ignore, logger): """ Do the initialization and setup for a Rubin/LSST deep coadd image. @@ -33,9 +33,9 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.nobjects = min(self.nobjects, input_nobj) logger.info('image %d: nobj = %d', image_num, self.nobjects) - req = { 'tract': int, 'patch': int } + req = { 'tract': int, 'patch': int, 'band': str } opt = { 'nbatch': int, 'nsubbatch': int, 'nbatch_fft': int, - 'npbatch_per_checkpoint': int} + 'nbatch_per_checkpoint': int} extra_ignore = ['image_pos', 'world_pos', 'stamp_size', 'stamp_xsize', 'stamp_ysize', 'nobjects' ] params = GetAllParams(config, base, req=req, opt=opt, @@ -43,6 +43,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.tract = params['tract'] self.patch = params['patch'] + self.band = params['band'] self.add_noise = True @@ -52,7 +53,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.nbatch_fft = params.get('nbatch_fft', 1) try: self.checkpoint = galsim.config.GetInputObj( - 'checkpoint', config, base, 'LSST_CoaddImageBuilder') + 'checkpoint', config, base, 'RubinCoaddImageBuilder') except galsim.config.GalSimConfigError: self.checkpoint = None @@ -80,4 +81,4 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, AddNoise(base, image, current_var, logger) -RegisterImageType('LSST_DeepCoadd', LSST_CoaddImageBuilder()) +RegisterImageType('RubinDeepCoadd', RubinCoaddImageBuilder()) diff --git a/imsim/deep_coadd_output.py b/imsim/deep_coadd_output.py new file mode 100644 index 00000000..6c5a2a03 --- /dev/null +++ b/imsim/deep_coadd_output.py @@ -0,0 +1,36 @@ +from galsim.config import (OutputBuilder, RegisterOutputType, ParseValue, + GetAllParams, GetInputObj, BuildImage, + eval_base_variables) + +__all__ = ['RubinDeepCoaddOutputBuilder'] + + +class RubinDeepCoaddOutputBuilder(OutputBuilder): + + _added_eval_base_variables = False + + def setup(self, config, base, file_num, logger): + coadd_num = ParseValue(config, 'coadd_num', base, int)[0] + if not self._added_eval_base_variables: + eval_base_variables.append('coadd_num') + self._added_eval_base_variables = True + base['coadd_num'] = coadd_num + + def getNFiles(self, config, base, logger=None): + if 'nfiles' in config: + return ParseValue(config, 'nfiles', base, int)[0] + else: + deep_coadds = GetInputObj('deep_coadd', config, base, + "RubinDeepCoaddOutputBuilder") + return len(deep_coadds.data_ids) + + def buildImages(self, config, base, file_num, image_num, obj_num, + ignore, logger): + ignore += [ 'file_name', 'dir', 'nfiles', 'coadd_num' ] + params, _ = GetAllParams(config, base, ignore=ignore) + + image = BuildImage(base, image_num, obj_num, logger=logger) + return [ image ] + + +RegisterOutputType('RubinDeepCoaddOutput', RubinDeepCoaddOutputBuilder()) diff --git a/imsim/deep_coadd_wcs.py b/imsim/deep_coadd_wcs.py index 004eafae..dee01aed 100644 --- a/imsim/deep_coadd_wcs.py +++ b/imsim/deep_coadd_wcs.py @@ -1,13 +1,13 @@ """ -WCS builder for LSST deep_coadds +WCS builder for Rubin deep_coadds """ from galsim.config import WCSBuilder, GetAllParams, GetInputObj, RegisterWCSType -__all__ = ["LSSTDeepCoaddWCSBuilder"] +__all__ = ["RubinDeepCoaddWCSBuilder"] -class LSSTDeepCoaddWCSBuilder(WCSBuilder): +class RubinDeepCoaddWCSBuilder(WCSBuilder): def buildWCS(self, config, base, logger): """Build an AstropyWCS from a deep_coadd @@ -28,8 +28,8 @@ def buildWCS(self, config, base, logger): params, _ = GetAllParams(config, base, req=req) data_id = { key: params[key] for key in req } deep_coadd = GetInputObj("deep_coadd", config, base, - "LSSTDeepCoaddWCSBuilder") + "RubinDeepCoaddWCSBuilder") return deep_coadd.getWcs(data_id) -RegisterWCSType("LSSTDeepCoaddWcs", LSSTDeepCoaddWCSBuilder()) +RegisterWCSType("RubinDeepCoaddWcs", RubinDeepCoaddWCSBuilder()) diff --git a/imsim/skycat.py b/imsim/skycat.py index dbc2931c..670a8a2f 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -15,7 +15,7 @@ class SkyCatalogInterface: """Interface to skyCatalogs package.""" - def __init__(self, file_name, wcs, band, mjd, fiducial_wcs=None, + def __init__(self, file_name, wcs, band, mjd=None, fiducial_wcs=None, fiducial_det_name=None, xsize=4096, ysize=4096, obj_types=None, skycatalog_root=None, edge_pix=200, camera=None, det_name=None, diff --git a/imsim/templates.py b/imsim/templates.py index 9cec8731..71d865fb 100644 --- a/imsim/templates.py +++ b/imsim/templates.py @@ -15,3 +15,4 @@ galsim.config.RegisterTemplate('imsim-config-photon-pooling', os.path.join(config_dir, 'imsim-config-photon-pooling.yaml')) galsim.config.RegisterTemplate('imsim-config-instcat-comcam', os.path.join(config_dir, 'imsim-config-instcat-comcam.yaml')) galsim.config.RegisterTemplate('imsim-config-skycat-comcam', os.path.join(config_dir, 'imsim-config-skycat-comcam.yaml')) +galsim.config.RegisterTemplate('imsim-config-deep-coadd', os.path.join(config_dir, 'imsim-config-deep-coadd.yaml')) From ac0a3da69714b69fd411f206cc60778686e3a446 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Wed, 9 Sep 2026 18:00:23 -0700 Subject: [PATCH 07/30] add deep-coadd template yaml --- config/imsim-config-deep-coadd.yaml | 57 +++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 config/imsim-config-deep-coadd.yaml diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml new file mode 100644 index 00000000..13b54142 --- /dev/null +++ b/config/imsim-config-deep-coadd.yaml @@ -0,0 +1,57 @@ +modules: + - imsim + +eval_variables: + itract: { type: DeepCoaddData, field: tract } + ipatch: { type: DeepCoaddData, field: patch } + sband: { type: DeepCoaddData, field: band } + +input: + sky_catalog: + file_name: default_sky_cat.yaml + band: $band + + checkpoint: + dir: checkpoint + file_name: + type: FormattedStr + format : checkpoint_%08d-%s-%s.hdf + items: + - $tract + - $patch + - $band + +image: + type: RubinDeepCoadd + random_seed: 42 + nproc: 1 + nbatch: 100 + nsubbatch: 50 + nbatch_fft: 1 + nbatch_per_checkpoint: 1 + band: $band + tract: $tract + patch: $patch + + wcs: + type: RubinDeepCoaddWcs + tract: $tract + patch: $patch + band: $band + +psf: + type: RubinCoaddPSF + +output: + type: RubinDeepCoaddOutput + coadd_num: 0 + nproc: 1 + nfiles: 1 + dir: fits + file_name: + type: FormattedStr + format: deep_coadd_%s_%d_%d.fits + items: + - $band + - $tract + - $patch From fc753cf3f17e7064bdab8b27f342385b7d420d33 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Sat, 12 Sep 2026 21:55:11 -0700 Subject: [PATCH 08/30] initial stamp and noise code --- config/imsim-config-deep-coadd.yaml | 18 ++++++++-- imsim/__init__.py | 1 + imsim/deep_coadd_image.py | 52 +++++++++++++++++++++++++++-- imsim/deep_coadd_output.py | 13 +++++++- imsim/deep_coadd_stamp.py | 29 ++++++++++++++++ imsim/lsst_image.py | 2 -- 6 files changed, 106 insertions(+), 9 deletions(-) create mode 100644 imsim/deep_coadd_stamp.py diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 13b54142..a6a979fa 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -15,11 +15,19 @@ input: dir: checkpoint file_name: type: FormattedStr - format : checkpoint_%08d-%s-%s.hdf + format : checkpoint_%s-%d-%d.hdf items: + - $band - $tract - $patch - - $band + +gal: + type: SkyCatObj + +stamp: + type: RubinDeepCoadd + world_pos: + type: SkyCatWorldPos image: type: RubinDeepCoadd @@ -35,8 +43,12 @@ image: wcs: type: RubinDeepCoaddWcs + band: $band tract: $tract patch: $patch + + bandpass: + type: RubinBandpass band: $band psf: @@ -50,7 +62,7 @@ output: dir: fits file_name: type: FormattedStr - format: deep_coadd_%s_%d_%d.fits + format: deep_coadd_%s-%d-%d.fits items: - $band - $tract diff --git a/imsim/__init__.py b/imsim/__init__.py index bfdb02a1..ec7924e2 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -45,3 +45,4 @@ from .deep_coadd_image import * from .deep_coadd_wcs import * from .deep_coadd_output import * +from .deep_coadd_stamp import * diff --git a/imsim/deep_coadd_image.py b/imsim/deep_coadd_image.py index 2f505a29..5b2196c5 100644 --- a/imsim/deep_coadd_image.py +++ b/imsim/deep_coadd_image.py @@ -1,5 +1,8 @@ +from itertools import pairwise +import numpy as np import galsim -from galsim.config import GetAllParams, AddNoise, RegisterImageType +from galsim.config import (GetAllParams, AddNoise, RegisterImageType, + GetInputObj) from .lsst_image import LSST_ImageBuilderBase @@ -52,16 +55,22 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.nsubbatch = params.get('nsubbatch', 50) self.nbatch_fft = params.get('nbatch_fft', 1) try: - self.checkpoint = galsim.config.GetInputObj( + self.checkpoint = GetInputObj( 'checkpoint', config, base, 'RubinCoaddImageBuilder') except galsim.config.GalSimConfigError: self.checkpoint = None - deep_coadd = galsim.config.GetInputObj( + deep_coadd = GetInputObj( 'deep_coadd', config, base, 'LSSTC_CoaddImageBuilder') bbox = deep_coadd.skymap[self.tract][self.patch].getOuterBBox() xsize, ysize = bbox.width, bbox.height + data_id = {'band': self.band, + 'tract': self.tract, + 'patch': self.patch} + self.noise_level = np.std( + deep_coadd.get(data_id=data_id).noise_realizations[0].array) + return xsize, ysize def addNoise(self, image, config, base, image_num, obj_num, current_var, @@ -77,8 +86,45 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, current_var: The current noise variance in each postage stamps. logger: A logger object to log progress. """ + if 'noise' not in base['image']: + base['image']['noise'] = { "type": "Gaussian", + "sigma": self.noise_level } if self.add_noise: AddNoise(base, image, current_var, logger) + def buildImage(self, config, base, image_num, obj_num, logger): + + all_stamps = [] + all_vars = [] + full_image = self._create_full_image(config, base) + min_batch_size = 10 + if self.nobjects//self.nbatch < min_batch_size: + num_batches = 1 + else: + num_batches = self.nbatch + bounds = np.linspace(0, self.nobjects, num_batches + 1, dtype=int) + for batch, (start_obj_num, end_obj_num) in enumerate(pairwise(bounds)): + nobj_batch = end_obj_num - start_obj_num + stamps, current_vars = galsim.config.BuildStamps( + nobj_batch, base, logger=logger, obj_num=start_obj_num, + do_noise=False) + base['index_key'] = 'image_num' + + for k in range(nobj_batch): + if stamps[k] is None: + continue + bounds = stamps[k].bounds & full_image.bounds + if not bounds.isDefined(): + continue + full_image[bounds] += stamps[k][bounds] + + nz_var = np.nonzero(current_vars)[0] + all_stamps.extend([stamps[k] for k in nz_var]) + all_vars.extend([current_vars[k] for k in nz_var]) + + current_var = galsim.config.FlattenNoiseVariance( + base, full_image, all_stamps, all_vars, logger) + return full_image, current_var + RegisterImageType('RubinDeepCoadd', RubinCoaddImageBuilder()) diff --git a/imsim/deep_coadd_output.py b/imsim/deep_coadd_output.py index 6c5a2a03..a32029ab 100644 --- a/imsim/deep_coadd_output.py +++ b/imsim/deep_coadd_output.py @@ -1,7 +1,9 @@ +import galsim from galsim.config import (OutputBuilder, RegisterOutputType, ParseValue, GetAllParams, GetInputObj, BuildImage, eval_base_variables) + __all__ = ['RubinDeepCoaddOutputBuilder'] @@ -15,13 +17,20 @@ def setup(self, config, base, file_num, logger): eval_base_variables.append('coadd_num') self._added_eval_base_variables = True base['coadd_num'] = coadd_num + deep_coadds = GetInputObj('deep_coadd', config, base, + 'RubinDeepCoaddOutputBuilder') + data_id = deep_coadds.data_ids[coadd_num] + patch = deep_coadds.skymap[data_id['tract']][data_id['patch']] + bbox = patch.getOuterBBox() + base['det_xsize'] = bbox.width + base['det_ysize'] = bbox.height def getNFiles(self, config, base, logger=None): if 'nfiles' in config: return ParseValue(config, 'nfiles', base, int)[0] else: deep_coadds = GetInputObj('deep_coadd', config, base, - "RubinDeepCoaddOutputBuilder") + 'RubinDeepCoaddOutputBuilder') return len(deep_coadds.data_ids) def buildImages(self, config, base, file_num, image_num, obj_num, @@ -30,6 +39,8 @@ def buildImages(self, config, base, file_num, image_num, obj_num, params, _ = GetAllParams(config, base, ignore=ignore) image = BuildImage(base, image_num, obj_num, logger=logger) + image.header = galsim.FitsHeader() + image.wcs.header = {} return [ image ] diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py new file mode 100644 index 00000000..7ff8a2bd --- /dev/null +++ b/imsim/deep_coadd_stamp.py @@ -0,0 +1,29 @@ +import galsim +from galsim.config import (ParseValue, ParseWorldPos, + StampBuilder, RegisterStampType, GetAllParams) + + +__all__ = ["RubinDeepCoaddStampBuilder"] + + +class RubinDeepCoaddStampBuilder(StampBuilder): + + def setup(self, config, base, xsize, ysize, ignore, logger): + + xsize, ysize, image_pos, world_pos \ + = super().setup(config, base, xsize, ysize, ignore, logger) + return xsize, ysize, image_pos, world_pos + +# if 'image_pos' in config: +# image_pos = ParseValue(config, 'image_pos', base, +# galsim.PositionD)[0] +# else: +# image_pos = None +# +# if 'world_pos' in config: +# world_pos = ParseWorldPos(config, 'world_pos', base, logger) +# else: +# world_pos = None + + +RegisterStampType('RubinDeepCoadd', RubinDeepCoaddStampBuilder()) diff --git a/imsim/lsst_image.py b/imsim/lsst_image.py index b936f266..43c0393d 100644 --- a/imsim/lsst_image.py +++ b/imsim/lsst_image.py @@ -243,7 +243,6 @@ def _set_config_image_pos(config, base): config: The configuration dictionary for the image field. base: The base configuration dictionary. """ - if 'image_pos' in config and 'world_pos' in config: raise galsim.config.GalSimConfigValueError( "Both image_pos and world_pos specified for LSST_Image.", @@ -267,7 +266,6 @@ def _set_config_image_pos(config, base): class LSST_ImageBuilder(LSST_ImageBuilderBase): """This is based on the GalSim "Scattered" image type. - Additional features supported here: * batching the stamp generation * checkpointing if input.checkpoint is given From 8d15b6c45db7b613a25fffd35338845eed1cdc43 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 14:24:24 -0700 Subject: [PATCH 09/30] add RubinDeepCoaddBandpass --- config/imsim-config-deep-coadd.yaml | 2 +- imsim/__init__.py | 1 + imsim/deep_coadd.py | 7 +++++++ imsim/deep_coadd_bandpass.py | 27 +++++++++++++++++++++++++++ 4 files changed, 36 insertions(+), 1 deletion(-) create mode 100644 imsim/deep_coadd_bandpass.py diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index a6a979fa..f42cdd1b 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -48,7 +48,7 @@ image: patch: $patch bandpass: - type: RubinBandpass + type: RubinDeepCoaddBandpass band: $band psf: diff --git a/imsim/__init__.py b/imsim/__init__.py index ec7924e2..916676f0 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -46,3 +46,4 @@ from .deep_coadd_wcs import * from .deep_coadd_output import * from .deep_coadd_stamp import * +from .deep_coadd_bandpass import * diff --git a/imsim/deep_coadd.py b/imsim/deep_coadd.py index 3ebc5675..6771dd25 100644 --- a/imsim/deep_coadd.py +++ b/imsim/deep_coadd.py @@ -48,6 +48,13 @@ def get(self, index=None, data_id=None): assert hasattr(deep_coadd, 'to_legacy') return deep_coadd + def getBandpass(self, band): + passband = self.butler.get('standard_passband', band=band) + lut = galsim.LookupTable(passband['wavelength'], + passband['throughput'], + interpolant='linear') + return galsim.Bandpass(lut, wave_type='nm').thin() + def getWcs(self, data_id): deep_coadd = self.get(data_id=data_id) return galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) diff --git a/imsim/deep_coadd_bandpass.py b/imsim/deep_coadd_bandpass.py new file mode 100644 index 00000000..c3df43f5 --- /dev/null +++ b/imsim/deep_coadd_bandpass.py @@ -0,0 +1,27 @@ +import galsim +from galsim.config import BandpassBuilder, RegisterBandpassType + + +__all__ = ['RubinDeepCoaddBandpassBuilder'] + + +class RubinDeepCoaddBandpassBuilder(BandpassBuilder): + """Build the Rubin bandpass to be used with deep coadd sims + from standard passbands.""" + def buildBandpass(self, config, base, logger): + """ + Parameters: + config: The configuration dict for the bandpass type. + base: The base configuration dict. + logger: If provided, a logger for logging debug statements. + Returns: + the constructed Bandpass object. + """ + req = { 'band': str } + params, safe = galsim.config.GetAllParams(config, base, req=req) + deep_coadds = galsim.config.GetInputObj('deep_coadd', config, base, + 'DeepCoaddBandpass') + return deep_coadds.getBandpass(params['band']), safe + + +RegisterBandpassType('RubinDeepCoaddBandpass', RubinDeepCoaddBandpassBuilder()) From d95c8172f0686d77054e0de1090afc00ad3f58bb Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 14:25:27 -0700 Subject: [PATCH 10/30] implement .draw method for RubinDeepCoaddStampBuilder --- imsim/deep_coadd_stamp.py | 167 +++++++++++++++++++++++++++++++++++--- 1 file changed, 154 insertions(+), 13 deletions(-) diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 7ff8a2bd..18686d0b 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -1,29 +1,170 @@ +import numpy as np import galsim -from galsim.config import (ParseValue, ParseWorldPos, - StampBuilder, RegisterStampType, GetAllParams) +from galsim.config import StampBuilder, RegisterStampType +from .stamp_utils import get_stamp_size __all__ = ["RubinDeepCoaddStampBuilder"] class RubinDeepCoaddStampBuilder(StampBuilder): + _pixel_scale = 0.2 + _tiny_flux = 10.0 + _Nmax = 4096 def setup(self, config, base, xsize, ysize, ignore, logger): - + # Use base class setup to find default xsize, ysize (and object + # position). xsize, ysize, image_pos, world_pos \ = super().setup(config, base, xsize, ysize, ignore, logger) + + obj = galsim.config.BuildGSObject(base, 'gal', logger=logger)[0] + if obj is None: + raise galsim.config.SkipThisObject( + 'gal is None (invalid parameters)') + self.obj = obj + + self.rng = galsim.config.GetRNG(config, base, logger, "RubinDeepCoadd") + self.image = base['current_image'] + bandpass = base['bandpass'] + if not hasattr(obj, 'flux'): + obj.flux = obj.calculateFlux(bandpass) + self.nominal_flux = obj.flux + self.phot_flux = galsim.PoissonDeviate(self.rng, mean=obj.flux)() + + # Save values for output truth catalog + base['nominal_flux'] = self.nominal_flux + base['phot_flux'] = self.phot_flux + base['realized_flux'] = 0 # This will be updated by .drawImage + + if self.phot_flux == 0: + raise galsim.config.SkipThisObject('phot_flux==0') + + if xsize > 0 and ysize > 0: + # Use the stamp size from the config. + pass + elif 'size' in config: + # Get the stamp size from the size config entry + xsize = ysize = galsim.config.ParseValue( + config, 'size', base, int)[0] + elif self.nominal_flux < self._tiny_flux: + xsize = ysize = 32 + else: + # Determine the stamp size from the object flux. + base['current_noise_image'] = base['current_image'] + noise_var = galsim.config.CalculateNoiseVariance(base) + obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) + stamp_size = get_stamp_size( + obj_achrom=obj_achrom, + nominal_flux=self.nominal_flux, + noise_var=noise_var, + Nmax=self._Nmax, + pixel_scale=self._pixel_scale, + logger=logger + ) + xsize = ysize = stamp_size + + logger.info('Object %d will use stamp size = %s, %s and nominal flux %s', + base.get('obj_num',0), xsize, ysize, self.nominal_flux) return xsize, ysize, image_pos, world_pos -# if 'image_pos' in config: -# image_pos = ParseValue(config, 'image_pos', base, -# galsim.PositionD)[0] -# else: -# image_pos = None -# -# if 'world_pos' in config: -# world_pos = ParseWorldPos(config, 'world_pos', base, logger) -# else: -# world_pos = None + def draw(self, prof, image, method, offset, config, base, logger): + if prof is None: + image.photons = galsim.PhotonArray(0) + return image + + # Prof is normally a convolution here with obj_list being + # [gal, psf1, psf2,...] for some number of component PSFs. + gal, *psfs = prof.obj_list if hasattr(prof, 'obj_list') else [prof] + obj_num = base.get('obj_num',0) + bandpass = base['bandpass'] + + max_flux_simple = config.get('max_flux_simple', 100) + faint = self.nominal_flux < max_flux_simple + + if faint: + logger.info("Flux = %.0f Using trivial sed", self.obj.flux) + for profile_wl in (bandpass.effective_wavelength, + bandpass.red_limit, + bandpass.blue_limit): + sed_value = gal.sed(profile_wl) + if sed_value != 0: + break + if sed_value == 0: + # We can't evaluate the profile for this object, so skip it. + obj_num = base.get('obj_num') + object_id = base.get('object_id') + logger.warning("Zero-valued SED for faint object %d, " + "object_id %s. Skipping.", obj_num, object_id) + return image + gal = gal.evaluateAtWavelength(profile_wl) + gal = gal * self._trivial_sed + else: + self._fix_seds(gal, bandpass, logger) + + image.wcs = base['wcs'] + + gal = gal.withFlux(self.phot_flux, bandpass) + + if not faint and 'photon_ops' in config: + photon_ops = galsim.config.BuildPhotonOps( + config, 'photon_ops', base, logger) + else: + photon_ops = [] + + # Put the psfs at the start of the photon_ops. + # Probably a little better to put them a bit later than the + # start in some cases (e.g. after TimeSampler, + # PupilAnnulusSampler), but leave that as a todo for now. + photon_ops = psfs + photon_ops + sensor = None + image = gal.drawImage(bandpass, + method='phot', + offset=offset, + rng=self.rng, + n_photons=self.phot_flux, + image=image, + sensor=None, + photon_ops=photon_ops, + add_to_image=True, + poisson_flux=False) + base['realized_flux'] = image.added_flux + logger.debug('After .drawImage(...), image.added_flux %s', image.added_flux) + return image + + @classmethod + def _fix_seds(cls, prof, bandpass, logger): + # If any SEDs are not currently using a LookupTable for the + # function or if they are using spline interpolation, then the + # codepath is quite slow. Better to fix them before doing + # WavelengthSampler. + if (isinstance(prof, galsim.SimpleChromaticTransformation) and + (not isinstance(prof._flux_ratio._spec, galsim.LookupTable) + or prof._flux_ratio._spec.interpolant != 'linear')): + if not cls._sed_logged: + logger.warning( + "Warning: Chromatic drawing is most efficient when SEDs have " + "interpolant='linear'. Switching LookupTables to use 'linear'." + ) + cls._sed_logged = True + sed = prof._flux_ratio + wave_list, _, _ = galsim.utilities.combine_wave_list(sed, bandpass) + f = np.broadcast_to(sed(wave_list), wave_list.shape) + new_spec = galsim.LookupTable(wave_list, f, interpolant='linear') + new_sed = galsim.SED( + new_spec, + 'nm', + 'fphotons' if sed.spectral else '1' + ) + prof._flux_ratio = new_sed + + # Also recurse onto any components. + if isinstance(prof, galsim.ChromaticObject): + if hasattr(prof, 'obj_list'): + for obj in prof.obj_list: + cls._fix_seds(obj, bandpass, logger) + if hasattr(prof, 'original'): + cls._fix_seds(prof.original, bandpass, logger) RegisterStampType('RubinDeepCoadd', RubinDeepCoaddStampBuilder()) From bcd33c2350bdf7d20137e60cccb1a2f7af489e60 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 16:53:51 -0700 Subject: [PATCH 11/30] use LSST_SiliconBuilder as base class; enable fft rendering --- config/imsim-config-deep-coadd.yaml | 2 + imsim/deep_coadd_stamp.py | 167 ++++++++++++++++------------ 2 files changed, 98 insertions(+), 71 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index f42cdd1b..1df19422 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -26,6 +26,8 @@ gal: stamp: type: RubinDeepCoadd + fft_sb_thresh: 2.0e5 # When to switch to fft rendering. + max_flux_simple: 100 # When to switch to simple SED. world_pos: type: SkyCatWorldPos diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 18686d0b..0d5a8adc 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -1,27 +1,23 @@ -import numpy as np import galsim from galsim.config import StampBuilder, RegisterStampType +from .stamp import LSST_SiliconBuilder from .stamp_utils import get_stamp_size __all__ = ["RubinDeepCoaddStampBuilder"] -class RubinDeepCoaddStampBuilder(StampBuilder): - _pixel_scale = 0.2 - _tiny_flux = 10.0 - _Nmax = 4096 +class RubinDeepCoaddStampBuilder(LSST_SiliconBuilder): def setup(self, config, base, xsize, ysize, ignore, logger): - # Use base class setup to find default xsize, ysize (and object - # position). - xsize, ysize, image_pos, world_pos \ - = super().setup(config, base, xsize, ysize, ignore, logger) + # Use StampBuilder.setup(...) to set the object position. + ignore = ignore + ['fft_sb_thresh', 'max_flux_simple'] + _, _, image_pos, world_pos = StampBuilder.setup( + self, config, base, xsize, ysize, ignore, logger) obj = galsim.config.BuildGSObject(base, 'gal', logger=logger)[0] if obj is None: - raise galsim.config.SkipThisObject( - 'gal is None (invalid parameters)') + raise galsim.config.SkipThisObject('obj is None') self.obj = obj self.rng = galsim.config.GetRNG(config, base, logger, "RubinDeepCoadd") @@ -64,10 +60,25 @@ def setup(self, config, base, xsize, ysize, ignore, logger): ) xsize = ysize = stamp_size - logger.info('Object %d will use stamp size = %s, %s and nominal flux %s', + logger.info('Object %d will use stamp size %s, %s and nominal flux %s', base.get('obj_num',0), xsize, ysize, self.nominal_flux) return xsize, ysize, image_pos, world_pos + def buildPSF(self, config, base, gsparams, logger): + psf = galsim.config.BuildGSObject( + base, 'psf', gsparams=gsparams, logger=logger)[0] + + # Check if fft rendering should be used. + if 'fft_sb_thresh' in config: + fft_sb_thresh = galsim.config.ParseValue( + config, 'fft_sb_thresh', base, float)[0] + else: + fft_sb_thresh = None + + self.use_fft = (fft_sb_thresh is not None + and self.nominal_flux > fft_sb_thresh) + return psf + def draw(self, prof, image, method, offset, config, base, logger): if prof is None: image.photons = galsim.PhotonArray(0) @@ -104,67 +115,81 @@ def draw(self, prof, image, method, offset, config, base, logger): image.wcs = base['wcs'] - gal = gal.withFlux(self.phot_flux, bandpass) - - if not faint and 'photon_ops' in config: - photon_ops = galsim.config.BuildPhotonOps( - config, 'photon_ops', base, logger) + if method == 'fft': + gal = gal.withFlux(self.nominal_flux, bandpass) + fft_image = image.copy() + fft_offset = offset + kwargs = dict( + method='fft', + offset=fft_offset, + image=fft_image + ) + if not faint and config.get('fft_photon_ops'): + fft_photon_ops = galsim.config.BuildPhotonOps( + config, 'fft_photon_ops', base, logger) + kwargs.update({ + "photon_ops": fft_photon_ops, + "rng": self.rng, + "n_subsample": 1, + }) + + # Go back to a combined convolution for fft drawing. + prof = galsim.Convolve([gal] + psfs) + try: + fft_image = prof.drawImage(bandpass, **kwargs) + except galsim.errors.GalSimFFTSizeError as e: + # I think this shouldn't happen with the updates I + # made to how the image size is calculated, even for + # extremely bright things. So it should be ok to just + # report what happened, give some extra information to + # diagonose the problem and raise the error. + logger.error('Caught error trying to draw using FFT:') + logger.error('%s', e) + logger.error('You may need to add a gsparams field with ' + 'maximum_fft_size to either') + logger.error('the psf or gal field to allow larger FFTs.') + logger.info('prof = %r', prof) + logger.info('fft_image = %s', fft_image) + logger.info('offset = %r', offset) + raise + # Some pixels can end up negative from FFT numerics. Just + # set them to 0. + fft_image.array[fft_image.array < 0] = 0. + fft_image.addNoise(galsim.PoissonNoise(rng=self.rng)) + # In case we had to make a bigger image, just copy the + # part we need. + image += fft_image[image.bounds] + base['realized_flux'] = fft_image.added_flux + logger.debug('After .drawImage(...), fft_image.added_flux %s', + fft_image.added_flux) else: - photon_ops = [] - - # Put the psfs at the start of the photon_ops. - # Probably a little better to put them a bit later than the - # start in some cases (e.g. after TimeSampler, - # PupilAnnulusSampler), but leave that as a todo for now. - photon_ops = psfs + photon_ops - sensor = None - image = gal.drawImage(bandpass, - method='phot', - offset=offset, - rng=self.rng, - n_photons=self.phot_flux, - image=image, - sensor=None, - photon_ops=photon_ops, - add_to_image=True, - poisson_flux=False) - base['realized_flux'] = image.added_flux - logger.debug('After .drawImage(...), image.added_flux %s', image.added_flux) - return image + gal = gal.withFlux(self.phot_flux, bandpass) + if not faint and 'photon_ops' in config: + photon_ops = galsim.config.BuildPhotonOps( + config, 'photon_ops', base, logger) + else: + photon_ops = [] + + # Put the psfs at the start of the photon_ops. + # Probably a little better to put them a bit later than the + # start in some cases (e.g. after TimeSampler, + # PupilAnnulusSampler), but leave that as a todo for now. + photon_ops = psfs + photon_ops + image = gal.drawImage(bandpass, + method='phot', + offset=offset, + rng=self.rng, + n_photons=self.phot_flux, + image=image, + sensor=None, + photon_ops=photon_ops, + add_to_image=True, + poisson_flux=False) + base['realized_flux'] = image.added_flux + logger.debug('After .drawImage(...), image.added_flux %s', + image.added_flux) - @classmethod - def _fix_seds(cls, prof, bandpass, logger): - # If any SEDs are not currently using a LookupTable for the - # function or if they are using spline interpolation, then the - # codepath is quite slow. Better to fix them before doing - # WavelengthSampler. - if (isinstance(prof, galsim.SimpleChromaticTransformation) and - (not isinstance(prof._flux_ratio._spec, galsim.LookupTable) - or prof._flux_ratio._spec.interpolant != 'linear')): - if not cls._sed_logged: - logger.warning( - "Warning: Chromatic drawing is most efficient when SEDs have " - "interpolant='linear'. Switching LookupTables to use 'linear'." - ) - cls._sed_logged = True - sed = prof._flux_ratio - wave_list, _, _ = galsim.utilities.combine_wave_list(sed, bandpass) - f = np.broadcast_to(sed(wave_list), wave_list.shape) - new_spec = galsim.LookupTable(wave_list, f, interpolant='linear') - new_sed = galsim.SED( - new_spec, - 'nm', - 'fphotons' if sed.spectral else '1' - ) - prof._flux_ratio = new_sed - - # Also recurse onto any components. - if isinstance(prof, galsim.ChromaticObject): - if hasattr(prof, 'obj_list'): - for obj in prof.obj_list: - cls._fix_seds(obj, bandpass, logger) - if hasattr(prof, 'original'): - cls._fix_seds(prof.original, bandpass, logger) + return image RegisterStampType('RubinDeepCoadd', RubinDeepCoaddStampBuilder()) From bedafcc7b7cd4b9ede1ec7788043788a5dc72bd2 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 17:19:27 -0700 Subject: [PATCH 12/30] enable iterating over coadds --- config/imsim-config-deep-coadd.yaml | 6 +++++- imsim/deep_coadd.py | 6 +++++- 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 1df19422..70bc33ee 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -5,6 +5,7 @@ eval_variables: itract: { type: DeepCoaddData, field: tract } ipatch: { type: DeepCoaddData, field: patch } sband: { type: DeepCoaddData, field: band } + incoadds: { type: DeepCoaddData, field: ncoadds } input: sky_catalog: @@ -58,7 +59,10 @@ psf: output: type: RubinDeepCoaddOutput - coadd_num: 0 + coadd_num: + type: Sequence + nitems: $ncoadds + first: 0 nproc: 1 nfiles: 1 dir: fits diff --git a/imsim/deep_coadd.py b/imsim/deep_coadd.py index 6771dd25..eeebcf53 100644 --- a/imsim/deep_coadd.py +++ b/imsim/deep_coadd.py @@ -191,12 +191,16 @@ def BuildRubinCoaddPSF(config, base, ignore, gsparams, logger): def DeepCoaddData(config, base, value_type): deep_coadd = GetInputObj('deep_coadd', config, base, 'DeepCoaddData') + num_coadds = len(deep_coadd.data_ids) + req = { 'field': str } params, safe = GetAllParams(config, base, req=req) field = params['field'] + if field == 'ncoadds': + return num_coadds, safe coadd_num = base['coadd_num'] - assert (coadd_num >= 0 and coadd_num < len(deep_coadd.data_ids)) + assert (coadd_num >= 0 and coadd_num < num_coadds) data_id = deep_coadd.data_ids[coadd_num] val = value_type(data_id.get(field, None)) From 5917ee6215d6ca27a37ed60d34e2999188979014 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 17:26:29 -0700 Subject: [PATCH 13/30] revert inadvertant changes --- imsim/lsst_image.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/imsim/lsst_image.py b/imsim/lsst_image.py index 43c0393d..6cbca210 100644 --- a/imsim/lsst_image.py +++ b/imsim/lsst_image.py @@ -243,6 +243,7 @@ def _set_config_image_pos(config, base): config: The configuration dictionary for the image field. base: The base configuration dictionary. """ + if 'image_pos' in config and 'world_pos' in config: raise galsim.config.GalSimConfigValueError( "Both image_pos and world_pos specified for LSST_Image.", @@ -266,6 +267,7 @@ def _set_config_image_pos(config, base): class LSST_ImageBuilder(LSST_ImageBuilderBase): """This is based on the GalSim "Scattered" image type. + Additional features supported here: * batching the stamp generation * checkpointing if input.checkpoint is given From 909073bbb84b7629b99ae6d6fb9c0ecb0275ed96 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 17:27:33 -0700 Subject: [PATCH 14/30] revert inadvertant whitespace change --- imsim/lsst_image.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/imsim/lsst_image.py b/imsim/lsst_image.py index 6cbca210..b936f266 100644 --- a/imsim/lsst_image.py +++ b/imsim/lsst_image.py @@ -267,7 +267,7 @@ def _set_config_image_pos(config, base): class LSST_ImageBuilder(LSST_ImageBuilderBase): """This is based on the GalSim "Scattered" image type. - + Additional features supported here: * batching the stamp generation * checkpointing if input.checkpoint is given From 27fc58f464afa71cde42073b289469d17e5dc342 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 15 Sep 2026 17:56:19 -0700 Subject: [PATCH 15/30] move psf builder to a separate file; try to use consistent naming throughout --- config/imsim-config-deep-coadd.yaml | 6 +-- imsim/__init__.py | 3 +- imsim/deep_coadd_bandpass.py | 4 +- imsim/deep_coadd_image.py | 14 +++---- imsim/deep_coadd_output.py | 4 +- imsim/deep_coadd_psf.py | 35 +++++++++++++++++ imsim/deep_coadd_stamp.py | 2 +- imsim/deep_coadd_wcs.py | 6 +-- imsim/{deep_coadd.py => deep_coadds.py} | 51 ++++++------------------- 9 files changed, 67 insertions(+), 58 deletions(-) create mode 100644 imsim/deep_coadd_psf.py rename imsim/{deep_coadd.py => deep_coadds.py} (80%) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 70bc33ee..d012770b 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -26,14 +26,14 @@ gal: type: SkyCatObj stamp: - type: RubinDeepCoadd + type: RubinDeepCoaddStamp fft_sb_thresh: 2.0e5 # When to switch to fft rendering. max_flux_simple: 100 # When to switch to simple SED. world_pos: type: SkyCatWorldPos image: - type: RubinDeepCoadd + type: RubinDeepCoaddImage random_seed: 42 nproc: 1 nbatch: 100 @@ -55,7 +55,7 @@ image: band: $band psf: - type: RubinCoaddPSF + type: RubinDeepCoaddPSF output: type: RubinDeepCoaddOutput diff --git a/imsim/__init__.py b/imsim/__init__.py index 916676f0..11c4ef47 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -41,9 +41,10 @@ from .table_row import * from .photon_pooling import * from .full_focal_plane import * -from .deep_coadd import * +from .deep_coadds import * from .deep_coadd_image import * from .deep_coadd_wcs import * from .deep_coadd_output import * from .deep_coadd_stamp import * from .deep_coadd_bandpass import * +from .deep_coadd_psf import * diff --git a/imsim/deep_coadd_bandpass.py b/imsim/deep_coadd_bandpass.py index c3df43f5..00e9f131 100644 --- a/imsim/deep_coadd_bandpass.py +++ b/imsim/deep_coadd_bandpass.py @@ -19,8 +19,8 @@ def buildBandpass(self, config, base, logger): """ req = { 'band': str } params, safe = galsim.config.GetAllParams(config, base, req=req) - deep_coadds = galsim.config.GetInputObj('deep_coadd', config, base, - 'DeepCoaddBandpass') + deep_coadds = galsim.config.GetInputObj('deep_coadds', config, base, + 'RubinDeepCoaddBandpass') return deep_coadds.getBandpass(params['band']), safe diff --git a/imsim/deep_coadd_image.py b/imsim/deep_coadd_image.py index 5b2196c5..1ab5fc51 100644 --- a/imsim/deep_coadd_image.py +++ b/imsim/deep_coadd_image.py @@ -6,10 +6,10 @@ from .lsst_image import LSST_ImageBuilderBase -__all__ = ['RubinCoaddImageBuilder'] +__all__ = ['RubinDeepCoaddImageBuilder'] -class RubinCoaddImageBuilder(LSST_ImageBuilderBase): # noqa: N801 +class RubinDeepCoaddImageBuilder(LSST_ImageBuilderBase): # noqa: N801 def setup(self, config, base, image_num, obj_num, ignore, logger): """ Do the initialization and setup for a Rubin/LSST deep coadd image. @@ -60,16 +60,16 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): except galsim.config.GalSimConfigError: self.checkpoint = None - deep_coadd = GetInputObj( - 'deep_coadd', config, base, 'LSSTC_CoaddImageBuilder') - bbox = deep_coadd.skymap[self.tract][self.patch].getOuterBBox() + deep_coadds = GetInputObj( + 'deep_coadds', config, base, 'LSSTC_CoaddImageBuilder') + bbox = deep_coadds.skymap[self.tract][self.patch].getOuterBBox() xsize, ysize = bbox.width, bbox.height data_id = {'band': self.band, 'tract': self.tract, 'patch': self.patch} self.noise_level = np.std( - deep_coadd.get(data_id=data_id).noise_realizations[0].array) + deep_coadds.get(data_id=data_id).noise_realizations[0].array) return xsize, ysize @@ -127,4 +127,4 @@ def buildImage(self, config, base, image_num, obj_num, logger): return full_image, current_var -RegisterImageType('RubinDeepCoadd', RubinCoaddImageBuilder()) +RegisterImageType('RubinDeepCoaddImage', RubinDeepCoaddImageBuilder()) diff --git a/imsim/deep_coadd_output.py b/imsim/deep_coadd_output.py index a32029ab..7741ce2b 100644 --- a/imsim/deep_coadd_output.py +++ b/imsim/deep_coadd_output.py @@ -17,7 +17,7 @@ def setup(self, config, base, file_num, logger): eval_base_variables.append('coadd_num') self._added_eval_base_variables = True base['coadd_num'] = coadd_num - deep_coadds = GetInputObj('deep_coadd', config, base, + deep_coadds = GetInputObj('deep_coadds', config, base, 'RubinDeepCoaddOutputBuilder') data_id = deep_coadds.data_ids[coadd_num] patch = deep_coadds.skymap[data_id['tract']][data_id['patch']] @@ -29,7 +29,7 @@ def getNFiles(self, config, base, logger=None): if 'nfiles' in config: return ParseValue(config, 'nfiles', base, int)[0] else: - deep_coadds = GetInputObj('deep_coadd', config, base, + deep_coadds = GetInputObj('deep_coadds', config, base, 'RubinDeepCoaddOutputBuilder') return len(deep_coadds.data_ids) diff --git a/imsim/deep_coadd_psf.py b/imsim/deep_coadd_psf.py new file mode 100644 index 00000000..d40142f5 --- /dev/null +++ b/imsim/deep_coadd_psf.py @@ -0,0 +1,35 @@ +import galsim +from galsim.config import GetInputObj, RegisterObjectType + + +__all__ = ["BuildRubinDeepCoaddPSF"] + + +def BuildRubinDeepCoaddPSF(config, base, ignore, gsparams, logger): + """ + Build PSFs from Rubin deep_coadds. The deep_coadd image will be + retrieved from the data repository using the butler for the + tract-patch combination corresponding to the object's sky position. + + Assuming the input.deep_coadd object is defined, to use this + PSF, add the following to the config yaml: + + input.atm_psf: "" # disable the atmospheric PSF + psf: + type: RubinDeepCoaddPSF + """ + deep_coadds = GetInputObj('deep_coadds', config, base, 'RubinDeepCoaddPSF') + coadd_num = base['coadd_num'] + assert (coadd_num >= 0 and coadd_num < len(deep_coadds.data_ids)) + data_id = deep_coadds.data_ids[coadd_num] + image_pos = base['image_pos'] + celestial_coord = base['wcs'].toWorld(image_pos) + ra = celestial_coord.ra / galsim.degrees + dec = celestial_coord.dec / galsim.degrees + band = data_id['band'] + safe = False + return deep_coadds.getPSF(ra, dec, band), safe + + +RegisterObjectType('RubinDeepCoaddPSF', BuildRubinDeepCoaddPSF, + input_type='deep_coadds') diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 0d5a8adc..090c7b72 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -192,4 +192,4 @@ def draw(self, prof, image, method, offset, config, base, logger): return image -RegisterStampType('RubinDeepCoadd', RubinDeepCoaddStampBuilder()) +RegisterStampType('RubinDeepCoaddStamp', RubinDeepCoaddStampBuilder()) diff --git a/imsim/deep_coadd_wcs.py b/imsim/deep_coadd_wcs.py index dee01aed..a12bbf4e 100644 --- a/imsim/deep_coadd_wcs.py +++ b/imsim/deep_coadd_wcs.py @@ -27,9 +27,9 @@ def buildWCS(self, config, base, logger): } params, _ = GetAllParams(config, base, req=req) data_id = { key: params[key] for key in req } - deep_coadd = GetInputObj("deep_coadd", config, base, - "RubinDeepCoaddWCSBuilder") - return deep_coadd.getWcs(data_id) + deep_coadds = GetInputObj("deep_coadds", config, base, + "RubinDeepCoaddWCSBuilder") + return deep_coadds.getWcs(data_id) RegisterWCSType("RubinDeepCoaddWcs", RubinDeepCoaddWCSBuilder()) diff --git a/imsim/deep_coadd.py b/imsim/deep_coadds.py similarity index 80% rename from imsim/deep_coadd.py rename to imsim/deep_coadds.py index eeebcf53..1395ea38 100644 --- a/imsim/deep_coadd.py +++ b/imsim/deep_coadds.py @@ -3,13 +3,13 @@ """ import pandas as pd import galsim -from galsim.config import (InputLoader, RegisterInputType, RegisterObjectType, - RegisterValueType, GetAllParams, GetInputObj) +from galsim.config import (InputLoader, RegisterInputType, RegisterValueType, + GetAllParams, GetInputObj) import lsst.daf.butler as daf_butler import lsst.geom -class DeepCoadd: +class DeepCoadds: def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """ Parameters @@ -117,20 +117,20 @@ def getPSF(self, ra, dec, band): class DeepCoaddLoader(InputLoader): """ - Load the deep_coadd input object. Here's an example yaml entry: + Load the deep_coadds input object. Here's an example yaml entry: - input.deep_coadd: + input.deep_coadds: repo: dp2 collection: dp2 """ def __init__(self): - super().__init__(init_func=DeepCoadd, takes_logger=True, + super().__init__(init_func=DeepCoadds, takes_logger=True, use_proxy=False) self.butler = None self.deep_coadd_list = None def getKwargs(self, config, base, logger): - logger.debug("Get kwargs for DeepCoadd") + logger.debug("Get kwargs for DeepCoadds") req = { "repo": str, "collection": str, @@ -163,35 +163,9 @@ def getKwargs(self, config, base, logger): return kwargs, safe -def BuildRubinCoaddPSF(config, base, ignore, gsparams, logger): - """ - Build PSFs from Rubin deep_coadds. The deep_coadd image will be - retrieved from the data repository using the butler for the - tract-patch combination corresponding to the object's sky position. - - Assuming the input.deep_coadd object is defined, to use this - PSF, add the following to the config yaml: - - input.atm_psf: "" # disable the atmospheric PSF - psf: - type: RubinCoaddPSF - """ - deep_coadd = GetInputObj('deep_coadd', config, base, 'RubinCoaddPSF') - coadd_num = base['coadd_num'] - assert (coadd_num >= 0 and coadd_num < len(deep_coadd.data_ids)) - data_id = deep_coadd.data_ids[coadd_num] - image_pos = base['image_pos'] - celestial_coord = base['wcs'].toWorld(image_pos) - ra = celestial_coord.ra / galsim.degrees - dec = celestial_coord.dec / galsim.degrees - band = data_id['band'] - safe = False - return deep_coadd.getPSF(ra, dec, band), safe - - def DeepCoaddData(config, base, value_type): - deep_coadd = GetInputObj('deep_coadd', config, base, 'DeepCoaddData') - num_coadds = len(deep_coadd.data_ids) + deep_coadds = GetInputObj('deep_coadds', config, base, 'DeepCoaddData') + num_coadds = len(deep_coadds.data_ids) req = { 'field': str } params, safe = GetAllParams(config, base, req=req) @@ -201,14 +175,13 @@ def DeepCoaddData(config, base, value_type): coadd_num = base['coadd_num'] assert (coadd_num >= 0 and coadd_num < num_coadds) - data_id = deep_coadd.data_ids[coadd_num] + data_id = deep_coadds.data_ids[coadd_num] val = value_type(data_id.get(field, None)) return val, safe -RegisterInputType('deep_coadd', DeepCoaddLoader()) -RegisterObjectType('RubinCoaddPSF', BuildRubinCoaddPSF, input_type='deep_coadd') +RegisterInputType('deep_coadds', DeepCoaddLoader()) RegisterValueType('DeepCoaddData', DeepCoaddData, [int, str], - input_type='deep_coadd') + input_type='deep_coadds') From a991c70615d03b8eb381eac8515f5d889e8c2e7c Mon Sep 17 00:00:00 2001 From: James Chiang Date: Thu, 17 Sep 2026 13:52:40 -0700 Subject: [PATCH 16/30] switch to using lsst.images version of deep_coadd throughout; also consider band-dependence in the internal caches; fix pixel cell center calculation --- imsim/deep_coadds.py | 80 +++++++++++++++++++++----------------------- 1 file changed, 39 insertions(+), 41 deletions(-) diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 1395ea38..390d4714 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -9,6 +9,9 @@ import lsst.geom +__all__ = ["DeepCoadds", "DeepCoaddLoader", "DeepCoaddData"] + + class DeepCoadds: def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """ @@ -29,7 +32,8 @@ def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): self.data_ids = data_ids self.dstype = dstype self._psf_cache = {} - self._butler_cache = {} + self._grid_cache = {} + self._wcs_cache = {} def get(self, index=None, data_id=None): """Return the deep_coadd using the index of the self.data_ids list, @@ -41,8 +45,9 @@ def get(self, index=None, data_id=None): data_id = self.data_ids[index] elif data_id not in self.data_ids: return None - data_id['skymap'] = self.skymap_name - deep_coadd = self.butler.get(self.dstype, **data_id) + my_data_id = data_id.copy() + my_data_id['skymap'] = self.skymap_name + deep_coadd = self.butler.get(self.dstype, **my_data_id) # Ensure this is a lsst.images version: assert hasattr(deep_coadd, 'to_legacy') @@ -56,8 +61,11 @@ def getBandpass(self, band): return galsim.Bandpass(lut, wave_type='nm').thin() def getWcs(self, data_id): - deep_coadd = self.get(data_id=data_id) - return galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) + key = data_id['tract'], data_id['patch'] + if key not in self._wcs_cache: + deep_coadd = self.get(data_id=data_id) + self._wcs_cache[key] = galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) + return self._wcs_cache[key] def getPSF(self, ra, dec, band): """Return the cell coadd PSF, evaluated at the center of @@ -71,47 +79,37 @@ def getPSF(self, ra, dec, band): tract_info = self.skymap.findTract(sky_coords) tract = tract_info.getId() patch = tract_info.findPatch(sky_coords).getSequentialIndex() - - # Cache info from the butler for the corresponding deep_coadd. - butler_key = tract, patch - if butler_key not in self._butler_cache: - dataId = dict(skymap=self.skymap_name, tract=tract, patch=patch, - band=band) - deep_coadd = self.butler.get(self.dstype, **dataId) - if hasattr(deep_coadd, "to_legacy"): - # This is a `lsst.images` object. Use the `.to_legacy` - # function to convert to a `lsst.afw.image.Exposure` so - # that we can access the functions to generate a PSF - # image array to pass to galsim. - deep_coadd = deep_coadd.to_legacy() - wcs = deep_coadd.getWcs() - psf_grid = deep_coadd.getPsf() # This is the grid of coadd psfs. - self._butler_cache[butler_key] = wcs, psf_grid + data_id = dict(tract=tract, patch=patch, band=band) + wcs = self.getWcs(data_id) + + # Cache the cell grid data for the corresponding deep_coadd, + # keyed by (tract, patch, band). + grid_key = tract, patch, band + if grid_key not in self._grid_cache: + deep_coadd = self.get(data_id=data_id) + grid = deep_coadd.grid + psf = deep_coadd.psf + self._grid_cache[grid_key] = grid, psf else: - wcs, psf_grid = self._butler_cache[butler_key] + grid, psf = self._grid_cache[grid_key] - # Cache cell PSFs based on tract, patch, psf_grid index. - pixel_coords = wcs.skyToPixel(sky_coords) - psf_grid_index = psf_grid.grid.index(lsst.geom.Point2I(pixel_coords)) - psf_key = tract, patch, psf_grid_index + # xy offsets for the current grid. + x_offset = grid.bbox.start.x + y_offset = grid.bbox.start.y + + # Cache cell PSFs, keyed by (tract, patch, band, grid_index). + x, y = wcs.toImage(ra, dec, units='degrees') + grid_index = grid.index_of(x=x + x_offset, y=y + y_offset) + psf_key = tract, patch, band, grid_index if psf_key not in self._psf_cache: # Evaluate PSF at cell center. - x0 = ((pixel_coords.x // psf_grid.grid.cell_size.x) - * psf_grid.grid.cell_size.x + psf_grid.grid.shape.x/2) - y0 = ((pixel_coords.y // psf_grid.grid.cell_size.y) - * psf_grid.grid.cell_size.y + psf_grid.grid.shape.y/2) - pixel_cell_center = lsst.geom.Point2D(x0, y0) - sky_cell_center = wcs.pixelToSky(pixel_cell_center) - # The following is based on the lsst/source_injection - # implementation. See https://github.com/lsst/source_injection/blob/w.2026.31/python/lsst/source/injection/inject_engine.py#L705 - mat = wcs.linearizePixelToSky( - sky_cell_center, lsst.geom.arcseconds).getMatrix() - galsim_wcs = galsim.JacobianWCS(mat[0, 0], mat[0, 1], - mat[1, 0], mat[1, 1]) - psf_array = psf_grid.computeKernelImage(pixel_cell_center).array + x0 = ( (x // grid.cell_shape.x) * grid.cell_shape.x + + grid.cell_shape.x/2 + x_offset ) + y0 = ( (y // grid.cell_shape.y) * grid.cell_shape.y + + grid.cell_shape.y/2 + y_offset ) + psf_array = psf.compute_kernel_image(x=x0, y=y0).array self._psf_cache[psf_key] = galsim.InterpolatedImage( - galsim.Image(psf_array), wcs=galsim_wcs) - + galsim.Image(psf_array), wcs=wcs) return self._psf_cache[psf_key] From cc1e198cba9a3d56c08c88a74e6d464cca3369a1 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 18 Sep 2026 13:24:37 -0700 Subject: [PATCH 17/30] add access to number of visits per cell; refactor caching --- imsim/deep_coadds.py | 93 ++++++++++++++++++++++++++++++-------------- 1 file changed, 64 insertions(+), 29 deletions(-) diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 390d4714..ce97b1e7 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -1,6 +1,7 @@ """ Interface to Rubin deep_coadds """ +from collections import namedtuple import pandas as pd import galsim from galsim.config import (InputLoader, RegisterInputType, RegisterValueType, @@ -12,6 +13,9 @@ __all__ = ["DeepCoadds", "DeepCoaddLoader", "DeepCoaddData"] +GridKey = namedtuple('GridKey', ['tract', 'patch', 'band']) + + class DeepCoadds: def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """ @@ -34,6 +38,7 @@ def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): self._psf_cache = {} self._grid_cache = {} self._wcs_cache = {} + self._num_visits = {} def get(self, index=None, data_id=None): """Return the deep_coadd using the index of the self.data_ids list, @@ -49,8 +54,16 @@ def get(self, index=None, data_id=None): my_data_id['skymap'] = self.skymap_name deep_coadd = self.butler.get(self.dstype, **my_data_id) - # Ensure this is a lsst.images version: - assert hasattr(deep_coadd, 'to_legacy') + # Fill in various cached quantities to avoid repeated butler + # gets. + wcs_key = data_id['tract'], data_id['patch'] + if wcs_key not in self._wcs_cache: + self._wcs_cache[wcs_key] = galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) + + grid_key = GridKey(data_id['tract'], data_id['patch'], data_id['band']) + self._num_visits[grid_key] = self.num_visits_per_cell(deep_coadd) + self._grid_cache[grid_key] = deep_coadd.grid, deep_coadd.psf + return deep_coadd def getBandpass(self, band): @@ -63,44 +76,24 @@ def getBandpass(self, band): def getWcs(self, data_id): key = data_id['tract'], data_id['patch'] if key not in self._wcs_cache: - deep_coadd = self.get(data_id=data_id) - self._wcs_cache[key] = galsim.AstropyWCS(wcs=deep_coadd.fits_wcs) + self.get(data_id=data_id) # This will fill the cache. return self._wcs_cache[key] def getPSF(self, ra, dec, band): """Return the cell coadd PSF, evaluated at the center of the cell containing this sky position. """ - # Find the tract, patch for this location. - sky_coords = lsst.geom.SpherePoint( - lsst.geom.Angle(ra*lsst.geom.degrees), - lsst.geom.Angle(dec*lsst.geom.degrees) - ) - tract_info = self.skymap.findTract(sky_coords) - tract = tract_info.getId() - patch = tract_info.findPatch(sky_coords).getSequentialIndex() - data_id = dict(tract=tract, patch=patch, band=band) - wcs = self.getWcs(data_id) - - # Cache the cell grid data for the corresponding deep_coadd, - # keyed by (tract, patch, band). - grid_key = tract, patch, band - if grid_key not in self._grid_cache: - deep_coadd = self.get(data_id=data_id) - grid = deep_coadd.grid - psf = deep_coadd.psf - self._grid_cache[grid_key] = grid, psf - else: - grid, psf = self._grid_cache[grid_key] + grid_key, cell_index, x, y = self._get_cache_keys(ra, dec, band) + data_id = dict(tract=grid_key.tract, patch=grid_key.patch, + band=grid_key.band) + grid, psf = self._grid_cache[grid_key] # xy offsets for the current grid. x_offset = grid.bbox.start.x y_offset = grid.bbox.start.y - # Cache cell PSFs, keyed by (tract, patch, band, grid_index). - x, y = wcs.toImage(ra, dec, units='degrees') - grid_index = grid.index_of(x=x + x_offset, y=y + y_offset) - psf_key = tract, patch, band, grid_index + # Access the cached cell PSFs. + psf_key = grid_key, cell_index if psf_key not in self._psf_cache: # Evaluate PSF at cell center. x0 = ( (x // grid.cell_shape.x) * grid.cell_shape.x @@ -108,10 +101,52 @@ def getPSF(self, ra, dec, band): y0 = ( (y // grid.cell_shape.y) * grid.cell_shape.y + grid.cell_shape.y/2 + y_offset ) psf_array = psf.compute_kernel_image(x=x0, y=y0).array + wcs = self.getWcs(data_id) self._psf_cache[psf_key] = galsim.InterpolatedImage( galsim.Image(psf_array), wcs=wcs) return self._psf_cache[psf_key] + def getNumVisits(self, ra, dec, band): + grid_key, cell_index, _, _ = self._get_cache_keys(ra, dec, band) + i = int(cell_index.i) + j = int(cell_index.j) + return self._num_visits[grid_key][(i, j)] + + def _get_cache_keys(self, ra, dec, band): + # Find the tract, patch for this location. + sky_coords = lsst.geom.SpherePoint( + lsst.geom.Angle(ra*lsst.geom.degrees), + lsst.geom.Angle(dec*lsst.geom.degrees) + ) + tract_info = self.skymap.findTract(sky_coords) + tract = tract_info.getId() + patch = tract_info.findPatch(sky_coords).getSequentialIndex() + grid_key = GridKey(tract, patch, band) + data_id = dict(tract=grid_key.tract, patch=grid_key.patch, + band=grid_key.band) + if grid_key not in self._grid_cache: + self.get(data_id=data_id) + grid, _ = self._grid_cache[grid_key] + + # xy offsets for the current grid. + x_offset = grid.bbox.start.x + y_offset = grid.bbox.start.y + + # Get cell_index for the requested location. + wcs = self.getWcs(data_id) + x, y = wcs.toImage(ra, dec, units='degrees') + cell_index = grid.index_of(x=x + x_offset, y=y + y_offset) + return grid_key, cell_index, x, y + + @staticmethod + def num_visits_per_cell(deep_coadd): + df0 = deep_coadd.provenance.contributions.to_pandas() + num_visits = {} + for i, j in sorted(set(zip(df0['cell_i'], df0['cell_j']))): + df = df0.query(f"cell_i == {i} and cell_j == {j}") + num_visits[(i, j)] = sum(df['overlap_fraction']) + return num_visits + class DeepCoaddLoader(InputLoader): """ From 86c42e3557d29592833059e1d05a9a025584caa4 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 18 Sep 2026 21:49:24 -0700 Subject: [PATCH 18/30] use photon ops weighting to scale object fluxes so that output units are nJy --- config/imsim-config-deep-coadd.yaml | 8 +++++ imsim/deep_coadd_stamp.py | 56 ++++++++++++++++++++++++----- imsim/deep_coadds.py | 3 +- 3 files changed, 58 insertions(+), 9 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index d012770b..26ac5461 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -6,6 +6,11 @@ eval_variables: ipatch: { type: DeepCoaddData, field: patch } sband: { type: DeepCoaddData, field: band } incoadds: { type: DeepCoaddData, field: ncoadds } + fpupil_R_outer: 4.18 + fpupil_R_inner: 2.55 # M1 inner diameter is 2.558, but we need a bit of slack for off-axis rays + fpupil_area: + type: Eval + str: "np.pi * (pupil_R_outer**2 - pupil_R_inner**2) * 100**2" input: sky_catalog: @@ -31,6 +36,9 @@ stamp: max_flux_simple: 100 # When to switch to simple SED. world_pos: type: SkyCatWorldPos + band: $band + pupil_area: $pupil_area + exptime: 30.0 # Exposure time per visit image: type: RubinDeepCoaddImage diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 090c7b72..ca41cc6b 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -1,3 +1,4 @@ +import astropy.units as u import galsim from galsim.config import StampBuilder, RegisterStampType from .stamp import LSST_SiliconBuilder @@ -7,11 +8,36 @@ __all__ = ["RubinDeepCoaddStampBuilder"] +class DeepCoaddWeighting(galsim.PhotonOp): + """ + Weight photons so that final object fluxes are in nJy. + """ + def __init__(self, mag, nphotons): + """ + Parameters + ---------- + mag: float + The bandpass magnitude for the object being rendered. + nphotons: float + The number of photons to simulate. This should be + the total number over all visits contributing to the + object exposure. Nominally, this is given by + nphotons = photon_flux * num_visits * exptime_per_visit * pupil_area + """ + # Convert magnitude to nJy: + fnu = 1e9 * (mag*u.ABmag).to(u.Jy).to_value() + self.weight = fnu / nphotons + + def applyTo(self, photon_array, local_wcs=None, rng=None): + photon_array.flux *= self.weight + + class RubinDeepCoaddStampBuilder(LSST_SiliconBuilder): def setup(self, config, base, xsize, ysize, ignore, logger): # Use StampBuilder.setup(...) to set the object position. - ignore = ignore + ['fft_sb_thresh', 'max_flux_simple'] + ignore = ignore + ['fft_sb_thresh', 'max_flux_simple', 'band', + 'pupil_area', 'exptime'] _, _, image_pos, world_pos = StampBuilder.setup( self, config, base, xsize, ysize, ignore, logger) @@ -22,11 +48,24 @@ def setup(self, config, base, xsize, ysize, ignore, logger): self.rng = galsim.config.GetRNG(config, base, logger, "RubinDeepCoadd") self.image = base['current_image'] + + deep_coadds = galsim.config.GetInputObj('deep_coadds', config, base, + 'RubinDeepCoaddStampBuilder') + band = galsim.config.ParseValue(config, 'band', base, str)[0] + num_visits = deep_coadds.getNumVisits(world_pos.ra / galsim.degrees, + world_pos.dec / galsim.degrees, + band) bandpass = base['bandpass'] if not hasattr(obj, 'flux'): obj.flux = obj.calculateFlux(bandpass) - self.nominal_flux = obj.flux - self.phot_flux = galsim.PoissonDeviate(self.rng, mean=obj.flux)() + self.nominal_flux = obj.flux * num_visits + self.phot_flux = galsim.PoissonDeviate(self.rng, mean=self.nominal_flux)() + + pupil_area = galsim.config.ParseValue(config, 'pupil_area', base, float)[0] + exptime = galsim.config.ParseValue(config, 'exptime', base, float)[0] + sed = obj.sed / pupil_area / exptime + mag = sed.calculateMagnitude(bandpass) + self.deep_coadd_phot_op = DeepCoaddWeighting(mag, self.nominal_flux) # Save values for output truth catalog base['nominal_flux'] = self.nominal_flux @@ -115,6 +154,8 @@ def draw(self, prof, image, method, offset, config, base, logger): image.wcs = base['wcs'] + photon_ops = [self.deep_coadd_phot_op] + if method == 'fft': gal = gal.withFlux(self.nominal_flux, bandpass) fft_image = image.copy() @@ -122,10 +163,11 @@ def draw(self, prof, image, method, offset, config, base, logger): kwargs = dict( method='fft', offset=fft_offset, - image=fft_image + image=fft_image, + photon_ops=photon_ops, ) if not faint and config.get('fft_photon_ops'): - fft_photon_ops = galsim.config.BuildPhotonOps( + fft_photon_ops = photon_ops + galsim.config.BuildPhotonOps( config, 'fft_photon_ops', base, logger) kwargs.update({ "photon_ops": fft_photon_ops, @@ -165,10 +207,8 @@ def draw(self, prof, image, method, offset, config, base, logger): else: gal = gal.withFlux(self.phot_flux, bandpass) if not faint and 'photon_ops' in config: - photon_ops = galsim.config.BuildPhotonOps( + photon_ops = photon_ops + galsim.config.BuildPhotonOps( config, 'photon_ops', base, logger) - else: - photon_ops = [] # Put the psfs at the start of the photon_ops. # Probably a little better to put them a bit later than the diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index ce97b1e7..273d0525 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -71,7 +71,8 @@ def getBandpass(self, band): lut = galsim.LookupTable(passband['wavelength'], passband['throughput'], interpolant='linear') - return galsim.Bandpass(lut, wave_type='nm').thin() + bp = galsim.Bandpass(lut, wave_type='nm').withZeropoint('AB').thin() + return bp def getWcs(self, data_id): key = data_id['tract'], data_id['patch'] From 523ccbe81924eaff0fcd995d1b4e8e628d12b2e1 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Sun, 20 Sep 2026 10:09:16 -0700 Subject: [PATCH 19/30] add default input.deep_coadds config --- config/imsim-config-deep-coadd.yaml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 26ac5461..bb57e03c 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -17,6 +17,11 @@ input: file_name: default_sky_cat.yaml band: $band + deep_coadds: + repo: dp2 + collection: dp2 + skymap_name: lsst_cells_v2 + checkpoint: dir: checkpoint file_name: From 5fddd54dd301ca14ca0986b3f6fe978c2832f9c9 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Sun, 20 Sep 2026 20:29:22 -0700 Subject: [PATCH 20/30] set maximum number of simulated photons to avoid excess memory usage --- config/imsim-config-deep-coadd.yaml | 1 + imsim/deep_coadd_stamp.py | 42 ++++++++++++++++++----------- 2 files changed, 27 insertions(+), 16 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index bb57e03c..d3ab0bcd 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -39,6 +39,7 @@ stamp: type: RubinDeepCoaddStamp fft_sb_thresh: 2.0e5 # When to switch to fft rendering. max_flux_simple: 100 # When to switch to simple SED. + max_simulated_flux: 3e8 # Maximum simulated flux to avoid excessive memory world_pos: type: SkyCatWorldPos band: $band diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index ca41cc6b..ee31f993 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -21,7 +21,7 @@ def __init__(self, mag, nphotons): nphotons: float The number of photons to simulate. This should be the total number over all visits contributing to the - object exposure. Nominally, this is given by + object's counts. Nominally, this is given by nphotons = photon_flux * num_visits * exptime_per_visit * pupil_area """ # Convert magnitude to nJy: @@ -36,7 +36,8 @@ class RubinDeepCoaddStampBuilder(LSST_SiliconBuilder): def setup(self, config, base, xsize, ysize, ignore, logger): # Use StampBuilder.setup(...) to set the object position. - ignore = ignore + ['fft_sb_thresh', 'max_flux_simple', 'band', + ignore = ignore + ['fft_sb_thresh', 'max_flux_simple', + 'max_simulated_flux', 'band', 'pupil_area', 'exptime'] _, _, image_pos, world_pos = StampBuilder.setup( self, config, base, xsize, ysize, ignore, logger) @@ -51,24 +52,33 @@ def setup(self, config, base, xsize, ysize, ignore, logger): deep_coadds = galsim.config.GetInputObj('deep_coadds', config, base, 'RubinDeepCoaddStampBuilder') + ra = world_pos.ra / galsim.degrees + dec = world_pos.dec / galsim.degrees band = galsim.config.ParseValue(config, 'band', base, str)[0] - num_visits = deep_coadds.getNumVisits(world_pos.ra / galsim.degrees, - world_pos.dec / galsim.degrees, - band) + num_visits = deep_coadds.getNumVisits(ra, dec, band) + bandpass = base['bandpass'] if not hasattr(obj, 'flux'): obj.flux = obj.calculateFlux(bandpass) - self.nominal_flux = obj.flux * num_visits - self.phot_flux = galsim.PoissonDeviate(self.rng, mean=self.nominal_flux)() pupil_area = galsim.config.ParseValue(config, 'pupil_area', base, float)[0] exptime = galsim.config.ParseValue(config, 'exptime', base, float)[0] sed = obj.sed / pupil_area / exptime mag = sed.calculateMagnitude(bandpass) - self.deep_coadd_phot_op = DeepCoaddWeighting(mag, self.nominal_flux) + + # Scale the single-visit object flux by number of visits in + # order to get correct signal variance, but apply a ceiling on + # simulated flux to avoid excessive memory usage. The deep + # coadd weighting ensure the final flux has the right value in + # nJy. + max_simulated_flux = galsim.config.ParseValue( + config, 'max_simulated_flux', base, float)[0] + self.simulated_flux = min(obj.flux * num_visits, max_simulated_flux) + self.deep_coadd_phot_op = DeepCoaddWeighting(mag, self.simulated_flux) + self.phot_flux = galsim.PoissonDeviate(self.rng, mean=self.simulated_flux)() # Save values for output truth catalog - base['nominal_flux'] = self.nominal_flux + base['simulated_flux'] = self.simulated_flux base['phot_flux'] = self.phot_flux base['realized_flux'] = 0 # This will be updated by .drawImage @@ -82,7 +92,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # Get the stamp size from the size config entry xsize = ysize = galsim.config.ParseValue( config, 'size', base, int)[0] - elif self.nominal_flux < self._tiny_flux: + elif self.simulated_flux < self._tiny_flux: xsize = ysize = 32 else: # Determine the stamp size from the object flux. @@ -91,7 +101,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) stamp_size = get_stamp_size( obj_achrom=obj_achrom, - nominal_flux=self.nominal_flux, + nominal_flux=self.simulated_flux, noise_var=noise_var, Nmax=self._Nmax, pixel_scale=self._pixel_scale, @@ -99,8 +109,8 @@ def setup(self, config, base, xsize, ysize, ignore, logger): ) xsize = ysize = stamp_size - logger.info('Object %d will use stamp size %s, %s and nominal flux %s', - base.get('obj_num',0), xsize, ysize, self.nominal_flux) + logger.info('Object %d will use stamp size %s, %s and simulated flux %s', + base.get('obj_num',0), xsize, ysize, self.simulated_flux) return xsize, ysize, image_pos, world_pos def buildPSF(self, config, base, gsparams, logger): @@ -115,7 +125,7 @@ def buildPSF(self, config, base, gsparams, logger): fft_sb_thresh = None self.use_fft = (fft_sb_thresh is not None - and self.nominal_flux > fft_sb_thresh) + and self.simulated_flux > fft_sb_thresh) return psf def draw(self, prof, image, method, offset, config, base, logger): @@ -130,7 +140,7 @@ def draw(self, prof, image, method, offset, config, base, logger): bandpass = base['bandpass'] max_flux_simple = config.get('max_flux_simple', 100) - faint = self.nominal_flux < max_flux_simple + faint = self.simulated_flux < max_flux_simple if faint: logger.info("Flux = %.0f Using trivial sed", self.obj.flux) @@ -157,7 +167,7 @@ def draw(self, prof, image, method, offset, config, base, logger): photon_ops = [self.deep_coadd_phot_op] if method == 'fft': - gal = gal.withFlux(self.nominal_flux, bandpass) + gal = gal.withFlux(self.simulated_flux, bandpass) fft_image = image.copy() fft_offset = offset kwargs = dict( From 3beeafa66c66e1db6c2a33699f6d667bf79353f5 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Sun, 20 Sep 2026 20:31:38 -0700 Subject: [PATCH 21/30] use the noise realization to estimate the noise variance for stamp size calculation --- imsim/deep_coadd_stamp.py | 4 ++-- imsim/deep_coadds.py | 8 ++++++++ 2 files changed, 10 insertions(+), 2 deletions(-) diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index ee31f993..22499392 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -96,8 +96,8 @@ def setup(self, config, base, xsize, ysize, ignore, logger): xsize = ysize = 32 else: # Determine the stamp size from the object flux. - base['current_noise_image'] = base['current_image'] - noise_var = galsim.config.CalculateNoiseVariance(base) + noise_var = (deep_coadds.getNoiseVar(ra, dec, band) + / self.deep_coadd_phot_op.weight) obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) stamp_size = get_stamp_size( obj_achrom=obj_achrom, diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 273d0525..a501b2e1 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -2,6 +2,7 @@ Interface to Rubin deep_coadds """ from collections import namedtuple +import numpy as np import pandas as pd import galsim from galsim.config import (InputLoader, RegisterInputType, RegisterValueType, @@ -39,6 +40,7 @@ def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): self._grid_cache = {} self._wcs_cache = {} self._num_visits = {} + self._noise_var = {} def get(self, index=None, data_id=None): """Return the deep_coadd using the index of the self.data_ids list, @@ -63,6 +65,8 @@ def get(self, index=None, data_id=None): grid_key = GridKey(data_id['tract'], data_id['patch'], data_id['band']) self._num_visits[grid_key] = self.num_visits_per_cell(deep_coadd) self._grid_cache[grid_key] = deep_coadd.grid, deep_coadd.psf + self._noise_var[grid_key] \ + = np.var(deep_coadd.noise_realizations[0].array) return deep_coadd @@ -113,6 +117,10 @@ def getNumVisits(self, ra, dec, band): j = int(cell_index.j) return self._num_visits[grid_key][(i, j)] + def getNoiseVar(self, ra, dec, band): + grid_key, cell_index, _, _ = self._get_cache_keys(ra, dec, band) + return self._noise_var[grid_key] + def _get_cache_keys(self, ra, dec, band): # Find the tract, patch for this location. sky_coords = lsst.geom.SpherePoint( From 2edcc271b19aee2283b5a4774302dbb9b1dc6da8 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Wed, 23 Sep 2026 21:21:55 -0700 Subject: [PATCH 22/30] use maxN in drawImage to limit the memory usage for bright objects --- config/imsim-config-deep-coadd.yaml | 1 - imsim/deep_coadd_stamp.py | 26 ++++++++++++++++---------- 2 files changed, 16 insertions(+), 11 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index d3ab0bcd..bb57e03c 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -39,7 +39,6 @@ stamp: type: RubinDeepCoaddStamp fft_sb_thresh: 2.0e5 # When to switch to fft rendering. max_flux_simple: 100 # When to switch to simple SED. - max_simulated_flux: 3e8 # Maximum simulated flux to avoid excessive memory world_pos: type: SkyCatWorldPos band: $band diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 22499392..49315b9b 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -37,8 +37,7 @@ class RubinDeepCoaddStampBuilder(LSST_SiliconBuilder): def setup(self, config, base, xsize, ysize, ignore, logger): # Use StampBuilder.setup(...) to set the object position. ignore = ignore + ['fft_sb_thresh', 'max_flux_simple', - 'max_simulated_flux', 'band', - 'pupil_area', 'exptime'] + 'band', 'pupil_area', 'exptime'] _, _, image_pos, world_pos = StampBuilder.setup( self, config, base, xsize, ysize, ignore, logger) @@ -61,19 +60,18 @@ def setup(self, config, base, xsize, ysize, ignore, logger): if not hasattr(obj, 'flux'): obj.flux = obj.calculateFlux(bandpass) + # To get the AB magnitude, we need to divide out the pupil + # area and single visit exposure time from the object SED as + # provided from skyCatalogs. pupil_area = galsim.config.ParseValue(config, 'pupil_area', base, float)[0] exptime = galsim.config.ParseValue(config, 'exptime', base, float)[0] sed = obj.sed / pupil_area / exptime mag = sed.calculateMagnitude(bandpass) - # Scale the single-visit object flux by number of visits in - # order to get correct signal variance, but apply a ceiling on - # simulated flux to avoid excessive memory usage. The deep - # coadd weighting ensure the final flux has the right value in - # nJy. - max_simulated_flux = galsim.config.ParseValue( - config, 'max_simulated_flux', base, float)[0] - self.simulated_flux = min(obj.flux * num_visits, max_simulated_flux) + # Scale the single-visit object flux by the number of visits + # in order to get the correct signal variance. The deep coadd + # weighting ensures that the final flux is in units of nJy. + self.simulated_flux = obj.flux * num_visits self.deep_coadd_phot_op = DeepCoaddWeighting(mag, self.simulated_flux) self.phot_flux = galsim.PoissonDeviate(self.rng, mean=self.simulated_flux)() @@ -166,6 +164,12 @@ def draw(self, prof, image, method, offset, config, base, logger): photon_ops = [self.deep_coadd_phot_op] + # Set limit on the size of photons batches to consider when + # calling gsobject.drawImage. + maxN = int(1e6) + if 'maxN' in config: + maxN = galsim.config.ParseValue(config, 'maxN', base, int)[0] + if method == 'fft': gal = gal.withFlux(self.simulated_flux, bandpass) fft_image = image.copy() @@ -181,6 +185,7 @@ def draw(self, prof, image, method, offset, config, base, logger): config, 'fft_photon_ops', base, logger) kwargs.update({ "photon_ops": fft_photon_ops, + "maxN": maxN, "rng": self.rng, "n_subsample": 1, }) @@ -229,6 +234,7 @@ def draw(self, prof, image, method, offset, config, base, logger): method='phot', offset=offset, rng=self.rng, + maxN=maxN, n_photons=self.phot_flux, image=image, sensor=None, From 29b56430d729cf8c28b3ead285dab1c2d6a34ae9 Mon Sep 17 00:00:00 2001 From: Jim Chiang Date: Wed, 23 Sep 2026 21:22:23 -0700 Subject: [PATCH 23/30] add class docstring and some comments --- imsim/deep_coadds.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index a501b2e1..15d76686 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -18,6 +18,12 @@ class DeepCoadds: + """ + This class serves as the access point for deep_coadd information + obtained via requiring butler queries. It caches PSF, wcs, visit, + and noise information in order to minimize the number of + butler.get calls. + """ def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """ Parameters @@ -100,7 +106,8 @@ def getPSF(self, ra, dec, band): # Access the cached cell PSFs. psf_key = grid_key, cell_index if psf_key not in self._psf_cache: - # Evaluate PSF at cell center. + # Evaluate PSF at cell center. This is based on + # https://github.com/lsst/source_injection/blob/w.2026.38/python/lsst/source/injection/inject_engine.py#L705 x0 = ( (x // grid.cell_shape.x) * grid.cell_shape.x + grid.cell_shape.x/2 + x_offset ) y0 = ( (y // grid.cell_shape.y) * grid.cell_shape.y From 729d09a93ae1be3430515c9f6ebf4eb63dfc7ee3 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Mon, 12 Oct 2026 14:14:48 -0700 Subject: [PATCH 24/30] disable fft option for deep_coadd sims --- config/imsim-config-deep-coadd.yaml | 2 - imsim/deep_coadd_image.py | 4 +- imsim/deep_coadd_stamp.py | 114 +++++++--------------------- 3 files changed, 29 insertions(+), 91 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index bb57e03c..210a92be 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -37,7 +37,6 @@ gal: stamp: type: RubinDeepCoaddStamp - fft_sb_thresh: 2.0e5 # When to switch to fft rendering. max_flux_simple: 100 # When to switch to simple SED. world_pos: type: SkyCatWorldPos @@ -51,7 +50,6 @@ image: nproc: 1 nbatch: 100 nsubbatch: 50 - nbatch_fft: 1 nbatch_per_checkpoint: 1 band: $band tract: $tract diff --git a/imsim/deep_coadd_image.py b/imsim/deep_coadd_image.py index 1ab5fc51..bfb94456 100644 --- a/imsim/deep_coadd_image.py +++ b/imsim/deep_coadd_image.py @@ -37,8 +37,7 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): logger.info('image %d: nobj = %d', image_num, self.nobjects) req = { 'tract': int, 'patch': int, 'band': str } - opt = { 'nbatch': int, 'nsubbatch': int, 'nbatch_fft': int, - 'nbatch_per_checkpoint': int} + opt = { 'nbatch': int, 'nsubbatch': int, 'nbatch_per_checkpoint': int } extra_ignore = ['image_pos', 'world_pos', 'stamp_size', 'stamp_xsize', 'stamp_ysize', 'nobjects' ] params = GetAllParams(config, base, req=req, opt=opt, @@ -53,7 +52,6 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.nbatch = params.get('nbatch', 10) self.nbatch_per_checkpoint = params.get('nbatch_per_checkpoint', 1) self.nsubbatch = params.get('nsubbatch', 50) - self.nbatch_fft = params.get('nbatch_fft', 1) try: self.checkpoint = GetInputObj( 'checkpoint', config, base, 'RubinCoaddImageBuilder') diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 49315b9b..87af6ac0 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -35,9 +35,11 @@ def applyTo(self, photon_array, local_wcs=None, rng=None): class RubinDeepCoaddStampBuilder(LSST_SiliconBuilder): def setup(self, config, base, xsize, ysize, ignore, logger): + # Always disable fft rendering: + self.use_fft = False + # Use StampBuilder.setup(...) to set the object position. - ignore = ignore + ['fft_sb_thresh', 'max_flux_simple', - 'band', 'pupil_area', 'exptime'] + ignore = ignore + ['max_flux_simple', 'band', 'pupil_area', 'exptime'] _, _, image_pos, world_pos = StampBuilder.setup( self, config, base, xsize, ysize, ignore, logger) @@ -114,16 +116,6 @@ def setup(self, config, base, xsize, ysize, ignore, logger): def buildPSF(self, config, base, gsparams, logger): psf = galsim.config.BuildGSObject( base, 'psf', gsparams=gsparams, logger=logger)[0] - - # Check if fft rendering should be used. - if 'fft_sb_thresh' in config: - fft_sb_thresh = galsim.config.ParseValue( - config, 'fft_sb_thresh', base, float)[0] - else: - fft_sb_thresh = None - - self.use_fft = (fft_sb_thresh is not None - and self.simulated_flux > fft_sb_thresh) return psf def draw(self, prof, image, method, offset, config, base, logger): @@ -170,80 +162,30 @@ def draw(self, prof, image, method, offset, config, base, logger): if 'maxN' in config: maxN = galsim.config.ParseValue(config, 'maxN', base, int)[0] - if method == 'fft': - gal = gal.withFlux(self.simulated_flux, bandpass) - fft_image = image.copy() - fft_offset = offset - kwargs = dict( - method='fft', - offset=fft_offset, - image=fft_image, - photon_ops=photon_ops, - ) - if not faint and config.get('fft_photon_ops'): - fft_photon_ops = photon_ops + galsim.config.BuildPhotonOps( - config, 'fft_photon_ops', base, logger) - kwargs.update({ - "photon_ops": fft_photon_ops, - "maxN": maxN, - "rng": self.rng, - "n_subsample": 1, - }) - - # Go back to a combined convolution for fft drawing. - prof = galsim.Convolve([gal] + psfs) - try: - fft_image = prof.drawImage(bandpass, **kwargs) - except galsim.errors.GalSimFFTSizeError as e: - # I think this shouldn't happen with the updates I - # made to how the image size is calculated, even for - # extremely bright things. So it should be ok to just - # report what happened, give some extra information to - # diagonose the problem and raise the error. - logger.error('Caught error trying to draw using FFT:') - logger.error('%s', e) - logger.error('You may need to add a gsparams field with ' - 'maximum_fft_size to either') - logger.error('the psf or gal field to allow larger FFTs.') - logger.info('prof = %r', prof) - logger.info('fft_image = %s', fft_image) - logger.info('offset = %r', offset) - raise - # Some pixels can end up negative from FFT numerics. Just - # set them to 0. - fft_image.array[fft_image.array < 0] = 0. - fft_image.addNoise(galsim.PoissonNoise(rng=self.rng)) - # In case we had to make a bigger image, just copy the - # part we need. - image += fft_image[image.bounds] - base['realized_flux'] = fft_image.added_flux - logger.debug('After .drawImage(...), fft_image.added_flux %s', - fft_image.added_flux) - else: - gal = gal.withFlux(self.phot_flux, bandpass) - if not faint and 'photon_ops' in config: - photon_ops = photon_ops + galsim.config.BuildPhotonOps( - config, 'photon_ops', base, logger) - - # Put the psfs at the start of the photon_ops. - # Probably a little better to put them a bit later than the - # start in some cases (e.g. after TimeSampler, - # PupilAnnulusSampler), but leave that as a todo for now. - photon_ops = psfs + photon_ops - image = gal.drawImage(bandpass, - method='phot', - offset=offset, - rng=self.rng, - maxN=maxN, - n_photons=self.phot_flux, - image=image, - sensor=None, - photon_ops=photon_ops, - add_to_image=True, - poisson_flux=False) - base['realized_flux'] = image.added_flux - logger.debug('After .drawImage(...), image.added_flux %s', - image.added_flux) + gal = gal.withFlux(self.phot_flux, bandpass) + if not faint and 'photon_ops' in config: + photon_ops = photon_ops + galsim.config.BuildPhotonOps( + config, 'photon_ops', base, logger) + + # Put the psfs at the start of the photon_ops. + # Probably a little better to put them a bit later than the + # start in some cases (e.g. after TimeSampler, + # PupilAnnulusSampler), but leave that as a todo for now. + photon_ops = psfs + photon_ops + image = gal.drawImage(bandpass, + method='phot', + offset=offset, + rng=self.rng, + maxN=maxN, + n_photons=self.phot_flux, + image=image, + sensor=None, + photon_ops=photon_ops, + add_to_image=True, + poisson_flux=False) + base['realized_flux'] = image.added_flux + logger.debug('After .drawImage(...), image.added_flux %s', + image.added_flux) return image From 02eda92a898db29679612622455f34f5dd78ef7d Mon Sep 17 00:00:00 2001 From: James Chiang Date: Tue, 13 Oct 2026 16:57:07 -0700 Subject: [PATCH 25/30] set noise_var to avoid apparent stamp edges; add output.truth section to deep coadd config --- config/imsim-config-deep-coadd.yaml | 20 ++++++++++++++++++++ imsim/deep_coadd_stamp.py | 4 ++-- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 210a92be..5a4f32c9 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -84,3 +84,23 @@ output: - $band - $tract - $patch + + truth: + dir: '@output.dir' + file_name: + type: FormattedStr + format: centroid_%s-%d-%d.txt.gz + items: + - $band + - $tract + - $patch + columns: + object_id: "@object_id" + ra: $sky_pos.ra.deg + dec: $sky_pos.dec.deg + x: $image_pos.x + y: $image_pos.y + nominal_flux: "@simulated_flux" + phot_flux: "@phot_flux" + realized_flux: "@realized_flux" + stamp_size: "@stamp_size" diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 87af6ac0..56b5e5d8 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -96,8 +96,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): xsize = ysize = 32 else: # Determine the stamp size from the object flux. - noise_var = (deep_coadds.getNoiseVar(ra, dec, band) - / self.deep_coadd_phot_op.weight) + noise_var = deep_coadds.getNoiseVar(ra, dec, band) / 3.0 obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) stamp_size = get_stamp_size( obj_achrom=obj_achrom, @@ -109,6 +108,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): ) xsize = ysize = stamp_size + base['stamp_size'] = xsize logger.info('Object %d will use stamp size %s, %s and simulated flux %s', base.get('obj_num',0), xsize, ysize, self.simulated_flux) return xsize, ysize, image_pos, world_pos From 58610d58c5830630454e5358a362eb9b4b4602c0 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 16 Oct 2026 08:41:57 -0700 Subject: [PATCH 26/30] update DeepCoadds.getPSF and .getNumVisits interfaces to pass the coordinate object instead of ra, dec --- imsim/deep_coadd_psf.py | 4 +--- imsim/deep_coadd_stamp.py | 4 +--- imsim/deep_coadds.py | 8 ++++++-- 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/imsim/deep_coadd_psf.py b/imsim/deep_coadd_psf.py index d40142f5..b5df6cdf 100644 --- a/imsim/deep_coadd_psf.py +++ b/imsim/deep_coadd_psf.py @@ -24,11 +24,9 @@ def BuildRubinDeepCoaddPSF(config, base, ignore, gsparams, logger): data_id = deep_coadds.data_ids[coadd_num] image_pos = base['image_pos'] celestial_coord = base['wcs'].toWorld(image_pos) - ra = celestial_coord.ra / galsim.degrees - dec = celestial_coord.dec / galsim.degrees band = data_id['band'] safe = False - return deep_coadds.getPSF(ra, dec, band), safe + return deep_coadds.getPSF(celestial_coord, band), safe RegisterObjectType('RubinDeepCoaddPSF', BuildRubinDeepCoaddPSF, diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 56b5e5d8..62b146eb 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -53,10 +53,8 @@ def setup(self, config, base, xsize, ysize, ignore, logger): deep_coadds = galsim.config.GetInputObj('deep_coadds', config, base, 'RubinDeepCoaddStampBuilder') - ra = world_pos.ra / galsim.degrees - dec = world_pos.dec / galsim.degrees band = galsim.config.ParseValue(config, 'band', base, str)[0] - num_visits = deep_coadds.getNumVisits(ra, dec, band) + num_visits = deep_coadds.getNumVisits(world_pos, band) bandpass = base['bandpass'] if not hasattr(obj, 'flux'): diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 15d76686..424f1e3b 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -90,10 +90,12 @@ def getWcs(self, data_id): self.get(data_id=data_id) # This will fill the cache. return self._wcs_cache[key] - def getPSF(self, ra, dec, band): + def getPSF(self, celestial_coord, band): """Return the cell coadd PSF, evaluated at the center of the cell containing this sky position. """ + ra = celestial_coord.ra / galsim.degrees + dec = celestial_coord.dec / galsim.degrees grid_key, cell_index, x, y = self._get_cache_keys(ra, dec, band) data_id = dict(tract=grid_key.tract, patch=grid_key.patch, band=grid_key.band) @@ -118,7 +120,9 @@ def getPSF(self, ra, dec, band): galsim.Image(psf_array), wcs=wcs) return self._psf_cache[psf_key] - def getNumVisits(self, ra, dec, band): + def getNumVisits(self, celestial_coord, band): + ra = celestial_coord.ra / galsim.degrees + dec = celestial_coord.dec / galsim.degrees grid_key, cell_index, _, _ = self._get_cache_keys(ra, dec, band) i = int(cell_index.i) j = int(cell_index.j) From dd3fb79f2cebf717209b0324593f1c4897ff3b94 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 16 Oct 2026 11:12:14 -0700 Subject: [PATCH 27/30] add RubinDeepCoaddGaussianNoise type and require user to provide a noise type explicitly in the config --- config/imsim-config-deep-coadd.yaml | 6 +++++ imsim/deep_coadd_image.py | 34 ++++++++++++++++++++++------- imsim/deep_coadd_stamp.py | 2 +- imsim/deep_coadds.py | 4 ++-- 4 files changed, 35 insertions(+), 11 deletions(-) diff --git a/config/imsim-config-deep-coadd.yaml b/config/imsim-config-deep-coadd.yaml index 5a4f32c9..b786f970 100644 --- a/config/imsim-config-deep-coadd.yaml +++ b/config/imsim-config-deep-coadd.yaml @@ -65,6 +65,12 @@ image: type: RubinDeepCoaddBandpass band: $band + noise: + type: RubinDeepCoaddGaussianNoise + band: $band + tract: $tract + patch: $patch + psf: type: RubinDeepCoaddPSF diff --git a/imsim/deep_coadd_image.py b/imsim/deep_coadd_image.py index bfb94456..d7f7e8cf 100644 --- a/imsim/deep_coadd_image.py +++ b/imsim/deep_coadd_image.py @@ -2,7 +2,7 @@ import numpy as np import galsim from galsim.config import (GetAllParams, AddNoise, RegisterImageType, - GetInputObj) + RegisterNoiseType, GetInputObj, GaussianNoiseBuilder) from .lsst_image import LSST_ImageBuilderBase @@ -26,6 +26,13 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): Returns: xsize, ysize """ + if 'noise' not in base['image']: + raise RuntimeError( + "A noise type needs to be specified for rendering the image. " + "The ``RubinDeepCoaddGaussianNoise`` type models the noise " + "level in the original coadd image and could be used here." + ) + self.nobjects = self.getNObj(config, base, image_num, logger=logger) if 'nobjects' in config: # User specified nobjects. @@ -47,8 +54,6 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): self.patch = params['patch'] self.band = params['band'] - self.add_noise = True - self.nbatch = params.get('nbatch', 10) self.nbatch_per_checkpoint = params.get('nbatch_per_checkpoint', 1) self.nsubbatch = params.get('nsubbatch', 50) @@ -84,11 +89,7 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, current_var: The current noise variance in each postage stamps. logger: A logger object to log progress. """ - if 'noise' not in base['image']: - base['image']['noise'] = { "type": "Gaussian", - "sigma": self.noise_level } - if self.add_noise: - AddNoise(base, image, current_var, logger) + AddNoise(base, image, current_var, logger) def buildImage(self, config, base, image_num, obj_num, logger): @@ -125,4 +126,21 @@ def buildImage(self, config, base, image_num, obj_num, logger): return full_image, current_var +class RubinDeepCoaddGaussianNoiseBuilder(GaussianNoiseBuilder): + """ + Treat the variance of the coadd noise realization plane + as Gaussian noise for a specified deep_coadd image. + """ + + def getNoiseVariance(self, config, base, full=False): + req = { 'tract': int, 'patch': int, 'band': str } + params = GetAllParams(config, base, req=req)[0] + deep_coadds = galsim.config.GetInputObj('deep_coadds', config, base, + 'RubinDeepCoaddGaussianNoiseBuilder') + return deep_coadds.getNoiseVar(params['tract'], params['patch'], + params['band']) + + RegisterImageType('RubinDeepCoaddImage', RubinDeepCoaddImageBuilder()) +RegisterNoiseType('RubinDeepCoaddGaussianNoise', + RubinDeepCoaddGaussianNoiseBuilder()) diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 62b146eb..b099340f 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -94,7 +94,7 @@ def setup(self, config, base, xsize, ysize, ignore, logger): xsize = ysize = 32 else: # Determine the stamp size from the object flux. - noise_var = deep_coadds.getNoiseVar(ra, dec, band) / 3.0 + noise_var = galsim.config.CalculateNoiseVariance(base) obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) stamp_size = get_stamp_size( obj_achrom=obj_achrom, diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 424f1e3b..0939de6b 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -128,8 +128,8 @@ def getNumVisits(self, celestial_coord, band): j = int(cell_index.j) return self._num_visits[grid_key][(i, j)] - def getNoiseVar(self, ra, dec, band): - grid_key, cell_index, _, _ = self._get_cache_keys(ra, dec, band) + def getNoiseVar(self, tract, patch, band): + grid_key = GridKey(tract, patch, band) return self._noise_var[grid_key] def _get_cache_keys(self, ra, dec, band): From 729e56a6493ce54531248515f22ce39418ec4d80 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 16 Oct 2026 11:18:50 -0700 Subject: [PATCH 28/30] divide noise_var value passed to get_stamp_size by 3 to suppress apparent stamp edges for long, thin galaxies --- imsim/deep_coadd_stamp.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index b099340f..85b3f858 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -96,10 +96,13 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # Determine the stamp size from the object flux. noise_var = galsim.config.CalculateNoiseVariance(base) obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) + # We're dividing the nominal noise variance by 3.0 when + # passing to get_stamp_size to suppress apparent stamp + # edges for some high aspect ratio galaxies. stamp_size = get_stamp_size( obj_achrom=obj_achrom, nominal_flux=self.simulated_flux, - noise_var=noise_var, + noise_var=noise_var/3.0, Nmax=self._Nmax, pixel_scale=self._pixel_scale, logger=logger From 31daf282438aec8b99ba1f382baef164e59e070b Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 16 Oct 2026 14:17:20 -0700 Subject: [PATCH 29/30] let GalSim handle the psfs instead of extracting them and adding them as photon_ops explicitly --- imsim/deep_coadd_stamp.py | 44 ++++++++++++++++----------------------- 1 file changed, 18 insertions(+), 26 deletions(-) diff --git a/imsim/deep_coadd_stamp.py b/imsim/deep_coadd_stamp.py index 85b3f858..5c300bbe 100644 --- a/imsim/deep_coadd_stamp.py +++ b/imsim/deep_coadd_stamp.py @@ -97,8 +97,8 @@ def setup(self, config, base, xsize, ysize, ignore, logger): noise_var = galsim.config.CalculateNoiseVariance(base) obj_achrom = obj.evaluateAtWavelength(bandpass.effective_wavelength) # We're dividing the nominal noise variance by 3.0 when - # passing to get_stamp_size to suppress apparent stamp - # edges for some high aspect ratio galaxies. + # passing the value to get_stamp_size in order to suppress + # apparent stamp edges for high aspect ratio galaxies. stamp_size = get_stamp_size( obj_achrom=obj_achrom, nominal_flux=self.simulated_flux, @@ -124,9 +124,6 @@ def draw(self, prof, image, method, offset, config, base, logger): image.photons = galsim.PhotonArray(0) return image - # Prof is normally a convolution here with obj_list being - # [gal, psf1, psf2,...] for some number of component PSFs. - gal, *psfs = prof.obj_list if hasattr(prof, 'obj_list') else [prof] obj_num = base.get('obj_num',0) bandpass = base['bandpass'] @@ -138,7 +135,7 @@ def draw(self, prof, image, method, offset, config, base, logger): for profile_wl in (bandpass.effective_wavelength, bandpass.red_limit, bandpass.blue_limit): - sed_value = gal.sed(profile_wl) + sed_value = prof.sed(profile_wl) if sed_value != 0: break if sed_value == 0: @@ -148,10 +145,10 @@ def draw(self, prof, image, method, offset, config, base, logger): logger.warning("Zero-valued SED for faint object %d, " "object_id %s. Skipping.", obj_num, object_id) return image - gal = gal.evaluateAtWavelength(profile_wl) - gal = gal * self._trivial_sed + prof = prof.evaluateAtWavelength(profile_wl) + prof = prof * self._trivial_sed else: - self._fix_seds(gal, bandpass, logger) + self._fix_seds(prof, bandpass, logger) image.wcs = base['wcs'] @@ -163,27 +160,22 @@ def draw(self, prof, image, method, offset, config, base, logger): if 'maxN' in config: maxN = galsim.config.ParseValue(config, 'maxN', base, int)[0] - gal = gal.withFlux(self.phot_flux, bandpass) + prof = prof.withFlux(self.phot_flux, bandpass) if not faint and 'photon_ops' in config: photon_ops = photon_ops + galsim.config.BuildPhotonOps( config, 'photon_ops', base, logger) - # Put the psfs at the start of the photon_ops. - # Probably a little better to put them a bit later than the - # start in some cases (e.g. after TimeSampler, - # PupilAnnulusSampler), but leave that as a todo for now. - photon_ops = psfs + photon_ops - image = gal.drawImage(bandpass, - method='phot', - offset=offset, - rng=self.rng, - maxN=maxN, - n_photons=self.phot_flux, - image=image, - sensor=None, - photon_ops=photon_ops, - add_to_image=True, - poisson_flux=False) + image = prof.drawImage(bandpass, + method='phot', + offset=offset, + rng=self.rng, + maxN=maxN, + n_photons=self.phot_flux, + image=image, + sensor=None, + photon_ops=photon_ops, + add_to_image=True, + poisson_flux=False) base['realized_flux'] = image.added_flux logger.debug('After .drawImage(...), image.added_flux %s', image.added_flux) From 3fa99a19ef392885befeaa397425c81e2111b2b8 Mon Sep 17 00:00:00 2001 From: James Chiang Date: Fri, 16 Oct 2026 15:56:25 -0700 Subject: [PATCH 30/30] tweak class docstring --- imsim/deep_coadds.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/imsim/deep_coadds.py b/imsim/deep_coadds.py index 0939de6b..5063b5dd 100644 --- a/imsim/deep_coadds.py +++ b/imsim/deep_coadds.py @@ -20,9 +20,8 @@ class DeepCoadds: """ This class serves as the access point for deep_coadd information - obtained via requiring butler queries. It caches PSF, wcs, visit, - and noise information in order to minimize the number of - butler.get calls. + obtained via butler queries. It caches PSF, wcs, visit, and noise + information in order to minimize the number of butler.get calls. """ def __init__(self, butler, skymap_name, data_ids=None, dstype="deep_coadd"): """