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 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.4418965816732188),
  'semiminor_sigma': np.float64(1.4314624444996362),
  'orientation': np.float64(-61.69453409913115)},
 'r-band': {'semimajor_sigma': np.float64(1.4418592935767625),
  'semiminor_sigma': np.float64(1.4314336530407448),
  'orientation': np.float64(-61.555816532699936)},
 'i-band': {'semimajor_sigma': np.float64(1.4419935689691024),
  'semiminor_sigma': np.float64(1.4314931321136346),
  'orientation': np.float64(-64.70440808182138)}}

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	437.844648	10.379990	99.996495	9.038971	60309.999797
1	440.881519	10.465685	99.445609	9.132651	60309.999809
2	434.890966	11.982308	99.865158	10.849745	60309.999820
3	437.757840	11.089946	100.465640	9.842834	60309.999832
4	438.789814	12.323376	100.538294	11.214403	60309.999843
5	436.318764	11.281866	100.015615	10.067294	60309.999855
6	437.089504	11.516709	100.865563	10.323189	60309.999866
7	437.848098	11.828927	100.547841	10.669860	60309.999878
8	440.728774	11.171752	99.742155	9.930760	60309.999890
9	444.063361	11.580025	99.858273	10.380200	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
    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 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,
        error=error,
        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	437.927508	9.786209	99.913635	8.350371	60309.999797
1	441.965180	11.001050	98.361948	9.741553	60309.999809
2	433.558514	12.660819	101.197610	11.594723	60309.999820
3	437.703572	10.901601	100.519908	9.630130	60309.999832
4	438.619174	11.924683	100.708934	10.774753	60309.999843
5	436.780879	10.284865	99.553500	8.935792	60309.999855
6	437.575242	11.528630	100.379825	10.336485	60309.999866
7	437.471952	12.744118	100.923987	11.676256	60309.999878
8	440.264885	11.707050	100.206044	10.529338	60309.999890
9	444.294972	10.604625	99.626662	9.279529	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 error [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.