Optimizing Pointing

In this notebook, we will explore the jwpoint functionality that can help us find the optimal pointing when planning an upcoming JWST observations.

The goal here, for a given instrument, detector and subaperture, is to select the optimal pointing position that minimizes the amount of bad pixels and respects a series of other criteria.

This is mostly done via do_region_search(). We will explain its arguments in the following sections.

Downloading a reference observation

Before we start optimizing the pointing, we need to download a reference observation. We need an observation with the same instrument, detector, subaperture and filter as our planned observation. Ideally, there should be at least once easy to find and isolated point source in the observation so we can use it as a reference PSF.

Here we use an data from the 1205 program which observed high-z quasars with NIRCam We use an observation from the long-wavelength (LW) channel in the F430M filter and the SUB400P subarray.

from pathlib import Path
from astroquery.mast import Observations

data_dir = Path("data")
program_dir = "01205"

filename = Path("jw01205002001_02104_00002_nrcblong_cal.fits")
uri = f"mast:JWST/product/{filename}"

filepath = data_dir / program_dir / filename
filepath.parent.mkdir(exist_ok=True, parents=True)
_ = Observations.download_file(uri, local_path=filepath)

INFO: Found cached file data/01205/jw01205002001_02104_00002_nrcblong_cal.fits with expected size 4613760. [astroquery.query]

import numpy as np
from astropy.io import fits

with fits.open(filepath) as hdul:
    hdr = hdul[0].header
    img = hdul[1].data
dq_mask = np.isnan(img)

from matplotlib import rcParams
import matplotlib.pyplot as plt

rcParams["image.origin"] = "lower"

plt.imshow(img, norm="symlog")
plt.xlabel("X [pixel]")
plt.ylabel("Y [pixel]")
plt.show()

It is pretty clear which is our point source target here. Let us extract it so it can serve as a reference PSF.

Extracting the reference PSF

First, we can use get_pointing_position() to get the position of our target.

from jwpoint import pointing
x_target, y_target = (int(round(off)) for off in pointing.get_pointing_position(hdr["XOFFSET"], hdr["YOFFSET"], filepath))

from matplotlib import rcParams
import matplotlib.pyplot as plt

rcParams["image.origin"] = "lower"

plt.imshow(img, norm="symlog")
plt.plot(x_target, y_target, "r*")
plt.xlabel("X [pixel]")
plt.ylabel("Y [pixel]")
plt.show()

from jwpoint.plot import plot_dithers, zoom_plot

zoom_plot(img, x_target, y_target, size=64, show_mask=False)
plt.xlabel("X [pixel]")
plt.ylabel("Y [pixel]")
plt.show()

This looks good! We can crop the image. The PSF must have the same shape as the region size we will be optimizing later. Here we set this to 70 pixels, which is a reasonable value if what we care about is the close surroundings of the source.

crop_size = 70
hs = crop_size // 2
psf_with_badpix = img[y_target-hs:y_target+hs, x_target-hs:x_target+hs]

plt.imshow(psf_with_badpix, norm="symlog")
plt.xlabel("X [pixel]")
plt.ylabel("Y [pixel]")
plt.show()

This cropped PSF has bad pixels set to NaN. To avoid propagating NaNs throughout our calculations, we can replace them with a median filter.

import jwpoint.utils as ut

psf = ut.filter_nans(psf_with_badpix)

plt.imshow(psf, norm="symlog")
plt.xlabel("X [pixel]")
plt.ylabel("Y [pixel]")
plt.show()

Optimal region search

The pointing.do_region_search() function uses the DQ (bad pixel) mask from the reference image to find regions in the data with as little bad pixels as possible. It does a naive grid search over all possible regions in the image.

By default, for each position, it defines a region of size crop_size and calculates the number of bad pixels. The n_top best regions are then returned. Let us test this for a single region.

from jwpoint.search import do_region_search

x_naive, y_naive = do_region_search(
    dq_mask, img, crop_size, psf=psf, n_top=1
)

Let us break down what is shown here:

On the left of the first plot, we se the bad pixel mask for the entire subaperture. In the middle, we see the bad pixel which is color-coded by the number of bad-pixels for a 70x70 region centered on this pixels. In short, the darker a point is, the better. On the right, we see the science image for reference. On each panel the “1” indicates the best pointing position

In the second plot, we see a zoom on this optimal region. The top panel shows the actual region in the image while the bottom shows the reference PSF with the bad pixel mask from this region. The latter is useful to get a sense of where bad pixels would fall on a point source.

Searching for multiple regions

If we wanted to get a few options, we could simply increase n_top.

x_few, y_few = do_region_search(
    dq_mask, img, crop_size, psf=psf, n_top=5
)
x_few, y_few

(array([ 81, 244, 310, 164, 236]), array([207, 260, 104, 249, 127]))

The x and Y positions are each returned in a NumPy array.

Weighting bad pixels with a PSF

The searches above simply summed the number of bad pixels everywhere in the image uniformly. While this can work, it may be desirable to penalize bad pixels near the core of the PSF more than on the edges. This can be done with the kernel="weighted" argument:

x_weighted, y_weighted = do_region_search(
    dq_mask, img, crop_size, psf=psf, n_top=1, kernel="weighted"
)

As we can see, in this case, it barely changes the result, but it does change how pixels are weighted (we see PPSFs instead of squares).

Rejecting bad pixels at the core of the PSF

Sometimes, the weighting from the previous section may not be enough. There are cases where any region with bad pixels in the core of the PSF should be rejected. This can be done with the forbidden_size argument, which defines the size of the PSF core where no bad pixels are accepted.

x_forbid, y_forbid = do_region_search(
    dq_mask, img, crop_size, psf=psf, n_top=1, forbidden_size=25,
)

As we can see, this changes the results quite a bit: many pixels are now rejected (NaN in the middle panel) and the pointing goes from middle left to bottom right.

However, be careful not to push this option too far: you will reach a regime where no region satisfies the requirement.

x_forbid, y_forbid = do_region_search(
    dq_mask,
    img,
    crop_size,
    psf=psf,
    n_top=1,
    forbidden_size=40,
)

WARNING: No optimal region was was found. Try relaxing the constraints.

Requiring a minimal edge distance

If we need our target to be centered in the subarray, or if we don’t want to take too many risks by putting it near the edge, we can request a minimum edge distance for the optimal region(s). By default, the search window (region size of crop_size = 70 here) is used as the minimum edge distance.

x_edge, y_edge = do_region_search(
    dq_mask,
    img,
    crop_size,
    psf=psf,
    n_top=1,
    min_edge_distance=150,
)

Accounting for the subarray in the short-wavelength channel

In all the searches above, we considered only the long-wavelength (LW) channel with the F430M filter. However, NIRCam also observed in the short-wavelength (SW) channel simultaneously.

Naive LW transformed to SW

The simplest thing we can do is just take our optimal pointing for the LW channel and see where it falls in SW. Let us first extract the LW filter and subarray for later use.

filt = hdr["FILTER"]
subarray = hdr["SUBARRAY"]
print(f"Filter: {filt}")
print(f"Subarray: {subarray}")

Filter: F430M
Subarray: SUB400P

We then infer the SW detector from the initial “naive” optimal position and download the file, all with the download_sw_file() function which is discussed in more detail in the Understanding JWST Pointing notebook.

filepath_sw = ut.download_sw_file(filepath, x_naive, y_naive)

INFO: Found cached file data/01205/jw01205002001_02104_00002_nrcb1_cal.fits with expected size 4610880. [astroquery.query]

with fits.open(filepath_sw) as hdul_sw:
    hdr_sw = hdul_sw[0].header
    img_sw = hdul_sw[1].data

filt_sw = hdr_sw["FILTER"]

x_naive_sw, y_naive_sw = pointing.long_to_short(x_naive, y_naive, filepath, filepath_sw)

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
ax_lw, ax_sw = axs
ax_lw.imshow(img, norm="symlog")
ax_lw.plot(x_naive, y_naive, "r*", label="naive position")
ax_lw.set_xlabel("X [pixel]")
ax_lw.set_ylabel("Y [pixel]")
ax_lw.set_title(f"LW image in filter {filt}")

