Merge pull request #81 from LSSTDESC/u/troxel/diffsky

U/troxel/diffsky

.github/workflows/ci.yml

@@ -39,12 +39,11 @@ jobs:
 
             - name: Install conda deps
               run: |
+                conda update -n base conda
                 conda info
                 conda list
-                conda install -y mamba
-                mamba install -y --file etc/conda_requirements.txt
+                conda install -y --file etc/conda_requirements.txt
                 conda info
-                mamba list
 
             - name: Install rubin_sim_data
               run: |

cfg/object_types.yaml

@@ -34,6 +34,39 @@ object_types :
       internal_extinction : CCM
       MW_extinction : F19
       spatial_model :  knots
+  diffsky_galaxy :
+      file_template : 'galaxy_(?P<healpix>\d+).parquet'
+      flux_file_template : 'galaxy_flux_(?P<healpix>\d+).parquet'
+      sed_file_template : 'galaxy_sed_(?P<healpix>\d+).hdf5'
+      data_file_type : parquet
+      area_partition :
+        { type: healpix, ordering : ring, nside : 32}
+      composite :
+        bulge : required
+        disk  : required
+        knots : optional
+  diffsky_bulge :
+      subtype : bulge
+      parent : diffsky_galaxy
+      sed_model : TBD
+      internal_extinction : CCM
+      MW_extinction : F19
+      spatial_model : sersic2D
+  diffsky_disk :
+      subtype : disk
+      parent : diffsky_galaxy
+      sed_model : TBD
+      internal_extinction : CCM
+      MW_extinction : F19
+      spatial_model :  sersic2D
+  diffsky_knots :
+      subtype : knots
+      parent : diffsky_galaxy
+      sed_model : TBD
+      internal_extinction : CCM
+      MW_extinction : F19
+      spatial_model :  knots
+
   star :
       file_template : 'pointsource_(?P<healpix>\d+).parquet'
       flux_file_template : 'pointsource_flux_(?P<healpix>\d+).parquet'

skycatalogs/_version.py

@@ -1 +1 @@
-__version__ = "1.7.0-rc1"
+__version__ = "1.7.0-rc2"

skycatalogs/objects/base_object.py

@@ -4,6 +4,9 @@
 import itertools
 import numpy as np
 import galsim
+from galsim.roman import longwave_bands as roman_longwave_bands
+from galsim.roman import shortwave_bands as roman_shortwave_bands
+from galsim.roman import getBandpasses as roman_getBandpasses
 
 from skycatalogs.utils.translate_utils import form_object_string
 from skycatalogs.utils.config_utils import Config
@@ -15,9 +18,11 @@
 '''
 
 __all__ = ['BaseObject', 'ObjectCollection', 'ObjectList',
-           'LSST_BANDS', 'load_lsst_bandpasses']
+           'LSST_BANDS', 'ROMAN_BANDS',
+           'load_lsst_bandpasses', 'load_roman_bandpasses']
 
 LSST_BANDS = ('ugrizy')
+ROMAN_BANDS = roman_shortwave_bands+roman_longwave_bands
 
 # global for easy access for code run within mp
 
@@ -68,6 +73,16 @@ def load_lsst_bandpasses():
     return lsst_bandpasses
 
 
+def load_roman_bandpasses():
+    '''
+    Read in Roman bandpasses from standard place, trim, and store in global dict
+    Returns: The dict
+    '''
+    global roman_bandpasses
+    roman_bandpasses = roman_getBandpasses()
+    return roman_bandpasses
+
+
 class BaseObject(object):
     '''
     Abstract base class for static (in position coordinates) objects.
