Injecting stream mocks into a survey#
This short tutorial shows how to:
Load a survey (e.g., LSST yr5) and inspect completeness/error models
Create a
StreamInjectorand inject photometric effects + detection flagsPlot against the footprint and visualize injected magnitudes
This is useful for people who want to convert their dynamical simulation results into realistic survey data.
import sys
import os
import pandas as pd
import yaml
import numpy as np
import matplotlib.pyplot as plt
%load_ext autoreload
%autoreload 2
from streamobs import surveys, observed
from streamobs.model import StreamModel
1) Load the survey (cached)#
We load LSST yr5 once; the data is cached for efficiency (bands, depth maps, extinction, completeness, and photometric error models).
lsst_yr4 = surveys.Survey.load(survey="lsst", release="yr4")
Loading survey data for 'lsst_yr4'...
Loading config from: lsst_yr4.yaml
======================================================================
LOADING SURVEY DATA FILES
======================================================================
Survey data directory: /Users/pelissier/Documents/Codes/packages/streamobs/streamobs/../data/surveys/lsst_yr4
Fallback directory for shared data files: /Users/pelissier/Documents/Codes/packages/streamobs/streamobs/../data/others
Available bands: g, i, r, u, y, z
Loading survey properties...
Loading magnitude limit maps...
✓ Success for g-band magnitude limit
⚠ Warning: 'maglim_map_i' not specified in config (skipping i-band)
✓ Success for r-band magnitude limit
⚠ Warning: 'maglim_map_u' not specified in config (skipping u-band)
⚠ Warning: 'maglim_map_y' not specified in config (skipping y-band)
⚠ Warning: 'maglim_map_z' not specified in config (skipping z-band)
Loading completeness/efficiency function...
Loading Completeness/efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading Detection efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading Classification efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading photometric error model...
Loading Photometric error model (catalog / reported)...
File: lsst_photoerror_r.csv
✓ Success
Loading band-independent maps...
Loading E(B-V) extinction map...
File: ebv_sfd98_lowres_nside_512_ring_equatorial.fits
✓ Success
⚠ Warning: 'coverage' not specified in config (skipping)
Building coverage map from magnitude limit maps...
✓ Built coverage map (nside=128, 133236 pixels covered)
Survey properties summary:
g-band:
Extinction coefficient: 3.661
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
i-band:
Extinction coefficient: 2.054
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
r-band:
Extinction coefficient: 2.701
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
u-band:
Extinction coefficient: 4.757
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
y-band:
Extinction coefficient: 1.308
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
z-band:
Extinction coefficient: 1.590
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
======================================================================
SURVEY DATA LOADED SUCCESSFULLY
======================================================================
✓ Survey 'lsst_yr4' loaded and cached successfully
Completeness and photometric error#
Both depend on magnitude relative to the local limit (delta_mag = mag - maglim). We can visualize them for the r-band at a chosen magnitude limit.
mag_r_test = np.linspace(15, 31, 1000)
maglim_r_dc2 = 26.8
completeness_val = lsst_yr4.get_completeness("r", mag_r_test, maglim_r_dc2)
mag_error_val = lsst_yr4.get_photo_error("r", mag_r_test, maglim_r_dc2)
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
ax[0].plot(mag_r_test, completeness_val, label="Completeness")
ax[0].axvline(maglim_r_dc2, color="red", linestyle="--", label="Magnitude Limit")
ax[0].set_xlabel("Apparent r-band magnitude")
ax[0].set_ylabel("Completeness")
ax[0].set_title("LSST DC2 Completeness Function")
ax[0].grid()
ax[0].legend()
ax[1].plot(mag_r_test, mag_error_val, label="Photometric Error", color="orange")
ax[1].axvline(maglim_r_dc2, color="red", linestyle="--", label="Magnitude Limit")
ax[1].set_yscale("log")
ax[1].set_xlabel("Apparent r-band magnitude")
ax[1].set_ylabel("Photometric Error (mag)")
ax[1].set_title("LSST DC2 Photometric Error Model")
ax[1].grid()
ax[1].legend()
plt.show()
Cache manipulation#
Once a survey has been loaded, it’s stored in the cache. You can use commands to check which surveys are stored. To force reloading of a survey, you must clear the cache first.
# Once a survey has been loaded once, it can be reused directly from the cache
lsst_yr4 = surveys.Survey.load(
survey="lsst", release="yr4"
) # reuses the cached survey object
# Check the cache content
print("Current cache content: ", surveys.SurveyFactory.list_cached_surveys())
surveys.SurveyFactory.clear_cache(survey="lsst", release="yr4")
# surveys.SurveyFactory.clear_cache(survey = 'lsst') # to clear all releases of 'lsst'
# surveys.SurveyFactory.clear_cache() # to clear the entire cache
print(
"Cache content after clearing 'lsst': ", surveys.SurveyFactory.list_cached_surveys()
)
✓ Using cached survey data for 'lsst_yr4'
Current cache content: ['lsst_yr4']
✓ Cleared cached survey 'lsst_yr4'
Cache content after clearing 'lsst': []
2) Create an injector#
The StreamInjector class wraps Survey objects and uses their properties to inject observational effects into a dataset.
# Create an injector using the Survey object
stream_injector = observed.StreamInjector(lsst_yr4)
# Or using survey name and release (this will create the survey object internally)
stream_injector = observed.StreamInjector(survey="lsst", release="yr4")
Loading survey data for 'lsst_yr4'...
Loading config from: lsst_yr4.yaml
======================================================================
LOADING SURVEY DATA FILES
======================================================================
Survey data directory: /Users/pelissier/Documents/Codes/packages/streamobs/streamobs/../data/surveys/lsst_yr4
Fallback directory for shared data files: /Users/pelissier/Documents/Codes/packages/streamobs/streamobs/../data/others
Available bands: g, i, r, u, y, z
Loading survey properties...
Loading magnitude limit maps...
✓ Success for g-band magnitude limit
⚠ Warning: 'maglim_map_i' not specified in config (skipping i-band)
✓ Success for r-band magnitude limit
⚠ Warning: 'maglim_map_u' not specified in config (skipping u-band)
⚠ Warning: 'maglim_map_y' not specified in config (skipping y-band)
⚠ Warning: 'maglim_map_z' not specified in config (skipping z-band)
Loading completeness/efficiency function...
Loading Completeness/efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading Detection efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading Classification efficiency function...
File: lsst_stellar_efficiency_cutr.csv
✓ Success
Loading photometric error model...
Loading Photometric error model (catalog / reported)...
File: lsst_photoerror_r.csv
✓ Success
Loading band-independent maps...
Loading E(B-V) extinction map...
File: ebv_sfd98_lowres_nside_512_ring_equatorial.fits
✓ Success
⚠ Warning: 'coverage' not specified in config (skipping)
Building coverage map from magnitude limit maps...
✓ Built coverage map (nside=128, 133236 pixels covered)
Survey properties summary:
g-band:
Extinction coefficient: 3.661
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
i-band:
Extinction coefficient: 2.054
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
r-band:
Extinction coefficient: 2.701
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
u-band:
Extinction coefficient: 4.757
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
y-band:
Extinction coefficient: 1.308
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
z-band:
Extinction coefficient: 1.590
Saturation limit: 16.0 mag
Systematic error: 0.0050 mag
======================================================================
SURVEY DATA LOADED SUCCESSFULLY
======================================================================
✓ Survey 'lsst_yr4' loaded and cached successfully
3) Inject different dataset types#
To convert a stream catalog into realistic survey data, you need:
Sky coordinates (ra, dec)
True apparent magnitudes of the stars
However, Streamobs can automatically sample these quantities if they’re missing, allowing incomplete datasets to be converted into realistic survey data.
Streamobs uses the magnitude limit at each star’s location to estimate:
The photometric error
The observed magnitude
Whether the object is detected and classified as a star
The output of the inject() method is a DataFrame containing:
ra,dec: positions of the stream’s stars on the skymag_{band}: true apparent magnitude in a given bandmag_{band}_obs: observed magnitude in a given bandmagerr_{band}: photometric error on the observed magnitude (in mag). Includes both systematic and statistical errors.flag_observed: 1 if the star has been detected by the survey and classified as a star. By default includes an SNR cut (SNR = 1/magerr > 5) in g and r bands.
The output can optionally include:
phi1,phi2: position of stars in the stream framedist: distance modulus of the stars
Dataset with only (phi1, phi2) coordinates#
If you only provide stream coordinates, Streamobs will:
1. Sample apparent magnitudes:
This is done using an isochrone and distance modulus model, which must be provided in a configuration dictionary. See tutorial_generate_datamocks.ipynb for more details.
2. Convert (phi1, phi2) → (ra, dec):
If you don’t provide (ra, dec) columns, Streamobs will randomly place the stream on the sky. You can restrict the placement using the mask_type argument to ensure the stream falls within:
The survey’s footprint
Regions with low dust extinction
Areas meeting a minimum magnitude limit in a given band
Any combination of the above
Streamobs searches uniformly across the sky for a reference frame where a specified fraction of the stream (percentile_threshold) lies inside the mask. For more details, see the documentation of the StreamInjector._find_gc_frame() method.
Alternatively, you can provide your own gala.coordinates.GreatCircleICRSFrame object via the gc_frame argument.
N = 1000
seed = 42
rng = np.random.default_rng(seed)
# Replace with your actual data
data_test = pd.DataFrame(
{
"phi1": rng.uniform(-5, 5, N),
"phi2": rng.uniform(-1, 1, N),
}
)
# Since our data set does not have magnitudes, we need to provide an isochrone and distance model
# - cf tutorial_generate_datamocks.ipynb
# This step can be skipped if the input data already has magnitudes in the desired bands
isochrone_config = {
"isochrone": {
"name": "Marigo2017", # isochrone set name
"survey": "lsst", # survey for filter set
"age": 12.0, # Age in Gyr of the population
"z": 0.0006, # Metallicity of the population
"band_1": "g", # first band for color-magnitude
"band_2": "r", # second band for color-magnitude
"band_1_detection": True,
},
}
distance_modulus_config = {
"distance_modulus": {
"center": {"type": "Constant", "value": 16.5},
"spread": {"type": "Constant", "value": 0.0},
}
}
stream_config = {**isochrone_config, **distance_modulus_config}
mask_type = [
"footprint",
"ebv",
] # Restrict to survey footprint and apply extinction mask
# By default the bands used for detection are (g and r)
# You can also add perfect_galstarsep=True to get bolean flag related to
# detection only (no classification).
injected_data_full = stream_injector.inject(
data_test,
seed=seed,
mask_type=mask_type,
stream_config=stream_config,
verbose=True,
perfect_galstarsep=True,
)
Filled 1000 dist values.
Filled ['lsst_g_true', 'lsst_r_true'] (missing rows only).
Building new mask for ['ebv', 'footprint']...
Resampling ebv from nside=512 to nside=128
✓ Mask created: valid pixels fraction = 0.5
Cached with key: ('lsst', ('ebv', 'footprint'), 0.2)
Found suitable great circle frame after 2 trials with 100.00% points inside the mask.
Applying dust correction for r-band on observed magnitudes.
Applying dust correction for g-band on observed magnitudes.
Applying detection cut on g-band with SNR >= 5.0
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
print(injected_data_full.columns)
injected_data_full.head()
Index(['phi1', 'phi2', 'dist', 'lsst_g_true', 'lsst_r_true', 'ra', 'dec',
'lsst_yr4_r_obs', 'lsst_yr4_r_err', 'lsst_yr4_g_obs', 'lsst_yr4_g_err',
'lsst_yr4_flag_observed', 'lsst_yr4_flag_perfect_galstarsep'],
dtype='object')
| phi1 | phi2 | dist | lsst_g_true | lsst_r_true | ra | dec | lsst_yr4_r_obs | lsst_yr4_r_err | lsst_yr4_g_obs | lsst_yr4_g_err | lsst_yr4_flag_observed | lsst_yr4_flag_perfect_galstarsep | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 2.739560 | -0.875874 | 16.5 | 30.303281 | 28.747045 | 309.295394 | -10.749590 | 27.944501 | 2.839885 | 27.617721 | 10.000001 | False | False |
| 1 | -0.611216 | -0.083476 | 16.5 | 29.728133 | 28.237061 | 310.268379 | -14.058956 | 27.607724 | 1.408137 | 28.568965 | 10.000001 | False | False |
| 2 | 3.585979 | -0.741940 | 16.5 | 28.225967 | 26.902620 | 308.757294 | -10.075699 | 26.892841 | 0.341235 | 28.845171 | 1.318412 | False | False |
| 3 | 1.973680 | -0.695347 | 16.5 | 29.385597 | 27.935658 | 309.516004 | -11.506027 | BAD_MAG | 1.061085 | 26.726983 | 10.000001 | False | False |
| 4 | -4.058227 | 0.264566 | 16.5 | 28.291092 | 26.960775 | 311.724304 | -17.227233 | 26.483797 | 0.409048 | 27.505773 | 1.932585 | False | False |
# If you want to work only with detected and classified stars
detected_data = injected_data_full[injected_data_full["lsst_yr4_flag_observed"]]
print(len(detected_data), "stars detected out of", len(injected_data_full))
# When applying the perfect galstar separation, all detected stars are
# classified as stars, so we have:
classified_data = injected_data_full[injected_data_full["lsst_yr4_flag_observed"]]
print(len(classified_data), "stars classified as stars out of", len(injected_data_full))
160 stars detected out of 1000
160 stars classified as stars out of 1000
Streamobs has automatically sampled the missing columns (positions and apparent magnitudes) and estimated observational quantities.
You can verify the injection behavior with some useful plots:
# Plot the injected stream in the survey mask
stream_injector.plot_stream_in_mask(injected_data_full, mask_type=mask_type)
(<Figure size 800x500 with 2 Axes>, <Axes: >)
The stream lies entirely within the chosen mask.
from streamobs.plotting import plot_inject
# Plot results
fig, ax = plot_inject(injected_data_full, lsst_yr4, bands=["g", "r"], save=False)
plt.show()
Dataset containing (ra, dec) coordinates and distance modulus#
Streamobs only needs to sample the missing magnitudes, so you must provide a configuration dictionary for the isochrone model.
# Let's build a data set containing only (ra, dec) coordinates
data_set_radec = injected_data_full.drop(
columns=[
col for col in injected_data_full.columns if col not in ["ra", "dec", "dist"]
]
)
data_set_radec.head()
| dist | ra | dec | |
|---|---|---|---|
| 0 | 16.5 | 309.295394 | -10.749590 |
| 1 | 16.5 | 310.268379 | -14.058956 |
| 2 | 16.5 | 308.757294 | -10.075699 |
| 3 | 16.5 | 309.516004 | -11.506027 |
| 4 | 16.5 | 311.724304 | -17.227233 |
injected_data = stream_injector.inject(
data_set_radec,
seed=seed,
mask_type=None,
stream_config=isochrone_config,
verbose=True,
)
injected_data.head()
Distance modulus model not defined; skipping distances.
Filled ['lsst_g_true', 'lsst_r_true'] (missing rows only).
Applying dust correction for r-band on observed magnitudes.
Applying dust correction for g-band on observed magnitudes.
Applying detection cut on g-band with SNR >= 5.0
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
| dist | ra | dec | lsst_g_true | lsst_r_true | lsst_yr4_r_obs | lsst_yr4_r_err | lsst_yr4_g_obs | lsst_yr4_g_err | lsst_yr4_flag_observed | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 16.5 | 309.295394 | -10.749590 | 22.307507 | 21.958448 | 21.955709 | 0.006555 | 22.315017 | 0.007065 | True |
| 1 | 16.5 | 310.268379 | -14.058956 | 31.246381 | 29.588484 | 27.543322 | 10.000001 | 30.050405 | 10.000001 | False |
| 2 | 16.5 | 308.757294 | -10.075699 | 28.129792 | 26.820850 | 26.811599 | 0.322715 | 28.015569 | 1.148001 | False |
| 3 | 16.5 | 309.516004 | -11.506027 | 29.071968 | 27.661313 | 29.617694 | 0.836162 | 28.525326 | 4.297743 | False |
| 4 | 16.5 | 311.724304 | -17.227233 | 30.192328 | 28.648052 | 26.836684 | 3.190956 | 26.780333 | 10.000001 | False |
Dataset containing (ra, dec) and true apparent magnitudes#
The dataset already contains all the required quantities, so you don’t need to provide a configuration dictionary or mask type.
# Let's build a dataset containing (ra, dec) coordinates and magnitudes
data_set_radecmag = injected_data_full.drop(
columns=[
col
for col in injected_data_full.columns
if col not in ["ra", "dec", "lsst_g_true", "lsst_r_true"]
]
)
data_set_radecmag.head()
| lsst_g_true | lsst_r_true | ra | dec | |
|---|---|---|---|---|
| 0 | 30.303281 | 28.747045 | 309.295394 | -10.749590 |
| 1 | 29.728133 | 28.237061 | 310.268379 | -14.058956 |
| 2 | 28.225967 | 26.902620 | 308.757294 | -10.075699 |
| 3 | 29.385597 | 27.935658 | 309.516004 | -11.506027 |
| 4 | 28.291092 | 26.960775 | 311.724304 | -17.227233 |
injected_data = stream_injector.inject(
data_set_radecmag, seed=seed, mask_type=None, stream_config=None, verbose=True
)
injected_data.head()
Applying dust correction for r-band on observed magnitudes.
Applying dust correction for g-band on observed magnitudes.
Applying detection cut on g-band with SNR >= 5.0
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
/opt/anaconda3/envs/streamsim_dev/lib/python3.11/site-packages/pandas/core/arraylike.py:399: RuntimeWarning: invalid value encountered in log10
result = getattr(ufunc, method)(*inputs, **kwargs)
| lsst_g_true | lsst_r_true | ra | dec | lsst_yr4_r_obs | lsst_yr4_r_err | lsst_yr4_g_obs | lsst_yr4_g_err | lsst_yr4_flag_observed | |
|---|---|---|---|---|---|---|---|---|---|
| 0 | 30.303281 | 28.747045 | 309.295394 | -10.749590 | 27.944501 | 2.839885 | BAD_MAG | 10.000001 | False |
| 1 | 29.728133 | 28.237061 | 310.268379 | -14.058956 | 27.607724 | 1.408137 | 28.532156 | 10.000001 | False |
| 2 | 28.225967 | 26.902620 | 308.757294 | -10.075699 | 26.892841 | 0.341235 | 28.095768 | 1.318412 | False |
| 3 | 29.385597 | 27.935658 | 309.516004 | -11.506027 | BAD_MAG | 1.061085 | 28.380889 | 10.000001 | False |
| 4 | 28.291092 | 26.960775 | 311.724304 | -17.227233 | 26.483797 | 0.409048 | 26.483845 | 1.932585 | False |
4) Usage example#
This can be useful for comparing analysis results between the ideal case of dynamic simulations and the more realistic case of observed data.
For example, we can calculate the 1D density along the stream and its power spectrum.
def plot_power_spectrum(data, fig=None, ax=None, phi1_bin_edges=None, label=None):
"""Plot 1D density along phi1 and its power spectrum."""
if fig is None or ax is None:
fig, ax = plt.subplots(1, 2, figsize=(12, 5))
# Compute 1D density along phi1
phi1 = np.array(data["phi1"])
if phi1_bin_edges is None:
phi1_bin_edges = np.linspace(np.min(phi1), np.max(phi1), 50)
counts, _ = np.histogram(phi1, bins=phi1_bin_edges)
bin_centers = 0.5 * (phi1_bin_edges[:-1] + phi1_bin_edges[1:])
# Compute power spectrum using cross-spectral density method from scipy
from scipy.signal import csd
fs = 1.0 / (bin_centers[1] - bin_centers[0]) # Sampling frequency
k, Pxx = csd(counts, counts, nperseg=len(counts), fs=fs)
# Plot density
ax[0].plot(bin_centers, counts, drawstyle="steps-mid", label=label)
ax[0].set_xlabel("phi1 (deg)")
ax[0].set_ylabel("Star counts")
ax[0].set_title("1D Density along phi1")
# Plot power spectrum (1/k, P(k)) in log-log scale
ax[1].loglog(1 / k, Pxx, label=label)
ax[1].set_xlabel("1/k (deg)")
ax[1].set_ylabel("P(k)")
ax[1].set_title("Power Spectrum of 1D Density")
for a in ax.flatten():
a.grid()
if label is not None:
a.legend()
return fig, ax
# Compare the power spectrum of all stars vs. detected stars only
phi1_bin_edges = np.linspace(-7, 7, 50)
fig, ax = plot_power_spectrum(
injected_data_full, label="All data", phi1_bin_edges=phi1_bin_edges
)
fig, ax = plot_power_spectrum(
injected_data_full[injected_data_full["lsst_yr4_flag_observed"]],
fig=fig,
ax=ax,
label="Detected data",
phi1_bin_edges=phi1_bin_edges,
)
fig.tight_layout()
/var/folders/9j/zk4h50g505d635j69p9tgwbc0000gp/T/ipykernel_10726/3483219967.py:25: RuntimeWarning: divide by zero encountered in divide
ax[1].loglog(1/k, Pxx, label=label)
5) Conclusion#
Streamobs provides a method for converting data mocks into realistic observed data for a given survey, with minimal computational cost. This can be easily achieved with just a few commands, helping to bridge the gap between simulations and observational challenges in the field of stellar streams.
You can find more informations in the full documentation.