Reduction

This notebook will demonstrate basic OPTICAM data reduction using opticam.

Performance Tip

Before starting, it’s worth noting that opticam uses the multiprocessing Python module to parallelize data reduction across multiple CPU cores. For the best scaling, it’s therefore a good idea to disable automatic parallelization by any underlying libraries (like `numpy <https://numpy.org>`__) using some combination of the following:

import os

# limit underlying math libraries to a single thread for better multiprocessing performance
os.environ['OMP_NUM_THREADS'] = '1'  # OpenMP
# os.environ['OPENBLAS_NUM_THREADS'] = '1'  # OpenBLAS
# os.environ['MKL_NUM_THREADS'] = '1'  # Intel Math Kernel Library
# os.environ['VECLIB_MAXIMUM_THREADS'] = '1'  # Apple Accelerate vector library

The relevant libraries will depend on your build and the specific system you’re using; in my case, I only need to include os.environ['OMP_NUM_THREADS'] = '1'.

Warning: if you do not disable automatic parallelization as shown above, cosmic ray removal performance will be extremely degraded and/or may hang.

Generating Data

Before we can beging reducing data, we need some data to reduce. opticam provides a routine for generating some dummy observations via the generate_observations() function:

from pathlib import Path

import opticam


out_dir = Path('out')

opticam.generate_observations(
    out_directory=out_dir / 'data',  # path to the directory where the generated data will be saved
    circular_aperture=False,  # disable circular aperture shadow
    n_images=20,
    )

[OPTICAM] variable source is at (131, 115)
[OPTICAM] variability RMS: 0.02 %
[OPTICAM] variability frequency: 0.135 Hz
[OPTICAM] variability phase lags:
    [OPTICAM] g-band: 0.000 radians
    [OPTICAM] r-band: 1.571 radians
    [OPTICAM] i-band: 3.142 radians

Generating observations: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:05<00:00]

These dummy observations come in three filters: \(g\), \(r\), and \(i\), and will be used in many of the guided tutorials. We can also see that there is a variable source at (131, 115), which we will use as our source of interest. We’re also given the fractional RMS amplitude of the variability, its frequency, the phase lags between the different filters. These values will be useful for sanity checking our results in the Timing Methods Tutorial.

There are a number of parameters that can be tweaked when calling generate_observations(), such as the number of images, the binning scale of the images, and whether a circular aperture shadow is applied to the images. For this example, I’ve only generated a small number of images to minimise storage and ensure this notebook runs quickly.

Defining an instrument

In opticam, Instruments are used to interface with data products. An Instrument defines the keywords and necessary routines required for opticam to interpret the headers of images from a particular instrument. Since opticam was developed primarily for use with OPTICAM-MX, the data generated by the generate_observations() function are compatible with the OPTICAM_MX instrument:

instrument = opticam.OPTICAM_MX()

All of opticam’s reduction routines require an Instrument instance, and we have tried to make it straightforward for users to define custom Instrument instances for their instruments (see the instruments tutorial for more details).

Pre-reduction checks

Before we reduce any data, it’s a good idea to run opticam.scan_data(), especially if you’re unfamiliar with the data. This scans the image headers and summarises the relevant information for you:

opticam.scan_data(
    out_directory=out_dir / 'reduced',
    data_directory=out_dir / 'data',
    instrument=instrument,
)

[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]

[OPTICAM] Checking instrument OPTICAM_MX.
[OPTICAM] OPTICAM_MX sucessfully passed all checks.

[OPTICAM] Parsing file headers: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:02<00:00]

Binning: 4x4
Filters: 1:g, 2:r, 3:i
20 1:g images.
20 2:r images.
20 3:i images.