@@ -187,7 +202,6 @@ def _get_sed(self, component=None, resolution=None, mjd=None):
 
         Must be implemented by subclass
         '''
-
         raise NotImplementedError('Must be implemented by BaseObject subclass if needed')
 
     def write_sed(self, sed_file_path, component=None, resolution=None,
@@ -340,6 +354,45 @@ def get_LSST_fluxes(self, cache=True, as_dict=True, mjd=None):
         else:
             return list(fluxes.values())
 
+    def get_roman_flux(self, band, sed=None, cache=True, mjd=None):
+        if band not in ROMAN_BANDS:
+            return None
+        att = f'roman_flux_{band}'
+
+        # Check if it's already an attribute
+        val = getattr(self, att, None)
+        if val is not None:
+            return val
+
+        if att in self.native_columns:
+            return self.get_native_attribute(att)
+
+        val = self.get_flux(roman_bandpasses[band], sed=sed, mjd=mjd)
+
+        if cache:
+            setattr(self, att, val)
+        return val
+
+    def get_roman_fluxes(self, cache=True, as_dict=True, mjd=None):
+        '''
+        Return a dict of fluxes for Roman bandpasses or, if as_dict is False,
+        just a list in the same order as ROMAN_BANDS
+        '''
+        fluxes = dict()
+        sed = self.get_total_observer_sed(mjd=mjd)
+
+        if sed is None:
+            for band in ROMAN_BANDS:
+                fluxes[band] = 0.0
+        else:
+            for band in ROMAN_BANDS:
+                fluxes[band] = self.get_roman_flux(band, sed=sed,
+                                                   cache=cache, mjd=mjd)
+        if as_dict:
+            return fluxes
+        else:
+            return list(fluxes.values())
+
 
 class ObjectCollection(Sequence):
     '''
@@ -436,6 +489,8 @@ def subcomponents(self):
             for s in ['bulge', 'disk', 'knots']:
                 if f'sed_val_{s}' in self.native_columns:
                     subs.append(s)
+        elif self._object_type_unique == 'diffsky_galaxy':
+            subs = ['bulge', 'disk', 'knots']
         else:
             return ['this_object']
         return subs

skycatalogs/objects/diffsky_object.py

