Generate Spatially Varying Magnitude Errors According to Observing Conditions
last run successfully: April 26, 2023
Note: If you’re planning to run this in a notebook, you may want to
use interactive mode instead. See
Spatial_Variability.ipynb
in the interactive_examples/creation_examples/ folder for a version
of this notebook in interactive mode.
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: IfnVisYris provided inmap_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 toTrue.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 inhealpixformat:m5,nVisYr,airmass,gamma,msky,theta,km, andtvis. Notice that exceptairmassandtvis, for all other arguements, numbers/paths for specific bands should be passed. OtherLsstErrorModelparameters can also be passed in this dictionary (e.g. a necessary one may benYrObsfor 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.
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
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
import os
import pandas as pd
import tables_io
import rail
from rail.core.stage import RailStage
from rail.utils.path_utils import find_rail_file
Let’s generate some fake data.
# Fake data with same magnitude in each band
i = np.ones(50_000)*23.
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 | 1.642107 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 |
| 1 | 1.350961 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 |
| 2 | 0.504739 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 |
| 3 | 0.061243 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 |
| 4 | 1.072184 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 | 23.0 |
Now let’s import the ObsCondition from rail.
from rail.creation.degraders import observing_condition_degrader
from rail.creation.degraders.observing_condition_degrader import ObsCondition
# First, let's use default arguments:
obs_cond_degrader = ObsCondition.make_stage()
# You can see what arguments have been entered by printing the degrader:
print(obs_cond_degrader)
Loaded observing conditions from configuration file:
nside = 128,
mask file: /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/../../examples_data/creation_data/data/survey_conditions/DC2-mask-neg-nside-128.fits,
weight file: /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/../../examples_data/creation_data/data/survey_conditions/DC2-dr6-galcounts-i20-i25.3-nside-128.fits,
tot_nVis_flag = True,
random_seed = 42,
map_dict contains the following items:
{'m5': {'i': '/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/../../examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_fiveSigmaDepth_i_and_nightlt1825_HEAL.fits'}, 'nYrObs': 5.0}
Let’s run the code and see how long it takes:
%%time
data_degraded = obs_cond_degrader(data)
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
CPU times: user 3.38 s, sys: 36.7 ms, total: 3.42 s
Wall time: 3.42 s
data_degraded.data.head()
| redshift | u | u_err | g | g_err | r | r_err | i | i_err | z | z_err | y | y_err | ra | decl | pixel | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.642107 | 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 | 1.350961 | 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.504739 | 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.061243 | 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.072184 | 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.data).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)
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).
airmass_degrader = ObsCondition.make_stage(
map_dict={"airmass": find_rail_file("examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_airmass_i_and_nightlt1825_HEAL.fits"),
"nYrObs": 5.0}
)
print(airmass_degrader)
Loaded observing conditions from configuration file:
nside = 128,
mask file: /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/../../examples_data/creation_data/data/survey_conditions/DC2-mask-neg-nside-128.fits,
weight file: /opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/creation/degraders/../../examples_data/creation_data/data/survey_conditions/DC2-dr6-galcounts-i20-i25.3-nside-128.fits,
tot_nVis_flag = True,
random_seed = 42,
map_dict contains the following items:
{'airmass': '/opt/hostedtoolcache/Python/3.10.20/x64/lib/python3.10/site-packages/rail/examples_data/creation_data/data/survey_conditions/minion_1016_dc2_Median_airmass_i_and_nightlt1825_HEAL.fits', 'nYrObs': 5.0}
data_degraded_airmass = airmass_degrader(data)
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.data.head()
| redshift | u | u_err | g | g_err | r | r_err | i | i_err | z | z_err | y | y_err | ra | decl | pixel | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1.642107 | 22.990381 | 0.026308 | 23.011222 | 0.010121 | 22.991259 | 0.008986 | 23.000944 | 0.013288 | 22.946682 | 0.023207 | 22.952134 | 0.052242 | 61.171875 | -40.620185 | 162135 |
| 1 | 1.350961 | 22.963601 | 0.025237 | 23.012513 | 0.010062 | 23.024069 | 0.009136 | 22.990413 | 0.013135 | 22.975912 | 0.023731 | 23.026800 | 0.055413 | 63.632812 | -34.953865 | 154458 |
| 2 | 0.504739 | 23.050276 | 0.027242 | 22.992825 | 0.009945 | 23.012555 | 0.009077 | 22.986365 | 0.013102 | 23.019244 | 0.024645 | 23.020765 | 0.055185 | 52.795276 | -42.210370 | 164170 |
| 3 | 0.061243 | 23.014777 | 0.026365 | 22.984848 | 0.009884 | 22.995184 | 0.008976 | 22.987306 | 0.013105 | 23.011070 | 0.024462 | 23.026709 | 0.055420 | 53.789062 | -39.450895 | 160588 |
| 4 | 1.072184 | 23.005223 | 0.025983 | 23.010704 | 0.010029 | 22.991407 | 0.008944 | 23.000856 | 0.013226 | 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.data).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)
avg_SNR_i_airmass
array([0., 0., 0., ..., 0., 0., 0.], shape=(196608,))
In both cases, we see a negative correlation between the airmass and the SNR in \(i\) and \(r\) bands, as expected.