Search code examples
pythonpython-3.xastronomyastropy

How to plot a star from TESScut to the right place on the sky using Python (and not with a ~10 degree error)?

I would like to plot TESS measurements of Beta Dor, using data from TESScut. In the Gaia DR2 archive I look up RA and DEC values of Beta Dor, I get:

RA: 83.40630967798376 DEC:  -62.48977125108528

I go to TESScut, put these two values in, select Sector 1, then Download FFI Cutout. Unzip downloaded file, then in Python I do:

import matplotlib.pyplot as plt
from lightkurve import TessTargetPixelFile
import numpy as np

sector1ffi_cutout='tess-s0001-4-4_83.40630967798376_-62.48977125108528_64x64_astrocut.fits'
tpf_s1 = TessTargetPixelFile(sector1ffi_cutout)

fig = plt.figure(figsize=(5,5))
fig.add_subplot(111, projection=tpf_s1.wcs)
plt.pcolormesh(np.log(tpf_s1.flux[0]))
plt.show()

Giving me:

enter image description here

Which, looking at the coordinates (~ -53°, ~6h30min) is clearly wrong, the star is not in this position according to Gaia (and other sources).

What am I doing wrongly, and how can I plot this star to the place it should be?

For completeness, when I plot the star, I get these warnings:

> /home/szabopal/.local/lib/python3.5/site-packages/ipykernel_launcher.py:10:
> RuntimeWarning: divide by zero encountered in log   # Remove the CWD
> from sys.path while we load stuff.
> /home/szabopal/.local/lib/python3.5/site-packages/ipykernel_launcher.py:10:
> RuntimeWarning: invalid value encountered in log   # Remove the CWD
> from sys.path while we load stuff.
> /home/szabopal/.local/lib/python3.5/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73:
> RuntimeWarning: invalid value encountered in greater   discontinuous =
> step[1:] > DISCONT_FACTOR * step[:-1]
> /home/szabopal/.local/lib/python3.5/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73:
> RuntimeWarning: invalid value encountered in greater   discontinuous =
> step[1:] > DISCONT_FACTOR * step[:-1]
> /home/szabopal/.local/lib/python3.5/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73:
> RuntimeWarning: invalid value encountered in greater   discontinuous =
> step[1:] > DISCONT_FACTOR * step[:-1]
> /home/szabopal/.local/lib/python3.5/site-packages/astropy/visualization/wcsaxes/grid_paths.py:73:
> RuntimeWarning: invalid value encountered in greater   discontinuous =
> step[1:] > DISCONT_FACTOR * step[:-1]

Further developments

I believe the issue above is caused by the wcs refering to the whole FFI, not the cutout.

print(tpf_s1.wcs)

gives:

WCS Keywords

Number of WCS axes: 2
CTYPE : 'RA---TAN-SIP'  'DEC--TAN-SIP'  
CRVAL : 90.634460449219  -57.666290283203  
CRPIX : 250.0  -984.0  
PC1_1 PC1_2  : 1.0  1.0  
PC2_1 PC2_2  : 1.0  1.0  
CDELT : 0.00571299832697903  0.005705604460241471  
NAXIS : 81986  1282
Solution

  • This issue stems from a subtlety in the way matplotlib handles WCS projections.

    WCSAxes, which creates the WCS projection for maplotlib plots does not take into account SIP distortions in the WCS (because it uses wcs.wcs_world2pix rather than wcs.all_world2pix).

    Often this doesn’t matter that much, however there are two factors that can make it matter a lot when displaying TESScut cutouts. The first is that the TPFs that TESScut produces come with the original WCS information from the full-frame image adjusted for the position of the cutout. This means that unlike the WCS information in TESS pipeline TPFs, the WCS information in TESScut TPFS includes a full complement of SIP distortions, which are being ignored by matplotlib. The second factor in how accurate the displayed coordinates will be is where on the TESS CCD the source falls. Because the TESS field of view is large, at the edges of a full-frame image the difference made by ignoring the SIP distortions can be on the order of 15 pixels.

    This particular source is very much on the edge of that FFI so the SIP distortions will make quite a large difference:

    TESS CCD with source shown near edge

    There are two options for handling this problem:

    1. Don't use the matplotlib projection option directly, and instead move everything into pixel space, manually calling wcs.all_world2pix or wcs.allpix2world as required.

    2. Create a new WCS for the cutout that does not include SIP distortions (which will be fine because the cutout is small enough to not really need them). Here is a Jupyter notebook that describes one way to make this new WCS for the cutout: https://github.com/ceb8/tessworkshop_wcs_hack/blob/master/tesscut_wcs_hack.ipynb.