Add simplified IMCOM coadd simulation

roman_imsim/__init__.py

@@ -23,6 +23,7 @@
 from .psf import *
 from .sca import *
 from .skycat import *
+from .roman_coadd import *
 from .stamp import *
 from .wcs import *
 

roman_imsim/config/sim_coadd.yaml

@@ -0,0 +1,189 @@
+# Default settings for roman simulation
+# Includes creation of noisless oversampled images (including PSF)
+#  -- processing of other detector and instrument effects are still handled in the
+#     python postprocessing layer to enable things not currently in galsim.roman
+
+modules:
+
+    # Including galsim.roman in the list of modules to import will add a number of Roman-specific
+    # functions and classes that we will use here.
+    - roman_imsim
+    - galsim.roman
+
+    # We need this for one of our Eval items.  GalSim does not by default import datetime into
+    # the globals dict it uses when evaluating Eval items, so we can tell it to import it here.
+    - datetime
+
+# Define some other information about the images
+image:
+
+    # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc.
+    # It also by default applies a number of detector effects, but these can be turned
+    # off if desired by setting some parameters (given below) to False.
+    type: roman_coadd
+    coadd_file: '/hpc/group/cosmology/yuedong/H158_Row00/prod_H_00_00_map.fits'
+    xsize: 2688
+    ysize: 2688
+    white_noise_weight: 0.7
+    pink_noise_weight: 0.3
+    pixel_scale: 0.0390625
+
+    wcs:
+        type: RomanWCS
+        # SCA: '@image.SCA'
+        # ra: { type: ObSeqData, field: ra }
+        # dec: { type: ObSeqData, field: dec }
+        # pa: { type: ObSeqData, field: pa }
+        # mjd: { type: ObSeqData, field: mjd }
+        # coadd_file: '/cwork/cmh215/OU24-Roman-Coadds/H158/prod_H_00_00_map.fits.gz'
+
+    bandpass:
+        type: RomanBandpass
+        # name: { type: ObSeqData, field: filter }
+        name: "H158"
+
+    # When you want to have multiple images generate the same random galaxies, then
+    # you can set up multiple random number generators with different update cadences
+    # by making random_seed a list.
+    # The default behavior is just to have the random seeds for each object go in order by
+    # object number across all images, but this shows how to set it up so we use two separate
+    # cadences.
+    # The first one behaves normally, which will be used for things like noise on the image.
+    # The second one sets the initial seed for each object to repeat to the same starting value
+    # at the start of each filter.  If we were doing more than 3 total files, it would then
+    # move on to another sequence for the next 3 and so on.
+    random_seed:
+        # Used for noise and nobjects.
+        - { type: ObSeqData, field: visit }
+
+        # Used for objects.  Repeats sequence for each filter
+        # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then
+        # treated the same way as an integer (i.e. making a regular sequence starting from that
+        # value).  Using an explicit dict with an Eval type means GalSim will leave it alone and
+        # evaluate it as is for each object.
+
+
+    # We're just doing one SCA here.
+    # If you wanted to do all of them in each of three filters (given below), you could use:
+    #
+    # SCA:
+    #     type: Sequence
+    #     first: 1
+    #     last: 18
+    #     repeat: 3  # repeat each SCA num 3 times before moving on, for the 3 filters.
+    #
+
+    SCA: 10
+    mjd: { type: ObSeqData, field: mjd }
+    filter: { type: ObSeqData, field: filter }
+    exptime: { type: ObSeqData, field: exptime }
+
+    draw_method: 'fft'
+    # Photon shooting is way faster for chromatic objects than fft, especially when most of them
+    # are fairly faint.  The cross-over point for achromatic objects is generally of order
+    # flux=1.e6 or so (depending on the profile).  Most of our objects here are much fainter than
+    # that.  The fft rendering for chromatic is a factor of 10 or so slower still, whereas
+    # chromatic photon shooting is only slighly slower than achromatic, so the difference
+    # is even more pronounced in this case.
+    use_fft_bright: True
+
+    # These are all by default turned on, but you can turn any of them off if desired:
+    # ignore_noise: True
+    # stray_light: False
+    # thermal_background: False
+    # reciprocity_failure: False
+    # dark_current: False
+    # nonlinearity: False
+    # ipc: False
+    # read_noise: False
+    # sky_subtract: False
+
+    # ignore_noise: False
+    # stray_light: True
+    # thermal_background: True
+    # reciprocity_failure: True
+    # dark_current: True
+    # nonlinearity: True
+    # ipc: True
+    # read_noise: True
+    # sky_subtract: False
+
+    # nobjects: 500
+
+stamp:
+    type: Roman_stamp
+    world_pos:
+        type: SkyCatWorldPos
+    exptime: { type: ObSeqData, field: exptime }
+    skip_failures: True
+    photon_ops:
+        -
+            type: ChargeDiff
+
+# psf:
+#     type: roman_psf
+#     # If omitted, it would figure this out automatically, because we are using the RomanSCA image
+#     # type.  But if we weren't, you'd have to tell it which SCA to build the PSF for.
+#     SCA: '@image.SCA'
+#     # n_waves defines how finely to sample the PSF profile over the bandpass.
+#     # Using 10 wavelengths usually gives decent accuracy.
+#     n_waves: 10
+
+psf:
+    type: RomanPSF
+    interpolator:
+        type: RomanPSFInterpolator
+
+# Define the galaxy type and positions to use
+gal:
+    type: SkyCatObj
+
+input:
+    obseq_data:
+        file_name: /hpc/group/cosmology/OpenUniverse2024/RomanWAS/Roman_WAS_obseq_11_1_23.fits
+        visit: 4907
+        SCA: '@image.SCA'
+    # roman_psf:
+    #     SCA: '@image.SCA'
+    #     n_waves: 5
+    RomanPSFInterpolator:
+        kind: corners
+        n_waves: 5
+    sky_catalog:
+        file_name: /hpc/home/yf194/Work/test_coadd_sims/roman_imsim/config/skyCatalog.yaml
+        edge_pix: 512
+        mjd: { type: ObSeqData, field: mjd }
+        exptime: { type: ObSeqData, field: exptime }
+        obj_types: ['diffsky_galaxy','star','snana']
+
+output:
+
+    nfiles: 1
+    dir: /hpc/group/cosmology/yuedong/roman_output/coadd_sim_new/images/truth
+    file_name:
+        type: FormattedStr
+        format: "Roman_WAS_truth_%s_%i_%i.fits.gz"
+        items:
+            - { type: ObSeqData, field: filter }
+            - { type: ObSeqData, field: visit }
+            - '@image.SCA'
+
+    truth:
+        dir: /hpc/group/cosmology/yuedong/roman_output/coadd_sim_new/truth
+        file_name:
+            type: FormattedStr
+            format: "Roman_WAS_index_%s_%i_%i.txt"
+            items:
+                - { type: ObSeqData, field: filter }
+                - { type: ObSeqData, field: visit }
+                - '@image.SCA'
+        columns:
+            object_id: "@object_id"
+            ra: "$sky_pos.ra.deg"
+            dec: "$sky_pos.dec.deg"
+            x: "$image_pos.x"
+            y: "$image_pos.y"
+            realized_flux: "@realized_flux"
+            flux: "@flux"
+            mag: "@mag"
+            obj_type: "@object_type"