ax_sw.imshow(img_sw, norm="symlog")
ax_sw.plot(x_naive_sw, y_naive_sw, "r*")
ax_sw.set_xlabel("X [pixel]")
ax_sw.set_ylabel("Y [pixel]")
ax_sw.set_title(f"SW Image in filter {filt_sw}")
fig.legend()
plt.show()

Oups! As we can see above, the pointing falls outside the subarray in the SW filter. This can happen because some subarrays (such as SUB400P used here) have different FOVs in the LW and SW channels. If we don’t care about the SW data, this is not a big deal, but if we do we must make sure our optimal position is usable in both channels.

Using the subarray option

Currently, jwpoint does not support joint optimization in the LW and SW filters since most of the science cases we used it for had a strong emphasis on LW. However, it should not be too hard to add, so feel free to open an issue or a pull request on GitHub!

For now, what we do have is a subarray argument that will restrict the search to regions that are visible in both channels for a given subarray, but the optimization takes place entirely in the provided LW image.

x_naive_sub, y_naive_sub = do_region_search(
    dq_mask, img, crop_size, psf=psf, n_top=1, subarray=subarray,
)

Again, we can propagate this to the SW filter and see where it falls on both detectors.

x_naive_sub_sw, y_naive_sub_sw = pointing.long_to_short(x_naive_sub, y_naive_sub, filepath, filepath_sw)

fig, axs = plt.subplots(1, 2, figsize=(10, 5))
ax_lw, ax_sw = axs
ax_lw.imshow(img, norm="symlog")
ax_lw.plot(x_naive_sub, y_naive_sub, "r*", label="Subarray-constrained position")
ax_lw.set_xlabel("X [pixel]")
ax_lw.set_ylabel("Y [pixel]")
ax_lw.set_title(f"LW image in filter {filt}")

ax_sw.imshow(img_sw, norm="symlog")
ax_sw.plot(x_naive_sub_sw, y_naive_sub_sw, "r*")
ax_sw.set_xlabel("X [pixel]")
ax_sw.set_ylabel("Y [pixel]")
ax_sw.set_title(f"SW Image in filter {filt_sw}")

fig.legend()
plt.show()

Looks like the contraint worked and the optimal position found is now within both detectors.

Optimizing pointing for multiple dithers

In Understanding Pointing, we discuss how to calculate the pointing for multiple dithers. Similarly, we can optimize the pointing by jointly optimizing the number of bad pixels across dithers instead of on a region.

This is done by specifying the joint_offsets argument. It needs to be a list of tuples, each with the x and y coordinates of the dithers, or a mapping to x and y offsets in pixels.

Let us say we wanted to apply the INTRAMODULEBOX pattern. We first need to fetch the dither offsets. By passing a detector argument, the offsets are converted from arcsec to pixels.

from jwpoint.dithers import get_dither_info

dither_pattern = "INTRAMODULEBOX"
n_dithers = 5
detector = hdr["DETECTOR"]

Then we can pass the offsets so that they are taken into account for the region search. Note that the plots currently only show the reference pointing position and do not show the dithers, the dithers are however considered when counting bad pixels.

x_dithers, y_dithers = do_region_search(
    dq_mask,
    img,
    crop_size,
    psf=psf,
    n_top=1,
    dither_pattern=dither_pattern,
    n_dithers=n_dithers,
    detector=detector
)

Now, instead of being 1D with shape (n_top,), the arrays will have shape (n_top, n_dithers).

Finding multiple pointings for multiple dithers

We can also do the multi-dither optimization while also optimizing for the n_top best regions while accounting for dithers

x_dithers_multi, y_dithers_multi = do_region_search(
    dq_mask,
    img,
    crop_size,
    psf=psf,
    n_top=2,
    dither_pattern=dither_pattern,
    n_dithers=n_dithers,
    detector=detector
)