Generate Spatially Varying Magnitude Errors According to Observing Conditions

last run successfully: Feb 9, 2026

The ObsCondition degrader can be used to generate spatially-varying photometric errors using input survey condition maps in healpix format, such as survey coadd depth, airmass, sky brightness etc. The photometric error is computed by photerr.LsstErrorModel, based on the LSST Overview Paper: https://arxiv.org/abs/0805.2366.

The degrader assigns each object in the input catalogue with a pixel within the survey footprint and computes the magnitude error (SNR) on each pixel. The degrader takes the following arguments:

nside: nside used for the HEALPIX maps.
mask: Path to the mask covering the survey footprint in HEALPIX format. Notice that all negative values will be set to zero.
weight: Path to the weights HEALPIX format, used to assign sample galaxies in pixels. Default is weight=““, which uniform weighting.
tot_nVis_flag: If nVisYr is provided in map_dict (see below), this flag indicates whether the map shows the total number of visits in nYrObs (tot_nVis_flag=True), or the average number of visits per year (tot_nVis_flag=False). The default is set to True.
random_seed: A random seed for reproducibility.
map_dict: A dictionary that contains the paths to the survey condition maps in HEALPIX format. This dictionary uses the same arguments as LSSTErrorModel. The following arguements, if supplied, may contain either a single number (as in the case of LSSTErrorModel), or a path to the corresponding survey condition map in healpix format:m5, nVisYr, airmass, gamma, msky, theta, km, and tvis. Notice that except airmass and tvis, for all other arguements, numbers/paths for specific bands should be passed. Other LsstErrorModel parameters can also be passed in this dictionary (e.g. a necessary one may be nYrObs for the survey condition maps; the default value is 10 years, although most may be interested in early data releases). If any arguement is not passed, the default value in https://arxiv.org/abs/0805.2366 is adopted. Example:

{
   "m5": {"u": "path", ...},
   "theta": {"u": "path", ...},
}

Argument defaults are determined by the defaults of the LsstErrorModel in PhotErr.

Note: If you’re interested in running this in pipeline mode, see 09_Spatial_Variability.ipynb in the pipeline_examples/creation_examples/ folder.

In this quick notebook we’ll generate the photometric error based on the DC2 Y5 LSST median \(5\sigma\) depth in \(i\)-band generated by OpSim minion_1016 database using the Rubin Observatory Metrics Analysis Framework (MAF).

import healpy as hp
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import rail.interactive as ri
from rail.utils.path_utils import find_rail_file

Install FSPS with the following commands:
pip uninstall fsps
git clone --recursive https://github.com/dfm/python-fsps.git
cd python-fsps
python -m pip install .
export SPS_HOME=$(pwd)/src/fsps/libfsps

LEPHAREDIR is being set to the default cache directory:
/home/runner/.cache/lephare/data
More than 1Gb may be written there.
LEPHAREWORK is being set to the default cache directory:
/home/runner/.cache/lephare/work
Default work cache is already linked.
This is linked to the run directory:
/home/runner/.cache/lephare/runs/20260601T134643

A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.

If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.