roman_imsim/psf.py

@@ -83,6 +83,18 @@ def _parse_pupil_bin(self, pupil_bin):
         else:
             return pupil_bin
 
+    def _psf_call_coadd(self, bpass, n_waves, logger):
+        # Currently only implementing a Gaussian PSF for each band
+        fwhm_dict = {
+            "Y106": 0.220,
+            "J129": 0.231,
+            "H158": 0.242,
+            "F184": 0.253,
+            "K213": 0.264,
+        }
+        psf = galsim.Gaussian(fwhm=fwhm_dict[bpass.name])
+        return psf.withGSParams(maximum_fft_size=16384)
+
     def _psf_call(self, SCA, bpass, SCA_pos, WCS, pupil_bin, n_waves, logger, extra_aberrations):
 
         if pupil_bin == 8:
@@ -227,7 +239,9 @@ def initPSF(
                 SCA, bandpass, cc, WCS, pupil_bin, n_waves, logger, self._extra_aberrations
             )
 
-    def getPSF(self, pupil_bin, pos):
+        self.PSF_coadd = self._psf_call_coadd(bandpass, n_waves, logger)
+
+    def getPSF(self, pupil_bin, pos, is_coadd=False):
         """
         Return a PSF to be convolved with sources.
 
@@ -256,6 +270,9 @@ def getPSF(self, pupil_bin, pos):
         #     psf = self.PSF[pupil_bin]['cc']
         # return psf
 
+        if is_coadd:
+            return self.PSF_coadd
+
         psf = self.PSF[pupil_bin]
         if pupil_bin != 8:
             return psf