Instruments

While phoptic was primarily designed for use with the OPTICAM-MX instrument, it includes an interface for compatibility with other instruments. In this notebook, we will explore this interface.

Generating observations from a new instrument

First, let’s generate some observations that are not compatible with the instruments currently supported by phoptic. For the data, we’ll generate a 1 s \(U\)-band exposure from an instrument that has a gain of 100 electrons/ADU. We’ll set the time of the observation to 00:00 on January 1st, 2026, with a pointing of 0 RA, 0 DEC:

from pathlib import Path

from photutils.datasets import make_100gaussians_image
from astropy.io import fits
import numpy as np


out_dir = Path('out')

image = make_100gaussians_image()

hdu = fits.PrimaryHDU(image)  # image data

hdu.header['FILTER'] = 'U'  # the observation filter
hdu.header['BINNING'] = '1 1'  # the image binning (i.e., 1x1)
hdu.header['GAIN'] = 100.  # the detector gain in adu/photon
hdu.header['EXPTIME'] = 1.  # the exposure time in seconds
hdu.header['RA'] = '0'  # the RA of the image in degrees
hdu.header['DEC'] = '0'  # the DEC of the image in degrees
hdu.header["DATE-OBS"] = f"2026-01-01T00:00:00.000"  # observation timestamp
hdu.header['INSTRUME'] = 'SYNTHETICAM'  # name of the instrument
hdu.header['RDNOISE'] = 0.  # the instrument's read noise in electrons/pixel
hdu.header['AIRMASS'] = 1.  # the observation's airmass

save_path = out_dir / 'data' / 'simple'  # directory to which image will be saved

# create save directory if it does not already exist
if not save_path.is_dir():
    save_path.mkdir(parents=True)

# save image
file_path = save_path / 'image.fits.gz'
hdu.writeto(
    file_path,
    overwrite=True,
    )

Let’s take a look at one of these images:

from astropy.visualization import simple_norm
from matplotlib import pyplot as plt


# get image data from file
with fits.open(file_path) as hdul:
    data = np.asarray(hdul[0].data, dtype=np.float64)

fig, ax = plt.subplots()

im = ax.imshow(
    data,
    norm=simple_norm(data, stretch='log'),
    cmap='Greys',
    origin='lower',
    )

plt.show()

The image looks good.

Let’s see what happens when we try to use the OPTICAM-MX instrument to reduce these data. Before performing reduction with a particular instrument, it’s a good idea to call the instrument’s run_checks() method. This method takes a path to an image and runs a series of checks to ensure that the data are consistent with those produced by the instrument:

import phoptic


instrument = phoptic.OPTICAM_MX()
instrument.run_checks(file_path)

In this case, we can see that we have several errors, and one warning:

• The first error is caused by OPTICAM_MX using a different header keyword for the image’s exposure time. Typically, the keyword “EXPTIME” is used to identify the image’s exposure time; this is the case for these data. However, OPTICAM_MX uses the “EXPOSURE” keyword instead.

• The second error is again caused by a header keyword mismatch, this time for the timestamp keyword. Typically, the keyword “DATE-OBS” is used to identify the image’s timestamp; this is the case for these data. However, OPTICAM_MX uses the “UT” keyword instead.

• The third error is due to the data using the \(U\)-band filter, which is not supported by OPTICAM_MX and therefore has no corresponding pixel scale.

• The final error is related to the second error: the Modified Julian Date (MJD) of the image cannot be inferred because of the timestamp keyword mismatch between these data and OPTICAM_MX.

• The warning is due to these data not containing a “DARKCURR” keyword in the header. phoptic can use the dark current value (represented by the header keyword “DARKCURR”) to infer and subtract the dark noise from the image. If this keyword does not exist, then dark noise can be removed using a DarkNoiseCorrector instance instead (see the corrections tutorial for more details). This warning can therefore be ignored in this case.

There are also a number of other issues with using OPTICAM_MX that didn’t result in any errors being raised. Let’s look at OPTICAM_MX to see what the instrument represents:

print(instrument)