@@ -0,0 +1,153 @@
+# import numpy as np
+import galsim
+import numpy as np
+from .base_object import BaseObject
+
+__all__ = ['DiffskyObject']
+
+
+class DiffskyObject(BaseObject):
+    _type_name = 'diffsky_galaxy'
+
+    # sersic values are constant for diffsky galaxies
+    _sersic_disk = 1
+    _sersic_bulge = 4
+
+    def _get_sed(self, component=None, resolution=None):
+        '''
+        Return sed and mag_norm for a galaxy component or for a star
+        Parameters
+        ----------
+        component    one of 'bulge', 'disk', 'knots' for now. Other components
+                     may be supported.  Ignored for stars
+        resolution   desired resolution of lambda in nanometers. Ignored
+                     for stars.
+
+        Returns
+        -------
+        A pair (sed, mag_norm)
+        '''
+
+        if component not in ['disk', 'bulge', 'knots']:
+            raise ValueError(f'Cannot fetch SED for component type {component}')
+
+        if not hasattr(self, '_seds'):
+            z_h = self.get_native_attribute('redshiftHubble')
+            z = self.get_native_attribute('redshift')
+            pixel = self.partition_id
+
+            sky_cat = self._belongs_to._sky_catalog
+            self._seds = sky_cat.observed_sed_factory.create(pixel, self.id, z_h, z)
+
+        magnorm = None
+        return self._seds[component], magnorm
+
+    def get_knot_size(self,z):
+        """
+        Return the angular knot size. Knots are modelled as the same physical size
+        """
+        # Deceleration paramameter
+        q  = -0.5
+        # Angular diameter scaling approximation in pc
+        dA = (3e9/q**2)*(z*q+(q-1)*(np.sqrt(2*q*z+1)-1))/(1+z)**2*(1.4-0.53*z)
+        # Using typical knot size 250pc, convert to sigma in arcmin
+        if z<0.6:
+            return 206264.8*250/dA/2.355
+        else:
+            # Above z=0.6, fractional contribution to post-convolved size 
+            # is <20% for smallest Roman PSF size, so can treat as point source
+            return None
+
+    def get_knot_n(self,rng=None):
+        """
+        Return random value for number of knots based on galaxy sm.
+        """
+        if rng is not None:
+            ud = galsim.UniformDeviate(rng)
+        else:
+            ud = galsim.UniformDeviate(int(self.id))
+        sm = np.log10(self.get_native_attribute('um_source_galaxy_obs_sm'))
+        m  = (50-3)/(12-6) # (knot_n range)/(logsm range)
+        n_knot_max = m*(sm-6)+3
+        n_knot = int(ud()*n_knot_max) # random n up to n_knot_max
+        if n_knot==0:
+            n_knot+=1 # need at least 1 knot
+        return n_knot
+
+    def get_wl_params(self):
+        """Return the weak lensing parameters, g1, g2, mu."""
+        gamma1 = self.get_native_attribute('shear1')
+        gamma2 = self.get_native_attribute('shear2')
+        kappa = self.get_native_attribute('convergence')
+        # Compute reduced shears and magnification.
+        g1 = gamma1/(1. - kappa)    # real part of reduced shear
+        g2 = gamma2/(1. - kappa)    # imaginary part of reduced shear
+        mu = 1./((1. - kappa)**2 - (gamma1**2 + gamma2**2))  # magnification
+        return g1, g2, mu
+
+    def get_total_observer_sed(self, mjd=None):
+        """
+        Return the SED summed over SEDs for the individual SkyCatalog
+        components.
+        """
+        sed = super().get_total_observer_sed()
+
+        if sed is None:
+            return sed
+
+        _, _, mu = self.get_wl_params()
+        sed *= mu
+        return sed
+
+    def get_gsobject_components(self, gsparams=None, rng=None):
+
+        if gsparams is not None:
+            gsparams = galsim.GSParams(**gsparams)
+
+        if rng is None:
+            rng = galsim.BaseDeviate(int(self.id))
+
+        obj_dict = {}
+        for component in self.subcomponents:
+            # knots use the same major/minor axes as the disk component.
+            my_cmp = 'disk' if component != 'bulge' else 'spheroid'
+            hlr = self.get_native_attribute(f'{my_cmp}HalfLightRadiusArcsec')
+
+            # Get ellipticities saved in catalog. Not sure they're what
+            # we need
+            e1 = self.get_native_attribute(f'{my_cmp}Ellipticity1')
+            e2 = self.get_native_attribute(f'{my_cmp}Ellipticity2')
+            shear = galsim.Shear(g1=e1, g2=e2)
+
+            if component == 'knots':
+                npoints = self.get_knot_n()
+                assert npoints > 0
+                knot_profile = galsim.Sersic(n=self._sersic_disk, 
+                                            half_light_radius=hlr/2.,
+                                            gsparams=gsparams)
+                knot_profile = knot_profile._shear(shear)
+                obj = galsim.RandomKnots(npoints=npoints,
+                                         profile=knot_profile, rng=rng,
+                                         gsparams=gsparams)
+                z = self.get_native_attribute('redshift')
+                size = self.get_knot_size(z) # get knot size
+                if size is not None:
+                    obj = galsim.Convolve(obj, galsim.Gaussian(sigma=size))
+                obj_dict[component] = obj
+            else:
+                n = self._sersic_disk if component == 'disk' else self._sersic_bulge
+                obj = galsim.Sersic(n=n, half_light_radius=hlr,
+                                    gsparams=gsparams)
+                obj_dict[component] = obj._shear(shear)
+
+        # Apply lensing
+        g1, g2, mu = self.get_wl_params()
+        for component in self.subcomponents:
+            obj_dict[component] = obj_dict[component]._lens(g1, g2, mu)
+        return obj_dict
+
+    def get_observer_sed_component(self, component, mjd=None, resolution=None):
+        sed, _ = self._get_sed(component)
+        if sed is not None:
+            sed = self._apply_component_extinction(sed)
+        return sed