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:04<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.413218686650065),
  'semiminor_sigma': np.float64(1.3790703013931664),
  'orientation': np.float64(-75.23359473318958)},
 'r-band': {'semimajor_sigma': np.float64(1.4132191514395016),
  'semiminor_sigma': np.float64(1.3790698107870591),
  'orientation': np.float64(-75.23520649261539)},
 'i-band': {'semimajor_sigma': np.float64(1.4132373104484879),
  'semiminor_sigma': np.float64(1.3790481577857459),
  'orientation': np.float64(-75.23624808916581)}}

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

	BMJD	flux	flux_err	background	background_err
0	60309.999797	259.655900	4.780290	101.069653	10.799642
1	60309.999809	263.028460	4.676428	100.614186	9.383994
2	60309.999820	251.555261	4.680888	100.500853	10.390129
3	60309.999832	261.706478	4.655758	102.183199	8.896059
4	60309.999843	252.587626	4.647972	101.334062	9.896461
5	60309.999855	239.812840	4.607177	100.130503	10.254058
6	60309.999866	253.302115	4.641585	99.434177	9.829064
7	60309.999878	256.523696	4.625194	99.544523	9.452757
8	60309.999890	248.272053	4.573641	99.356804	9.532862
9	60309.999901	263.470737	4.633110	99.240122	9.023487

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

	BMJD	flux	flux_err	background	background_err
0	60309.999797	259.961588	4.838659	100.763965	11.418520
1	60309.999809	262.492278	4.676677	101.150368	9.387024
2	60309.999820	250.787204	4.680005	101.268910	10.380381
3	60309.999832	260.873759	4.711841	103.015917	9.592005
4	60309.999843	251.363319	4.614560	102.558369	9.505810
5	60309.999855	240.030602	4.583176	99.912741	9.987173
6	60309.999866	253.752504	4.683014	98.983788	10.299108
7	60309.999878	256.814488	4.518520	99.253731	8.091121
8	60309.999890	248.345301	4.571694	99.283556	9.509957
9	60309.999901	263.119422	4.707714	99.591437	9.924216

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.