created a skeleton function opsim_time_series_images_data

.github/workflows/tests.yml

@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.8]
+        python-version: [3.9]
 
     steps:
       - uses: actions/checkout@v4

.gitignore

@@ -35,6 +35,9 @@ share/python-wheels/
 .installed.cfg
 *.egg
 MANIFEST
+data/OpSim_database/
+data/Scotch/
+data/SALT3.NIR_WAVEEXT/
 
 # PyInstaller
 #  Usually these files are written by a python script from a template

slsim/LsstSciencePipeline/opsim_pipeline.py

@@ -0,0 +1,219 @@
+import numpy as np
+from astropy.table import Table
+from astropy.table import Column
+from slsim.image_simulation import lens_image
+import lenstronomy.Util.util as util
+import lenstronomy.Util.kernel_util as kernel_util
+import lenstronomy.Util.data_util as data_util
+
+
+def opsim_time_series_images_data(
+    ra_list,
+    dec_list,
+    obs_strategy,
+    MJD_min=60000,
+    MJD_max=64000,
+    size=101,
+    moffat_beta=3.1,
+    readout_noise=10,
+    delta_pix=0.2,
+    print_warning=True,
+    opsim_path=None,
+):
+    """Creates time series data from opsim database.
+
+    :param ra_list: a list of ra points (in degrees) from objects we want to collect
+        observations for
+    :param dec_list: a list of dec points (in degrees) from objects we want to collect
+        observations for
+    :param obs_strategy: version of observing strategy corresponding to opsim database.
+        for example "baseline_v3.0_10yrs" (string)
+    :param MJD_min: minimum MJD for the observations
+    :param MJD_max: maximum MJD for the observations
+    :param size: cutout size of images (in pixels)
+    :param moffat_beta: power index of the moffat psf kernel
+    :param readout_noise: noise added per readout
+    :param delta_pix: size of pixel in units arcseonds
+    :param print_warning: if True, prints a warning of coordinates outside of the LSST
+        footprint
+    :param opsim_path: optional: provide a path to the opsim database.
+        if None: use "../data/OpSim_database/" + obs_strategy + ".db" as default path.
+    :return: a list of astropy tables containing observation information for each
+        coordinate
+    """
+
+    # Import OpSimSummaryV2
+    try:
+        import opsimsummaryv2 as op
+    except ImportError:
+        raise ImportError(
+            "Users need to have OpSimSummaryV2 installed (https://github.com/bastiencarreres/OpSimSummaryV2)"
+        )
+
+    # Initialise OpSimSummaryV2 with opsim database
+    if opsim_path is None:
+        opsim_path = "../data/OpSim_database/" + obs_strategy + ".db"
+    try:
+        OpSimSurv = op.OpSimSurvey(opsim_path)
+    except FileNotFoundError:
+        raise FileNotFoundError(
+            "File not found: "
+            + opsim_path
+            + ". Input variable 'obs_strategy' should correspond to the name of an opsim database saved in the folder ../data/OpSim_database"
+        )
+
+    # Collect observations that cover the coordinates in ra_list and dec_list
+    gen = OpSimSurv.get_obs_from_coords(
+        ra_list,
+        dec_list,
+        is_deg=True,
+        formatobs=True,
+        keep_keys=["visitExposureTime", "seeingFwhmGeom", "fieldRA", "fieldDec"],
+    )
+
+    table_data_list = []
+
+    # Loop through all coordinates and compute the table_data
+    for i in range(len(ra_list)):
+
+        # Collect the next observation sequence from the opsim generator
+        seq = next(gen)
+        seq = seq.sort_values(by=["expMJD"])
+
+        # Check if the coordinates are in the opsim LSST footprint
+        opsim_ra = np.mean(seq["fieldRA"])
+        opsim_dec = np.mean(seq["fieldDec"])
+
+        if np.isnan(opsim_ra) or np.isnan(opsim_dec):
+            if print_warning:
+                print(
+                    f"Coordinate ({ra_list[i]}, {dec_list[i]}) is not in the LSST footprint. This entry is skipped."
+                )
+            continue
+
+        # Get the relevant properties from opsim
+        obs_time = np.array(seq["expMJD"])
+
+        # Only give the observations between MJD_min and MJD_max
+        mask = (obs_time > MJD_min) & (obs_time < MJD_max)
+        obs_time = obs_time[mask]
+
+        expo_time = np.array(seq["visitExposureTime"])[mask]
+        bandpass = np.array(seq["BAND"])[mask]
+        zero_point_mag = np.array(seq["ZPT"])[mask]
+        sky_brightness = np.array(seq["SKYSIG"])[mask] ** 2 / (delta_pix**2 * expo_time)
+        psf_fwhm = np.array(seq["seeingFwhmGeom"])[mask]
+        # Question: use 'FWHMeff' or 'seeingFwhmGeom' for the psf?
+
+        radec_list = [(ra_list[i], dec_list[i])] * len(obs_time)
+
+        # Create a Moffat psf kernel for each epoch
+
+        psf_kernels = []
+
+        for psf in psf_fwhm:
+            psf_kernel = kernel_util.kernel_moffat(
+                num_pix=size, delta_pix=delta_pix, fwhm=psf, moffat_beta=moffat_beta
+            )
+            psf_kernel = util.array2image(psf_kernel)
+
+            psf_kernels.append(psf_kernel)
+
+        psf_kernels = np.array(psf_kernels)
+
+        # Calculate background noise
+        bkg_noise = data_util.bkg_noise(
+            readout_noise, expo_time, sky_brightness, delta_pix, num_exposures=1
+        )
+
+        table_data = Table(
+            [
+                bkg_noise,
+                psf_kernels,
+                obs_time,
+                expo_time,
+                zero_point_mag,
+                radec_list,
+                bandpass,
+            ],
+            names=(
+                "bkg_noise",
+                "psf_kernel",
+                "obs_time",
+                "expo_time",
+                "zero_point",
+                "calexp_center",
+                "band",
+            ),
+        )
+
+        table_data_list.append(table_data)
+    return table_data_list
+
+
+def opsim_variable_lens_injection(
+    lens_class, bands, num_pix, transform_pix2angle, exposure_data
+):
+    """Injects variable lens to the OpSim time series data (1 object).
+
+    :param lens_class: Lens() object
+    :param bands: list of imaging bands of interest
+    :param num_pix: number of pixels per axis
+    :param transform_pix2angle: transformation matrix (2x2) of pixels into coordinate
+        displacements
+    :param exposure_data: An astropy table of exposure data. One entry of
+        table_list_data generated from the opsim_time_series_images_data function. It
+        must contain the rms of background noise fluctuations (column name should be
+        "bkg_noise"), psf kernel for each exposure (column name should be "psf_kernel",
+        these are pixel psf kernel for each single exposure images in time series
+        image), observation time (column name should be "obs_time", these are
+        observation time in days for each single exposure images in time series images),
+        exposure time (column name should be "expo_time", these are exposure time for
+        each single exposure images in time series images), magnitude zero point (column
+        name should be "zero_point", these are zero point magnitudes for each single
+        exposure images in time series image), coordinates of the object (column name
+        should be "calexp_center"), these are the coordinates in (ra, dec), and the band
+        in which the observation is taken (column name should be "band").
+    :return: Astropy table of injected lenses and exposure information of dp0 data
+    """
+
+    final_image = []
+
+    for obs in range(len(exposure_data["obs_time"])):
+
+        exposure_data_obs = exposure_data[obs]
+
+        if exposure_data_obs["band"] not in bands:
+            continue
+
+        if "bkg_noise" in exposure_data_obs.keys():
+            std_gaussian_noise = exposure_data_obs["bkg_noise"]
+        else:
+            std_gaussian_noise = None
+
+        lens_images = lens_image(
+            lens_class,
+            band=exposure_data_obs["band"],
+            mag_zero_point=exposure_data_obs["zero_point"],
+            num_pix=num_pix,
+            psf_kernel=exposure_data_obs["psf_kernel"],
+            transform_pix2angle=transform_pix2angle,
+            exposure_time=exposure_data_obs["expo_time"],
+            t_obs=exposure_data_obs["obs_time"],
+            std_gaussian_noise=std_gaussian_noise,
+        )
+
+        final_image.append(lens_images)
+
+    lens_col = Column(name="lens", data=final_image)
+    final_image_col = Column(name="injected_lens", data=final_image)
+
+    # Create a new Table with only the bands of interest
+    expo_bands = np.array([b for b in exposure_data["band"]])
+    mask = np.isin(expo_bands, bands)
+    exposure_data_new = exposure_data[mask]
+
+    # if len(exposure_data_new) > 0:
+    exposure_data_new.add_columns([lens_col, final_image_col])
+
+    return exposure_data_new