We can see that OPTICAM_MX has a predefined position, a dictionary of pixel scales for each camera, a read noise value, and a number of header keywords. Many of these parameters are likely not be representative of our new instrument. Notably, location and read noise values. The location of the instrument is used when applying Barycentric corrections, and the read noise is used to calculate the total photometric error. Our new instrument therefore needs to redefine many of these parameters for accurate results.

There are two ways to define a new instrument. The simplest is to generate a configuration file template, edit the necessary values, and then use that configuration file to define a new instrument. This works well provided the instrument adheres to the conventional FITS header keywords (discussed in more detail below) and represents values in FITS format. If this is not the case, then you will need to define your instrument as a class that inherits from the phoptic.Instrument base class.

Defining an instrument from a configuration file

In this case, our data use conventional FITS header keywords and all values are given in FITS format. Explicitly:

• The corresponding filter is given by the “FILTER” keyword.

• The image binning is given by the “BINNING” keyword.

• The detector’s gain value is given by the “GAIN” keyword in electrons/ADU.

• The corresponding exposure time is given by the “EXPTIME” keyword in seconds.

• The corresponding RA is given by the “RA” keyword.

• The corresponding DEC is given by the “DEC” keyword.

• The corresponding timestamp is given by the “DATE-OBS” keyword in FITS format (YYYY-MM-DDTHH:MM:SS.sss).

If any of these values were given by different keywords or were in different units/formats, then it would be necessary to define the instrument using a custom class instead of from configuration file.

To create a configuration file, it is suggested to generate and then edit a configuration template. Configuration templates can be generated using the generate_instrument_json_template() function. This function takes a directory path and writes an instrument_template.json file to that directory:

phoptic.generate_instrument_json_template(out_dir)

Let’s take a look at the contents of this template:

import json


with open(out_dir / 'instrument_template.json', 'r') as file:
    config_dict = json.load(file)

for k, v in config_dict.items():
    print(k)
    print(v)
    print()

As we can see, the template lists a number of parameters that can be configured, including descriptions for each parameter. In practise, one would probably edit this template in their preferred text editor. For demonstrative purposes, however, I will edit the dictionary in this notebook.

For this example, I’ll ignore the location parameters (longitude, latitude, and height), as well as the read noise, since these values vary on an instrument-by-instrument basis. The only parameter that needs changing to prevent errors is pixel_scales, since the exptime_kw and timestamp_kw parameters of the template are already set to the conventional “EXPTIME” and “DATE-OBS” values, respectively:

config_dict['pixel_scales'] = {
    'SYNTHETICAM': 1.0,  # 1.0 arcsec/pixel
    }

Let’s save this edited config file under a new name:

with open(out_dir / 'new_instrument_config.json', 'w') as file:
    json.dump(
        config_dict,
        file,
        indent=4,
        )

Now we can use this file to create our instrument using the from_json() method of phoptic.Instrument. This method takes either a path to a config JSON file, or a config can be passed directly:

# using the saved config file
new_instrument = phoptic.Instrument.from_json(file_path=out_dir / 'new_instrument_config.json')

# using the config dictionary
new_instrument_alt = phoptic.Instrument.from_json(config=config_dict)

new_instrument == new_instrument_alt

assert new_instrument == new_instrument_alt

Let’s take a look at our new instrument and check it’s using the correct values:

print(new_instrument)

Looks good!

New let’s run the instrument checks for this new instrument:

new_instrument.run_checks(file_path)

No more errors! We still have a warning about the dark current keyword, but we can safely ignore this for now.

As noted above, creating instruments from a configuration template works well provided the instrument follows the convetional FITS header structure. Now let’s take a look at what to do if this is not the case:

from astropy.time import Time

image = make_100gaussians_image()

hdu = fits.PrimaryHDU(image)  # image data

