Local background estimation

phoptic 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 phoptic, as well as explain phoptic’s default local background estimator.

Generating Data

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

from pathlib import Path

import phoptic


out_dir = Path('out')

phoptic.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:01<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 = phoptic.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/tmp_gahbtm5/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:02<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/tmp_gahbtm5/out/reduced/cat/catalogs.pdf.

[OPTICAM] Plot saved to /tmp/tmp_gahbtm5/out/reduced/diag/systematics.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_axis': np.float64(1.413218686650065),
  'semiminor_axis': np.float64(1.3790703013931664),
  'orientation': np.float64(276.00161285522074)},
 '2:r': {'semimajor_axis': np.float64(1.4132191514395016),
  'semiminor_axis': np.float64(1.3790698107870591),
  'orientation': np.float64(275.9819920134871)},
 '3:i': {'semimajor_axis': np.float64(1.4132373104484879),
  'semiminor_axis': np.float64(1.3790481577857459),
  'orientation': np.float64(276.0491613402856)}}

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 = phoptic.OptimalPhotometer(
    local_background_estimator=phoptic.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:03<00:00]
[OPTICAM] Performing photometry on 2:r images: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:03<00:00]
[OPTICAM] Performing photometry on 3:i images: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:03<00:00]

[OPTICAM] Plot saved to /tmp/tmp_gahbtm5/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	524.170037	8.450782	100.975446	10.847552
1	60309.999814	529.932059	8.150545	100.583883	9.409263
2	60309.999826	513.412788	8.239396	100.554593	10.304351
3	60309.999837	516.973537	7.986058	102.256224	8.902649
4	60309.999849	513.618897	8.158841	101.243097	9.932992
5	60309.999861	486.294126	8.052939	99.703673	10.229240
6	60309.999872	512.419160	8.135108	99.289159	9.816574
7	60309.999884	520.007481	8.107673	99.366936	9.506616
8	60309.999895	508.965617	8.006183	99.385057	9.449908
9	60309.999907	535.980101	8.095239	99.211028	9.012000

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 out/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 out/custom_local_background.py

Due to compatibility issues, routines defined inside of IPython notebooks cannot be used by multiprocessing (used by phoptic). 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 in our out directory. 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 out.custom_local_background import custom_local_background


custom_phot = phoptic.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:03<00:00]
[OPTICAM] Performing photometry on 2:r images: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:03<00:00]
[OPTICAM] Performing photometry on 3:i images: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████|[00:03<00:00]

[OPTICAM] Plot saved to /tmp/tmp_gahbtm5/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.681022	8.674337	100.676065	11.760272
1	60309.999814	528.653084	8.070145	101.083609	9.071487
2	60309.999826	511.291899	8.223679	101.372644	10.280293
3	60309.999837	515.206137	8.131329	103.015917	9.592005
4	60309.999849	510.117002	8.140133	102.558369	9.924983
5	60309.999861	487.331359	7.862913	99.527012	9.383600
6	60309.999872	513.351212	8.225331	98.891797	10.191924
7	60309.999884	520.744030	7.993161	99.173442	8.968707
8	60309.999895	508.551660	7.998921	99.548126	9.426302
9	60309.999907	534.882355	8.277368	99.425698	9.864661

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 phoptic! 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.