Traceback (most recent call last):  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/runpy.py", line 196, in _run_module_as_main
    return _run_code(code, main_globals, None,
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/runpy.py", line 86, in _run_code
    exec(code, run_globals)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel_launcher.py", line 18, in <module>
    app.launch_new_instance()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/traitlets/config/application.py", line 1082, in launch_instance
    app.start()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/kernelapp.py", line 758, in start
    self.io_loop.start()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/tornado/platform/asyncio.py", line 211, in start
    self.asyncio_loop.run_forever()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/asyncio/base_events.py", line 603, in run_forever
    self._run_once()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/asyncio/base_events.py", line 1909, in _run_once
    handle._run()
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/asyncio/events.py", line 80, in _run
    self._context.run(self._callback, *self._args)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/utils.py", line 71, in preserve_context
    return await f(*args, **kwargs)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/kernelbase.py", line 621, in shell_main
    await self.dispatch_shell(msg, subshell_id=subshell_id)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/kernelbase.py", line 478, in dispatch_shell
    await result
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/ipkernel.py", line 372, in execute_request
    await super().execute_request(stream, ident, parent)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/kernelbase.py", line 834, in execute_request
    reply_content = await reply_content
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/ipkernel.py", line 464, in do_execute
    res = shell.run_cell(
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/ipykernel/zmqshell.py", line 663, in run_cell
    return super().run_cell(*args, **kwargs)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3077, in run_cell
    result = self._run_cell(
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3132, in _run_cell
    result = runner(coro)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/async_helpers.py", line 128, in _pseudo_sync_runner
    coro.send(None)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3336, in run_cell_async
    has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3519, in run_ast_nodes
    if await self.run_code(code, result, async_=asy):
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/IPython/core/interactiveshell.py", line 3579, in run_code
    exec(code_obj, self.user_global_ns, self.user_ns)
  File "/tmp/ipykernel_5342/2407287589.py", line 5, in <module>
    import rail.interactive as ri
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/interactive/__init__.py", line 3, in <module>
    from . import calib, creation, estimation, evaluation, tools
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/interactive/calib/__init__.py", line 3, in <module>
    from rail.utils.interactive.initialize_utils import _initialize_interactive_module
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/utils/interactive/initialize_utils.py", line 17, in <module>
    from rail.utils.interactive.base_utils import (
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/utils/interactive/base_utils.py", line 10, in <module>
    rail.stages.import_and_attach_all(silent=True)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/stages/__init__.py", line 74, in import_and_attach_all
    RailEnv.import_all_packages(silent=silent)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/core/introspection.py", line 541, in import_all_packages
    _imported_module = importlib.import_module(pkg)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/importlib/__init__.py", line 126, in import_module
    return _bootstrap._gcd_import(name[level:], package, level)
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/som/__init__.py", line 1, in <module>
    from rail.creation.degraders.specz_som import *
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/specz_som.py", line 15, in <module>
    from somoclu import Somoclu
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/somoclu/__init__.py", line 11, in <module>
    from .train import Somoclu
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/somoclu/train.py", line 25, in <module>
    from .somoclu_wrap import train as wrap_train
  File "/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/somoclu/somoclu_wrap.py", line 11, in <module>
    import _somoclu_wrap

---------------------------------------------------------------------------

ImportError                               Traceback (most recent call last)

File /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/numpy/core/_multiarray_umath.py:44, in __getattr__(attr_name)
     39     # Also print the message (with traceback).  This is because old versions
     40     # of NumPy unfortunately set up the import to replace (and hide) the
     41     # error.  The traceback shouldn't be needed, but e.g. pytest plugins
     42     # seem to swallow it and we should be failing anyway...
     43     sys.stderr.write(msg + tb_msg)
---> 44     raise ImportError(msg)
     46 ret = getattr(_multiarray_umath, attr_name, None)
     47 if ret is None:


ImportError:
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.

If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.

Warning: the binary library cannot be imported. You cannot train maps, but you can load and analyze ones that you have already saved.
The problem occurs because either compilation failed when you installed Somoclu or a path is missing from the dependencies when you are trying to import it. Please refer to the documentation to see your options.

Let’s generate some fake data.

# Fake data with same magnitude in each band
i = np.ones(50_000) * 23.0
u = np.full_like(i, 23.0, dtype=np.double)
g = np.full_like(i, 23.0, dtype=np.double)
r = np.full_like(i, 23.0, dtype=np.double)
y = np.full_like(i, 23.0, dtype=np.double)
z = np.full_like(i, 23.0, dtype=np.double)
redshift = np.random.uniform(size=len(i)) * 2

mockdict = {}
for label, item in zip(
    ["redshift", "u", "g", "r", "i", "z", "y"], [redshift, u, g, r, i, z, y]
):
    mockdict[f"{label}"] = item

data = pd.DataFrame(mockdict)
data.head()

	redshift	u	g	r	i	z	y
0	0.828264	23.0	23.0	23.0	23.0	23.0	23.0
1	0.975312	23.0	23.0	23.0	23.0	23.0	23.0
2	0.746123	23.0	23.0	23.0	23.0	23.0	23.0
3	0.058103	23.0	23.0	23.0	23.0	23.0	23.0
4	1.045653	23.0	23.0	23.0	23.0	23.0	23.0

data_degraded = ri.creation.degraders.observing_condition_degrader.obs_condition(
    sample=data,
    map_dict={
        "m5": {
            "i": find_rail_file(
                "examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_fiveSigmaDepth_i_and_nightlt1825_HEAL.fits"
            ),
        },
        "nYrObs": 5.0,
    },
)

Inserting handle into data store.  input: None, ObsCondition
Assigning pixels.
No ra, dec found in catalogue, randomly assign pixels with weights.
Warning: objects found outside given mask, pixel assigned=-99. These objects will be assigned with defualt error from LSST error model!

Inserting handle into data store.  output: inprogress_output.pq, ObsCondition

data_degraded["output"].head()

	redshift	u	u_err	g	g_err	r	r_err	i	i_err	z	z_err	y	y_err	ra	decl	pixel
0	0.828264	22.990404	0.026246	23.011212	0.010113	22.991270	0.008975	23.001056	0.014845	22.946670	0.023212	22.952105	0.052272	61.171875	-40.620185	162135
1	0.975312	22.962996	0.025645	23.012586	0.010122	23.024129	0.009159	22.990276	0.013320	22.975839	0.023801	23.027019	0.055860	63.632812	-34.953865	154458
2	0.746123	23.050979	0.027628	22.992791	0.009991	23.012578	0.009093	22.985111	0.014289	23.019294	0.024710	23.020908	0.055559	52.795276	-42.210370	164170
3	0.058103	23.015021	0.026798	22.984763	0.009939	22.995173	0.008997	22.985953	0.014481	23.011103	0.024536	23.026922	0.055856	53.789062	-39.450895	160588
4	1.045653	23.005344	0.026579	23.010790	0.010110	22.991376	0.008976	23.000950	0.014682	23.005917	0.024426	22.946816	0.052028	69.609375	-28.971532	145763

We see that extra columns containing the magnitude errors: u_err, g_err… have been added to the catalogue. Notice that since we have only provided the limiting magnitude for \(i\)-band, the errors in all other bands except \(i\) are computed using the default parameters in LsstErrorModel (see: https://github.com/jfcrenshaw/photerr/blob/main/photerr/lsst.py).

The last column shows the pixel of the survey condition map that is assigned to each object.

We can check if the spatial dependence has been implemented by looking at the SNR at different area of the sky, and compare that with the \(i\)-band depth:

mask = hp.read_map(
    find_rail_file(
        "examples_data/creation_data/data/survey_conditions/DC2-mask-neg-nside-128.fits"
    )
)
weight = hp.read_map(
    find_rail_file(
        "examples_data/creation_data/data/survey_conditions/DC2-dr6-galcounts-i20-i25.3-nside-128.fits"
    )
)
Med_5sd_i = hp.read_map(
    find_rail_file(
        "examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_fiveSigmaDepth_i_and_nightlt1825_HEAL.fits"
    )
)

# Set negative values in mask to zero
mask[mask < 0] = 0

# Compute the average SNR in each pixel
avg_SNR_i = np.zeros(len(mask))
for pix, pix_cat in (data_degraded["output"]).groupby("pixel"):
    avg_SNR_i[pix] = np.mean((pix_cat["i"] / pix_cat["i_err"]).to_numpy())

# View the healpix map

fig, axarr = plt.subplots(1, 3, figsize=[12, 6])

plt.sca(axarr[0])
hp.gnomview(
    weight * mask / sum(weight),
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="weight",
    hold=True,
)
plt.sca(axarr[1])
hp.gnomview(
    Med_5sd_i * mask,
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="5sigmadepth i",
    hold=True,
)
plt.sca(axarr[2])
hp.gnomview(
    avg_SNR_i,
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="avg SNR i",
    min=1400,
    max=1750,
    hold=True,
)

../../../_images/Spatial_Variability_15_0.png

Now if we want to change any of the default settings, we can supply them in ObsCondition.make_stage(). In this example, instead of supplying the median \(5\sigma\) depth in \(i\)-band, we supply the median airmass in \(i\)-band. In this case, the \(i\)-band limiting magnitude m5 will be computed explicitly (notice that if m5 is also supplied, then it will overwrite the explicitly computed m5).

data_degraded_airmass = ri.creation.degraders.observing_condition_degrader.obs_condition(
    sample=data,
    map_dict={
        "airmass": find_rail_file(
            "examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_airmass_i_and_nightlt1825_HEAL.fits"
        ),
        "m5": {
            "i": find_rail_file(
                "examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_fiveSigmaDepth_i_and_nightlt1825_HEAL.fits"
            ),
        },
        "nYrObs": 5.0,
    },
)

Inserting handle into data store.  input: None, ObsCondition
Assigning pixels.
No ra, dec found in catalogue, randomly assign pixels with weights.
Warning: objects found outside given mask, pixel assigned=-99. These objects will be assigned with defualt error from LSST error model!

Inserting handle into data store.  output: inprogress_output.pq, ObsCondition

data_degraded_airmass["output"].head()

	redshift	u	u_err	g	g_err	r	r_err	i	i_err	z	z_err	y	y_err	ra	decl	pixel
0	0.828264	22.990381	0.026308	23.011222	0.010121	22.991259	0.008986	23.001056	0.014845	22.946682	0.023207	22.952134	0.052242	61.171875	-40.620185	162135
1	0.975312	22.963601	0.025237	23.012513	0.010062	23.024069	0.009136	22.990276	0.013320	22.975912	0.023731	23.026800	0.055413	63.632812	-34.953865	154458
2	0.746123	23.050276	0.027242	22.992825	0.009945	23.012555	0.009077	22.985111	0.014289	23.019244	0.024645	23.020765	0.055185	52.795276	-42.210370	164170
3	0.058103	23.014777	0.026365	22.984848	0.009884	22.995184	0.008976	22.985953	0.014481	23.011070	0.024462	23.026709	0.055420	53.789062	-39.450895	160588
4	1.045653	23.005223	0.025983	23.010704	0.010029	22.991407	0.008944	23.000950	0.014682	23.005894	0.024330	22.947353	0.051527	69.609375	-28.971532	145763

Again, we can examine whether the spatial dependence is indeed applied. Here, LsstErrorModel does not have band-dependent airmass, so it affects all bands. The default airmass is \(X=1.2\), but the input median airmass is more optimistic, thus reducing the magnitude errors.

Med_airmass_i = hp.read_map(
    find_rail_file(
        "examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_airmass_i_and_nightlt1825_HEAL.fits"
    )
)

Compute the average SNR in each pixel for i and r bands:

avg_SNR_i_airmass = np.zeros(len(mask))
avg_SNR_r_airmass = np.zeros(len(mask))
for pix, pix_cat in (data_degraded_airmass["output"]).groupby("pixel"):
    avg_SNR_i_airmass[pix] = np.mean((pix_cat["i"] / pix_cat["i_err"]).to_numpy())
    avg_SNR_r_airmass[pix] = np.mean((pix_cat["r"] / pix_cat["r_err"]).to_numpy())

View the healpix map:

fig, axarr = plt.subplots(1, 3, figsize=[12, 6])

plt.sca(axarr[0])
hp.gnomview(
    Med_airmass_i * mask,
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="airmass i",
    hold=True,
)
plt.sca(axarr[1])
hp.gnomview(
    avg_SNR_i_airmass,
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="avg SNR i",
    min=2240,
    max=2280,
    hold=True,
)

plt.sca(axarr[2])
hp.gnomview(
    avg_SNR_r_airmass,
    rot=(62, -36.5, 0),
    xsize=100,
    ysize=100,
    reso=16,
    title="avg SNR r",
    min=2930,
    max=2970,
    hold=True,
)

../../../_images/Spatial_Variability_24_0.png

In both cases, we see a negative correlation between the airmass and the SNR in \(i\) and \(r\) bands, as expected.