Initial imports and data download if you need the KPF input files
[ ]:
# Jupyter auto-reload source code
%load_ext autoreload
%autoreload 2
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from astropy.io import fits
from rvdata.core.models.level2 import RV2
#To download the input files for this example, set download=1 to run the code block below
download=1
if download == 1:
import requests
def download_file(url, filename):
response = requests.get(url)
response.raise_for_status() # Check if the request was successful
with open(filename, "wb") as file:
file.write(response.content)
#Download the native, input FITS files for this example if not already on your computer
file_urls = {
"KPF": [
"http://grinnell.as.arizona.edu/~rvdata/kpf/KP.20250208.17485.59.fits",
"http://grinnell.as.arizona.edu/~rvdata/kpf/KP.20250208.17485.59_L1.fits",
]
}
for file in file_urls["KPF"]:
filename = file.split("/")[-1].split("?")[0]
if not os.path.exists(filename):
print(f"Downloading {filename}...")
download_file(file, filename)
else:
print(f"{filename} already exists. Skipping download.")
KP.20250208.17485.59.fits already exists. Skipping download.
KP.20250208.17485.59_L1.fits already exists. Skipping download.
Using the L2 translator to go from native KPF data format to the EPRV Standard Format
[ ]:
#Translating a KPF File [change file destination to where you have placed the L0 and L1 KPF files from Zenodo]
l1file = 'KP.20250208.17485.59_L1.fits'
l0file = 'KP.20250208.17485.59.fits'
with fits.open(l1file) as hdul:
datetime = hdul[0].header['DATE-BEG']
timestamp, seconds = datetime.split('.')
filetime = timestamp.replace(":","").replace("-","")
kpf_l2_filename = './docs/tutorials/KPFL2_'+filetime.replace("-:","")+'.fits'
kpf_l2 = RV2.from_fits(l1file, l0file=l0file, instrument="KPF")
kpf_l2.to_fits(kpf_l2_filename)
Examining What’s Inside the Community Standard L2 File
Once you have created the community standard L2 fits file, you can read it into python either as a traditional FITS object using astropy’s fits.open() function, or as a python object using the from_fits() function within RV2:
from astropy.io import fits l2 = fits.open(kpf_l2_filename)
OR
from rvdata.core.models.level2 import RV2 l2_obj = RV2.from_fits(l2_standard)
[ ]:
#In this example we'll use the astropy fits.open() approach
l2 = fits.open(kpf_l2_filename)
#The primary header contains useful info about what was observed, what's in each trace, etc. The keywords will be the same regardless
#of what facility was used to take the data, so you can always query the file in the same way.
primary_header = l2[0].header
print('What telescope produced this observation? ',primary_header['TELESCOP'])
print('What instrument was used? ',primary_header['INSTRUME'])
print('What object was observed? ',primary_header['OBJECT'],)
print('How many traces are produced by the spectrograph? ',primary_header['NUMTRACE'] )
print('What is in each of the traces? ')
for i in range(1,primary_header['NUMTRACE']+1):
print('Trace',i,primary_header['TRACE'+str(i)])
The Extension Description table, located in the “EXT_DESCRIPT” extension, lists all of the FITS extensions in the L2 file
[ ]:
EXT_DESCRIPT = pd.DataFrame(l2['EXT_DESCRIPT'].data)
EXT_DESCRIPT
And the Order Table, located in the “ORDER_TABLE” extension, details the start and end wavelengths of each echelle order
[ ]:
ORDER_TABLE = pd.DataFrame(l2['ORDER_TABLE'].data)
ORDER_TABLE
Plotting Examples of L2 Spectra
From the Extension Description table, we know that the stellar data is contained in Traces 2, 3, and 4 – so we can take a look at those
[ ]:
#specify what order you'd like to plot
order = 20
#standardize the y-axis range based on the order with this highest flux measurements so it's easier to compare
flux_max = np.max([np.max(l2['TRACE2_FLUX'].data[order]), np.max(l2['TRACE3_FLUX'].data[order]), np.max(l2['TRACE4_FLUX'].data[order])])
ymax = flux_max*1.05
#Calculate scaling factors between the blaze extensions and the flux extensions for easier visual comparison
blz_scale2 = np.nanmax(l2['TRACE2_FLUX'].data[order]/np.nanmax(l2['TRACE2_BLAZE'].data[order]))
blz_scale3 = np.nanmax(l2['TRACE3_FLUX'].data[order]/np.nanmax(l2['TRACE3_BLAZE'].data[order]))
blz_scale4 = np.nanmax(l2['TRACE4_FLUX'].data[order]/np.nanmax(l2['TRACE4_BLAZE'].data[order]))
fig = plt.figure(1,[10,8])
ax1 = fig.add_subplot(3,1,1)
ax1.scatter(l2['TRACE2_WAVE'].data[order],l2['TRACE2_FLUX'].data[order],marker='.',s=5)
ax1.plot(l2['TRACE2_WAVE'].data[order],l2['TRACE2_BLAZE'].data[order]*blz_scale2,color='orange')
ax1.set_ylim(0,ymax)
ax1.set_ylabel('Counts',fontsize=15,fontweight='bold')
ax2 = fig.add_subplot(3,1,2)
ax2.scatter(l2['TRACE3_WAVE'].data[order],l2['TRACE3_FLUX'].data[order],marker='.',s=5)
ax2.plot(l2['TRACE3_WAVE'].data[order],l2['TRACE3_BLAZE'].data[order]*blz_scale3,color='orange')
ax2.set_ylim(0,ymax)
ax2.set_ylabel('Counts',fontsize=15,fontweight='bold')
ax3 = fig.add_subplot(3,1,3)
ax3.scatter(l2['TRACE4_WAVE'].data[order],l2['TRACE4_FLUX'].data[order],marker='.',s=5)
ax3.plot(l2['TRACE4_WAVE'].data[order],l2['TRACE4_BLAZE'].data[order]*blz_scale4,color='orange')
ax3.set_ylim(0,ymax)
ax3.set_ylabel('Counts',fontsize=15,fontweight='bold')
ax3.set_xlabel('Wavelength (A)',fontsize=15,fontweight='bold')
filename_pieces = kpf_l2_filename.split('/')
fig.suptitle('Science Traces for Order '+str(order)+' of file '+filename_pieces[-1],fontsize=15,fontweight='bold')
plt.show()
Read in data from multiple instruments on the same star
Generate an EPRV Standard L2 file for NEID & KPF
[ ]:
from rvdata.instruments.neid.level2 import RV2, NEIDRV2
neid_L2_Filename = 'NEIDL2_20240921T052641.fits'
neid_L2_obj = NEIDRV2.from_fits('neidL2_20240921T052641.fits', instrument="NEID")
neid_L2_obj.to_fits(neid_L2_Filename)
kpf_l1file = 'KP.20240920.37616.85_L1.fits'
kpf_l0file = 'KP.20240920.37616.85.fits'
kpf_l2_filename = 'KPFL2_20240920T102656.fits'
kpf_l2 = RV2.from_fits(l1file, l0file=l0file, instrument="KPF")
kpf_l2.to_fits(kpf_l2_filename)
Plot an overlapping order, first with the blaze function and then after dividing it out
[ ]:
kpf_obs = fits.open('./docs/tutorials/KPFL2_20240920T102656.fits')
neid_obs = fits.open('./docs/tutorials/NEIDL2_20240921T052641.fits')
kpf_order = 19
neid_order1 = 55
fig = plt.figure(1,[12,6])
ax1 = fig.add_subplot(1,1,1)
ax1.plot(kpf_obs['TRACE2_WAVE'].data[kpf_order],kpf_obs['TRACE2_FLUX'].data[kpf_order],color='blue')
ax1.plot(neid_obs['TRACE1_WAVE'].data[neid_order1],neid_obs['TRACE1_FLUX'].data[neid_order1],color='orange')
ax1.set_ylabel('Counts',fontsize=15,fontweight='bold')
ax1.set_xlabel('Wavelength (A)',fontsize=15,fontweight='bold')
ax1.text(5145,90000,'KPF obs of HD 4268',color='blue',fontsize=20,fontweight='bold')
ax1.text(5145,80000,'NEID obs of HD 4268',color='darkorange',fontsize=20,fontweight='bold')
plt.show()
[ ]:
kpf_obs = fits.open('./docs/tutorials/KPFL2_20240920T102656.fits')
neid_obs = fits.open('./docs/tutorials/NEIDL2_20240921T052641.fits')
kpf_order = 19
neid_order = 55
#Calculate scaling factors between the blaze extensions and the flux extensions for easier visual comparison
blz_scale_kpf = np.nanmax(kpf_obs['TRACE2_FLUX'].data[kpf_order]/np.nanmax(kpf_obs['TRACE2_BLAZE'].data[kpf_order]))
blz_scale_neid = np.nanmax(neid_obs['TRACE1_FLUX'].data[neid_order]/np.nanmax(neid_obs['TRACE1_BLAZE'].data[neid_order]))
fig = plt.figure(1,[12,6])
ax1 = fig.add_subplot(1,1,1)
ax1.plot(neid_obs['TRACE1_WAVE'].data[neid_order],neid_obs['TRACE1_FLUX'].data[neid_order] /
(neid_obs['TRACE1_BLAZE'].data[neid_order] * blz_scale_neid) ,color='orange')
ax1.plot(kpf_obs['TRACE2_WAVE'].data[kpf_order],kpf_obs['TRACE2_FLUX'].data[kpf_order] /
(kpf_obs['TRACE2_BLAZE'].data[kpf_order] * blz_scale_kpf) ,color='blue',alpha=.6)
ax1.set_ylabel('Counts',fontsize=15,fontweight='bold')
ax1.set_xlabel('Wavelength (A)',fontsize=15,fontweight='bold')
ax1.set_ylim(0,1.3)
ax1.text(5145,1.2,'KPF obs of HD 4268 - deblazed',color='blue',fontsize=16,fontweight='bold')
ax1.text(5145,1.12,'NEID obs of HD 4268 - deblazed',color='darkorange',fontsize=16,fontweight='bold')
plt.show()