diff --git a/config/imsim-config-instcat.yaml b/config/imsim-config-instcat.yaml index 82e44fc6..3e28ff11 100644 --- a/config/imsim-config-instcat.yaml +++ b/config/imsim-config-instcat.yaml @@ -25,6 +25,33 @@ input.instance_catalog: file_name: default_catalog_file.txt sed_dir: $os.environ.get('SIMS_SED_LIBRARY_DIR') pupil_area: $pupil_area + # A single WCS calculated used when determining which objects are to be drawn + # by which detector. We could in principle use any WCS. Here we start from the + # same Batoid WCS as in the image. In the code, any SIP components of the + # fiducial WCS, once it's been calculated, are removed. + fiducial_wcs: + type: Batoid + + # These are required: + camera: $camera + boresight: $boresight + + obstime: + type: Eval + str: "astropy.time.Time(mjd_val, format='mjd', scale='tai')" + fmjd_val: { type: OpsimData, field: mjd } + + # The fiducial WCS will use this detector. We default here to the + # central detector, but it can be overrdden if needed. + det_name: "R22_S11" + wavelength: $(@image.bandpass).effective_wavelength + + # The rest can be omitted, since these are the default values, but shown here + # for reference. + temperature: 280 # Kelvin + pressure: 72.7 # kPa + H2O_pressure: 1.0 # kPa + order: 3 # Order of the SIP polynomial input.opsim_data: # Read the visit meta data. By default, we use the same file as the above diff --git a/config/imsim-config-photon-pooling.yaml b/config/imsim-config-photon-pooling.yaml index 541813e6..58a56a1c 100644 --- a/config/imsim-config-photon-pooling.yaml +++ b/config/imsim-config-photon-pooling.yaml @@ -65,4 +65,4 @@ output.photon_pooling_truth: nominal_flux: "@nominal_flux" # The nominal "expectation value" flux. phot_flux: "@phot_flux" # The realized flux for photon shooting. fft_flux: "@fft_flux" # If FFT rendering, then the flux used, including vignetting. - incident_flux: "@incident_flux" # The sum of the fluxes of all photons from the object. \ No newline at end of file + incident_flux: "@incident_flux" # The sum of the fluxes of all photons from the object. diff --git a/config/imsim-config-skycat.yaml b/config/imsim-config-skycat.yaml index 57ba1ab3..f13d8a52 100644 --- a/config/imsim-config-skycat.yaml +++ b/config/imsim-config-skycat.yaml @@ -26,6 +26,33 @@ input.sky_catalog: band: $band mjd: { type: OpsimData, field: mjd } pupil_area: $pupil_area + # A single WCS calculated used when determining which objects are to be drawn + # by which detector. We could in principle use any WCS. Here we start from the + # same Batoid WCS as in the image. In the code, any SIP components of the + # fiducial WCS, once it's been calculated, are removed. + fiducial_wcs: + type: Batoid + + # These are required: + camera: $camera + boresight: $boresight + + obstime: + type: Eval + str: "astropy.time.Time(mjd_val, format='mjd', scale='tai')" + fmjd_val: { type: OpsimData, field: mjd } + + # The fiducial WCS will use this detector. We default here to the + # central detector, but it can be overrdden if needed. + det_name: "R22_S11" + wavelength: $(@image.bandpass).effective_wavelength + + # The rest can be omitted, since these are the default values, but shown here + # for reference. + temperature: 280 # Kelvin + pressure: 72.7 # kPa + H2O_pressure: 1.0 # kPa + order: 3 # Order of the SIP polynomial input.opsim_data: file_name: default_opsim.db diff --git a/examples/example_instance_catalog_ffp.txt b/examples/example_instance_catalog_ffp.txt new file mode 100644 index 00000000..a2dadb84 --- /dev/null +++ b/examples/example_instance_catalog_ffp.txt @@ -0,0 +1,23 @@ +rightascension 60.49045502638662697 +declination -38.16437495898705379 +mjd 60143.42156961110595148 +altitude 53.16185928082865786 +azimuth 114.38574818197552929 +filter 2 +rotskypos 131.23506207988341998 +dist2moon 88.2886136 +moonalt -28.7317781 +moondec 24.6999775 +moonphase 0.5804871 +moonra 127.0562204 +nsnap 1 +obshistid 398414 +rottelpos 40.0386345 +seed 398414 +seeing 0.7501810 +sunalt -19.2243392 +vistime 30.0000000 +seqnum 0 +object 32166149840900 60.5808675 -38.07410111 13.8487927 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.0112685 3.1 +object 32166139536388 60.55213375 -38.04768278 13.2693929 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.01014165 3.1 +object 32166139536388 60.55989542 -38.038889 13.2693929 starSED/phoSimMLT/lte035-4.5-1.0a+0.4.BT-Settl.spec.gz 0 0 0 0 0 0 point none CCM 0.01014165 3.1 \ No newline at end of file diff --git a/examples/imsim-user-instcat-ffp-pass1.yaml b/examples/imsim-user-instcat-ffp-pass1.yaml new file mode 100644 index 00000000..e8860e30 --- /dev/null +++ b/examples/imsim-user-instcat-ffp-pass1.yaml @@ -0,0 +1,56 @@ +# Use imSim custom modules +modules: + - imsim + +# Get most of the configuration from the imSim config-template +# for instance catalogs. +template: imsim-config-instcat + +################################################################ +# Make your changes below. +################################################################ + +# Put your own commands that override the defaults below here. For example +# input.instance_catalog.file_name: ./imsim_cat_197356.txt +# input.instance_catalog.sort_mag: False +# input.tree_rings.only_dets: [R22_S11] +# image.nobjects: 5 + +input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' +input.instance_catalog.sort_mag: False + +# For now, disable tree rings to make the diffraction spikes +# easier to see in the images. +input.tree_rings: "" +image.sensor.treering_center: "" +image.sensor.treering_func: "" + +input.atm_psf: "" +psf: + type: Convolve + items: + - + type: Gaussian + fwhm: 0.3 + +image.nobjects: 10 + +stamp.fft_sb_thresh: 1.e30 + +output.det_num.first: 94 +output.nfiles: 2 +output.nproc: 1 + +output.dir: output + +output.off_detector_photons: + dir: output + file_name: + type: FormattedStr + format : off_det_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" diff --git a/examples/imsim-user-instcat-ffp-pass2.yaml b/examples/imsim-user-instcat-ffp-pass2.yaml new file mode 100644 index 00000000..6b136b39 --- /dev/null +++ b/examples/imsim-user-instcat-ffp-pass2.yaml @@ -0,0 +1,81 @@ +# Use imSim custom modules +modules: + - imsim + +# Get most of the configuration from the imSim config-template +# for instance catalogs. +template: imsim-config-instcat + +################################################################ +# Make your changes below. +################################################################ + +input.checkpoint: "" +# input.instance_catalog: "" +input.instance_catalog.file_name: $os.environ.get('IMSIM_HOME')+'/imSim/examples/example_instance_catalog_ffp.txt' +input.instance_catalog.sort_mag: False +# input.opsim_data: "" + +input.tree_rings: "" +input.atm_psf: "" + +# In the second pass, we add off-detector photons to the image generated in the +# first pass from the the on-detector photons. Give the image here. +input.initial_image: + dir: output + file_name: + type: FormattedStr + format : eimage_%08d-%1d-%s-%s-det%03d.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + - $det_name + - "@output.det_num" + +# Two methods can be used to read files containing off-detector photons. If a +# single file only is to be read, then give the file_name. Alternatively, we can +# give a file_pattern containing the fields DETNAME and DETNUM, which are +# wildcards for the detector name and number in an expected set of files +# matching the pattern. The loader will then attempt to read from all files it +# finds matching this patter, excluding the file with the same name/num as the +# current detector. +input.off_detector_photons: + dir: output + # file_name: + # type: FormattedStr + # format: off_det_%08d-%1d-%s-%s-det%03d.fits + # items: + # - { type: OpsimData, field: observationId } + # - { type: OpsimData, field: snap } + # - $band + # - $det_name + # - "@output.det_num" + file_pattern: + type: FormattedStr + format: off_det_%08d-%1d-%s-DETNAME-detDETNUM.fits + items: + - { type: OpsimData, field: observationId } + - { type: OpsimData, field: snap } + - $band + camera: $camera + det_name: $det_name + +psf: + type: Convolve + items: + - + type: Gaussian + fwhm: 0.3 + +stamp: "" + +image.type: LSST_FocalPlaneImage +image.sensor.treering_center: "" +image.sensor.treering_func: "" + +output.det_num.first: 95 +output.nfiles: 1 +output.file_name.format: "eimage_full_focal_%08d-%1d-%s-%s-det%03d.fits" + +output.dir: output diff --git a/imsim/__init__.py b/imsim/__init__.py index b982ffde..69965189 100644 --- a/imsim/__init__.py +++ b/imsim/__init__.py @@ -40,3 +40,4 @@ from .process_info import * from .table_row import * from .photon_pooling import * +from .full_focal_plane import * diff --git a/imsim/camera.py b/imsim/camera.py index d3cb45dc..0537ad21 100644 --- a/imsim/camera.py +++ b/imsim/camera.py @@ -1,6 +1,8 @@ import os +import re import json from collections import defaultdict +import numpy as np import galsim import lsst.utils from .meta_data import data_dir @@ -214,3 +216,40 @@ def update(self, other): def __getattr__(self, attr): """Provide access to the attributes of the underlying lsst_camera.""" return getattr(self.lsst_camera, attr) + + def get_adjacent_detectors(self, det_name): + """ + Given a science detector's name, return a list of the names of science + detectors within the 3x3 grid around it, including itself. + """ + if det_name not in self: + raise galsim.GalSimValueError("Detector is not in camera", det_name) + match = re.match(r"R(\d{2})_S(\d{2})", det_name) + if not match: + # If the given detector is in the list but couldn't be matched, then + # it's a wavefront or guide detector. + raise galsim.GalSimValueError( + "Arg must be a science detector, not wavefront or guide", det_name) + r_det = (int(match.group(1)[0]), int(match.group(1)[1])) + s_det = (int(match.group(2)[0]), int(match.group(2)[1])) + + # The following method of manipulating the indices in the detector name + # is about 5000 times faster than using lsst.afw.cameraGeom to find the + # detector's center and then the names of those shifted away by one + # detector width/height. + + # We'll always have a full set of the nine S index pairs, but they need + # to be permuted to match the position within the raft. + s_rows = np.roll(np.array([0, 1, 2]), 1 - s_det[0]) + s_cols = np.roll(np.array([0, 1, 2]), 1 - s_det[1]) + + # A little integer arithmetic is used to increment or decrement the raft + # indices if the detector is positioned on an edge of the raft. + r_rows = np.array([r_det[0]-s_rows[0]//2, r_det[0], r_det[0]+(2-s_rows[2])//2]) + r_cols = np.array([r_det[1]-s_cols[0]//2, r_det[1], r_det[1]+(2-s_cols[2])//2]) + + adjacent_detectors = [f"R{ri}{rj}_S{si}{sj}" + for ri, si in zip(r_rows, s_rows) + for rj, sj in zip(r_cols, s_cols) + if f"R{ri}{rj}_S{si}{sj}" in self] + return adjacent_detectors diff --git a/imsim/full_focal_plane.py b/imsim/full_focal_plane.py new file mode 100644 index 00000000..ddbb0b0b --- /dev/null +++ b/imsim/full_focal_plane.py @@ -0,0 +1,238 @@ +import os + +import galsim +from galsim.config.extra import ExtraOutputBuilder, RegisterExtraOutput +from galsim.config import RegisterImageType, InputLoader, RegisterInputType +from galsim.sensor import Sensor +from astropy.io.fits import BinTableHDU + +from .lsst_image import LSST_ImageBuilderBase +from .utils import pixel_to_focal, focal_to_pixel +from .camera import get_camera + +def gather_out_of_bounds_photons(image_bounds, photons): + """ Given image bounds and a list of photon arrays, gather the photons + that fall outside the bounds, copy them into a new PhotonArray, and return + them. + + Parameters: + image_bounds: galsim.Bounds + The bounds used to determine which photons fall outside. + photons: List(galsim.PhotonArray) + A list of PhotonArrays containing all the photons to check. + + Returns: + out_of_bounds_photons: galsim.PhotonArray + A new PhotonArray containing the photons falling outside the bounds. + If no such photons are found, the PhotonArray will be empty (N=0). + """ + out_of_bounds_indices = [i for i in range(len(photons)) + if not image_bounds.includes(photons.x[i], photons.y[i])] + if len(out_of_bounds_indices) > 0: + out_of_bounds_photons = galsim.PhotonArray(len(out_of_bounds_indices)) + out_of_bounds_photons.copyFrom(photons, + target_indices=slice(len(out_of_bounds_indices)), + source_indices=out_of_bounds_indices, + do_xy=True, + do_flux=True, + do_other=False) + else: + out_of_bounds_photons = galsim.PhotonArray(N=0) + return out_of_bounds_photons + +# An extra output type to enable two-pass drawing of off-detector photons. This +# output should be used in the first pass to write off-detector photons to file, +# while the second pass will read the off-detector photons then go through all +# the images again and draw the photons on top of the first pass image. +class OffDetectorPhotonsBuilder(ExtraOutputBuilder): + """Build photon arrays containing the off-detector photons found during an + image's construction and write them to file. This item will typically be + included in the config for the first pass of a full focal plane run. + """ + + def processImage(self, index, obj_nums, config, base, logger): + """After each image is drawn, concatenate the list of off-detector + photon arrays found from each batch/sub-batch and store in data ready + to be written to file later on. + """ + if 'off_detector_photons' in base and len(base['off_detector_photons']) > 0: + self.data[index] = galsim.PhotonArray.concatenate(base['off_detector_photons']) + # We need to store the photons using focal plane coordinates so they + # can be accumulated in the second pass on an arbitrary sensor. + detector = get_camera(base['output']['camera'])[base['det_name']] + self.data[index].x, self.data[index].y = pixel_to_focal(self.data[index].x, self.data[index].y, detector) + else: + # If we've been told to write off-detector photons but found none, + # let's at least write an empty PhotonArray. + self.data[index] = galsim.PhotonArray(N=0) + + def finalize(self, config, base, main_data, logger): + return self.data + + def writeFile(self, file_name, config, base, logger): + """Write the photon array to fits. + """ + for photon_array in self.final_data: + photon_array.write(file_name) + + def writeHdu(self, config, base, logger): + """We don't want to write the off-detector photons to FITS, so return an + empty BinTable in lieu of something else. + """ + return BinTableHDU(data=None) + + +class OffDetectorPhotons(object): + """A class to hold the photons which fell outside the sensor being drawn + by the task that createed them during the first pass. They were saved to file + then, and will now be read in this second pass to be accumulated on other sensors. + """ + + def __init__(self, camera, det_name, file_name=None, file_pattern=None, photons=None, dir=None, logger=None): + """ + Initialize the off-detector photons input class. + + Parameters: + camera: str + The name of the camera containing the detector. + det_name: str + The name of the detector to use for photon coordinate transformations. + file_name: str + The name of the file to read. (default: None) + file_pattern: str + The file pattern to use if up to 189 are to be read. (default: None) + dir: str + The directory where the file is. (default: None) + logger: logging.logger + A logger object. (default: None) + """ + logger = galsim.config.LoggerWrapper(logger) + camera = get_camera(camera) + self.det_name = det_name + self.det = camera[det_name] + # Lazy read the off-detector photons when they're accessed for the first time. + self.photons = None + # The Loader should have ensured we have one but not both of file_name + # or file_pattern from the config, but here we also allow the photons to + # be set directly during initialisation in case when instantiating in code. + # So, make sure again that two of these are None. + if sum([file_name is None, file_pattern is None, photons is None]) != 2: + raise galsim.GalSimIncompatibleValuesError("OffDetectorPhotons must be initialized with exactly one of file_name, file_pattern, or photons.",values={"file_name": file_name, "file_pattern": file_pattern, "photons": photons}) + + if file_name is not None: + # Only have a single file to read. + file_names = [file_name] + elif file_pattern is not None: + # We want to read up to 188 files: potentially photons from every + # detector except this one (as those photons were drawn normally in + # the first pass). We may want change this to allow det_num_start + # and det_num_end to be set in the config in a similar way to output + # where it's controlled with det_num.first, nfiles and only_dets. + det_num_start = 0 + det_num_end = 189 + file_names = [file_pattern.replace("DETNAME", camera[i].getName()).replace("DETNUM", "{:03d}".format(i)) + for i in range(det_num_start, det_num_end) + if camera[i].getName() != det_name + ] + else: + # An OffDetectorPhotons object can be be directly initialized in + # code with a photon array instead of reading from file. + file_names = [] + self.photons = photons + if dir is not None: + file_names = [os.path.join(dir, fname) for fname in file_names] + # Allow for the possibility that not all the 188 files exist. + self.file_names = [fname for fname in file_names if os.path.isfile(fname)] + if len(self.file_names) == 0: + logger.warning(f"Detector {det_name} did not find any expected off detector photon files: {file_names}. No photons can be read from file.") + + @property + def photons(self): + if self._photons is None: + self._photons = self._read_photons() + return self._photons + + @photons.setter + def photons(self, photon_array): + # Allow the photons to be set manually from outside the class if necessary. + self._photons = photon_array + + def _read_photons(self, logger=None): + logger = galsim.config.LoggerWrapper(logger) + if len(self.file_names) > 0: + photon_storage = [] + logger.info(f"Detector {self.det_name} reading {len(self.file_names)} files for off-detector photons.") + for file_name in self.file_names: + photon_storage.append(galsim.PhotonArray.read(file_name)) + photons = galsim.PhotonArray.concatenate(photon_storage) + photons.x, photons.y = focal_to_pixel(photons.x, photons.y, self.det) + else: + # If we've reached this point, we've initialised the class without + # providing the names of any usable files and have now tried to + # access the photons without setting them from outside. + # Raise an error. + raise galsim.GalSimConfigError(f"No readable files provided for OffDetectorPhotons for {self.det_name}; photons cannot be read.") + return photons + + +class OffDetectorPhotonsLoader(InputLoader): + """Class to load off-detector photons from file. + These options should be provided in the config for the second pass of a full + focal plane run. + """ + def getKwargs(self, config, base, logger): + req = {'camera': str, 'det_name': str} + opt = {'file_name': str, 'file_pattern': str, 'dir': str} + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, + opt=opt) + if ('file_name' in kwargs) == ('file_pattern' in kwargs): + raise galsim.GalSimConfigError('off_detector_photons must contain one of either file_name or file_pattern in config') + if 'file_pattern' in kwargs and ("DETNAME" not in kwargs['file_pattern'] or "DETNUM" not in kwargs['file_pattern']): + raise galsim.GalSimConfigError('file_pattern must contain both "DETNAME" and "DETNUM" placeholders to be replaced ' + 'with the detector name and number when reading off-detector photon files.') + kwargs['logger'] = logger + + # Base assumption (for now) is that we have a per-detector set of inputs + # that cannot be reused for other detectors. + safe = False + return kwargs, safe + + +# A class to build one of the images making up a full focal plane, used in a +# second pass which uses as input images from the first pass along with any +# photons which fell off-detector. +class LSST_FocalPlaneImageBuilder(LSST_ImageBuilderBase): + + def setup(self, config, base, image_num, obj_num, ignore, logger): + xsize, ysize = super().setup(config, base, image_num, obj_num, ignore, logger) + # Disable application of the addNoise function in + # LSST_ImageBuilderBase.addNoise since it was already applied + # in pass1 and not needed for the added off-detector photons + self.add_noise = False + return xsize, ysize + + def buildImage(self, config, base, image_num, obj_num, logger): + """Draw the off-detector photons to the image. + """ + # Make sure we have an input image and off-detector photons to draw. + # Without them, this type of image won't work. + image = base['current_image'] + if not isinstance(base['_input_objs']['off_detector_photons'][0], OffDetectorPhotons): + raise galsim.config.GalSimConfigError( + "When using LSST_FocalPlaneImage, you must provide an off_detector_photons input.",) + + off_detector_photons = base['_input_objs']['off_detector_photons'][0].photons + + if len(off_detector_photons) > 0: + # There are photons to be accumulated. They should already have been + # transformed to pixel coordinates when they were read from file, so + # go ahead and accumulate them. + sensor = base.get('sensor', Sensor()) + sensor.accumulate(off_detector_photons, image) + + return image, [] + + +RegisterExtraOutput('off_detector_photons', OffDetectorPhotonsBuilder()) +RegisterInputType('off_detector_photons', OffDetectorPhotonsLoader(OffDetectorPhotons)) +RegisterImageType('LSST_FocalPlaneImage', LSST_FocalPlaneImageBuilder()) diff --git a/imsim/instcat.py b/imsim/instcat.py index ff6b006e..34e66236 100644 --- a/imsim/instcat.py +++ b/imsim/instcat.py @@ -13,7 +13,11 @@ from galsim import CelestialCoord import galsim import pickle -from .utils import RUBIN_AREA +from .camera import Camera +from .utils import RUBIN_AREA, focal_to_pixel +import lsst.afw.cameraGeom as cameraGeom + +from copy import deepcopy def clarify_radec_limits( @@ -168,9 +172,12 @@ class InstCatalog(object): # cached SEDs. fnu = (0 * u.ABmag).to(u.erg/u.s/u.cm**2/u.Hz) _flux_density = fnu.to_value(u.ph/u.nm/u.s/u.cm**2, u.spectral_density(500*u.nm)) - def __init__(self, file_name, wcs, xsize=4096, ysize=4096, sed_dir=None, - edge_pix=100, sort_mag=True, flip_g2=True, + def __init__(self, file_name, wcs, fiducial_wcs=None, fiducial_det_name=None, + xsize=4096, ysize=4096, sed_dir=None, + edge_pix=200, camera=None, det_name=None, + sort_mag=True, flip_g2=True, pupil_area=RUBIN_AREA, min_source=None, skip_invalid=True, + full_focal_plane=False, logger=None): logger = galsim.config.LoggerWrapper(logger) self.file_name = file_name @@ -178,6 +185,16 @@ def __init__(self, file_name, wcs, xsize=4096, ysize=4096, sed_dir=None, self.xsize = xsize self.ysize = ysize self.edge_pix = edge_pix + if full_focal_plane: + if camera is None or det_name is None or fiducial_wcs is None: + raise galsim.GalSimConfigValueError("full_focal_plane is set to True: camera, det_name and fiducial_wcs must be provided.") + self.full_focal_plane = True + self.camera = Camera(camera_class=camera) + self.det_name = det_name + self.fiducial_wcs = fiducial_wcs + self.fiducial_det_name = fiducial_det_name + else: + self.full_focal_plane = False self.sort_mag = sort_mag self.flip_g2 = flip_g2 self.min_source = min_source @@ -227,6 +244,19 @@ def _load(self, logger=None): nuse = 0 ntot = 0 + if self.full_focal_plane: + # Get adjacent detectors; then their centers in pixel coordinates, then + # store after converting to RA and dec. + adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) + detector_centers = {} + transform_mapping = self.camera.lsst_camera[self.fiducial_det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() + + # Get the focal plane centers for each detector in the list. + for detector in adjacent_detectors: + focal_plane_center = self.camera.lsst_camera[detector].getCenter(cameraGeom.FOCAL_PLANE) + pixel_center = galsim.PositionD(transform_mapping.applyForward(focal_plane_center)) + detector_centers[detector] = self.fiducial_wcs.toWorld(pixel_center) + with fopen(self.file_name, mode='rt') as _input: for line in _input: if ' inf ' in line: continue @@ -244,8 +274,16 @@ def _load(self, logger=None): and min_dec <= dec <= max_dec ): continue + world_pos = galsim.CelestialCoord(ra, dec) #logger.debug('world_pos = %s',world_pos) + if self.full_focal_plane: + # Check whether the object is closer to this detector's center + # than any of the adjacent detectors. If not, then skip it. + obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} + closest_detector = min(obj_dist, key=obj_dist.get) + if closest_detector != self.det_name: + continue try: image_pos = self.wcs.toImage(world_pos) except RuntimeError as e: @@ -640,12 +678,29 @@ def getKwargs(self, config, base, logger): 'pupil_area' : float, 'min_source' : int, 'skip_invalid' : bool, + 'full_focal_plane' : bool, } - kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt) + ignore = ['fiducial_wcs'] # Handled separately. + kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt, ignore=ignore) + if config.get('fiducial_wcs', None) is not None: + fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) + if fiducial_wcs.wcs_type == "TAN-SIP": + # If it's a TAN-SIP WCS (such as from Batoid) then disable the + # SIP distortions, since we'll potentially use this WCS a long + # way from the detector where the polynomial was fitted. + fiducial_wcs.ab = None + kwargs['fiducial_wcs'] = fiducial_wcs + kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] + else: + kwargs['fiducial_wcs'] = None + # Must build the 'real' WCS after the fiducial WCS as BatoidWCSBuilder.buildWCS() + # stores base['_icrf_to_field'], and we want that to be the 'real' one. wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) kwargs['wcs'] = wcs kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) + kwargs['camera'] = base['output'].get('camera', None) + kwargs['det_name'] = base.get('det_name', None) kwargs['logger'] = logger return kwargs, False diff --git a/imsim/lsst_image.py b/imsim/lsst_image.py index da85cf99..b936f266 100644 --- a/imsim/lsst_image.py +++ b/imsim/lsst_image.py @@ -88,6 +88,8 @@ def setup(self, config, base, image_num, obj_num, ignore, logger): else: self.boresight = params.get('boresight') + self.add_noise = True + self.camera_name = params.get('camera', 'LsstCamSim') # If using a SiliconSensor and not overridden in config, determine sensor type to use. @@ -197,7 +199,8 @@ def addNoise(self, image, config, base, image_num, obj_num, current_var, logger) sky.array[:] *= fringing_map image += sky - AddNoise(base,image,current_var,logger) + if self.add_noise: + AddNoise(base,image,current_var,logger) @staticmethod def _create_full_image(config, base): diff --git a/imsim/photon_pooling.py b/imsim/photon_pooling.py index 04001169..62d0d467 100644 --- a/imsim/photon_pooling.py +++ b/imsim/photon_pooling.py @@ -13,6 +13,7 @@ from .stamp import ProcessingMode, ObjectInfo, build_obj from .lsst_image import LSST_ImageBuilderBase +from .full_focal_plane import gather_out_of_bounds_photons class LSST_PhotonPoolingImageBuilder(LSST_ImageBuilderBase): @@ -138,6 +139,8 @@ def buildImage(self, config, base, image_num, _obj_num, logger): local_wcs = base['wcs'].local(full_image.true_center) if sensor is None: sensor = Sensor() + if 'off_detector_photons' in base['output']: + base['off_detector_photons'] = [] # Initialize the off_detector_photons list for batch_num, batch in enumerate(phot_batches, start=current_photon_batch_num): if not batch: continue @@ -154,9 +157,15 @@ def buildImage(self, config, base, image_num, _obj_num, logger): for op in photon_ops: op.applyTo(photons, local_wcs, rng) + # Gather off-detector photons if we're going to be outputting them + # (likely to draw them on the other sensors in a second pass). + if 'off_detector_photons' in base['output']: + base['off_detector_photons'].append(gather_out_of_bounds_photons(full_image.bounds, photons)) + # Now accumulate the photons onto the sensor. Resume is true for all calls but the first. Recalculate the pixel # boundaries on the first subbatch of each full batch. self.accumulate_photons(photons, full_image, sensor, resume=(batch_num > current_photon_batch_num or subbatch_num > 0), recalc=(subbatch_num == 0)) + del photons # As with the stamps above, let the garbage collector know we don't need the photons anymore. # Note: in typical imsim usage, all current_vars will be 0. So this normally doesn't diff --git a/imsim/skycat.py b/imsim/skycat.py index 054af0b5..c8c511f4 100644 --- a/imsim/skycat.py +++ b/imsim/skycat.py @@ -8,15 +8,20 @@ RegisterObjectType from skycatalogs import skyCatalogs from skycatalogs.utils import PolygonalRegion +from .camera import Camera from .utils import RUBIN_AREA +import lsst.afw.cameraGeom as cameraGeom class SkyCatalogInterface: """Interface to skyCatalogs package.""" - def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, - obj_types=None, skycatalog_root=None, edge_pix=100, + def __init__(self, file_name, wcs, band, mjd, 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, pupil_area=RUBIN_AREA, max_flux=None, logger=None, - apply_dc2_dilation=False, approx_nobjects=None): + apply_dc2_dilation=False, approx_nobjects=None, + full_focal_plane=False): """ Parameters ---------- @@ -28,6 +33,11 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, LSST band, one of ugrizy mjd : float MJD of the midpoint of the exposure. + fiducial_wcs : galsim.WCS [None] + The fiducial WCS, used across all detectors to determine which + should handle which objects. + fiducial_det_name : str [None] + The name of the detector used to create the fiducial WCS. xsize : int Size in pixels of CCD in x-direction. ysize : int @@ -39,9 +49,13 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, Root directory containing skyCatalogs files. If None, then set skycatalog_root to the same directory containing the yaml config file. - edge_pix : float [100] + edge_pix : float [200] Size in pixels of the buffer region around nominal image to consider objects. + camera : str [None] + Name of the imSim camera to use for a full_focal_plane run. + det_name : str [None] + Name of the current detector. pupil_area : float [RUBIN_AREA] The area of the telescope pupil in cm**2. The default value uses R_outer=418 cm and R_inner=255 cm. @@ -60,6 +74,8 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, Approximate number of objects per CCD used by galsim to set up the image processing. If None, then the actual number of objects found by skyCatalogs, via .getNObjects, will be used. + full_focal_plane : bool [False] + Enables or disables two-pass full focal plane processing. """ self.file_name = file_name self.wcs = wcs @@ -73,6 +89,16 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, else: self.skycatalog_root = skycatalog_root self.edge_pix = edge_pix + if full_focal_plane: + if camera is None or det_name is None or fiducial_wcs is None: + raise galsim.GalSimConfigValueError("full_focal_plane is set to True: camera, det_name and fiducial_wcs must be provided.") + self.full_focal_plane = True + self.camera = Camera(camera_class=camera) + self.det_name = det_name + self.fiducial_wcs = fiducial_wcs + self.fiducial_det_name = fiducial_det_name + else: + self.full_focal_plane = False self.pupil_area = pupil_area self.max_flux = max_flux self.logger = galsim.config.LoggerWrapper(logger) @@ -83,10 +109,24 @@ def __init__(self, file_name, wcs, band, mjd, xsize=4096, ysize=4096, self.logger.info(f'Object types restricted to {obj_types}') self.ccd_center = wcs.toWorld(galsim.PositionD(xsize/2.0, ysize/2.0)) self._objects = None + self._objects_to_skip = set() @property def objects(self): if self._objects is None: + if self.full_focal_plane: + # Get adjacent detectors; then their centers in pixel coordinates, then + # store after converting to RA and dec. + adjacent_detectors = self.camera.get_adjacent_detectors(self.det_name) + detector_centers = {} + transform_mapping = self.camera.lsst_camera[self.fiducial_det_name].getTransform(cameraGeom.FOCAL_PLANE, cameraGeom.PIXELS).getMapping() + + # Get the focal plane centers for each detector in the list. + for detector in adjacent_detectors: + focal_plane_center = self.camera.lsst_camera[detector].getCenter(cameraGeom.FOCAL_PLANE) + pixel_center = galsim.PositionD(transform_mapping.applyForward(focal_plane_center)) + detector_centers[detector] = self.fiducial_wcs.toWorld(pixel_center) + # Select objects from polygonal region bounded by CCD edges corners = ((-self.edge_pix, -self.edge_pix), (self.xsize + self.edge_pix, -self.edge_pix), @@ -102,6 +142,20 @@ def objects(self): self.file_name, skycatalog_root=self.skycatalog_root) self._objects = sky_cat.get_objects_by_region( region, obj_type_set=self.obj_types, mjd=self.mjd) + if self.full_focal_plane: + # Work out the which of the adjacent detectors are closest + # to each object found in the region. Add the object to the + # set of objects to skip if the current detector isn't the + # closest one. + for obj in self._objects: + world_pos = galsim.CelestialCoord(obj.ra*galsim.degrees, + obj.dec*galsim.degrees) + obj_dist = {det: detector_centers[det].distanceTo(world_pos) for det in adjacent_detectors} + closest_detector = min(obj_dist, key=obj_dist.get) + if closest_detector != self.det_name: + self._objects_to_skip.add(obj.id) + self.logger.info(f"Skipping {len(self._objects_to_skip)} " + "objects for this detector.") if not self._objects: self.logger.warning("No objects found on image.") return self._objects @@ -165,6 +219,8 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): raise RuntimeError("Trying to get an object from an empty sky catalog") skycat_obj = self.objects[index] + if skycat_obj.id in self._objects_to_skip: + return None gsobjs = skycat_obj.get_gsobject_components(gsparams) if self.apply_dc2_dilation and skycat_obj.object_type == 'galaxy': @@ -217,13 +273,30 @@ def getKwargs(self, config, base, logger): 'approx_nobjects': int, 'mjd': float, 'pupil_area': float, + 'full_focal_plane': bool, } + ignore = ['fiducial_wcs'] # Handled separately. kwargs, safe = galsim.config.GetAllParams(config, base, req=req, - opt=opt) + opt=opt, ignore=ignore) + if config.get('fiducial_wcs', None) is not None: + fiducial_wcs = galsim.config.BuildWCS(config, 'fiducial_wcs', base, logger=logger) + if fiducial_wcs.wcs_type == "TAN-SIP": + # If it's a TAN-SIP WCS (such as from Batoid) then disable the + # SIP distortions, since we'll potentially use this WCS a long + # way from the detector where the polynomial was fitted. + fiducial_wcs.ab = None + kwargs['fiducial_wcs'] = fiducial_wcs + kwargs['fiducial_det_name'] = config['fiducial_wcs']['det_name'] + else: + kwargs['fiducial_wcs'] = None + # Must build the 'real' WCS after the fiducial WCS as BatoidWCSBuilder.buildWCS() + # stores base['_icrf_to_field'], and we want that to be the 'real' one. wcs = galsim.config.BuildWCS(base['image'], 'wcs', base, logger=logger) kwargs['wcs'] = wcs kwargs['xsize'] = base.get('det_xsize', 4096) kwargs['ysize'] = base.get('det_ysize', 4096) + kwargs['camera'] = base['output'].get('camera', None) + kwargs['det_name'] = base.get('det_name', None) kwargs['logger'] = logger # Sky catalog object lists are created per CCD, so they are diff --git a/tests/data/off_det_00398414-0-r-R22_S11-det094.fits b/tests/data/off_det_00398414-0-r-R22_S11-det094.fits new file mode 100644 index 00000000..28abc953 Binary files /dev/null and b/tests/data/off_det_00398414-0-r-R22_S11-det094.fits differ diff --git a/tests/data/off_det_00398414-0-r-R22_S12-det095.fits b/tests/data/off_det_00398414-0-r-R22_S12-det095.fits new file mode 100644 index 00000000..7438632c Binary files /dev/null and b/tests/data/off_det_00398414-0-r-R22_S12-det095.fits differ diff --git a/tests/data/off_det_multi_00398414-0-r-R22_S10-det093.fits b/tests/data/off_det_multi_00398414-0-r-R22_S10-det093.fits new file mode 100644 index 00000000..28abc953 Binary files /dev/null and b/tests/data/off_det_multi_00398414-0-r-R22_S10-det093.fits differ diff --git a/tests/data/off_det_multi_00398414-0-r-R22_S12-det095.fits b/tests/data/off_det_multi_00398414-0-r-R22_S12-det095.fits new file mode 100644 index 00000000..0f3d1779 Binary files /dev/null and b/tests/data/off_det_multi_00398414-0-r-R22_S12-det095.fits differ diff --git a/tests/test_camera.py b/tests/test_camera.py index 47a1a318..ba129fc8 100644 --- a/tests/test_camera.py +++ b/tests/test_camera.py @@ -6,6 +6,7 @@ import json from pathlib import Path import imsim +from galsim import GalSimValueError DATA_DIR = str(Path(__file__).parent.parent / 'data') @@ -43,6 +44,109 @@ def test_bias_levels(self): for amp_name, amp in ccd.items(): self.assertEqual(bias_level, amp.bias_level) + def check_adjacent_detectors_match(self, camera, det_name, expected_detectors): + camera = imsim.Camera(camera) + adjacent_detectors = camera.get_adjacent_detectors(det_name) + self.assertCountEqual(expected_detectors, adjacent_detectors, + f"Adjacent detectors for {det_name} do not match expected values.") + + def test_get_adjacent_detectors_central(self): + # Case of detectors adjacent to one central in a raft. + det_name = 'R22_S11' + expected_detectors = [ + 'R22_S00', 'R22_S01', 'R22_S02', + 'R22_S10', 'R22_S11', 'R22_S12', + 'R22_S20', 'R22_S21', 'R22_S22', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + + def test_get_adjacent_detectors_offset(self): + # Case of detectors adjacent to one offset to a corner in a raft. + # Top left corner. + det_name = 'R22_S20' + expected_detectors = [ + 'R31_S02', 'R32_S00', 'R32_S01', + 'R21_S22', 'R22_S20', 'R22_S21', + 'R21_S12', 'R22_S10', 'R22_S11', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + # Bottom right corner. + det_name = 'R22_S02' + expected_detectors = [ + 'R22_S11', 'R22_S12', 'R23_S10', + 'R22_S01', 'R22_S02', 'R23_S00', + 'R12_S21', 'R12_S22', 'R13_S20', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + + def test_get_adjacent_detectors_edge(self): + # Case of detectors adjacent to one on the edge of the camera. + # Bottom edge. + det_name = 'R02_S01' + expected_detectors = [ + 'R02_S10', 'R02_S11', 'R02_S12', + 'R02_S00', 'R02_S01', 'R02_S02', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + # Left edge. + det_name = 'R20_S10' + expected_detectors = [ + 'R20_S20', 'R20_S21', + 'R20_S10', 'R20_S11', + 'R20_S00', 'R20_S01', + ] + # Top edge. + det_name = 'R42_S21' + expected_detectors = [ + 'R42_S20', 'R42_S21', 'R42_S22', + 'R42_S10', 'R42_S11', 'R42_S12', + ] + # Right edge. + det_name = 'R24_S12' + expected_detectors = [ + 'R24_S21', 'R24_S22', + 'R24_S11', 'R24_S12', + 'R24_S01', 'R24_S02', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + + def test_get_adjacent_detectors_corner(self): + # Case of detectors adjacent to one in a corner of the camera. + # We don't expect wavefront or guider detectors to be returned. + # On the side of a corner cutout. + det_name = 'R01_S10' + expected_detectors = [ + 'R01_S20', 'R01_S21', + 'R01_S10', 'R01_S11', + 'R01_S00', 'R01_S01', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + # On the corner of a corner cutout. + det_name = 'R11_S00' + expected_detectors = [ + 'R10_S12', 'R11_S10', 'R11_S11', + 'R10_S02', 'R11_S00', 'R11_S01', + 'R01_S20', 'R01_S21', + ] + self.check_adjacent_detectors_match('LsstCamSim', det_name, expected_detectors) + + def test_get_adjacent_detectors_exceptions(self): + # Ensure that get_adjacent_detectors raises errors for invalid detecor + # names -- either ones which don't exist, or non science ones. + camera = imsim.Camera('LsstCamSim') + + # Non-existent detectors. + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R99_S99') + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_S00') + + # Non-science detectors. + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_SG0') + with self.assertRaises(GalSimValueError): + camera.get_adjacent_detectors('R00_SW0') + if __name__ == '__main__': unittest.main() diff --git a/tests/test_full_focal_plane.py b/tests/test_full_focal_plane.py new file mode 100644 index 00000000..ca6d68aa --- /dev/null +++ b/tests/test_full_focal_plane.py @@ -0,0 +1,275 @@ +import numpy as np +import galsim + +import os +from pathlib import Path +import pytest + +from imsim.camera import get_camera +from imsim.full_focal_plane import gather_out_of_bounds_photons, OffDetectorPhotons +from imsim.utils import focal_to_pixel + +DATA_DIR = Path(__file__).parent / 'data' +OUTPUT_DIR = Path(__file__).parent / 'output' + +def test_gather_out_of_bounds_photons(): + """Make sure that photons inside and outside the image bounds are gathered correctly.""" + xymin = 0 + xymax = 100 + nphotons = 10 + bounds = galsim.BoundsI(xymin, xymax, xymin, xymax) + + rng = np.random.default_rng(42) + + # Create test photon array + photons = galsim.PhotonArray(nphotons) + photons.flux = np.ones(nphotons) + + # Test all photons with random positions inside the bounds -- none should be gathered. + photons.x = rng.uniform(xymin, xymax, size=nphotons) + photons.y = rng.uniform(xymin, xymax, size=nphotons) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) + assert len(gathered_photons) == 0, "Expected 0 photons to be gathered when all are inside the bounds." + + # Test all photons with random positions outside the bounds -- all should be gathered. + photons.x = rng.uniform(xymax+100, xymax+200, size=nphotons) + photons.y = rng.uniform(xymin+100, xymax+200, size=nphotons) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) + assert len(gathered_photons) == nphotons, "Expected all {} photons to be gathered when they are outside the bounds.".format(nphotons) + + # Test a mix of photons inside and outside the bounds. + n_inside = nphotons // 2 + n_outside = nphotons - n_inside + photons.x[:n_inside] = rng.uniform(xymin, xymax, size=n_inside) + photons.x[n_inside:] = rng.uniform(xymax+100, xymax+200, size=n_outside) + photons.y[:n_inside] = rng.uniform(xymin, xymax, size=n_inside) + photons.y[n_inside:] = rng.uniform(xymin+100, xymax+200, size=n_outside) + gathered_photons = gather_out_of_bounds_photons(bounds, photons) + assert len(gathered_photons) == n_outside, "Of {} photons total, expected the {} outside the bounds to be gathered.".format(nphotons, n_outside) + + +def process_off_detector_photons_output(pa_list): + # Given a list of photon arrays (as would be generated by + # gather_out_of_bounds_photons across an image's photon pools), create a + # config and process it to write the photons to file. + base = { + "det_name": "R22_S11", + "output": { + "camera": "LsstCamSim", + "off_detector_photons": { + "dir": OUTPUT_DIR, + "file_name": "off_det_photons.fits", + }, + }, + "file_num": 0, + "off_detector_photons": pa_list, + } + + # This should now generate a concatenated photon array in + # output/off_det_photons.fits with positions transformed to focal plane + # coordinates. + galsim.config.SetupExtraOutput(base) + galsim.config.ProcessExtraOutputsForImage(base) + galsim.config.WriteExtraOutputs(base, []) + + # Check that the output file exists. + assert os.path.exists(OUTPUT_DIR / "off_det_photons.fits"), "Output file 'output/off_det_photons.fits' was not created." + + # Assert that the photons stored in the output match. They should need to be + # transformed from focal plane to detector coordinates first. + output_photons = galsim.PhotonArray.read(OUTPUT_DIR / "off_det_photons.fits") + if len(output_photons) > 0: + detector = get_camera(base['output']['camera'])[base['det_name']] + output_photons.x, output_photons.y = focal_to_pixel(output_photons.x, output_photons.y, detector) + focal_plane_photons = galsim.PhotonArray.concatenate(pa_list) + assert np.allclose(output_photons.x, focal_plane_photons.x), "X coordinates of output photons do not match original photons in focal plane coordinates." + assert np.allclose(output_photons.y, focal_plane_photons.y), "Y coordinates of output photons do not match original photons in focal plane coordinates." + assert np.allclose(output_photons.flux, focal_plane_photons.flux), "Flux of output photons does not match original photons." + + # Final cleanup. + os.remove(OUTPUT_DIR / "off_det_photons.fits") + + +def test_off_detector_photons_output(): + # Create a list of photon arrays which are all of > 0 length, corresponding + # to the case for which all photon pools found and stored off-detector + # photons. + nphotons_total = 100 + nbatch = 10 + nphot_per_batch = nphotons_total // nbatch + pa_list = [galsim.PhotonArray(nphot_per_batch, + x=np.arange(i*nphot_per_batch, (i+1)*nphot_per_batch), + y=np.arange(i*nphot_per_batch, (i+1)*nphot_per_batch), + flux=np.ones(nphot_per_batch)) + for i in range(nbatch)] + process_off_detector_photons_output(pa_list) + + +def test_off_detector_photons_output_with_none(): + # A list of empty photon arrays should produce output with an empty photon + # array. This is the case if no photons were found outside the image bounds + # for the entire image. + nbatch = 10 + pa_list = [galsim.PhotonArray(N=0) for _ in range(nbatch)] + process_off_detector_photons_output(pa_list) + + +def test_off_detector_photons_output_with_mix(): + # Create a list of photon arrays of different sizes, including some empty + # ones. This is the case when some photon pools contained off-detector + # photons and some didn't. + rng = np.random.default_rng(42) + n_off_det = [0, 10, 20, 15, 0, 10, 0, 15, 0, 20, 10, 0] + xymin = 0.0 + xymax = 100.0 + pa_list = [galsim.PhotonArray(nphot, + x=rng.uniform(xymin, xymax, size=nphot), + y=rng.uniform(xymin, xymax, size=nphot), + flux=np.ones(nphot)) + for nphot in n_off_det] + process_off_detector_photons_output(pa_list) + + +def create_base_off_detector_config(det_name="R22_S11"): + base = { + "input": { + "off_detector_photons": { + "camera": "LsstCamSim", + "det_name": det_name, + "dir": DATA_DIR, + } + }, + } + return base + + +def process_single_off_detector_load(det_name, file_name, expected_photons): + config = create_base_off_detector_config(det_name) + config["input"]["off_detector_photons"]["file_name"] = file_name + galsim.config.ProcessInput(config) + photons = config['_input_objs']['off_detector_photons'][0].photons + np.testing.assert_allclose(photons.x, expected_photons.x, rtol=1.e-10) + np.testing.assert_allclose(photons.y, expected_photons.y, rtol=1.e-10) + np.testing.assert_allclose(photons.flux, expected_photons.flux, rtol=1.e-10) + + +def test_load_photons_from_single_file(): + # Test that the OffDetectorPhotonsLoader can load photons using the 'file_name' + # config option, working as both R22_S11 and R22_S12 to ensure photons + # are transformed from focal plane to detector coordinates. + expected_photons = galsim.PhotonArray.fromArrays(x=np.arange(1000., 1010.), + y=np.arange(2000., 2010.), + flux=np.ones(10), + ) + process_single_off_detector_load("R22_S11", "off_det_00398414-0-r-R22_S11-det094.fits", expected_photons) + expected_photons = galsim.PhotonArray.fromArrays(x=np.arange(-3000., -3010., -1.), + y=np.arange(-4000., -4010., -1), + flux=np.ones(10), + ) + process_single_off_detector_load("R22_S12", "off_det_00398414-0-r-R22_S12-det095.fits", expected_photons) + + +def test_load_photons_from_multiple_files(): + # Test that the OffDetectorPhotonsLoader can load photons from multiple files + # given a file name pattern. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_pattern"] = "off_det_multi_00398414-0-r-DETNAME-detDETNUM.fits" + # Processing this config should read the photons from the two files + # data/off_det_multi_00398414-0-r-R22_S10-det093.fits and + # data/off_det_multi_00398414-0-r-R22_S12-det095.fits, concatenate them, and + # transform them to the R22_S11 detector coordinates. + galsim.config.ProcessInput(config) + photons = config['_input_objs']['off_detector_photons'][0].photons + + # We expect the photons from the two files to have been transformed to the + # following in the R22_S11 detector coordinates. + x = np.concatenate((np.arange(1000., 1010.), np.arange(5000., 5010.))) + y = np.concatenate((np.arange(2000., 2010.), np.arange(10000., 10010.))) + expected_photons = galsim.PhotonArray.fromArrays(x=x, + y=y, + flux=np.ones(20), + ) + np.testing.assert_allclose(photons.x, expected_photons.x, rtol=1.e-10) + np.testing.assert_allclose(photons.y, expected_photons.y, rtol=1.e-10) + np.testing.assert_allclose(photons.flux, expected_photons.flux, rtol=1.e-10) + + +def test_off_detector_loader_file_naming(): + # Make sure that OffDetectorPhotonsLoader raises an error if neither a file + # name nor a file pattern is provided. + config = create_base_off_detector_config() + with pytest.raises(galsim.GalSimConfigError, + match='off_detector_photons must contain one of either file_name or file_pattern in config'): + galsim.config.ProcessInput(config) + + # Make sure the same error is raised if both are provided. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_name"] = "some_file_name.fits" + config["input"]["off_detector_photons"]["file_pattern"] = "some_file_pattern.fits" + with pytest.raises(galsim.GalSimConfigError, + match='off_detector_photons must contain one of either file_name or file_pattern in config'): + galsim.config.ProcessInput(config) + + # Now ensure that, if a file pattern is provided, the config includes the + # DETNAME and DETNUM fields that are used to read the range of files. + config = create_base_off_detector_config() + config["input"]["off_detector_photons"]["file_pattern"] = "some_file_pattern.fits" + with pytest.raises(galsim.GalSimConfigError, + match='file_pattern must contain both "DETNAME" and "DETNUM" placeholders to be replaced ' + 'with the detector name and number when reading off-detector photon files.'): + galsim.config.ProcessInput(config) + + +def test_load_photons_raises_error_with_no_files(): + # Photons are read lazily, so check that an error is raised when we try to access them + # but the files they are to be read from don't exist. + + # Firstly check the case when a single non-existent file_name is provided. + with pytest.raises(galsim.GalSimConfigError): + off_detector_photons_input = OffDetectorPhotons("LsstCamSim", "R22_S11", file_name="non_existent_file.fits") + photons = off_detector_photons_input.photons + + # The same but with a file_pattern that doesn't match any existing files. + with pytest.raises(galsim.GalSimConfigError): + off_detector_photons_input = OffDetectorPhotons("LsstCamSim", "R22_S11", file_pattern="non_existent_file_pattern_DETNAME_DETNUM.fits") + photons = off_detector_photons_input.photons + + +def test_draw_off_detector_photons(): + # Test that the LSST_FocalPlaneImageBuilder can draw photons loaded by the + # OffDetectorPhotonsLoader on top of the image. + xsize = 10 + ysize = 10 + camera = "LsstCamSim" + det_name = "R22_S11" + # Manually create some off-detector photons in the base config to be drawn + # on the image. We expect a diagonal line in the resulting image. + photons = galsim.PhotonArray.fromArrays(x=np.arange(1., 11.), + y=np.arange(1., 11.), + flux=np.ones(10), + ) + # Assume that an existing image as well as off_detector_photons from the + # first pass have been read from file and are at this point available in the + # base config. + base = { + "image": { + "xsize": xsize, + "ysize": ysize, + "det_name": det_name, + "nobjects": 1, + "type": "LSST_FocalPlaneImage", + }, + "current_image": galsim.Image(xsize, ysize), + '_input_objs': { + 'off_detector_photons': [OffDetectorPhotons(camera, det_name, photons=photons)], + } + } + expected = np.diag(np.ones(10)) + image = galsim.config.BuildImage(base) + np.testing.assert_allclose(image.array, expected, rtol=1.e-10) + + +if __name__ == "__main__": + testfns = [v for k, v in vars().items() if k.startswith("test_") and callable(v)] + for testfn in testfns: + testfn() diff --git a/tests/test_lsst_image.py b/tests/test_lsst_image.py index 8f59a353..9fdf1c28 100644 --- a/tests/test_lsst_image.py +++ b/tests/test_lsst_image.py @@ -35,6 +35,7 @@ def test_image_nobjects(): config = {'modules': ['imsim'], 'template': template, 'input.instance_catalog.file_name': instcat_file, + 'input.instance_catalog.edge_pix': 100, 'input.opsim_data.file_name': instcat_file, 'input.tree_rings': '', 'input.atm_psf': '', diff --git a/tests/test_stamp.py b/tests/test_stamp.py index 9b691b5d..60a56236 100644 --- a/tests/test_stamp.py +++ b/tests/test_stamp.py @@ -204,6 +204,7 @@ def test_stamp_sizes(): input.sky_catalog.pupil_area: type: Eval str: "0.25 * np.pi * 649**2" + input.sky_catalog.edge_pix: 100 input.opsim_data.file_name: data/small_opsim_9683.db input.opsim_data.visit: 449053 input.tree_rings.only_dets: [R22_S11] @@ -426,6 +427,7 @@ def test_faint_high_redshift_stamp(): input.sky_catalog.file_name: data/sky_cat_9683.yaml input.sky_catalog.obj_types: [galaxy] input.sky_catalog.band: u + input.sky_catalog.edge_pix: 100 input.opsim_data.file_name: data/small_opsim_9683.db input.opsim_data.visit: 449053 input.tree_rings.only_dets: [R22_S11] diff --git a/tests/test_trimmer.py b/tests/test_trimmer.py index 698b9e23..13faf4e2 100644 --- a/tests/test_trimmer.py +++ b/tests/test_trimmer.py @@ -45,8 +45,8 @@ def test_InstCatTrimmer(self): instcat = imsim.InstCatalog(instcat_file, wcs, edge_pix=1000, skip_invalid=False) self.assertEqual(instcat.nobjects, 24) - # With the default edge_pix=100, only 17 make the cut. - instcat = imsim.InstCatalog(instcat_file, wcs, skip_invalid=False) + # With a smaller than default edge_pix=100, only 17 make the cut. + instcat = imsim.InstCatalog(instcat_file, wcs, edge_pix=100, skip_invalid=False) self.assertEqual(instcat.nobjects, 17) # Check the application of min_source.