Add random magnitude generation

doc/lsst.source.injection/reference/12_generate_injection_catalog.rst

@@ -28,6 +28,7 @@ Catalog Parameters
 
 The tools on this page assist in generating a synthetic source injection catalog with sources quasi-randomly distributed across the sky.
 Sources will be scattered between the right ascension and declination limits specified by the user (using, for example, the `-a` and `-d` arguments respectively on the command line).
+Optional magnitudes may also be generated using the same sequence (using the `-m` command line argument).
 The number of repeats for each unique combination of profile parameters is also specified by the user (using, for example, the `-n` argument on the command line).
 
 .. tip::
@@ -165,6 +166,7 @@ Generate an Injection Catalog on the Command Line
 
 The :doc:`generate_injection_catalog <../scripts/generate_injection_catalog>` command line script is used to generate a synthetic source injection catalog.
 This script takes a number of arguments, including the right ascension and declination limits of the quasi-randomly generated positions and the number of sources to inject.
+Optional magnitudes may also be generated using the same random sequence.
 More information on the arguments accepted by this script may be found by running:
 
 .. code-block:: shell

python/lsst/source/injection/bin/generate_injection_catalog.py

@@ -86,6 +86,15 @@ def build_argparser():
         metavar="VALUE",
         nargs=2,
     )
+    parser_general.add_argument(
+        "-m",
+        "--mag-lim",
+        type=float,
+        help="The magnitude limits of the catalog in magnitudes.",
+        required=False,
+        metavar="VALUE",
+        nargs=2,
+    )
     parser_general.add_argument(
         "-n",
         "--number",
@@ -265,11 +274,13 @@ def main():
         logger.info("Using WCS in %s for %s.", wcs_type_name, dataset_ref.dataId)
 
     # Generate the source injection catalog.
+    mag_lim = vars(args).get("mag_lim", None)
     density = vars(args).get("density", None)
     seed = vars(args).get("seed", None)
     table = generate_injection_catalog(
         ra_lim=args.ra_lim,
         dec_lim=args.dec_lim,
+        mag_lim=mag_lim,
         wcs=wcs,
         number=args.number,
         density=density,

python/lsst/source/injection/utils/generate_injection_catalog.py

@@ -38,6 +38,7 @@
 def generate_injection_catalog(
     ra_lim: Sequence[float],
     dec_lim: Sequence[float],
+    mag_lim: Sequence[float] | None = None,
     wcs: SkyWcs = None,
     number: int = 1,
     density: int | None = None,
@@ -51,6 +52,7 @@ def generate_injection_catalog(
     supplied input parameters. The catalog is returned as an astropy Table.
 
     On-sky source positions are generated using the quasi-random Halton
+    sequence. Optional magnitudes may also be generated using the same
     sequence. By default, the Halton sequence is seeded using the product of
     the right ascension and declination limit ranges. This ensures that the
     same sequence is always generated for the same limits. This seed may be
@@ -74,6 +76,8 @@ def generate_injection_catalog(
         The right ascension limits of the catalog in degrees.
     dec_lim : `Sequence` [`float`]
         The declination limits of the catalog in degrees.
+    mag_lim : `Sequence` [`float`], optional
+        The magnitude limits of the catalog in magnitudes.
     wcs : `lsst.afw.geom.SkyWcs`, optional
         The WCS associated with these data. If not given or ``None`` (default),
         the catalog will be generated using Cartesian geometry.
@@ -171,6 +175,13 @@ def generate_injection_catalog(
     rng = np.random.default_rng(hashed_seed)
     sky_coords = Table(rng.permutation(sky_coords))
 
+    # Generate random magnitudes if limits are specified
+    if mag_lim:
+        mag_sampler = qmc.Halton(d=1, seed=hashed_seed)
+        mag_sample = mag_sampler.random(n=len(param_table))
+        mags = Table(qmc.scale(mag_sample, mag_lim[0], mag_lim[1]), names=("mag",))
+        sky_coords = hstack([sky_coords, mags])
+
     # Generate the unique injection ID and construct the final table.
     source_id = np.concatenate([([i] * number) for i in range(int(len(param_table) / number))])
     injection_id = param_table["version_id"] + source_id * int(10 ** np.ceil(np.log10(number)))