FFIMachine Tutorial 2.0¶
import numpy as np
import pandas as pd
import lightkurve as lk
import psfmachine as pm
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.table import Table
from scipy import sparse
import warnings
warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=sparse.SparseEfficiencyWarning)
Kepler FFI¶
FFIMachine is a submodule of PSFMachine that allows users work with Kepler and K2 Full Frame Images (FFIs), TESS support will be added later.
It is meant to be used to build robust, detailed, and smooth PRF models of Kepler at different quarters and channels.
These PRF models can be saved and later used by TPFMachine to do PSF photometry on Kepler's TPFs instead of computing a
"sparse" PRF model from the TPF pixel data.
It also includes the option to perform PSF photometry of sources in the data and save them to a source catalog
in FITS format.
First, let's load a Kepler FFI image. This FITS file can be download from MAST archive: bulk downloads, Full Frame Images, Calibrated.
file_path = "/Volumes/ADAP-Kepler/Work/BAERI/data"
fname = f"{file_path}/kepler/ffi/kplr2009231194831_ffi-cal.fits"
FFI files are large (~400 MB) and contain 84 extension, one per channel.
hdu = fits.open(fname)
hdu.info()
Filename: /Volumes/ADAP-Kepler/Work/BAERI/data/kepler/ffi/kplr2009231194831_ffi-cal.fits No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 58 () 1 MOD.OUT 2.1 1 ImageHDU 100 (1132, 1070) float32 2 MOD.OUT 2.2 1 ImageHDU 100 (1132, 1070) float32 3 MOD.OUT 2.3 1 ImageHDU 100 (1132, 1070) float32 4 MOD.OUT 2.4 1 ImageHDU 100 (1132, 1070) float32 5 MOD.OUT 3.1 1 ImageHDU 100 (1132, 1070) float32 6 MOD.OUT 3.2 1 ImageHDU 100 (1132, 1070) float32 7 MOD.OUT 3.3 1 ImageHDU 100 (1132, 1070) float32 8 MOD.OUT 3.4 1 ImageHDU 100 (1132, 1070) float32 9 MOD.OUT 4.1 1 ImageHDU 100 (1132, 1070) float32 10 MOD.OUT 4.2 1 ImageHDU 100 (1132, 1070) float32 11 MOD.OUT 4.3 1 ImageHDU 100 (1132, 1070) float32 12 MOD.OUT 4.4 1 ImageHDU 100 (1132, 1070) float32 13 MOD.OUT 6.1 1 ImageHDU 100 (1132, 1070) float32 14 MOD.OUT 6.2 1 ImageHDU 100 (1132, 1070) float32 15 MOD.OUT 6.3 1 ImageHDU 100 (1132, 1070) float32 16 MOD.OUT 6.4 1 ImageHDU 100 (1132, 1070) float32 17 MOD.OUT 7.1 1 ImageHDU 100 (1132, 1070) float32 18 MOD.OUT 7.2 1 ImageHDU 100 (1132, 1070) float32 19 MOD.OUT 7.3 1 ImageHDU 100 (1132, 1070) float32 20 MOD.OUT 7.4 1 ImageHDU 100 (1132, 1070) float32 21 MOD.OUT 8.1 1 ImageHDU 100 (1132, 1070) float32 22 MOD.OUT 8.2 1 ImageHDU 100 (1132, 1070) float32 23 MOD.OUT 8.3 1 ImageHDU 100 (1132, 1070) float32 24 MOD.OUT 8.4 1 ImageHDU 100 (1132, 1070) float32 25 MOD.OUT 9.1 1 ImageHDU 100 (1132, 1070) float32 26 MOD.OUT 9.2 1 ImageHDU 100 (1132, 1070) float32 27 MOD.OUT 9.3 1 ImageHDU 100 (1132, 1070) float32 28 MOD.OUT 9.4 1 ImageHDU 100 (1132, 1070) float32 29 MOD.OUT 10.1 1 ImageHDU 100 (1132, 1070) float32 30 MOD.OUT 10.2 1 ImageHDU 100 (1132, 1070) float32 31 MOD.OUT 10.3 1 ImageHDU 100 (1132, 1070) float32 32 MOD.OUT 10.4 1 ImageHDU 100 (1132, 1070) float32 33 MOD.OUT 11.1 1 ImageHDU 100 (1132, 1070) float32 34 MOD.OUT 11.2 1 ImageHDU 100 (1132, 1070) float32 35 MOD.OUT 11.3 1 ImageHDU 100 (1132, 1070) float32 36 MOD.OUT 11.4 1 ImageHDU 100 (1132, 1070) float32 37 MOD.OUT 12.1 1 ImageHDU 100 (1132, 1070) float32 38 MOD.OUT 12.2 1 ImageHDU 100 (1132, 1070) float32 39 MOD.OUT 12.3 1 ImageHDU 100 (1132, 1070) float32 40 MOD.OUT 12.4 1 ImageHDU 100 (1132, 1070) float32 41 MOD.OUT 13.1 1 ImageHDU 100 (1132, 1070) float32 42 MOD.OUT 13.2 1 ImageHDU 100 (1132, 1070) float32 43 MOD.OUT 13.3 1 ImageHDU 100 (1132, 1070) float32 44 MOD.OUT 13.4 1 ImageHDU 100 (1132, 1070) float32 45 MOD.OUT 14.1 1 ImageHDU 100 (1132, 1070) float32 46 MOD.OUT 14.2 1 ImageHDU 100 (1132, 1070) float32 47 MOD.OUT 14.3 1 ImageHDU 100 (1132, 1070) float32 48 MOD.OUT 14.4 1 ImageHDU 100 (1132, 1070) float32 49 MOD.OUT 15.1 1 ImageHDU 100 (1132, 1070) float32 50 MOD.OUT 15.2 1 ImageHDU 100 (1132, 1070) float32 51 MOD.OUT 15.3 1 ImageHDU 100 (1132, 1070) float32 52 MOD.OUT 15.4 1 ImageHDU 100 (1132, 1070) float32 53 MOD.OUT 16.1 1 ImageHDU 100 (1132, 1070) float32 54 MOD.OUT 16.2 1 ImageHDU 100 (1132, 1070) float32 55 MOD.OUT 16.3 1 ImageHDU 100 (1132, 1070) float32 56 MOD.OUT 16.4 1 ImageHDU 100 (1132, 1070) float32 57 MOD.OUT 17.1 1 ImageHDU 100 (1132, 1070) float32 58 MOD.OUT 17.2 1 ImageHDU 100 (1132, 1070) float32 59 MOD.OUT 17.3 1 ImageHDU 100 (1132, 1070) float32 60 MOD.OUT 17.4 1 ImageHDU 100 (1132, 1070) float32 61 MOD.OUT 18.1 1 ImageHDU 100 (1132, 1070) float32 62 MOD.OUT 18.2 1 ImageHDU 100 (1132, 1070) float32 63 MOD.OUT 18.3 1 ImageHDU 100 (1132, 1070) float32 64 MOD.OUT 18.4 1 ImageHDU 100 (1132, 1070) float32 65 MOD.OUT 19.1 1 ImageHDU 100 (1132, 1070) float32 66 MOD.OUT 19.2 1 ImageHDU 100 (1132, 1070) float32 67 MOD.OUT 19.3 1 ImageHDU 100 (1132, 1070) float32 68 MOD.OUT 19.4 1 ImageHDU 100 (1132, 1070) float32 69 MOD.OUT 20.1 1 ImageHDU 100 (1132, 1070) float32 70 MOD.OUT 20.2 1 ImageHDU 100 (1132, 1070) float32 71 MOD.OUT 20.3 1 ImageHDU 100 (1132, 1070) float32 72 MOD.OUT 20.4 1 ImageHDU 100 (1132, 1070) float32 73 MOD.OUT 22.1 1 ImageHDU 100 (1132, 1070) float32 74 MOD.OUT 22.2 1 ImageHDU 100 (1132, 1070) float32 75 MOD.OUT 22.3 1 ImageHDU 100 (1132, 1070) float32 76 MOD.OUT 22.4 1 ImageHDU 100 (1132, 1070) float32 77 MOD.OUT 23.1 1 ImageHDU 100 (1132, 1070) float32 78 MOD.OUT 23.2 1 ImageHDU 100 (1132, 1070) float32 79 MOD.OUT 23.3 1 ImageHDU 100 (1132, 1070) float32 80 MOD.OUT 23.4 1 ImageHDU 100 (1132, 1070) float32 81 MOD.OUT 24.1 1 ImageHDU 100 (1132, 1070) float32 82 MOD.OUT 24.2 1 ImageHDU 100 (1132, 1070) float32 83 MOD.OUT 24.3 1 ImageHDU 100 (1132, 1070) float32 84 MOD.OUT 24.4 1 ImageHDU 100 (1132, 1070) float32
FFIMachine object¶
Now we create an FFIMachine object using the pm.FFIMachine.from_file() method that inputs the file path (it also accepts a list of FFIs) and the extension (channel) to be used.
FFI images are 1100 x 1024 pixels, users can also extract a cutout of the CCD by providing the cutout_size and cutout_origin
parameters that sets the size in pixels of a square stamp and the origin (in row, column indexing). By default the complete FFI is loaded.
Internally, FFIMachine will query Gaia EDR3 dr=3 wit a limiting magnitude 18 magnitude_limit=18.
Additionally, FFIMachine removes the background flux using photutils.Background2D, it masks saturated pixels and bleed columns,
and attempts to remove pixels where bright sources produce halos from light scatered though the instrument optics.
You'll also see a progress bar that says: Creating delta sparse, this is the creation of sparse matrices that contain
information such as $\delta$RA and $\delta$Dec of pixels around sources, and are cital for the use of PSFMachine.
Printing the ffi will tell us the total number of sources, and usable pixels in the CCD. In this example, there are 22,123 Gaia sources and more than 1.1M pixels.
# this will load the full FFI but will take time
# ffi = pm.FFIMachine.from_file(fname, extension=1)
# we can also do a cutout of the FFI by specifying the cutout size and origin
# this is convinient for faster development or testing
ffi = pm.FFIMachine.from_file(fname, extension=1, cutout_size=400, cutout_origin=[300, 300], magnitude_limit=18)
ffi
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 55062.804927 from DATE-OBS. Set MJD-END to 55062.825361 from DATE-END'. [astropy.wcs.wcs] 2023-06-09 18:41:55,455 - astropy - WARNING - FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 55062.804927 from DATE-OBS. Set MJD-END to 55062.825361 from DATE-END'. 2023-06-09 18:42:00,491 - astroquery - INFO - Query finished.
INFO: Query finished. [astroquery.utils.tap.core]
2023-06-09 18:42:05,104 - astroquery - INFO - Query finished.
INFO: Query finished. [astroquery.utils.tap.core]
2023-06-09 18:42:10,970 - astroquery - INFO - Query finished.
INFO: Query finished. [astroquery.utils.tap.core]
2023-06-09 18:42:15,387 - astroquery - INFO - Query finished.
INFO: Query finished. [astroquery.utils.tap.core]
Creating delta arrays: 100%|██████████████████████████████████████████████████████| 2811/2811 [00:06<00:00, 466.66it/s]
FFIMachine (N sources, N times, N pixels): (2811, 1, 160000)
There is an attribuite with metadata information gathered from the FITS header
ffi.meta
{'TELESCOP': 'Kepler',
'INSTRUME': 'Kepler Photometer',
'MISSION': 'Kepler',
'DATSETNM': 'kplr2009231194831',
'RADESYS': 'ICRS',
'EQUINOX': 2000.0,
'BACKAPP': True,
'EXTENSION': 1,
'CHANNEL': 1,
'CAMERA': 1,
'QUARTER': 2,
'CAMPAIGN': 2,
'DCT_TYPE': 'FFI'}
We can visualize the FFI using ffi.plot_image(), sources can be overplotted, but for a large FOV with so many sources the figure gets too busy.
ffi.plot_image(sources=True)
plt.show()
The mask with rejected pixels can be visualized with ffi.plot_pixel_mask(). Here, yellow patches are pixels masked bright stars that produce halos,
while red markers are saturated pixels.
_ = ffi.plot_pixel_masks()
PSF model¶
Let's now fit a PSF model using all the sources in the scene. This is done using ffi.build_shape_model() and can show a diagnostic plot
with the pixel data and model in both Cartesian and Polar coordinates
ffi.build_shape_model(plot=True)
plt.show()
WARNING: Input data contains invalid values (NaNs or infs), which were automatically clipped. [astropy.stats.sigma_clipping] 2023-06-09 18:42:26,047 - astroquery - WARNING - Input data contains invalid values (NaNs or infs), which were automatically clipped.
Such a funky looking PSF, this is because channel 1 is at one of the edges of Kepler's focal plane where distortion is large.
This awesome PSF model can be saved to disk (in FITS format) to be used later by PSFMachine to do PSF photometry
on Kepler's TPFs instead of using the low number of sources and pixels that the TPFs have.
We save the model with ffi.save_shape_model(), by default it is saved in the working directory with a name describing the data we use. It also accepts custom names with the output parameter.
ffi.save_shape_model()
We can check the new file
!ls -lh *.fits
-rw-r--r-- 1 jorgemarpa staff 8.4K Jun 9 18:42 Kepler_ffi_shape_model_ext1_q2.fits
!fitsheader Kepler_ffi_shape_model_ext1_q2.fits
# HDU 0 in Kepler_ffi_shape_model_ext1_q2.fits: SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T CHECKSUM= 'B4eFD4Z9B4dCB4Z9' / HDU checksum updated 2023-06-09T18:42:31 DATASUM = '0 ' / data unit checksum updated 2023-06-09T18:42:31 # HDU 1 in Kepler_ffi_shape_model_ext1_q2.fits: XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 8 / length of dimension 1 NAXIS2 = 169 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 1 / number of table fields TTYPE1 = 'psf_w ' TFORM1 = 'D ' OBJECT = 'PRF shape' / PRF shape parameters DATATYPE= 'FFI ' / Type of data used to fit shape model ORIGIN = 'PSFmachine.FFIMachine' / Software of origin VERSION = '1.1.3 ' / Software version TELESCOP= 'Kepler ' / Telescope name MISSION = 'Kepler ' / Mission name QUARTER = 2 / Quarter/Campaign/Sector of observations CHANNEL = 1 / Channel/Camera-CCD output MJD-OBS = 2455063.31514402 / MJD of observation N_RKNOTS= 10 / Number of knots for spline basis in radial axis N_PKNOTS= 15 / Number of knots for spline basis in angle axis RMIN = 1 / Minimum value for knot spacing RMAX = 16 / Maximum value for knot spacing CUT_R = 6 / Radial distance to remove angle dependency SPLN_DEG= 3 / Degree of the spline basis NORM = 'False ' / Normalized model CHECKSUM= 'fG2AiD28fD2AfD25' / HDU checksum updated 2023-06-09T18:42:31 DATASUM = '451725647' / data unit checksum updated 2023-06-09T18:42:31
Let's check that the model was saved correctly by loading it with ffi.load_shape_model() we can plot the diagnostic figure to double check it looks OK.
# this might take some time to run
ffi.load_shape_model(input="./Kepler_ffi_shape_model_ext1_q2.fits", plot=True)
plt.show()
Awesome! we recover our nice PSF model.
Do PSF photometry¶
Also we can get PSF photometry of every source with ffi.fit_model(fit_va=False), here we set the fit_va parameter to False,
because we are working with just one cadence, even if loading multiple FFI epochs, fitting a time model like in PSFMachine doesn't make
sense due to the sparse observations. The photometry is stored in the ffi.ws and ffi.werrs attributes.
The source catalog can be saved to a FITS file with ffi.save_flux_values() that has the output parameter for custom file names,
by default the file is saved in the working directory.
This function internally runs ffi.fit_model(fit_va=False), so no need to do both here.
# this will take lot of time if working with a full channel FFI
ffi.save_flux_values()
Fitting 2811 Sources (w. VA): 100%|██████████████████████████████████████████████████████| 1/1 [00:00<00:00, 1.14it/s]
K2 FFI¶
FFIMachine also supports K2 FFIs.
Let's see how this work!
This time lets take a cutout of channel 43.
file_path = "/Volumes/ADAP-Kepler/Work/BAERI/data"
fname = f"{file_path}/k2/ffi/ktwo2014070234206-c00_ffi-cal.fits"
k2ffi = pm.FFIMachine.from_file(fname, extension=43, cutout_size=300, cutout_origin=[600, 600])
k2ffi
WARNING: FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56727.967142 from DATE-OBS. Set MJD-END to 56727.987575 from DATE-END'. [astropy.wcs.wcs] 2023-06-09 19:01:17,461 - astroquery - WARNING - FITSFixedWarning: 'datfix' made the change 'Set MJD-OBS to 56727.967142 from DATE-OBS. Set MJD-END to 56727.987575 from DATE-END'. Creating delta arrays: 100%|██████████████████████████████████████████████████████| 1339/1339 [00:01<00:00, 692.60it/s]
FFIMachine (N sources, N times, N pixels): (1339, 1, 90000)
This time the object initialization was much faster, we are just loading 300 x 300 pixels.
We can check the metadata
k2ffi.meta
{'TELESCOP': 'Kepler',
'INSTRUME': 'Kepler Photometer',
'MISSION': 'K2',
'DATSETNM': 'ktwo2014070234206-c00',
'RADESYS': 'ICRS',
'EQUINOX': 2000.0,
'BACKAPP': True,
'EXTENSION': 43,
'CHANNEL': 43,
'CAMERA': 43,
'QUARTER': 0,
'CAMPAIGN': 0,
'DCT_TYPE': 'FFI'}
A visualization of the image, this time marking the sources with red circles
_ = k2ffi.plot_image(sources=True)
We build a shape model again
_ = k2ffi.build_shape_model(plot=True)
WARNING: Input data contains invalid values (NaNs or infs), which were automatically clipped. [astropy.stats.sigma_clipping] 2023-06-09 19:01:21,965 - astroquery - WARNING - Input data contains invalid values (NaNs or infs), which were automatically clipped.
This PSF model is more rounded because channel 43 is in the center of Kepler's focal plane.
Lets save the photometry to disk and see how the FITS file looks like:
k2ffi.save_flux_values()
Fitting 1339 Sources (w. VA): 100%|██████████████████████████████████████████████████████| 1/1 [00:00<00:00, 6.49it/s]
!ls -lh *.fits
-rw-r--r-- 1 jorgemarpa staff 87K Jun 9 19:01 K2_source_catalog_ext43_q0_mjd2456728.47735842.fits -rw-r--r-- 1 jorgemarpa staff 8.4K Jun 9 18:42 Kepler_ffi_shape_model_ext1_q2.fits -rw-r--r-- 1 jorgemarpa staff 174K Jun 9 18:42 Kepler_source_catalog_ext1_q2_mjd2455063.31514402.fits
!fitsheader K2_source_catalog_ext43_q0_mjd2456728.47735842.fits
# HDU 0 in K2_source_catalog_ext43_q0_mjd2456728.47735842.fits: SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T OBJECT = 'Photometric Catalog' / Photometry ORIGIN = 'PSFmachine.FFIMachine' / Software of origin VERSION = '1.1.3 ' / Software version TELESCOP= 'Kepler ' / Telescope MISSION = 'K2 ' / Mission name DCT_TYPE= 'FFI ' / Data type QUARTER = 0 / Quarter/Campaign/Sector of observations CHANNEL = 43 / Channel/Camera-CCD output APERTURE= 'PSF ' / Type of photometry N_OBS = 1 / Number of cadences DATSETNM= 'ktwo2014070234206-c00' / data set name RADESYS = 'ICRS ' / reference frame of celestial coordinates EQUINOX = 2000.0 / equinox of celestial coordinate system CHECKSUM= 'oY6SqW6RoW6RoW6R' / HDU checksum updated 2023-06-09T19:01:29 DATASUM = '0 ' / data unit checksum updated 2023-06-09T19:01:29 # HDU 1 in K2_source_catalog_ext43_q0_mjd2456728.47735842.fits: XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 61 / length of dimension 1 NAXIS2 = 1339 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 5 / number of table fields TTYPE1 = 'gaia_id ' TFORM1 = '29A ' TTYPE2 = 'ra ' TFORM2 = 'D ' TUNIT2 = 'deg ' TTYPE3 = 'dec ' TFORM3 = 'D ' TUNIT3 = 'deg ' TTYPE4 = 'psf_flux' TFORM4 = 'D ' TUNIT4 = '-e/s ' TTYPE5 = 'psf_flux_err' TFORM5 = 'D ' TUNIT5 = '-e/s ' EXTNAME = 'CATALOG ' MJD-OBS = 2456728.47735842 / MJD of observation CHECKSUM= '7XAp7U4o7UAo7U3o' / HDU checksum updated 2023-06-09T19:01:29 DATASUM = '1750058098' / data unit checksum updated 2023-06-09T19:01:29
table = Table.read("K2_source_catalog_ext43_q0_mjd2456728.47735842.fits")
table
WARNING: UnitsWarning: '-e/s' did not parse as fits unit: Invalid character at col 0 If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html [astropy.units.core] 2023-06-09 19:01:32,597 - astroquery - WARNING - UnitsWarning: '-e/s' did not parse as fits unit: Invalid character at col 0 If this is meant to be a custom unit, define it with 'u.def_unit'. To have it recognized inside a file reader or other code, enable it with 'u.add_enabled_units'. For details, see https://docs.astropy.org/en/latest/units/combining_and_defining.html
| gaia_id | ra | dec | psf_flux | psf_flux_err |
|---|---|---|---|---|
| deg | deg | -e/s | -e/s | |
| bytes29 | float64 | float64 | float64 | float64 |
| Gaia DR3 3373071429411856000 | 98.65870854646455 | 21.136655773508906 | 1111.6548987961082 | 50.67865893353805 |
| ... | ... | ... | ... | ... |
| Gaia DR3 3373119498684323456 | 98.84477190034355 | 21.492906885836497 | 3634.1502069045546 | 71.01945021888659 |