slsim/Plots/plot_functions.py

@@ -7,7 +7,13 @@
 
 
 def create_image_montage_from_image_list(
-    num_rows, num_cols, images, time=None, image_type="other", image_center=None
+    num_rows,
+    num_cols,
+    images,
+    time=None,
+    band=None,
+    image_type="other",
+    image_center=None,
 ):
     """Creates an image montage from an image list.
 
@@ -16,6 +22,8 @@ def create_image_montage_from_image_list(
     :param images: list of images
     :param time: array of observation time for point source images. If None, considers
         static case.
+    :param band: array of bands corresponding to the observations. If None, does not
+        display any information regarding the band.
     :param image_type: type of the provided image. It could be 'dp0' or any other name.
     :param image_center: center of the source images.
     :type image_center: array. eg: for two image, it should be like
@@ -32,6 +40,10 @@ def create_image_montage_from_image_list(
     global_min = min(all_min)
     global_max = max(all_max)
 
+    # If band is one string, extend to list
+    if isinstance(band, str):
+        band = [band] * len(images)
+
     fig, axes = plt.subplots(num_rows, num_cols, figsize=(num_cols * 3, num_rows * 3))
 
     for i in range(num_rows):
@@ -61,6 +73,17 @@ def create_image_montage_from_image_list(
                         horizontalalignment="left",
                         transform=axes[i, j].transAxes,
                     )
+                if band is not None:
+                    axes[i, j].text(
+                        0.05,
+                        0.10,
+                        f"Band: {band[i * num_cols + j]}",
+                        fontsize=10,
+                        color="white",
+                        verticalalignment="top",
+                        horizontalalignment="left",
+                        transform=axes[i, j].transAxes,
+                    )
                 if image_center is not None:
                     for k in range(len(image_center)):
                         axes[i, j].scatter(

slsim/image_simulation.py

@@ -605,8 +605,9 @@ def lens_image_series(
     to simulate time series images of a lens.
 
     :param lens_class: Lens() object
-    :param band: imaging band
-    :param mag_zero_point: list of magnitude zero point for sqeuence of exposure
+    :param band: imaging band (or list of bands).
+        if float: assumed to apply to the full image series.
+    :param mag_zero_point: list of magnitude zero point for sequence of exposure
     :param num_pix: number of pixels per axis
     :param psf_kernels: list of psf kernel for each exposure.
     :param transform_pix2angle: list of transformation matrix (2x2) of pixels into
@@ -621,13 +622,18 @@ def lens_image_series(
     :param with_deflector: If True, simulates image with deflector.
     :return: list of series of images of a lens
     """