hdu.header['FILTER'] = 'I'  # the observation filter
hdu.header['BINNING'] = '2x2'  # the image binning
hdu.header['GAIN'] = 500.  # the detector gain in adu/photon
hdu.header['EXPOSURE'] = 1.  # the exposure time in seconds
hdu.header['RA'] = '12 30 49.4'  # the RA of the image in degrees
hdu.header['DEC'] = '+12 23 28.0'  # the DEC of the image in degrees
hdu.header['INSTRUME'] = 'SYNTHETICAM'  # the instrument (synthetic camera)
hdu.header['RDNOISE'] = 0.  # the instrument's read noise in electrons/pixel
hdu.header['AIRMASS'] = 1.  # the observation's airmass


fits_time = f"2026-01-01T00:00:00.000"  # observation timestamp
mjd_time = float(np.asarray(Time(fits_time, format='fits').mjd))  # convert time from FITS format to MJD
hdu.header['DATE-OBS'] = mjd_time

save_path = out_dir / 'data' / 'custom'  # directory to which image will be saved

# create save directory if it does not already exist
if not save_path.is_dir():
    save_path.mkdir(parents=True)

# save image
new_file_path = save_path / 'image.fits.gz'
hdu.writeto(
    new_file_path,
    overwrite=True,
    )

These new data are very similar to the previous example, but replace the “EXPTIME” keyword with “EXPOSURE”, and store the image’s timestamp in MJD format instead of the conventional FITS format.

Let’s look at the header of this image:

header = fits.getheader(new_file_path)
print(repr(header))

Predictably, these new data will fail with both OPTICAM_MX and our newly created instrument:

instrument.run_checks(new_file_path)

new_instrument.run_checks(new_file_path)

Defining an instrument from the phoptic.Instrument base class

To reduce these data, we cannot create a new instrument from a configuration file because the timestamps are not in FITS format. Instead, we need to define a custom class that inherits from phoptic.Instrument and implement a custom get_mjd() method:

%%writefile out/custom_instrument.py

import phoptic


class CustomInstrument(phoptic.Instrument):

    def get_mjd(self, file=None, header=None):

        if header is None:
            header = file.get_header()

        return float(header[self.dateobs_kw])

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 CustomInstrument class to a temporary module called custom_instrument.py in our out directory. This is required in Python \(\geq\) 3.14.

The get_mjd() method of an Instrument instance takes two inputs: file and header. file will take a `phoptic.MEFSlice <https://phoptic.readthedocs.io/en/latest/autoapi/phoptic/mef_slice/index.html#phoptic.mef_slice.MEFSlice>`__ instance, while header will take an astropy.io.fits.Header instance. phoptic’s MEFSlice class is used internally to represent a single extension of a multi-extension FITS file (or the only extension of a standard FITS file). This allows phoptic to easily get the header of a FITS HDU by calling the get_header() method of MEFSlice. When called, phoptic will only define one of these parameters, depending on whether the image is already open or not, and so both parameters should assume default values (i.e., None). The value returned by get_mjd() should be the MJD of the image as a float.

In this example, the MJD is given by the header, so it can be returned directly. If the header gives the timestamp in an alternative format, it may be convenient to use the `astropy.time <https://docs.astropy.org/en/stable/time/index.html>`__ module to convert it to MJD.

To instance this CustomInstrument, we can again generate and edit a configuration template:

from out.custom_instrument import CustomInstrument

with open(out_dir / 'instrument_template.json', 'r') as file:
    custom_config = json.load(file)

custom_config['exptime_kw'] = 'EXPOSURE'  # change EXPTIME keyword
custom_config['pixel_scales'] = {
    'SYNTHETICAM': 1.0,  # 1.0 arcsec/pixel
    }

custom_instrument = CustomInstrument.from_json(config=custom_config)

print(custom_instrument)

custom_instrument.run_checks(new_file_path)

Alternatively, we could define the instrument by hardcoding the required parameters, which is a bit more convenient for reusing this instrument many times:

%%writefile out/custom_instrument_class.py
from astropy.coordinates import EarthLocation
from astropy.io import fits
from astropy import units as u

import phoptic


__all__ = ['CustomInstrumentClass']


