Local background estimation

opticam supports specification of local background estimators, which can improve the background estimation in cases where the background varies considerably across an image. In this notebook, I will demonstrate how to define a local background estimator for use with opticam, as well as explain opticam’s default local background estimator.

Generating Data

First thing’s first, we need to generate some dummy data:

import opticam

opticam.generate_observations(
    out_dir='local_backgrounds_tutorial/data',
    n_images=10,
    circular_aperture=False,
    )

[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:02<00:00]

Now we can use these dummy observations to compute and compare local backgrounds.

Generating Catalogs

Before we can look at the local background estimator, we need to create our source catalogs:

reducer = opticam.Reducer(
    data_directory='local_backgrounds_tutorial/data',
    out_directory='local_backgrounds_tutorial/reduced',
    remove_cosmic_rays=False,
    )

reducer.create_catalogs()

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

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

[OPTICAM] Binning: 4x4
[OPTICAM] Filters: g-band, r-band, i-band
[OPTICAM] 10 g-band images.
[OPTICAM] 10 r-band images.
[OPTICAM] 10 i-band images.

[OPTICAM] Creating source catalogs

[OPTICAM] Aligning g-band images: 100%|██████████|[00:00<00:00]

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

[OPTICAM] Aligning r-band images: 100%|██████████|[00:00<00:00]

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

[OPTICAM] Aligning i-band images: 100%|██████████|[00:00<00:00]

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

DefaultLocalBackground

Now we look at the default local background estimator, DefaultLocalBackground. The default local background estimator uses an elliptical annulus that scales the semimajor_sigma and semiminor_sigma parameters of the catalogue’s psf_params attribute, which is defined when the source catalogs are created (or read from file if the catalogs were created previously). Let’s see what this dictionary looks like:

reducer.psf_params

{'g-band': {'semimajor_sigma': np.float64(1.4116196327047383),
  'semiminor_sigma': np.float64(1.3723953838086822),
  'orientation': np.float64(-77.33308094038188)},
 'r-band': {'semimajor_sigma': np.float64(1.4116462274747996),
  'semiminor_sigma': np.float64(1.372388390108369),
  'orientation': np.float64(-77.31455733431797)},
 'i-band': {'semimajor_sigma': np.float64(1.4116806391924892),
  'semiminor_sigma': np.float64(1.3723754235701116),
  'orientation': np.float64(-77.31313781077361)}}

We can see that all three cameras have the same PSF parameters, since these simulated data are the same for each camera. In practise, the PSF parameters will generally be different between cameras.

We can see each camera has three PSF parameters:

• semimajor_sigma: the standard deviation of the PSF along the semi-major axis (assuming a 2D Gaussian PSF)

• semiminor_sigma: the standard deviation of the PSF along the semi-minor axis (assuming a 2D Gaussian PSF)

• orientation: the orientation of the PSF in degrees

When defining a custom local background estimator, it’s a good idea to use these values to define your annulus (more on this later). For now, let’s perform photometry using the default local background estimator:

phot = opticam.OptimalPhotometer(
    local_background_estimator=opticam.DefaultLocalBackground(),
)

reducer.photometry(phot)

[OPTICAM] Photometry results will be saved to lcs/optimal_annulus in local_backgrounds_tutorial/reduced.

[OPTICAM] Performing photometry on g-band images: 100%|██████████|[00:00<00:00]
[OPTICAM] Performing photometry on r-band images: 100%|██████████|[00:00<00:00]
[OPTICAM] Performing photometry on i-band images: 100%|██████████|[00:00<00:00]

Let’s take a look at the resulting light curve file for a random source:

import pandas as pd

default_g_band_source_2_df = pd.read_csv(
    'local_backgrounds_tutorial/reduced/lcs/optimal_annulus/g-band_source_2.csv',
    )

default_g_band_source_2_df

	flux	flux_err	background	background_err	BMJD
0	258.504068	4.326146	100.956145	10.831136	60309.999797
1	262.614531	4.200837	100.720807	9.233961	60309.999809
2	254.273346	4.234523	100.502657	10.302317	60309.999820
3	261.075775	4.176462	102.102402	8.892589	60309.999832
4	253.055360	4.225237	101.188853	10.259006	60309.999843
5	242.286897	4.176203	99.968004	10.426320	60309.999855
6	252.002279	4.194522	99.665602	10.000045	60309.999866
7	255.867419	4.174600	99.591541	9.564518	60309.999878
8	251.248432	4.128897	99.170680	9.489917	60309.999890
9	263.289786	4.175057	99.363499	9.015452	60309.999901

We can see that the light curve file has background and background_err columns, which are only included if a local background estimator was used to compute the light curve. The average background is about 100, which is the expected value, and the average background error is about \(\sqrt{100} = 10\), which is also the expected value, showing that the default local background estimator is working well.

Custom Local Background

The default local background estimator uses the mean and standard deviation in the sigma-clipped annulus to estimate the background and its error, respectively. We will therefore define a custom local background estimator that works slightly differently, using the median instead of the mean, and the median absolute deviation instead of the standard deviation.

To define a custom local background estimator, you should either inherit from the BaseLocalBackground base class, or define a class that implements a __call__() method with the following input parameters:

• data: the image data as a 2D numpy array

• error: the error in the image data as a 2D numpy array

• position: the position of a source as a tuple (x, y)

• semimajor_axis: the semi-major standard deviation of the PSF (assumed to be a 2D Gaussian)

