ENH add DES gal cat (#30)

* ENH add DES gal cat

* REF add cuts and start in on script for catalogs

* add use_all_gals option

* ENH add cli to make catalogs

* BUG need cli file at top-level

* BUG fix typo in import

* TST install ruff

* BUG make sure to use larger area

* BUG wrong field name

* ENH add code for v2 cosmos catalogs

* ENH v3 cats

* ENH v4 cats

* DOC add doc strings

* BUG this import is gone

* ENH v6 model

* Update tests.yml

* ENH v7

* TST added tests for new utils

---------

Co-authored-by: Sidney Mau <sidneymau@gmail.com>

Author: beckermr

.github/workflows/tests.yml

@@ -32,8 +32,8 @@ jobs:
       - name: configure conda and install code
         shell: bash -l {0}
         run: |
-          mamba install pytest flake8
-          python -m pip install -e .
+          conda install pytest ruff --yes --quiet
+          python -m pip install --no-deps --no-build-isolation -e .
 
       - name: test
         shell: bash -l {0}
@@ -43,4 +43,4 @@ jobs:
       - name: lint
         shell: bash -l {0}
         run: |
-          flake8 montara
+          ruff check montara

environment.yaml

@@ -2,12 +2,17 @@ channels:
   - conda-forge
 dependencies:
   # run
+  - click
   - galsim
   - fitsio
   - numpy
   - astropy
+  - hpgeom
+  - scipy
+  - esutil
   - pandas
   - des-eastlake
+  - des-y6utils
   - galsim_extra
   - importlib_metadata
   - cfitsio  # for fpack

montara/__init__.py

@@ -20,3 +20,5 @@
 from . import badpixfromfits  # noqa
 from . import eastlake_step  # noqa
 from . import eval_gsobject  # noqa
+from . import input_desgal  # noqa
+from . import make_input_desgal  # noqa

montara/des_tile.py

@@ -539,27 +539,41 @@ def setup(self, config, base, file_num, logger):
         if isinstance(nobjects, dict):
             # ^ this should return True for both dict and OrderedDict
             if base['image']['nobjects']['type'] == "MixedNObjects":
-                # First get the number of galaxies. Either this will be an int, in
-                # which case
-                # ParseValue will work straightaway, or a random variable, in which
-                # case we'll
-                # need to initalize the rng and then try ParseValue again.
-                try:
-                    ngalaxies = galsim.config.ParseValue(
-                        nobjects, 'ngalaxies', base, int)[0]
-                except TypeError:
-                    seed = galsim.config.ParseValue(
-                        base['image'], 'random_seed', base, int)[0]
+                if nobjects.get("use_all_gals", False):
+                    # if we use all of the galaxies, then we look at the input desgals
+                    # object to get the number of galaxies
+                    galsim.config.input.SetupInput(base, logger=logger)
+                    key = 'desgal'
+                    field = base['input'][key]
+                    loader = galsim.config.input.valid_input_types[key]
+                    gal_input = galsim.config.GetInputObj("desgal", config, base, "desgal")
+                    if gal_input is None:
+                        kwargs, safe = loader.getKwargs(field, base, logger)
+                        kwargs['_nobjects_only'] = True
+                        gal_input = loader.init_func(**kwargs)
+                    ngalaxies = gal_input.getNObjects()
+                else:
+                    # First get the number of galaxies. Either this will be an int, in
+                    # which case
+                    # ParseValue will work straightaway, or a random variable, in which
+                    # case we'll
+                    # need to initalize the rng and then try ParseValue again.
                     try:
-                        assert (isinstance(seed, int) and (seed != 0))
-                    except AssertionError as e:
-                        logger.critical(
-                            "image.random_seed must be set to a non-zero integer for "
-                            "output type DES_Tile")
-                        raise e
-                    base['rng'] = galsim.BaseDeviate(seed)
-                    ngalaxies = galsim.config.ParseValue(nobjects, 'ngalaxies',
-                                                         base, int)[0]
+                        ngalaxies = galsim.config.ParseValue(
+                            nobjects, 'ngalaxies', base, int)[0]
+                    except TypeError:
+                        seed = galsim.config.ParseValue(
+                            base['image'], 'random_seed', base, int)[0]
+                        try:
+                            assert (isinstance(seed, int) and (seed != 0))
+                        except AssertionError as e:
+                            logger.critical(
+                                "image.random_seed must be set to a non-zero integer for "
+                                "output type DES_Tile")
+                            raise e
+                        base['rng'] = galsim.BaseDeviate(seed)
+                        ngalaxies = galsim.config.ParseValue(nobjects, 'ngalaxies',
+                                                            base, int)[0]
                 logger.log(logging.CRITICAL, "simulating %d galaxies" % ngalaxies)
 
                 if nobjects.get("use_all_stars", True):

montara/input_desgal.py

@@ -0,0 +1,113 @@
+import logging
+
+import fitsio
+import galsim
+import numpy as np
+from galsim.config.input import InputLoader, RegisterInputType, GetInputObj
+from galsim.config.value import RegisterValueType
+
+from .utils import add_field
+
+logger = logging.getLogger("pipeline")
+
+
+class DESGalCatalog(object):
+    _req_params = {'file_name': str}
+    _opt_params = {'cuts': dict, 'verbose': bool}
+    _single_params = []
+    _takes_rng = False
+
+    def __init__(self, file_name, cuts=None, _nobjects_only=False, verbose=False):
+        self.gal_data = fitsio.read(file_name)
+        self.gal_data = add_field(
+            self.gal_data,
+            [("catalog_row", int)],
+            [np.arange(len(self.gal_data))],
+        )
+        self.cuts = cuts
+        if cuts is not None:
+            self.apply_cuts(cuts, verbose=verbose)
+        self.nobjects = len(self.gal_data)
+
+    def apply_cuts(self, cuts, verbose=False):
+        """Apply some cuts.
+
+         - `cuts` is a dictionary
+         - A key should be the same as a column in
+           the catalog.
+         - The value corresponding to that cut can either be
+           a single integer, in which case that value will be retained
+           e.g. "flags":0 would just keep objects with flags==0.
+           Or it can be a range e.g. "hlr":[a,b] in which case only
+           objects with a<=hlr<b will be retained.
+        """
+
+        use = np.ones(len(self.gal_data), dtype=bool)
+
+        for key, val in cuts.items():
+            col_data = self.gal_data[key]
+
+            if len(val) == 1:
+                mask = col_data == int(val[0])
+                cut_string = "{0} == {1}".format(key, val[0])
+            elif len(val) == 2:
+                mask = (col_data >= val[0]) * (col_data < val[1])
+                cut_string = "{0} <= {1} < {2}".format(val[0], key, val[1])
+            else:
+                raise ValueError("cut value should be length 1 or 2")
+
+            use[~mask] = False
+
+            if verbose:
+                print('applying {0} leaves fraction {1}'.format(
+                    cut_string, float(mask.sum())/len(mask))
+                )
+
+        if verbose:
+            print("%d/%d objects remaining after cuts" % (use.sum(), len(use)))
+
+        self.gal_data = self.gal_data[use]
+
+    def get(self, index, col):
+        return self.gal_data[index][col]
+
+    def getNObjects(self):
+        return self.nobjects
+
+
+def _GenerateFromDESGalCatalog(config, base, value_type):
+    """@brief Return a value read from an input catalog
+    """
+    gal_input = GetInputObj('desgal', config, base, 'DESGalValue')
+
+    # Setup the indexing sequence if it hasn't been specified.
+    # The normal thing with a Catalog is to just use each object in order,
+    # so we don't require the user to specify that by hand.  We can do it for them.
+    galsim.config.SetDefaultIndex(config, gal_input.getNObjects())
+
+    req = {'col': str, 'index': int}
+    opt = {'num': int}
+    kwargs, safe = galsim.config.GetAllParams(config, base, req=req, opt=opt)
+    col = kwargs['col']
+    index = kwargs['index']
+
+    logger.log(
+        logging.DEBUG,
+        "sampling fixed gal catalog band|index|col: %s %s %s" % (
+            base["eval_variables"]["sband"],
+            index,
+            col,
+        ),
+    )
+
+    val = gal_input.get(index, col)
+    return val, safe
+
+
+RegisterInputType('desgal', InputLoader(DESGalCatalog, has_nobj=True))
+RegisterValueType(
+    'DESGalValue',
+    _GenerateFromDESGalCatalog,
+    [float],
+    input_type='desgal',
+)