import os
import re
import subprocess
import astropy.io.fits as fits
from astropy import wcs
from astropy.modeling import models, fitting
import numpy as np
import scipy.ndimage as ndimage
import scipy.stats
from scipy.interpolate import interp1d
import time
from copy import copy
import multiprocessing as mp
#different imports depending on if python2.7 or python3
import sys
if sys.version_info < (3,0):
#python 2.7 behavior
import ConfigParser
else:
import configparser as ConfigParser
import pyklip
from pyklip.instruments.Instrument import Data
from pyklip.instruments.utils.nair import nMathar
from pyklip.parallelized import high_pass_filter_imgs
from pyklip.fakes import gaussfit2d
from pyklip.fakes import gaussfit2dLSQ
[docs]
class NIRC2Data(Data):
"""
A sequence of Keck NIRC2 ADI Data. Each NIRC2Data object has the following fields and functions
Args:
filepaths: list of filepaths to files
highpass: if True, run a Gaussian high pass filter (default size is sigma=imgsize/10)
can also be a number specifying FWHM of box in pixel units
find_star: (default) 'auto' will first try to get the star center coordinates from the FITS
header PSFCENTX & PSFCENTY keywords, and if that fails it will do a Radon transform to
locate the star via the diffraction spikes (and store the star center in the header for
future use). True will force the Radon transform; False will skip the Radon transform
even if no center is found in the header.
guess_star: (default) None. Otherwise an [x, y] coordinate guess for where the star is for find_star.
The code will automatially guess the star location based on the FPM otherwise.
Attributes:
input: Array of shape (N,y,x) for N images of shape (y,x)
centers: Array of shape (N,2) for N centers in the format [x_cent, y_cent]
filenums: Array of size N for the numerical index to map data to file that was passed in
filenames: Array of size N for the actual filepath of the file that corresponds to the data
PAs: Array of N for the parallactic angle rotation of the target (used for ADI) [in degrees]
wvs: Array of N wavelengths of the images (used for SDI) [in microns]. For polarization data, defaults to "None"
wcs: Array of N wcs astormetry headers for each image.
IWA: a floating point scalar (not array). Specifies to inner working angle in pixels
output: Array of shape (b, len(files), len(uniq_wvs), y, x) where b is the number of different KL basis cutoffs
creator: string for creator of the data (used to identify pipelines that call pyklip)
klipparams: a string that saves the most recent KLIP parameters
Methods:
readdata(): reread in the dadta
savedata(): save a specified data in the GPI datacube format (in the 1st extension header)
calibrate_output(): flux calibrate the output data
"""
##########################
###Class Initilization ###
##########################
#some static variables to define the GPI instrument
lenslet_scales = {} # in arcsec/pixel
centralwave = {} # in microns
fpm_diam = {} # in pixels
fpm_yx = {} # in pixels
flux_zeropt = {}
pupil_diam = {} # in meters
#spot_ratio = {} #w.r.t. central star
#ifs_rotation = 0.0 # degrees CCW from +x axis to zenith
observatory_latitude = None
## read in GPI configuration file and set these static variables
package_directory = os.path.dirname(os.path.abspath(__file__))
configfile = package_directory + "/" + "NIRC2.ini"
config = ConfigParser.ConfigParser()
try:
config.read(configfile)
#get list of pixel scales for different cameras and time periods
cameras = ['narrow_pre150413', 'narrow_post150413', 'medium', 'wide']
for cam in cameras:
lenslet_scales[cam] = float(config.get("pixel_scales", "pixel_scale_{0}".format(cam))) # arcsec/pix
#get IFS rotation
#ifs_rotation = float(config.get("instrument", "ifs_rotation")) #degrees
#get some information specific to each band
bands = ['J', 'H', 'K', 'Ks', 'Kp', 'Lp', 'Ms']
for band in bands:
centralwave[band] = float(config.get("instrument", "cen_wave_{0}".format(band)))
flux_zeropt[band] = float(config.get("instrument", "zero_pt_flux_{0}".format(band)))
#spot_ratio[band] = float(config.get("instrument", "APOD_{0}".format(band)))
fpms = ['corona100', 'corona150', 'corona200', 'corona300', 'corona400',
'corona600', 'corona800', 'corona1000', 'corona1500', 'corona2000']
for fpm in fpms:
fpm_diam[fpm] = float(config.get("instrument", "fpm_diam_{0}".format(fpm))) # arcsec
for fpm in fpms:
fpm_yx['narrow_' + fpm] = eval(config.get("instrument", "fpm_yx_narrow_{0}".format(fpm))) # (y,x) tuple [pixels]
fpm_yx['medium_' + fpm] = eval(config.get("instrument", "fpm_yx_medium_{0}".format(fpm))) # (y,x) tuple [pixels]
fpm_yx['wide_' + fpm] = eval(config.get("instrument", "fpm_yx_wide_{0}".format(fpm))) # (y,x) tuple [pixels]
pupils = ['incircle', 'largehex', 'smallhex', 'open']
for pupil in pupils:
pupil_diam[pupil] = float(config.get("instrument", "pupil_diam_{0}".format(pupil))) # meters
observatory_latitude = float(config.get("observatory", "observatory_lat"))
except ConfigParser.Error as e:
print("Error reading GPI configuration file: {0}".format(e.message))
raise e
####################
### Constructors ###
####################
def __init__(self, filepaths=None, highpass=False, find_star='auto', meas_star_flux=False, guess_star=None):
"""
Initialization code for NIRC2Data
Note:
see class docstring for argument details
"""
super(NIRC2Data, self).__init__()
self._output = None
self.flipx = True ### TODO: is this right?????
if filepaths is None:
self._input = None
self._centers = None
self._filenums = None
self._filenames = None
self._PAs = None
self._wvs = None
self._wcs = None
self._IWA = None
self.spot_flux = None
self.star_flux = None
self.contrast_scaling = None
self.prihdrs = None
self.exthdrs = None
self.lenslet_scale = None
self.pupil_diam = None
else:
self.readdata(filepaths, highpass=highpass, find_star=find_star, meas_star_flux=meas_star_flux, guess_star=guess_star)
################################
### Instance Required Fields ###
################################
@property
def input(self):
return self._input
@input.setter
def input(self, newval):
self._input = newval
@property
def centers(self):
return self._centers
@centers.setter
def centers(self, newval):
self._centers = newval
@property
def filenums(self):
return self._filenums
@filenums.setter
def filenums(self, newval):
self._filenums = newval
@property
def filenames(self):
return self._filenames
@filenames.setter
def filenames(self, newval):
self._filenames = newval
@property
def PAs(self):
return self._PAs
@PAs.setter
def PAs(self, newval):
self._PAs = newval
@property
def wvs(self):
return self._wvs
@wvs.setter
def wvs(self, newval):
self._wvs = newval
@property
def wcs(self):
return self._wcs
@wcs.setter
def wcs(self, newval):
self._wcs = newval
@property
def IWA(self):
return self._IWA
@IWA.setter
def IWA(self, newval):
self._IWA = newval
@property
def output(self):
return self._output
@output.setter
def output(self, newval):
self._output = newval
###############
### Methods ###
###############
[docs]
def readdata(self, filepaths, highpass=False, find_star='auto', meas_star_flux=False, guess_star=None):
"""
Method to open and read a list of NIRC2 data
Args:
filespaths: a list of filepaths
highpass: if True, run a Gaussian high pass filter (default size is sigma=imgsize/10)
can also be a number specifying FWHM of box in pixel units
find_star: (default) 'auto' will first try to get the star center coordinates from the FITS
header PSFCENTX & PSFCENTY keywords, and if that fails it will do a Radon transform to
locate the star via the diffraction spikes (and store the star center in the header for
future use). True will force the Radon transform; False will skip the Radon transform
even if no center is found in the header.
guess_star: (default) None. Otherwise an [x, y] coordinate guess for where the star is for find_star.
The code will automatially guess the star location based on the FPM otherwise.
Returns:
Technically none. It saves things to fields of the NIRC2Data object. See object doc string
"""
#check to see if user just inputted a single filename string
if isinstance(filepaths, str):
filepaths = [filepaths]
#make some lists for quick appending
data = []
filenums = []
filenames = []
rot_angles = []
wvs = []
centers = []
wcs_hdrs = []
star_fluxes = []
spot_fluxes = []
prihdrs = []
#Create a threadpool for high pass filter
pool = None
if highpass:
pool = mp.Pool()
#extract data from each file
for index, filepath in enumerate(filepaths):
cube, center, pa, wv, astr_hdrs, filt_band, fpm_band, pupil, star_flux, spot_flux, prihdr, exthdr, camera, obsdate =\
_nirc2_process_file(filepath, highpass=highpass, find_star=find_star, meas_star_flux=meas_star_flux, pool=pool, guess_star=guess_star)
data.append(cube)
centers.append(center)
star_fluxes.append(star_flux)
spot_fluxes.append(spot_flux)
rot_angles.append(pa)
wvs.append(wv)
filenums.append(np.ones(pa.shape[0]) * index)
wcs_hdrs.append(astr_hdrs)
prihdrs.append(prihdr)
#filename = np.chararray(pa.shape[0])
#filename[:] = filepath
filenames.append([filepath for i in range(pa.shape[0])])
# Close threadpool
if highpass:
pool.close()
pool.join()
#convert everything into numpy arrays
#reshape arrays so that we collapse all the files together (i.e. don't care about distinguishing files)
data = np.array(data)
dims = data.shape
data = data.reshape([dims[0] * dims[1], dims[2], dims[3]])
filenums = np.array(filenums).reshape([dims[0] * dims[1]])
filenames = np.array(filenames).reshape([dims[0] * dims[1]])
rot_angles = np.array(rot_angles).reshape([dims[0] * dims[1]]) # want North Up
wvs = np.array(wvs).reshape([dims[0] * dims[1]])
wcs_hdrs = np.array(wcs_hdrs).reshape([dims[0] * dims[1]])
centers = np.array(centers).reshape([dims[0] * dims[1], 2])
star_fluxes = np.array(star_fluxes).reshape([dims[0] * dims[1]])
spot_fluxes = np.array(spot_fluxes).reshape([dims[0] * dims[1]])
#select correct pixel scale
#narrow camera pixel scale changed after servicing on 2015-04-13
date_2015_4_13 = time.strptime("2015-4-13", "%Y-%m-%d")
if camera=='narrow':
if obsdate < date_2015_4_13:
lenslet_scale = NIRC2Data.lenslet_scales['narrow_pre150413']
elif obsdate >= date_2015_4_13:
lenslet_scale = NIRC2Data.lenslet_scales['narrow_post150413']
else:
lenslet_scale = NIRC2Data.lenslet_scales[camera]
self.lenslet_scale = lenslet_scale
#set these as the fields for the GPIData object
self._input = data
self._centers = centers
self._filenums = filenums
self._filenames = filenames
self._PAs = rot_angles
self._wvs = wvs
self._wcs = [None] #wcs_hdrs
self.spot_flux = spot_fluxes
if fpm_band in NIRC2Data.fpm_diam:
self._IWA = NIRC2Data.fpm_diam[fpm_band]/lenslet_scale/2.0
else:
# for unocculted data
self._IWA = 0
self.star_flux = star_fluxes
self.contrast_scaling = 1./star_fluxes #GPIData.spot_ratio[ppm_band]/np.tile(np.mean(spot_fluxes.reshape(dims[0], dims[1]), axis=0), dims[0])
self.prihdrs = prihdrs
self.pupil_diam = NIRC2Data.pupil_diam[pupil]
self.flipx = False
#self.exthdrs = exthdrs
[docs]
def savedata(self, filepath, data, klipparams = None, filetype = None, zaxis = None, center=None, astr_hdr=None,
fakePlparams = None, more_keywords=None):
"""
Save data in a GPI-like fashion. Aka, data and header are in the first extension header
Note: In principle, the function only works inside klip_dataset(). In order to use it outside of klip_dataset,
you need to define the following attribute:
dataset.output_centers = dataset.centers
Inputs:
filepath: path to file to output
data: 2D or 3D data to save
klipparams: a string of klip parameters
filetype: filetype of the object (e.g. "KL Mode Cube", "PSF Subtracted Spectral Cube")
zaxis: a list of values for the zaxis of the datacub (for KL mode cubes currently)
astr_hdr: wcs astrometry header (None for NIRC2)
center: center of the image to be saved in the header as the keywords PSFCENTX and PSFCENTY in pixels.
The first pixel has coordinates (0,0)
fakePlparams: fake planet params
more_keywords (dictionary) : a dictionary {key: value, key:value} of header keywords and values which will
written into the primary header
"""
hdulist = fits.HDUList()
hdulist.append(fits.PrimaryHDU(header=self.prihdrs[0]))
hdulist.append(fits.ImageHDU(data=data, name="Sci"))
# save all the files we used in the reduction
# we'll assume you used all the input files
# remove duplicates from list
filenames = np.unique(self.filenames)
nfiles = np.size(filenames)
hdulist[0].header["DRPNFILE"] = nfiles
for i, thispath in enumerate(filenames):
thispath = thispath.replace("\\", '/')
splited = thispath.split("/")
fname = splited[-1]
# matches = re.search('S20[0-9]{6}[SE][0-9]{4}', fname)
filename = fname#matches.group(0)
hdulist[0].header["FILE_{0}".format(i)] = filename
# write out psf subtraction parameters
# get pyKLIP revision number
pykliproot = os.path.dirname(os.path.dirname(os.path.realpath(__file__)))
# the universal_newline argument is just so python3 returns a string instead of bytes
# this will probably come to bite me later
try:
pyklipver = pyklip.__version__
except:
pyklipver = "unknown"
hdulist[0].header['PSFSUB'] = "pyKLIP"
hdulist[0].header.add_history("Reduced with pyKLIP using commit {0}".format(pyklipver))
#if self.creator is None:
# hdulist[0].header['CREATOR'] = "pyKLIP-{0}".format(pyklipver)
#else:
# hdulist[0].header['CREATOR'] = self.creator
# hdulist[0].header.add_history("Reduced by {0}".self.creator)
# store commit number for pyklip
hdulist[0].header['pyklipv'] = pyklipver
if klipparams is not None:
hdulist[0].header['PSFPARAM'] = klipparams
hdulist[0].header.add_history("pyKLIP reduction with parameters {0}".format(klipparams))
if fakePlparams is not None:
hdulist[0].header['FAKPLPAR'] = fakePlparams
hdulist[0].header.add_history("pyKLIP reduction with fake planet injection parameters {0}".format(fakePlparams))
if filetype is not None:
hdulist[0].header['FILETYPE'] = filetype
if zaxis is not None:
#Writing a KL mode Cube
if "KL Mode" in filetype:
hdulist[0].header['CTYPE3'] = 'KLMODES'
#write them individually
for i, klmode in enumerate(zaxis):
hdulist[0].header['KLMODE{0}'.format(i)] = klmode
#use the dataset astr hdr if none was passed in
#if astr_hdr is None:
# print(self.wcs[0])
# astr_hdr = self.wcs[0]
if astr_hdr is not None:
#update astro header
#I don't have a better way doing this so we'll just inject all the values by hand
astroheader = astr_hdr.to_header()
exthdr = hdulist[0].header
exthdr['PC1_1'] = astroheader['PC1_1']
exthdr['PC2_2'] = astroheader['PC2_2']
try:
exthdr['PC1_2'] = astroheader['PC1_2']
exthdr['PC2_1'] = astroheader['PC2_1']
except KeyError:
exthdr['PC1_2'] = 0.0
exthdr['PC2_1'] = 0.0
#remove CD values as those are confusing
exthdr.remove('CD1_1')
exthdr.remove('CD1_2')
exthdr.remove('CD2_1')
exthdr.remove('CD2_2')
exthdr['CDELT1'] = 1
exthdr['CDELT2'] = 1
#use the dataset center if none was passed in
if center is None:
center = self.output_centers[0]
if center is not None:
hdulist[0].header.update({'PSFCENTX':center[0],'PSFCENTY':center[1]})
hdulist[0].header.update({'CRPIX1':center[0],'CRPIX2':center[1]})
hdulist[0].header.add_history("Image recentered to {0}".format(str(center)))
# store extra keywords in header
if more_keywords is not None:
for hdr_key in more_keywords:
hdulist[0].header[hdr_key] = more_keywords[hdr_key]
try:
hdulist.writeto(filepath, overwrite=True)
except TypeError:
hdulist.writeto(filepath, clobber=True)
hdulist.close()
[docs]
def calibrate_data(self, units="contrast"):
"""
Calibrates the flux of the output of PSF subtracted data.
Args:
img: unclaibrated image.
If spectral is not set, this can either be a 2-D or 3-D broadband image
where the last two dimensions are [y,x]
If specetral is True, this is a 3-D spectral cube with shape [wv,y,x]
spectral: if True, this is a spectral datacube. Otherwise, it is a broadband image.
units: currently only support "contrast" w.r.t central star
Return:
img: calibrated image of the same shape (this is the same object as the input!!!)
"""
if units == "contrast":
for i in range(self.output.shape[0]):
self.output[i] *= self.contrast_scaling[:, None, None]
[docs]
def calibrate_output(self, img, spectral=False, units="contrast"):
"""
Calibrates the flux of the output of PSF subtracted data.
Assumes the broadband flux calibration is just multiplication by a single scalar number whereas spectral
datacubes may have a separate calibration value for each wavelength
Args:
img: unclaibrated image.
If spectral is not set, this can either be a 2-D or 3-D broadband image
where the last two dimensions are [y,x]
If specetral is True, this is a 3-D spectral cube with shape [wv,y,x]
spectral: if True, this is a spectral datacube. Otherwise, it is a broadband image.
units: currently only support "contrast" w.r.t central star
Return:
img: calibrated image of the same shape (this is the same object as the input!!!)
"""
if units == "contrast":
if spectral:
# spectral cube, each slice needs it's own calibration
numwvs = img.shape[0]
img *= self.contrast_scaling[:numwvs, None, None]
else:
# broadband image
img *= np.nanmean(self.contrast_scaling)
return img
######################
## Static Functions ##
######################
def _nirc2_process_file(filepath, highpass=False, find_star='auto', meas_star_flux=False, pool=None, guess_star=None):
"""
Method to open and parse a NIRC2 file
Args:
filepath: the file to open
highpass: if True, run a Gaussian high pass filter (default size is sigma=imgsize/10)
can also be a number specifying FWHM of box in pixel units
find_star: (default) 'auto' will first try to get the star center coordinates from the FITS
header PSFCENTX & PSFCENTY keywords, and if that fails it will do a Radon transform to
locate the star via the diffraction spikes (and store the star center in the header for
future use). True will force the Radon transform; False will skip the Radon transform
even if no center is found in the header.
meas_star_flux: if True, measures the stellar flux
pool: optional to pass along a threadpool (mainly for highpass filtering multiple images)
guess_star: (default) None. Otherwise an [x, y] coordinate guess for where the star is for find_star.
The code will automatially guess the star location based on the FPM otherwise.
Returns: (using z as size of 3rd dimension, z=1 for NIRC2)
cube: 3D data cube from the file. Shape is (z,256,256)
center: array of shape (z,2) giving each datacube slice a [xcenter,ycenter] in that order
parang: array of z of the parallactic angle of the target (same value just repeated z times)
wvs: array of z of the wavelength of each datacube slice. (For NIRC2, wvs = [None])
astr_hdrs: array of z of the WCS header for each datacube slice. (For NIRC2, wcs = [None])
filt_band: the band (Ks, Lp, Ms) used in the data (string)
fpm_band: For NIRC2, fpm_band = [None]
ppm_band: For NIRC2, ppm_band = [None]
spot_fluxes: For NIRC2, array of z containing 1.0 for each image
prihdr: primary header of the FITS file
exthdr: For NIRC2, None
"""
print("Reading File: {0}".format(filepath))
hdulist = fits.open(filepath, mode='update', ignore_missing_end=True) # some NIRC2 FITS files are missing END card
try:
#grab the data and headers
cube = hdulist[0].data
exthdr = None #hdulist[1].header
prihdr = hdulist[0].header
obsdate = time.strptime(prihdr['DATE-OBS'], "%Y-%m-%d") # observation date
#get some instrument configuration from the primary header
filt_band = prihdr['FILTER'].split('+')[0].strip()
fpm_band = prihdr['SLITNAME']
camera = prihdr['CAMNAME']
pupil = prihdr['PMSNAME']
rotmode = prihdr['ROTMODE'].lower()
#for NIRC2, we only have broadband but want to keep the GPI array shape to make processing easier
if prihdr['CURRINST'].strip() == 'NIRC2':
wvs = [1.0]
# Pull PA from header if already there, or try to calculate it from header.
if 'ROTNORTH' in prihdr.keys():
parang = prihdr['ROTNORTH']*np.ones(1)
else:
try:
parang = get_pa(hdulist, obsdate=obsdate, rotmode=rotmode, write_hdr=True)*np.ones(1)
except:
parang = np.nan*np.ones(1)
# manual guess for star centering?
if guess_star is None:
# automatic guess based on fpm
guess_star = NIRC2Data.fpm_yx[('_').join([camera, fpm_band])]
if find_star is True:
# Ignore header and use Radon transform to find star center.
try:
center = [get_star(hdulist, ctr=guess_star, obsdate=obsdate,
hp_size=0, im_smooth=0, sp_width=0, rad=100, radon_wdw=400, smooth=1,
write_hdr=True, pool=pool, silent=True)]
except:
center = [[np.nan, np.nan]]
elif find_star is False:
# Only try to get star center from headers (no Radon transform).
try:
center = [[prihdr['PSFCENTX'], prihdr['PSFCENTY']]]
except:
center = [[np.nan, np.nan]]
elif find_star=='auto':
# Pull star center from header if already there, or try to find it via Radon transform.
# Radon assumes star is near center of FPM and may fail if otherwise.
if 'PSFCENTX' and 'PSFCENTY' in prihdr.keys():
center = [[prihdr['PSFCENTX'], prihdr['PSFCENTY']]]
else:
try:
center = [get_star(hdulist, ctr=guess_star, obsdate=obsdate,
hp_size=0, im_smooth=0, sp_width=0, rad=100, radon_wdw=400, smooth=1,
write_hdr=True, pool=pool, silent=True)]
except:
center = [[np.nan, np.nan]]
else:
raise ValueError("Unsupported value for find_star; only 'auto', True, or False are accepted.")
# No cube flipping anymore!
flipped_cube = cube
star_flux = np.nan
if meas_star_flux:
star_flux = calc_starflux(flipped_cube, center)
cube = flipped_cube.reshape([1, flipped_cube.shape[0], flipped_cube.shape[1]]) #maintain 3d-ness
astr_hdrs = np.repeat(None, 1)
spot_fluxes = [[1]] #not suported currently
finally:
pass
#high pass filter
highpassed = False
if isinstance(highpass, bool):
if highpass:
cube = high_pass_filter_imgs(cube, pool=pool)
highpassed = True
else:
# should be a number
if isinstance(highpass, (float, int)):
highpass = float(highpass)
fourier_sigma_size = (cube.shape[1]/(highpass)) / (2*np.sqrt(2*np.log(2)))
cube = high_pass_filter_imgs(cube, filtersize=fourier_sigma_size, pool=pool)
highpassed = True
return cube, center, parang, wvs, astr_hdrs, filt_band, fpm_band, pupil, star_flux, spot_fluxes, prihdr, exthdr, camera, obsdate
[docs]
def calc_starflux(cube, center):
"""
Fits a 2D Gaussian to an image to calculate the peak pixel value of
the central star. The code assumes an unobscurated PSF.
Args:
cube: 2D image array. Shape is (256,256)
center: star center in image in (x,y)
Returns:
Amplitude: Best fit amplitude of the 2D Gaussian.
"""
dims = cube.shape
y, x = np.meshgrid( np.arange(dims[0]), np.arange(dims[1]) )
# Initializing Model. Fixing the rotation and the X, Y location of the star.
g_init = models.Gaussian2D(cube.max(), x_mean=center[0][0], y_mean=center[0][1], x_stddev=5, y_stddev=5, \
fixed={'x_mean':True,'y_mean':True,'theta':True})
# Initializing Levenburg-Marquart Least-Squares fitting routine.
fit_g = fitting.LevMarLSQFitter()
# Fitting the amplitude, x_stddev and y_stddev
g = fit_g(g_init, y, x, cube)
return [[g.amplitude]]
[docs]
def measure_star_flux(img, star_x, star_y):
"""
Measure star peak fluxes using a Gaussian matched filter
Args:
img: 2D frame with unobscured, unsaturated PSF
star_x, star_y: coordinates of the star
Return:
star_f: star flux
"""
flux, fwhm, xfit, yfit = gaussfit2d(img, star_x, star_y, refinefit=False)
if flux == np.inf: flux == np.nan
print(flux, fwhm, xfit, yfit)
return flux
[docs]
def get_pa(hdulist, obsdate=None, rotmode=None, mean_PA=True, write_hdr=True):
"""
Given a FITS data-header unit list (HDUList), returns the NIRC2 PA in [radians].
PA is angle of detector relative to sky; ROTMODE is rotator tracking mode;
PARANG is parallactic angle astrometric; INSTANGL is instrument angle;
ROTPOSN is rotator physical position.
Additional PA offset of -0.252 or -0.262 deg is applied for NIRC2 narrow cam
depending on observation date.
NOTE that the PA sign is flipped at the very end before output to conform to
pyKLIP's rotation convention.
Inputs:
hdulist: a FITS HDUList (NOT a single HDU).
obsdate: date of observation; will try to get from prihdr if not provided.
rotmode: 'vertical angle' for ADI mode with PA rotating on detector, or
'position angle' for mode with PA orientation fixed on detector.
mean_pa: if True (default), return the mean PA during the exposure.
If False, return the PA at the start of the exposure only.
Only applies to vertical angle mode.
write_hdr: if True (default), writes keys to file header and saves them.
"""
prihdr = hdulist[0].header
# Date of NIRC2 servicing that changed PA offset.
date_2015_4_13 = time.strptime("2015-4-13", "%Y-%m-%d")
# If don't have it, get observation date (UT) from header and make a time object.
if obsdate is None:
obsdate = time.strptime(prihdr['DATE-OBS'], "%Y-%m-%d")
# Additional offset to narrow cam PA not included in INSTANGL keyword.
# This offset changed after instrument servicing on April 13, 2015.
if prihdr['CAMNAME'].lower() == 'narrow':
if obsdate < date_2015_4_13:
zp_offset = -0.252 # [deg]; from Yelda et al. 2010
elif obsdate >= date_2015_4_13:
zp_offset = -0.262 # [deg]; from Service et al. 2016
else:
zp_offset = 0.
print("WARNING: No PA offset applied.")
if rotmode is None:
global _last_rotmode
try:
rotmode = prihdr['ROTMODE'].strip().lower()
_last_rotmode = rotmode
except:
rotmode = _last_rotmode.lower()
rotposn = prihdr['ROTPOSN'] # [deg]
instangl = prihdr['INSTANGL'] # [deg]
if rotmode.lower() == 'vertical angle':
parang = prihdr['PARANG'] # [deg]
pa_deg = parang + rotposn - instangl + zp_offset # [deg]
elif rotmode.lower() == 'position angle':
pa_deg = rotposn - instangl + zp_offset # [deg]
else:
raise NotImplementedError
if mean_PA and (rotmode.lower() == 'vertical angle'):
# Get info for PA smearing calculation.
epochobj = prihdr['DATE-OBS']
name = prihdr['targname']
expref = prihdr['itime']
coaddref = prihdr['coadds']
sampref = prihdr['sampmode']
msrref = prihdr['MULTISAM']
xdimref = prihdr['naxis1']
ydimref = prihdr['naxis2']
tel = prihdr['TELESCOP']
dec = prihdr['DEC'] + prihdr['DECOFF']
if tel.lower() == 'keck ii':
tel = 'keck2' # just cleaning up str
# Calculate total time of exposure (integration + readout).
if sampref == 2: totexp = ( expref + 0.18*(xdimref/1024.)**2) * coaddref
if sampref == 3: totexp = ( expref + (msrref-1)*0.18*(xdimref/1024.)**2)*coaddref
tinteg = totexp # [seconds]
totexp = totexp/3600. # [hours]
# Get hour angle at start of exposure.
tmpahinit = prihdr['HA'] # [deg]
ahobs = 24.*tmpahinit/360. # [hours]
# Estimate vertical position angle at each second of the exposure.
vp = [] #fltarr(round(3600.*totexp))
for j in range(0, int(round(3600.*totexp))-1):
ahtmp = ahobs + (j*1.+0.001)/3600. # [hours]
vp.append(par_angle(ahtmp, dec, NIRC2Data.observatory_latitude))
if j == 0: vpref = vp[0]
vp = np.array(vp)
# Handle case where PA crosses 0 <--> 360.
vp[vp < 0.] += 360.
vp[vp > 360.] -= 360.
if vpref < 0.: vpref += 360.
if vpref > 360.: vpref -= 360.
# Check that images near PA=0 are handled correctly.
if any(vp > 350) & any(vp < 10):
vp[vp > 350] -= 360
vpmean = np.nanmean(vp)
if (vpmean < 0) & (vpref > 350):
vpmean += 360.
pa_deg_mean = pa_deg + (vpmean - vpref)
else:
pa_deg_mean = pa_deg
vpmean = np.nan
vpref = np.nan
totexp = np.nan
if write_hdr:
prihdr['TOTEXP'] = (totexp, 'Total exposure time [hours]')
prihdr['PASTART'] = (pa_deg, "Position angle at exposure start [deg]")
prihdr['PASMEAR'] = ((vpmean - vpref), "Exposure's weighted-mean PA minus PASTART [deg]")
prihdr['ROTNORTH'] = (pa_deg_mean, "Mean PA of North during exposure [deg]")
hdulist.flush(output_verify='ignore')
# Flip sign to conform to pyKLIP rotation convention.
return pa_deg_mean
[docs]
def par_angle(HA, dec, lat):
"""
Compute the parallactic angle, given hour angle (HA [hours]),
declination (dec [deg]), and latitude (lat [deg]). Returns
parallactic angle in [deg].
"""
HA_rad = np.radians(HA*15.) # [hours] -> [rad]
dec_rad = np.radians(dec) # [deg] -> [rad]
lat_rad = np.radians(lat) # [deg] -> [rad]
parallang = -np.arctan2(-np.sin(HA_rad), # [rad]
np.cos(dec_rad)*np.tan(lat_rad) - np.sin(dec_rad)*np.cos(HA_rad))
return np.degrees(parallang) # [deg]
[docs]
def get_star(hdulist, ctr, obsdate, hp_size=0, im_smooth=0, sp_width=0, spike_angles=None,
r_mask='all', rad=100, rad_out=np.inf, radon_wdw=400, smooth=1, PAadj=0.,
write_hdr=True, pool=None, silent=True):
"""
Runs Radon transform star-finding algorithm on image and (by default) saves the results
in the original FITS header.
Inputs:
hdulist: a FITS HDUList (NOT a single HDU).
ctr: (y,x) coordinate pair estimate for star position for image [pix].
obsdate: date of observation; will try to get from prihdr if not provided.
hp_size: size of high-pass filter box (via Fourier transform) in [pix].
im_smooth: sigma of smoothing Gaussian function in [pix].
sp_width: width of diffraction spike mask in [pix]; default is 0 (no masking).
spike_angles: list of discrete angles from the assumed star positions along
which the radon transform will sum intensity to search for the true star
position (it picks the maximum sum). These should match the angles
of the strongest diffraction spikes [deg].
r_mask: 'all' to mask out circle around ctr coords; anything else to do no radial masking.
rad: r_mask=='all' will mask out all r <= rad [pix].
rad_out: r_mask=='all' will mask out all r > rad_out [pix].
radon_window: half width of the radon sampling region; size_window = image.shape[0]/2 is suggested.
m & M: The sampling region will be (-M*size_window, -m*size_window)U(m*size_window, M*size_window).
smooth: smooth the radon cost function; for one pixel, replace it by the
average of its +/- smooth neighbours; default = 2.
PAadj: optional angle by which to rotate diffraction spike angles in [radians].
write_hdr: (default) True will write the Radon transform star center to the original
FITS header in PSFCENTX & PSFCENTY keywords.
pool: multiprocessing pool for highpass filtering and other parallel uses.
silent: (default) True to suppress additional output to stdout.
Outputs:
Returns [X,Y] list of Radon transform star center. Default is to also write
the star coordinates to PSFCENTX & PSFCENTY in original FITS header.
"""
from pyklip.instruments.utils.radonCenter import searchCenter
from scipy.ndimage.filters import median_filter, gaussian_filter
if not silent: print("Finding star...")
data = hdulist[0].data.copy()
hdr = hdulist[0].header
# median filter data to replace NaN's with median values of neighbors.
wh = np.where(np.isnan(data))
data_medfilt = median_filter(np.nan_to_num(data), size=9)
# replace all NaN in f with median values of neighbors.
data[wh] = data_medfilt[wh]
data = np.ma.masked_invalid(data)
# Highpass filtering.
if hp_size != 0:
# should be a number
fourier_sigma_size = (data.shape[1]/float(hp_size)) / (2*np.sqrt(2*np.log(2)))
data_filt = high_pass_filter_imgs(np.array([data]), filtersize=fourier_sigma_size, pool=pool)[0]
else:
data_filt = data
# Image smoothing.
if im_smooth != 0:
data_filt = gaussian_filter(data_filt, im_smooth)
# Build cartesian coordinate grid and radius map.
yy, xx = np.mgrid[0:data_filt.shape[0]:1, 0:data_filt.shape[1]:1]
radii = np.sqrt((yy - ctr[0])**2 + (xx - ctr[1])**2)
# Additional offset to narrow cam PA not included in INSTANGL keyword.
# This offset changed after instrument servicing on April 13, 2015.
date_2015_4_13 = time.strptime("2015-4-13", "%Y-%m-%d")
if hdr['CAMNAME'].lower()=='narrow':
if obsdate < date_2015_4_13:
zp_offset = -0.252 # [deg]; from Yelda et al. 2010
elif obsdate >= date_2015_4_13:
zp_offset = -0.262 # [deg]; from Service et al. 2016
elif hdr['CAMNAME'].lower()=='wide':
zp_offset = 0.
# Select angles along which to perform radon based on rotator mode.
# Position angle (pa) of telescope optics only, from header info.
if hdr['ROTMODE'].lower() == 'vertical angle':
pa_tele = hdr['INSTANGL'] - zp_offset - hdr['ROTPOSN'] # [deg]
# Other case is either 'position angle' mode or unidentified.
else:
pa_tele = hdr['PARANG'] # [deg]
if spike_angles is None:
# Angles at which diffraction spikes occur in NIRC2 data [deg].
spike_angles = pa_tele + 30.0 + np.arange(3)*60.0 + PAadj # [deg]
# Boolean mask excluding everything except diffraction spikes.
spikemask = ~ make_spikemask(data_filt, hdr, ctr, spike_angles, yy, xx, sp_width) # ~ is inverse boolean
mask_total = spikemask.copy()
# Optionally mask stellar halo with additional circular mask centered on ctr.
# Masked regions are r <= rad and r > rad_out.
if r_mask=='all':
mask_total[(spikemask==False) & (radii <= rad)] = True
mask_total[(spikemask==False) & (radii > rad_out)] = True
else:
pass
data_masked = np.ma.array(data_filt, mask=mask_total)
# Perform radon transform search for star.
if not silent: print("Performing radon transform search...")
(x_radon, y_radon) = searchCenter(data_masked.filled(0.), ctr[1], ctr[0],
size_window=radon_wdw, size_cost=7, m=0.2, M=0.8,
smooth=smooth, theta=spike_angles)
if not silent: print("radon y,x = {0}, {1}".format(y_radon, x_radon))
if write_hdr:
# Update the original FITS header with new star coordinates.
hdr['PSFCENTX'] = (x_radon, 'Star X numpy coord')
hdr['PSFCENTY'] = (y_radon, 'Star Y numpy coord')
hdulist.flush(output_verify='ignore')
if not silent: print("Wrote new star coordinates to FITS header in PSFCENTX & PSFCENTY.")
return [x_radon, y_radon]
[docs]
def make_spikemask(data, hdr, ctr, spike_angles, yy, xx, width=31):
"""
Construct diffraction spike mask from FITS header information.
data: 2-D ndarray image to be masked (just to get size of array).
hdr: FITS header for image constructing mask for.
ctr: (y,x) coordinates for center of diffraction spike pattern
(usually the estimated star location).
spike_angles: position angles for diffraction spikes in image [radians].
yy: mgrid or indices 2-D array of pixel y-coordinates.
xx: mgrid or indices 2-D array of pixel x-coordinates.
width: int or float width of spike mask in [pixels]; 0 for no mask.
"""
mask = np.zeros(data.shape, dtype=np.bool)
if (width > 0) & (not np.isnan(width)):
yy_ctr = yy - ctr[0]
xx_ctr = xx - ctr[1]
for ang in np.radians(spike_angles):
# Slope times x for each spike.
mx = np.tan(ang)*xx_ctr
# Mask with span 0.5*width above and below spike.
band = np.abs(0.5*width/np.cos(ang))
# Change spikemask to True anywhere inside mask.
spikemask = (yy_ctr <= (mx + band)) & (yy_ctr >= (mx - band))
# Replace False with True in mask wherever either spikemask or mask is True.
mask = mask | spikemask
else:
mask = np.ones(data.shape, dtype=np.bool)
return mask