• semiminor_axis: the semi-minor standard deviation of the PSF (assumed to be a 2D Gaussian)

• orientation: the orientation of the PSF in degrees

For simplicity, I will inherit from the base class since this is how I’d recommend defining a custom local background estimator:

from photutils.aperture import EllipticalAnnulus, ApertureStats
from astropy.stats import SigmaClip

def custom_local_background(
    data,  # the image data as a 2D numpy array
    position,  # the position of a source as a tuple (x, y)
    semimajor_axis,  # the semimajor axis of the PSF in pixels (assuming a 2D Gaussian PSF)
    semiminor_axis,  # the semiminor axis of the PSF in pixels (assuming a 2D Gaussian PSF)
    orientation,  # the orientation of the PSF in degrees
    ):

    # set inner and outer radius scale of annulus
    r_in_scale = 5
    r_out_scale = 7.5

    # define the annulus using photutils
    annulus = EllipticalAnnulus(
        positions=position,
        a_in=r_in_scale * semimajor_axis,
        a_out=r_out_scale * semimajor_axis,
        b_in=r_in_scale * semiminor_axis,
        b_out=r_out_scale * semiminor_axis,
        theta=orientation,
    )

    # compute the statistics in the annulus
    stats = ApertureStats(
        data,
        annulus,
        sigma_clip=SigmaClip(sigma=3.5, maxiters=10),
    )

    median_local_background = float(stats.median)

    # median absolute deviation (MAD) of the local background
    mad_std_local_background = float(stats.mad_std)

    return median_local_background, mad_std_local_background

When inheriting from the BaseLocalBackground base class, there is no need to implement an __init__() method.

The order of the parameters in the __call__() method is important, since opticam’s photometry routines assume a local background estimator’s __call__() method takes the above parameters in the above order. These parameters are handled by the catalogue, and so you don’t need to worry about computing them. In particular, the semimajor_axis, semiminor_axis, and orientation parameters are passed from the catalogue’s psf_params dictionary using the appropriate values for the camera corresponding to the image, so you don’t need to worry about unpacking this dictionary.

Finally, there is no requirement to use photutils for estimating local backgrounds, but photutils does make it a lot easier.

Let’s now perform photometry using our custom local background estimator:

custom_phot = opticam.OptimalPhotometer(
    local_background_estimator=custom_local_background,
    )

reducer.photometry(custom_phot, overwrite=True)

[OPTICAM] Photometry results will be saved to lcs/optimal_annulus in local_backgrounds_tutorial/reduced.

[OPTICAM] Performing photometry on g-band images: 100%|██████████|[00:00<00:00]
[OPTICAM] Performing photometry on r-band images: 100%|██████████|[00:00<00:00]
[OPTICAM] Performing photometry on i-band images: 100%|██████████|[00:00<00:00]

custom_g_band_source_2_df = pd.read_csv(
    'local_backgrounds_tutorial/reduced/lcs/optimal_annulus/g-band_source_2.csv',
    )

custom_g_band_source_2_df

	flux	flux_err	background	background_err	BMJD
0	258.928104	4.446490	100.532109	11.941155	60309.999797
1	262.121385	4.200649	101.213953	9.258552	60309.999809
2	253.525722	4.231791	101.250281	10.311268	60309.999820
3	260.401429	4.217775	102.776748	9.391214	60309.999832
4	251.761243	4.160446	102.482970	9.659930	60309.999843
5	242.341960	4.135236	99.912942	10.018176	60309.999855
6	252.683782	4.232157	98.984098	10.346014	60309.999866
7	256.336001	4.045020	99.122959	8.067330	60309.999878
8	251.334000	4.121537	99.085112	9.407160	60309.999890
9	263.054171	4.254553	99.599114	9.890801	60309.999901

Again, the background and its error are both close to the expected values. Let’s therefore compare the local background estimations more clearly:

from matplotlib import pyplot as plt

fig, axes = plt.subplots(
    nrows=2,
    tight_layout=True,
    sharex=True,
    )

axes[0].plot(
    default_g_band_source_2_df['BMJD'],
    default_g_band_source_2_df['background'],
    label='Default',
)
axes[0].plot(
    custom_g_band_source_2_df['BMJD'],
    custom_g_band_source_2_df['background'],
    ls='--',
    label='Custom',
)
axes[0].axhline(
    100,
    color='k',
    zorder=0,
)

axes[1].plot(
    default_g_band_source_2_df['BMJD'],
    default_g_band_source_2_df['background_err'],
    label='Default',
)
axes[1].plot(
    custom_g_band_source_2_df['BMJD'],
    custom_g_band_source_2_df['background_err'],
    ls='--',
    label='Custom',
)
axes[1].axhline(
    10,
    color='k',
    zorder=0,
)

axes[1].set_xlabel('BMJD [days]')
axes[0].set_ylabel('Background [Counts]')
axes[1].set_ylabel('Background RMS [Counts]')

axes[0].legend()
axes[1].legend()

plt.show()

As we can see, both methods are around the expected values, and so either local background estimator would likely give pretty good results in practise.

That concludes the local background tutorial for opticam! The default local background estimator should be “good enough” in most cases, but we have seen how to implement custom local background estimators if required. Once again, I refer to the excellent photutils documentation for further details on implementing custom local background estimators: https://photutils.readthedocs.io/en/stable/user_guide/aperture.html#local-background-subtraction.