As we can see, the images use the 4x4 binning mode, we have \(g\)-, \(r\)-, and \(i\)-band images, and there are 20 images for each filter. If there were more than three filters, this would suggest that the data directory contains images from different observations (since OPTICAM-MX can only observe in three bands at once). If this were the case, it would be a good idea to organise the data directory to ensure that it contained only images from a single observation. Since there were no issues with the data, however, we can instance proceed with reduction by defining a Reducer instance.

Defining a Reducer

When defining a Reducer, there are two required parameters: data_directory and out_directory. data_directory should be a path to a directory that contains a series of FITS files. out_directory determines the directory to which any output files are written. If out_directory does not exist, opticam will attempt to create it. We can also apply flat-field corrections to our images by passing a FlatFieldCorrector, though we will not do this here (see the corrections tutorial for details on applying flat-field corrections). Another useful parameter is remove_cosmic_rays, which can be used to clip cosmic rays from your images. If remove_cosmic_rays=True, cosmic rays are clipped using `astroscrappy <https://astroscrappy.readthedocs.io/en/latest/index.html>`__’s L.A.Cosmic implementation. By default, remove_cosmic_rays is set to False.

Let’s initialise a Reducer using its default paramters:

reducer = opticam.Reducer(
    data_directory=out_dir / 'data',  # path to the data
    out_directory=out_dir / 'reduced',  # path to where the output will be saved
    instrument=instrument,  # the instrument (defaults to OPTICAM_MX())
    )

[OPTICAM] out/reduced not found, attempting to create ...
[OPTICAM] out/reduced created.

[OPTICAM] Scanning data directory: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:00<00:00]

[OPTICAM] Checking instrument OPTICAM_MX.
[OPTICAM] OPTICAM_MX sucessfully passed all checks.

[OPTICAM] Parsing file headers: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:02<00:00]

[OPTICAM] Binning: 4x4
[OPTICAM] Filters: 1:g, 2:r, 3:i
[OPTICAM] 20 1:g images.
[OPTICAM] 20 2:r images.
[OPTICAM] 20 3:i images.

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/diag/header_times.pdf.

After creating a Reducer instance, opticam.scan_data() is ran to ensure there are no issues with the data; as mentioned above, this may result in errors that will need to be resolved. In this case, however, we can see that there are no errors, and so we can proceed with reduction.

We can also see that our chosen instrument, OPTICAM_MX passed the internal checks. See the custom instruments tutorial for guidance on how to create custom instruments for use with opticam, and fix potential errors that might arise when doing so.

Create Source Catalogs

The next step is to create source catalogs for each camera. Creating source catalogs requires aligning each camera’s images to track sources over time and ensure consistent labelling. The way in which images are aligned can be customised. By default, images will be aligned using transform_type='affine' which uses the `astroalign <https://astroalign.quatrope.org/en/latest/>`__ package. Alternatively, passing transform_type='translation' will compute simple (x, y) translations between images. The number of reference sources can also be changed by passing the desired value to n_alignment_sources. When using transform_type='translation', n_alignment_sources must be \(\geq 1\) while transform_type='affine' requires n_alignment_sources be \(\geq 3\).

To create the source catalogs, we must call the create_catalogs() method:

reducer.create_catalogs()

[OPTICAM] Creating source catalogs.

[OPTICAM] Aligning 1:g images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]

[OPTICAM] [OPTICAM] Done.
[OPTICAM] [OPTICAM] 20 image(s) aligned.
[OPTICAM] [OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Aligning 2:r images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]

[OPTICAM] [OPTICAM] Done.
[OPTICAM] [OPTICAM] 20 image(s) aligned.
[OPTICAM] [OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Aligning 3:i images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]

[OPTICAM] [OPTICAM] Done.
[OPTICAM] [OPTICAM] 20 image(s) aligned.
[OPTICAM] [OPTICAM] 0 image(s) could not be aligned.

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/cat/catalogs.pdf.

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/diag/background.pdf.