class CustomInstrumentClass(phoptic.Instrument):

    def __init__(
        self,
        diameter=1.0 * u.m,  # telescope aperture
        location = EarthLocation.from_geodetic(
            lon=0.0,  # observatory longitude in degrees
            lat=0.0,  # observatory East latitude in degrees
            height=0.0,  # observatory height in meters
            ),
        pixel_scales = {
            'SYNTHETICAM': 1.0,
            },
        exptime_kw = 'EXPOSURE',
        ):

        super().__init__(
            diameter=diameter,
            location=location,
            pixel_scales=pixel_scales,
            exptime_kw=exptime_kw,
            )

    def get_mjd(
        self,
        file=None,
        header=None,
        ):

        if header is None:
            header = file.get_header()

        return float(header[self.dateobs_kw])

from out.custom_instrument_class import CustomInstrumentClass

custom_instrument_class = CustomInstrumentClass()

print(custom_instrument_class)

custom_instrument_class.run_checks(new_file_path)

While creating custom instrument classes is more advanced than editing a configuration template, it is far more configurable. Additionally, instrument classes only need to be created once, and can then be reused indefinitely (unless the instrument changes the format of its data products). If you write a class for an instrument, we encourage you to consider opening a pull request to merge your instrument into phoptic.

Gotchas

Of course, not all instruments operate in the same way and may not have some of the keywords assumed by phoptic. In these cases, values can be faked by creating methods that return fixed values, or infer values from other key words.

Example: missing keyword

As an example, let’s imagine we have an instrument that does not include a “BINNING” keyword in its data products. Assuming all images will have the same binning value, we can fix this by returning a fixed binning value for all images:

%%writefile out/instrument_without_binning.py
import phoptic


class InstrumentWithoutBinning(phoptic.Instrument):

    def get_binning(self, file=None, header=None):

        return '1x1'

from out.instrument_without_binning import InstrumentWithoutBinning

instrument_without_binning = InstrumentWithoutBinning.from_json(config=config_dict)

print(instrument_without_binning)

All methods of an Instrument take two parameters: file and header. When an Instument’s methods are called, only one of file or header are passed, so custom methods should set defaults for each parameter and define the method to return the same value regardless of which parameter is defined. In this example, however, neither parameter is needed.

Now let’s generate an image without a binning keyword:

image = make_100gaussians_image()

hdu = fits.PrimaryHDU(image)  # image data

hdu.header['FILTER'] = 'U'
hdu.header['GAIN'] = 1.
hdu.header['EXPTIME'] = 1.
hdu.header['RA'] = '12 30 49.4'
hdu.header['DEC'] = '+12 23 28.0'
hdu.header['DATE-OBS'] = f"2026-01-01T00:00:00.000"
hdu.header['INSTRUME'] = 'SYNTHETICAM'
hdu.header['RDNOISE'] = 0.
hdu.header['AIRMASS'] = 0.

new_file_path = out_dir / 'data' / 'no_binnning'
hdu.writeto(
    new_file_path,
    overwrite=True,
    )

When we check a previous instrument that would otherwise work:

new_instrument.run_checks(new_file_path)

We see that an error is raised. With our new instrument, however:

instrument_without_binning.run_checks(new_file_path)

We don’t get any errors.

However, this will cause issues if the binning changes between images. In this case, the binning information should be inferred from other keywords (e.g., NAXIS1 and NAXIS2).

For example, if your instrument has a resolution of 2048x2048, then you can infer the binning mode by dividing the instrument’s resolution by the image resolution:

%%writefile out/instrument_with_inferred_binning.py
import phoptic


class InstrumentWithInferredBinning(phoptic.Instrument):

    def get_binning(self, file=None, header=None):

        if header is None:
            header = file.get_header()

        x_size = int(header['NAXIS1'])
        y_size = int(header['NAXIS2'])

        # divide instrument resolution by image resolution
        x_binning = 2048 // x_size
        y_binning = 2048 // y_size

        return f'{x_binning}x{y_binning}'

from out.instrument_with_inferred_binning import InstrumentWithInferredBinning


instrument_with_inferred_binning = InstrumentWithInferredBinning.from_json(config=config_dict)
instrument_with_inferred_binning.run_checks(new_file_path)

The above can be extended to more-or-less any case in which an instrument does not include required keywords in its data products.