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:

from pathlib import Path

import opticam


out_dir = Path('out')

opticam.generate_observations(
    out_directory=out_dir / '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=out_dir / 'data',
    out_directory=out_dir / 'reduced',
    )

reducer.create_catalogs()

[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] 10 1:g images.
[OPTICAM] 10 2:r images.
[OPTICAM] 10 3:i images.

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

[OPTICAM] Creating source catalogs.

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

[OPTICAM] [OPTICAM] Done.
[OPTICAM] [OPTICAM] 10 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] 10 image(s) aligned.
[OPTICAM] [OPTICAM] 0 image(s) could not be aligned.

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

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

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

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

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

{'1:g': {'semimajor_sigma': np.float64(1.413218686650065),
  'semiminor_sigma': np.float64(1.3790703013931664),
  'orientation': np.float64(-75.23359473318958)},
 '2:r': {'semimajor_sigma': np.float64(1.4132191514395016),
  'semiminor_sigma': np.float64(1.3790698107870591),
  'orientation': np.float64(-75.23520649261539)},
 '3:i': {'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] [OPTICAM] Photometry results will be saved to lcs/optimal_annulus 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/tmptpynue59/out/reduced/diag/optimal_annulus_rms_vs_median.pdf.

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(
    out_dir / 'reduced' / 'lcs' / 'optimal_annulus' / '1:g_source_2.csv',
    )

default_g_band_source_2_df

	BMJD	flux	flux_err	bkg	bkg_err
0	60309.999803	523.848788	8.282679	101.069653	10.799642
1	60309.999814	529.822842	7.979934	100.614186	9.383994
2	60309.999826	513.390410	8.106454	100.500853	10.390129
3	60309.999837	517.183106	7.825227	102.183199	8.896059
4	60309.999849	513.369114	7.994656	101.334062	9.896461
5	60309.999861	485.271491	7.913833	100.130503	10.254058
6	60309.999872	512.076309	7.982431	99.434177	9.829064
7	60309.999884	519.516636	7.933739	99.544523	9.452757
8	60309.999895	508.916331	7.870173	99.356804	9.532862
9	60309.999907	536.017228	7.928559	99.240122	9.023487

We can see that the light curve file has bkg and bkg_err columns, which are only included if a local background estimator was used to compute the light curve. We can see that 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. The default local background estimator is therefore working well in this case.

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, any function with the following input parameters can be used:

• 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 this example, I’ll use the `photutils.aperture <https://photutils.readthedocs.io/en/stable/user_guide/aperture.html>`__ package to handle the heavy lifting, and define a custom local background function that returns the median local background and its median absolute deviation:

%%writefile custom_local_background.py
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
    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),  # clip outliers from the annulus
    )

    median_local_background = stats.median

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

    return median_local_background, mad_std_local_background

Writing custom_local_background.py

Due to compatibility issues, routines defined inside of IPython notebooks cannot be used by multiprocessing (used by opticam). To overcome this issue, we have used the %%writefile magic command to write our custom_local_background() function to a temporary module called custom_local_background.py, such that it can be imported by multiprocessing. This is required in Python \(\geq\) 3.14.

There is no requirement to use photutils for estimating local backgrounds, but it does make things easier.

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

from custom_local_background import custom_local_background


custom_phot = opticam.OptimalPhotometer(
    local_background_estimator=custom_local_background,
    )

reducer.photometry(custom_phot, overwrite=True)

[OPTICAM] [OPTICAM] Photometry results will be saved to lcs/optimal_annulus 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/tmptpynue59/out/reduced/diag/optimal_annulus_rms_vs_median.pdf.

custom_g_band_source_2_df = pd.read_csv(
    out_dir / 'reduced' / 'lcs' / 'optimal_annulus' / '1:g_source_2.csv',
    )

custom_g_band_source_2_df

	BMJD	flux	flux_err	bkg	bkg_err
0	60309.999803	524.477905	8.436872	100.763965	11.418520
1	60309.999814	528.340224	7.974216	101.150368	9.387024
2	60309.999826	511.393319	8.095579	101.268910	10.380381
3	60309.999837	515.207166	7.975145	103.015917	9.592005
4	60309.999849	510.318046	7.891368	102.558369	9.505810
5	60309.999861	486.010733	7.853347	99.912741	9.987173
6	60309.999872	513.119950	8.096745	98.983788	10.299108
7	60309.999884	520.800217	7.642016	99.253731	8.091121
8	60309.999895	509.130285	7.865915	99.283556	9.509957
9	60309.999907	534.516505	8.124746	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['bkg'],
    label='Default',
)
axes[0].plot(
    custom_g_band_source_2_df['BMJD'],
    custom_g_band_source_2_df['bkg'],
    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['bkg_err'],
    label='Default',
)
axes[1].plot(
    custom_g_band_source_2_df['BMJD'],
    custom_g_band_source_2_df['bkg_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 via the photutils.aperture subpackage: https://photutils.readthedocs.io/en/stable/user_guide/aperture.html#local-background-subtraction.