We can see that six sources have been identified in each of the three sets of images. In this example, the source labelling is consistent across the catalogs because there are no field-of-view or pixel-scale differences between the simulated cameras. In practise, the source labelling will not usually be consistent across the catalogs, and so care must be taken when performing differential photometry that the same sources are being used for each filter (more on this later).

It’s also worth discussing opticam’s labelling convention. For compatibility with multi-camera instruments, like OPTICAM-MX, opticam groups images according to both the camera and filter used to take the image. In the above catalogue, you can see the catalogues are titled “1:g”, “2:r”, and “3:i”. The value before the colon represents the camera, which may simply be the name of the instrument if only has one camera, while the value after the colon represents the filter. In the case of OPTICAM-MX, filters are unique to each camera, though this may not be the case for all instruments. See the custom instruments tutorial for more details.

Performing Photometry

With our catalogs defined, we can now perform photometry. Performing photometry in opticam requires a photometer object that inherits from opticam.photometers.BasePhotometer. Currently, opticam provides two photometers: AperturePhotometer, for performing simple aperture photometry, and OptimalPhotometer, for performing optimal photometry (as described in Naylor 1998). Both photometers can be configured to perform “forced photometry”, and can use either the Reducer’s 2D background estimator or estimate the local background around each source using an annulus. In this example, I’ll show a couple of different photometry configurations:

# aperture photometer with local background estimations
default_annulus_photometer = opticam.AperturePhotometer(
    local_background_estimator=opticam.DefaultLocalBackground(),  # use the default local background estimator
)

# optimal photometer
# implements the method described in Naylor 1998, MNRAS, 296, 339-346
optimal_photometer = opticam.OptimalPhotometer()

Once a photometer has been defined, it can be passed to the photometry() method of Catalog to compute the raw light curves:

reducer.photometry(default_annulus_photometer)  # using the aperture photometer with local background estimations

[OPTICAM] [OPTICAM] Photometry results will be saved to lcs/aperture_annulus in out/reduced.

[OPTICAM] Performing photometry on 1:g images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]
[OPTICAM] Performing photometry on 2:r images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]
[OPTICAM] Performing photometry on 3:i images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:01<00:00]

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/diag/aperture_annulus_rms_vs_median.pdf.

When you perform photometry on a catalog, the directory to which the light curves are saved will depend on how the photometer is configured:

1. If the photometer is defined with forced=True, the directory will have a “forced” prefix.

2. If a local background estimator is passed to the photometer’s local_background_estimator parameter, then the directory will also have an “annulus” suffix.

In the above case, we used the AperturePhotometer with a local background estimator and so we can see that the raw light curves have been saved to a lcs/aperture_annulus directory. Let’s compare this to the subdirectory for light curves produced by our optimal_photometer:

reducer.photometry(optimal_photometer)  # using the optimal photometer

[OPTICAM] [OPTICAM] Photometry results will be saved to lcs/optimal in out/reduced.

[OPTICAM] Performing photometry on 1:g images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:02<00:00]
[OPTICAM] Performing photometry on 2:r images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:02<00:00]
[OPTICAM] Performing photometry on 3:i images: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:02<00:00]

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/diag/optimal_rms_vs_median.pdf.

In this case there is no prefix nor suffix, and only the name of the photometer class is included.

Computing Relative Light Curves

We now have some raw light curves using a couple of different photometry configurations. However, raw light curves contain a lot of atmospheric and systematic variability, and so we often want to compute the relative light curve between our source of interest and some comparison sources to reduce this atmospheric/systematic variability.

In this example, let’s say that Source 6 is our target of interest (since we know that it’s varying), and use Source 5 for comparison. In this example, the choice of comparison source(s) is arbitraray (since there are no atmospheric/systematic variations in our simulated data). In practise, however, choosing suitable comparison sources is vital for obtaining quality light curves, though a discussion of this is beyond the scope of this tutorial.