+
+    # If band is one string, extend to list
+    if isinstance(band, str):
+        band = [band] * len(mag_zero_point)
+
     image_series = []
-    for time, psf_kern, mag_zero, transf_matrix, expo_time in zip(
-        t_obs, psf_kernel, mag_zero_point, transform_pix2angle, exposure_time
+    for time, psf_kern, mag_zero, transf_matrix, expo_time, band_obs in zip(
+        t_obs, psf_kernel, mag_zero_point, transform_pix2angle, exposure_time, band
     ):
         image = lens_image(
             lens_class=lens_class,
-            band=band,
+            band=band_obs,
             mag_zero_point=mag_zero,
             num_pix=num_pix,
             psf_kernel=psf_kern,

test_requirements.txt

@@ -1,3 +1,4 @@
 # git+https://github.com/aymgal/coolest@main#egg=coolest
+git+https://github.com/LSSTDESC/OpSimSummaryV2
 coolest==0.1.8
 galsim

tests/test_lsst_science_pipeline.py → tests/test_LsstSciencePipeline/test_lsst_science_pipeline.py

@@ -3,7 +3,7 @@
 from astropy.table import Table
 from astropy.cosmology import FlatLambdaCDM
 from slsim.lens import Lens
-from slsim.lsst_science_pipeline import (
+from slsim.LsstSciencePipeline.lsst_science_pipeline import (
     variable_lens_injection,
     multiple_variable_lens_injection,
 )
@@ -14,10 +14,10 @@
 def pes_lens_instance():
     path = os.path.dirname(__file__)
     source_dict = Table.read(
-        os.path.join(path, "TestData/source_dict_ps.fits"), format="fits"
+        os.path.join(path, "../TestData/source_dict_ps.fits"), format="fits"
     )
     deflector_dict = Table.read(
-        os.path.join(path, "TestData/deflector_dict_ps.fits"), format="fits"
+        os.path.join(path, "../TestData/deflector_dict_ps.fits"), format="fits"
     )
 
     cosmo = FlatLambdaCDM(H0=70, Om0=0.3)
@@ -41,7 +41,7 @@ def test_variable_lens_injection(pes_lens_instance):
     lens_class = pes_lens_instance
     path = os.path.dirname(__file__)
     expo_data = Table.read(
-        os.path.join(path, "TestData/expo_data_1.fits"), format="fits"
+        os.path.join(path, "../TestData/expo_data_1.fits"), format="fits"
     )
     transf_matrix_single = np.array([[0.2, 0], [0, 0.2]])
     transform_matrices = [transf_matrix_single.copy() for _ in range(len(expo_data))]
@@ -59,10 +59,10 @@ def test_multiple_variable_lens_injection(pes_lens_instance):
     lens_class = [pes_lens_instance, pes_lens_instance]
     path = os.path.dirname(__file__)
     expo_data_1 = Table.read(
-        os.path.join(path, "TestData/expo_data_1.fits"), format="fits"
+        os.path.join(path, "../TestData/expo_data_1.fits"), format="fits"
     )
     expo_data_2 = Table.read(
-        os.path.join(path, "TestData/expo_data_2.fits"), format="fits"
+        os.path.join(path, "../TestData/expo_data_2.fits"), format="fits"
     )
     expo_data = [expo_data_1, expo_data_2]
     transf_matrix_single = np.array([[0.2, 0], [0, 0.2]])