Update unit tests for the new code installed via setup.py

.github/workflows/ci.yml

@@ -0,0 +1,70 @@
+name: imSim CI
+
+on:
+    push:
+        branches:
+            - main
+            - releases/*
+
+    pull_request:
+        branches:
+            - main
+            - releases/*
+
+jobs:
+    build:
+        runs-on: ${{ matrix.os }}
+
+        strategy:
+            matrix:
+                # For now, just ubuntu, 3.8.  Can add more later.
+                os: [ ubuntu-latest ]
+                py: [ 3.8 ]
+                CC: [ gcc ]
+                CXX: [ g++ ]
+
+        defaults:
+            run:
+                # cf. https://github.com/conda-incubator/setup-miniconda#important
+                shell: bash -l {0}
+
+        steps:
+            - uses: actions/checkout@v2
+
+            - name: Setup conda
+              uses: conda-incubator/setup-miniconda@v2
+              with:
+                  activate-environment: stack
+                  python-version: 3.8
+                  condarc-file: etc/.condarc
+
+            - name: Install conda deps
+              run: |
+                conda info
+                conda list
+                conda install -y mamba
+                mamba install -y --file conda_requirements.txt
+                conda info
+                conda list
+
+            - name: Install pip deps
+              run: |
+                # We need to get batoid onto conda, but for now, this is a separate step.
+                pip install batoid
+                conda info
+                conda list
+
+            - name: Install imSim
+              run:
+                pip install .
+
+            - name: Install test deps
+              run:
+                conda install -y pytest nose
+
+            - name: Run tests
+              run: |
+                cd tests
+                # We're working towards getting all the tests working, but so far these are
+                # the ones that work withe the pip installation.
+                pytest test_FWHMgeom.py test_atmPSF.py test_batoid_wcs.py test_cosmic_rays.py test_fopen.py test_instcat_parser.py test_optical_zernikes.py test_psf.py test_tree_rings.py test_trimmer.py

conda_requirements.txt

@@ -0,0 +1,4 @@
+# conda install --file conda_requirements should install all required dependencies of imSim.
+
+stackvana>=0.2021.30
+galsim>=2.3

etc/.condarc

@@ -0,0 +1,5 @@
+channels:
+  - conda-forge
+  - defaults
+ssl_verify: true
+channel_priority: strict

imsim/atmPSF.py

@@ -6,11 +6,38 @@
 import numpy as np
 from scipy.optimize import bisect
 
+import pickle
 import galsim
 
 from .optical_system import OpticalZernikes, mock_deviations
 
 
+def save_psf(psf, outfile):
+    """
+    Save the psf as a pickle file.
+    """
+    # Set any logger attribute to None since loggers cannot be persisted.
+    if hasattr(psf, 'logger'):
+        psf.logger = None
+    with open(outfile, 'wb') as output:
+        with galsim.utilities.pickle_shared():
+            pickle.dump(psf, output)
+
+def load_psf(psf_file, log_level='INFO'):
+    """
+    Load a psf from a pickle file.
+    """
+    with open(psf_file, 'rb') as fd:
+        psf = pickle.load(fd)
+
+    # Since save_psf sets any logger attribute to None, restore
+    # it here.
+    if hasattr(psf, 'logger'):
+        psf.logger = get_logger(log_level, 'psf')
+
+    return psf
+
+
 class OptWF(object):
     def __init__(self, rng, wavelength, gsparams=None):
         u = galsim.UniformDeviate(rng)
@@ -302,6 +329,7 @@ def __init__(self, airmass, rawSeeing, band, boresight, rng,
         # Instead, the user can choose to convolve this by a Gaussian in the config file.
         self.atm = AtmosphericPSF(airmass, rawSeeing, band, rng,
                                   t0=t0, exptime=exptime, kcrit=kcrit, gaussianFWHM=0.,
+                                  screen_size=screen_size, screen_scale=screen_scale,
                                   doOpt=doOpt, logger=logger, nproc=nproc)
         # The other change is that we need the boresight to do image_pos -> field_pos
         # I think it makes sense to take that as an input here rather than in the
@@ -328,6 +356,102 @@ def BuildAtmosphericPSF(config, base, ignore, gsparams, logger):
     psf = atm.getPSF(field_pos, gsparams)
     return psf, False
 
+# These next two are approximations to the above atmospheric PSF, which might be
+# useful in contexts where accuracy of the PSF isn't so important.
+def BuildDoubleGaussianPSF(config, base, ignore, gsparams, logger):
+    """
+    This is an example implementation of a wavelength- and
+    position-independent Double Gaussian PSF.  See the documentation
+    in PSFbase to learn how it is used.
+
+    This specific PSF comes from equation(30) of the signal-to-noise
+    document (LSE-40), which can be found at
+
+    www.astro.washington.edu/users/ivezic/Astr511/LSST_SNRdoc.pdf
+
+    The required fwhm parameter is the Full Width at Half Max of the total
+    PSF.  This is given in arcseconds.
+    """
+    req = {'fwhm': float}
+    opt = {'pixel_scale': float}
+
+    params, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt)
+    fwhm = params['fwhm']
+    pixel_scale = params.get('pixel_scale', 0.2)
+    if gsparams: gsparams = GSParams(**gsparams)
+    else: gsparams = None
+
+    # the expression below is derived by solving equation (30) of
+    # the signal-to-noise document
+    # (www.astro.washington.edu/uses/ivezic/Astr511/LSST_SNRdoc.pdf)
+    # for r at half the maximum of the PSF
+    alpha = fwhm/2.3835
+
+    eff_pixel_sigma_sq = pixel_scale*pixel_scale/12.0
+
+    sigma = np.sqrt(alpha*alpha - eff_pixel_sigma_sq)
+    gaussian1 = galsim.Gaussian(sigma=sigma, gsparams=gsparams)
+
+    sigma = np.sqrt(4.0*alpha*alpha - eff_pixel_sigma_sq)
+    gaussian2 = galsim.Gaussian(sigma=sigma, gsparams=gsparams)
+
+    psf = 0.909*(gaussian1 + 0.1*gaussian2)
+
+    return psf, safe
+
+
+def BuildKolmogorovPSF(config, base, ignore, gsparams, logger):
+    """
+    This PSF class is based on David Kirkby's presentation to the DESC
+    Survey Simulations working group on 23 March 2017.
+
+    https://confluence.slac.stanford.edu/pages/viewpage.action?spaceKey=LSSTDESC&title=SSim+2017-03-23
+
+    (you will need a SLAC Confluence account to access that link)
+
+    Parameters
+    ----------
+    airmass
+
+    rawSeeing is the FWHM seeing at zenith at 500 nm in arc seconds
+    (provided by OpSim)
+
+    band is the bandpass of the observation [u,g,r,i,z,y]
+    """
+
+    req = {
+            'airmass': float,
+            'rawSeeing': float,
+            'band': str,
+          }
+
+    params, safe = galsim.config.GetAllParams(config, base, req=req)
+    airmass = params['airmass']
+    rawSeeing = params['rawSeeing']
+    band = params['band']
+    if gsparams: gsparams = GSParams(**gsparams)
+    else: gsparams = None
+
+    # This code was provided by David Kirkby in a private communication
+
+    wlen_eff = dict(u=365.49, g=480.03, r=622.20, i=754.06, z=868.21,
+                    y=991.66)[band]
+    # wlen_eff is from Table 2 of LSE-40 (y=y2)
+
+    FWHMatm = rawSeeing*(wlen_eff/500.)**-0.3*airmass**0.6
+    # From LSST-20160 eqn (4.1)
+
+    FWHMsys = np.sqrt(0.25**2 + 0.3**2 + 0.08**2)*airmass**0.6
+    # From LSST-20160 eqn (4.2)
+
+    atm = galsim.Kolmogorov(fwhm=FWHMatm, gsparams=gsparams)
+    sys = galsim.Gaussian(fwhm=FWHMsys, gsparams=gsparams)
+    psf = galsim.Convolve((atm, sys))
+
+    return psf, safe
+
 from galsim.config import InputLoader, RegisterInputType, RegisterObjectType
 RegisterInputType('atm_psf', InputLoader(AtmosphericPSFBuilder, takes_logger=True))
 RegisterObjectType('AtmosphericPSF', BuildAtmosphericPSF, input_type='atm_psf')
+RegisterObjectType('DoubleGaussianPSF', BuildDoubleGaussianPSF)
+RegisterObjectType('KolmogorovPSF', BuildKolmogorovPSF)

imsim/batoid_wcs.py

@@ -476,13 +476,14 @@ def pixel_to_ICRF(self, x, y, det):
 
 
 class BatoidWCSBuilder(WCSBuilder):
+
     def __init__(self):
-        self._camera = None  # It's slow to make a camera instance, so only make it once.
+        self._camera = None
 
     @property
     def camera(self):
         if self._camera is None:
-            self._camera = get_camera(self._camera_class)
+            self._camera = get_camera(self._camera_name)
         return self._camera
 
     def buildWCS(self, config, base, logger):
@@ -497,7 +498,6 @@ def buildWCS(self, config, base, logger):
             the constructed WCS object (a galsim.GSFitsWCS instance)
         """
         req = {
-                "camera": str,
                 "boresight": galsim.CelestialCoord,
                 "rotTelPos": galsim.Angle,
                 "obstime": None,  # Either str or astropy.time.Time instance
@@ -506,6 +506,7 @@ def buildWCS(self, config, base, logger):
                               # become optional, since other telescopes don't use it.
               }
         opt = {
+                "camera": str,
                 "telescope": str,
                 "temperature": float,
                 "pressure": float,
@@ -521,50 +522,89 @@ def buildWCS(self, config, base, logger):
             base['bandpass'] = bp
 
         kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt)
-        self._camera_class = kwargs.pop('camera')
         logger.info("Building Batoid WCS for %s", kwargs['det_name'])
+        kwargs['bandpass'] = base.get('bandpass', None)
+        return self.makeWCS(**kwargs)
+
+    def makeWCS(self, boresight, rotTelPos, obstime, det_name, band, camera='LsstCam',
+                telescope='LSST', temperature=None, pressure=None, H2O_pressure=None,
+                wavelength=None, bandpass=None, order=3):
+        """Make the WCS object given the parameters explicitly rather than via a config dict.
+
+        It mostly just calls the BatoidWCSFactory.getWCS function, but it has sensible defaults
+        for many parameters.
+
+        Parameters
+        ----------
+        boresight : galsim.CelestialCoord
+            The ICRF coordinate of light that reaches the boresight.  Note that this
+            is distinct from the spherical coordinates of the boresight with respect
+            to the ICRF axes.
+        rotTelPos : galsim.Angle
+            Camera rotator angle.
+        obstime : astropy.time.Time or str
+            Mean time of observation.  Note: if this is a string, it is assumed to be in TAI scale,
+            which seems to be standard in the Rubin project.
+        det_name : str
+            Detector name in the format e.g. R22_S11
+        band : str
+            The name of the bandpass
+        telescope : str
+            The name of the telescope. [default: 'LSST']  Currenly only 'LSST' is functional.
+        temperature : float
+            Ambient temperature in Kelvin [default: 280 K]
+        pressure : float
+            Ambient pressure in kPa [default: based on LSST heigh of 2715 meters]
+        H2O_pressure : float
+            Water vapor pressure in kPa [default: 1.0 kPa]
+        wavelength : float
+            Nanometers [default: None, which means use the bandpass effective wavelength]
+        bandpass : galsim.Bandpass
+            If wavelegnth is None, use this to get the effective wavelength. [default: None,
+            which means the default LSST bandpass will be used given the (required) band parameter.
+        order : int
+            The order of the SIP polynomial to use. [default: 3]
+
+        Returns:
+            the constructed WCS object (a galsim.GSFitsWCS instance)
+        """
 
         # If a string, convert it to astropy.time.Time.
-        # XXX Assumption is that string time is in scale 'TAI'.  Should make sure
-        # to make this consistent with OpSim.
-        if isinstance(kwargs['obstime'], str):
-            kwargs['obstime'] = astropy.time.Time(kwargs['obstime'], scale='tai')
+        if isinstance(obstime, str):
+            obstime = astropy.time.Time(obstime, scale='tai')
 
-        telescope = kwargs.pop('telescope', 'LSST')
         if telescope != 'LSST':
             raise NotImplementedError("Batoid WCS only valid for telescope='LSST' currently")
-        band = kwargs.pop('band')
-        kwargs['fiducial_telescope'] = batoid.Optic.fromYaml(f"{telescope}_{band}.yaml")
-        kwargs['camera'] = self.camera
+        fiducial_telescope = batoid.Optic.fromYaml(f"{telescope}_{band}.yaml")
+        self._camera_name = camera
 
         # Update optional kwargs
 
-        if 'wavelength' not in kwargs:
-            kwargs['wavelength'] = base['bandpass'].effective_wavelength
+        if wavelength is None:
+            if bandpass is None:
+                bandpass = galsim.Bandpass('LSST_%s.dat'%band, wave_type='nm')
+            wavelength = bandpass.effective_wavelength
 
-        if 'temperature' not in kwargs:
+        if temperature is None:
             # cf. https://www.meteoblue.com/en/weather/historyclimate/climatemodelled/Cerro+Pachon
             # Average minimum temp is around 45 F = 7 C, but obviously varies a lot.
-            kwargs['temperature'] = 280 # Kelvin
+            temperature = 280 # Kelvin
 
-        if 'pressure' not in kwargs:
+        if pressure is None:
             # cf. https://www.engineeringtoolbox.com/air-altitude-pressure-d_462.html
             # p = 101.325 kPa (1 - 2.25577e-5 (h / 1 m))**5.25588
             # Cerro Pachon  altitude = 2715 m
             h = 2715
-            kwargs['pressure'] = 101.325 * (1-2.25577e-5*h)**5.25588
+            pressure = 101.325 * (1-2.25577e-5*h)**5.25588
 
-        if 'H2O_pressure' not in kwargs:
+        if H2O_pressure is None:
             # I have no idea what a good default is, but this seems like a minor enough effect
             # that we should not require the user to pick something.
-            kwargs['H2O_pressure'] = 1.0 # kPa
+            H2O_pressure = 1.0 # kPa
 
         # Finally, build the WCS.
-        det_name = kwargs.pop('det_name')
-        order = kwargs.pop('order', 3)
-        factory = BatoidWCSFactory(**kwargs)
-        wcs = factory.getWCS(self.camera[det_name], order=order)
-
-        return wcs
+        factory = BatoidWCSFactory(boresight, rotTelPos, obstime, fiducial_telescope,
+                                   wavelength, self.camera, temperature, pressure, H2O_pressure)
+        return factory.getWCS(self.camera[det_name], order=order)
 
 RegisterWCSType('Batoid', BatoidWCSBuilder())