Let’s now produce a relative light curve for Source 6 using the “aperture_annulus” light curves created by default_annulus_photometer. First, however, we need to initialise a DifferentialPhotometer object. When initialising a DifferentialPhotometer object, we need to pass the directory path to the reduced data created by Reducer:

dphot = opticam.DifferentialPhotometer(
    out_directory=out_dir / 'reduced',  # same as the catalog's out_directory
    show_plots=True,  # show plots (useful for diagnosis and debugging)
    )

[OPTICAM] Keys: 1:g, 2:r, 3:i

When initialising a DifferentialPhotometer object, the source catalogs are output for convenience (unless show_plots=False). We can now create our relative light curves:

target = 2  # source of interest
comparisons = [1]  # comparison sources

annulus_analyzer = dphot.get_relative_light_curve(
    '1:g',  # filter for which to compute the relative light curve
    target,  # source of interest
    comparisons,  # comparison sources
    phot_label='aperture_annulus',  # label for the photometry results
    prefix='test',  # prefix for the output files (e.g., the name of the target source)
    match_other_cameras=True,  # match sources across cameras
    show_diagnostics=True,  # show diagnostic plots (useful for diagnosis and debugging)
    )

[OPTICAM] 1:g target ID 2 was matched to 2:r target ID 2
[OPTICAM] 2:r comparison ID 1 was matched to 2:r comparison ID 1
[OPTICAM] 1:g target ID 2 was matched to 3:i target ID 2
[OPTICAM] 3:i comparison ID 1 was matched to 3:i comparison ID 1

While the relative light curves are generated, a number of diagnostic plots are rendered. These diagnostic plots can be suppressed by passing show_diagnostics=False. Regardless of the value of show_diagnostics, these diagnostic plots are saved to out_directory/relative_light_curves/diag. The resulting relative light curves are also shown for convenience. After relative light curves have been computed, an Analyzer object is returned (see the timing methods tutorial for more details on this).

Above, we used the get_relative_light_curve() method to automatically match sources across filters by setting match_other_cameras=True. However, this can misidentify sources, so care should be taken to check the correct sources are identified. In this case, we can see that the identified sources are correct, and so we don’t have to manually create relative light curves for each filter. For the optimal light curves, however, I will show how to merge Analyzer instances using the join() method should you ever have to match the sources manually. First, let’s define the individual Analyzer instances:

g_band_optimal_analyzer = dphot.get_relative_light_curve(
    '1:g',
    target,
    comparisons,
    phot_label='optimal',
    prefix='test',
    match_other_cameras=False,
    show_diagnostics=False,
    )

r_band_optimal_analyzer = dphot.get_relative_light_curve(
    '2:r',
    target,
    comparisons,
    phot_label='optimal',
    prefix='test',
    match_other_cameras=False,
    show_diagnostics=False,
    )

i_band_optimal_analyzer = dphot.get_relative_light_curve(
    '3:i',
    target,
    comparisons,
    phot_label='optimal',
    prefix='test',
    match_other_cameras=False,
    show_diagnostics=False,
    )

Now let’s combine them all into a single instance:

# join g-band with r-band
optimal_analyzer = g_band_optimal_analyzer.join(r_band_optimal_analyzer)

# join g-band and r-band with i-band
optimal_analyzer = optimal_analyzer.join(i_band_optimal_analyzer)

To check the Analyzer instances have merged properly, let’s try plotting the light curves:

fig = optimal_analyzer.plot()

[OPTICAM] Plot saved to /tmp/tmp_k5k9d11/out/reduced/plots/test_optimal_light_curves.pdf.

As we can see, all three light curves have been combined successfully combined into a single Analyzer instance.

That concludes the reduction tutorial for opticam! If you’re using opticam to reduce data from OPTICAM-MX, I recommend checking out the applying corrections and photometry tutorials next. To use opticam with another instrument, check out the custom instruments tutorial. To explore the quick-look timing analysis routines available in opticam, see the timing methods tutorial.