pyklip package
Subpackages
- pyklip.fmlib package
- pyklip.instruments package
- Subpackages
- pyklip.instruments.P1640_support package
- Submodules
- pyklip.instruments.P1640_support.P1640_cube_checker module
- pyklip.instruments.P1640_support.P1640_cube_checker_interactive module
- pyklip.instruments.P1640_support.P1640_spot_checker module
- pyklip.instruments.P1640_support.P1640contrast module
- pyklip.instruments.P1640_support.P1640cores module
- pyklip.instruments.P1640_support.P1640spots module
- pyklip.instruments.P1640_support.P1640utils module
- Module contents
- pyklip.instruments.utils package
- pyklip.instruments.P1640_support package
- Submodules
- pyklip.instruments.CHARIS module
CHARISData
CHARISData.input
CHARISData.centers
CHARISData.filenums
CHARISData.filenames
CHARISData.PAs
CHARISData.wvs
CHARISData.wcs
CHARISData.IWA
CHARISData.OWA
CHARISData.output
CHARISData.wv_indices
CHARISData.spot_flux
CHARISData.dn_per_contrast
CHARISData.flux_units
CHARISData.prihdrs
CHARISData.exthdrs
CHARISData.bad_sat_spots
CHARISData.readdata()
CHARISData.savedata()
CHARISData.calibrate_output()
CHARISData.IWA
CHARISData.OWA
CHARISData.PAs
CHARISData.calibrate_output()
CHARISData.centers
CHARISData.centralwave
CHARISData.filenames
CHARISData.filenums
CHARISData.flux_zeropt
CHARISData.fpm_diam
CHARISData.generate_host_star_psfs()
CHARISData.generate_psfs()
CHARISData.ifs_rotation
CHARISData.input
CHARISData.ivars
CHARISData.lenslet_scale
CHARISData.lenslet_scale_x
CHARISData.lenslet_scale_x_err
CHARISData.lenslet_scale_y
CHARISData.lenslet_scale_y_err
CHARISData.measure_spot_over_star_ratio()
CHARISData.obs_latitude
CHARISData.obs_longitude
CHARISData.output
CHARISData.readdata()
CHARISData.ref_spot_contrast
CHARISData.ref_spot_contrast_error_fraction
CHARISData.ref_spot_contrast_gridamp
CHARISData.ref_spot_contrast_wv
CHARISData.save_platecal_cubes()
CHARISData.savedata()
CHARISData.spot_ratio
CHARISData.update_input_cubes()
CHARISData.wcs
CHARISData.wvs
generate_psf()
- pyklip.instruments.GPI module
GPIData
GPIData.input
GPIData.centers
GPIData.filenums
GPIData.filenames
GPIData.PAs
GPIData.wvs
GPIData.wcs
GPIData.IWA
GPIData.output
GPIData.wv_indices
GPIData.spot_flux
GPIData.dn_per_contrast
GPIData.flux_units
GPIData.prihdrs
GPIData.exthdrs
GPIData.bad_sat_spots
GPIData.readdata()
GPIData.savedata()
GPIData.calibrate_output()
GPIData.IWA
GPIData.PAs
GPIData.band
GPIData.bands
GPIData.calibrate_output()
GPIData.centers
GPIData.centralwave
GPIData.config
GPIData.configfile
GPIData.coronagraph_throughputs
GPIData.filenames
GPIData.filenums
GPIData.flux_zeropt
GPIData.fpm_diam
GPIData.generate_psf_cube()
GPIData.generate_psfs()
GPIData.get_radial_psf()
GPIData.hdulist
GPIData.ifs_rotation
GPIData.input
GPIData.lenslet_scale
GPIData.observatory_latitude
GPIData.output
GPIData.package_directory
GPIData.profile_dir
GPIData.profile_filename
GPIData.readdata()
GPIData.savedata()
GPIData.seps
GPIData.spectral_collapse()
GPIData.spot_ratio
GPIData.spot_ratio2
GPIData.spot_ratio_h
GPIData.spot_ratio_h2
GPIData.throughputs
GPIData.wcs
GPIData.wvs
as2pix()
calc_center()
calc_center_least_squares()
generate_psf()
get_gpi_wavelength_sampling()
measure_sat_spot_fluxes()
pix2as()
recalculate_sat_spot_fluxes()
rescale_wvs()
subtract_satspots()
- pyklip.instruments.Instrument module
Data
Data.input
Data.centers
Data.filenums
Data.filenames
Data.PAs
Data.wvs
Data.wcs
Data.IWA
Data.OWA
Data.output
Data.output_centers
Data.output_wcs
Data.creator
Data.klipparams
Data.flipx
Data.readdata()
Data.savedata()
Data.calibrate_output()
Data.IWA
Data.PAs
Data.calibrate_output()
Data.centers
Data.filenames
Data.filenums
Data.input
Data.numwvs
Data.output
Data.readdata()
Data.savedata()
Data.spectral_collapse()
Data.wcs
Data.wvs
GenericData
GenericData.input
GenericData.centers
GenericData.filenums
GenericData.filenames
GenericData.PAs
GenericData.wvs
GenericData.wcs
GenericData.IWA
GenericData.output
GenericData.IWA
GenericData.PAs
GenericData.calibrate_output()
GenericData.centers
GenericData.filenames
GenericData.filenums
GenericData.input
GenericData.output
GenericData.readdata()
GenericData.savedata()
GenericData.wcs
GenericData.wvs
- pyklip.instruments.MagAO module
MagAOData
MagAOData.input
MagAOData.centers
MagAOData.filenums
MagAOData.filenames
MagAOData.PAs
MagAOData.wvs
MagAOData.wcs
MagAOData.IWA
MagAOData.output
MagAOData.spot_flux
MagAOData.contrast_scaling
MagAOData.flux_units
MagAOData.prihdrs
MagAOData.exthdrs
MagAOData.readdata()
MagAOData.savedata()
MagAOData.calibrate_output()
MagAOData.IWA
MagAOData.OWA
MagAOData.PAs
MagAOData.band
MagAOData.bands
MagAOData.calibrate_output()
MagAOData.centers
MagAOData.centralwave
MagAOData.config
MagAOData.configfile
MagAOData.filenums
MagAOData.flipx
MagAOData.flux_zeropt
MagAOData.ghstpeak_ratio
MagAOData.ifs_rotation
MagAOData.input
MagAOData.lenslet_scale
MagAOData.observatory_latitude
MagAOData.output
MagAOData.package_directory
MagAOData.readdata()
MagAOData.savedata()
MagAOData.wcs
MagAOData.wvs
- pyklip.instruments.NIRC2 module
NIRC2Data
NIRC2Data.input
NIRC2Data.centers
NIRC2Data.filenums
NIRC2Data.filenames
NIRC2Data.PAs
NIRC2Data.wvs
NIRC2Data.wcs
NIRC2Data.IWA
NIRC2Data.output
NIRC2Data.creator
NIRC2Data.klipparams
NIRC2Data.readdata()
NIRC2Data.savedata()
NIRC2Data.calibrate_output()
NIRC2Data.IWA
NIRC2Data.PAs
NIRC2Data.band
NIRC2Data.bands
NIRC2Data.calibrate_data()
NIRC2Data.calibrate_output()
NIRC2Data.cam
NIRC2Data.cameras
NIRC2Data.centers
NIRC2Data.centralwave
NIRC2Data.config
NIRC2Data.configfile
NIRC2Data.filenames
NIRC2Data.filenums
NIRC2Data.flux_zeropt
NIRC2Data.fpm
NIRC2Data.fpm_diam
NIRC2Data.fpm_yx
NIRC2Data.fpms
NIRC2Data.input
NIRC2Data.lenslet_scales
NIRC2Data.observatory_latitude
NIRC2Data.output
NIRC2Data.package_directory
NIRC2Data.pupil
NIRC2Data.pupil_diam
NIRC2Data.pupils
NIRC2Data.readdata()
NIRC2Data.savedata()
NIRC2Data.wcs
NIRC2Data.wvs
calc_starflux()
get_pa()
get_star()
make_spikemask()
measure_star_flux()
par_angle()
- pyklip.instruments.P1640 module
- pyklip.instruments.SPHERE module
Ifs
Ifs.input
Ifs.centers
Ifs.filenums
Ifs.filenames
Ifs.ifs_rdp
Ifs.PAs
Ifs.wvs
Ifs.IWA
Ifs.output
Ifs.psfs
Ifs.psf_center
Ifs.flipx
Ifs.nfiles
Ifs.nwvs
Ifs.IWA
Ifs.PAs
Ifs.calibrate_output()
Ifs.centers
Ifs.filenames
Ifs.filenums
Ifs.ifs_rdp
Ifs.input
Ifs.output
Ifs.platescale
Ifs.readdata()
Ifs.savedata()
Ifs.wcs
Ifs.wvs
Irdis
Irdis.input
Irdis.centers
Irdis.filenums
Irdis.filenames
Irdis.PAs
Irdis.wvs
Irdis.IWA
Irdis.OWA
Irdis.output
Irdis.psfs
Irdis.psf_center
Irdis.flipx
Irdis.nfiles
Irdis.prihdrs
Irdis.nwvs
Irdis.IWA
Irdis.OWA
Irdis.PAs
Irdis.calibrate_output()
Irdis.centers
Irdis.filenames
Irdis.filenums
Irdis.input
Irdis.output
Irdis.platescale
Irdis.readdata()
Irdis.savedata()
Irdis.wavelengths
Irdis.wcs
Irdis.wvs
- pyklip.instruments.osiris module
- Module contents
- Subpackages
- pyklip.kpp package
- Subpackages
- pyklip.kpp.detection package
- pyklip.kpp.metrics package
- pyklip.kpp.stat package
- pyklip.kpp.utils package
- Module contents
- Subpackages
Submodules
pyklip.covars module
- pyklip.covars.delta(x, y, sigmas, *args)[source]
Generates a diagonal covariance matrix
C_ij = sigma_i sigma_j delta_{ij}
- Parameters:
x (np.array) – 1-D array of x coordinates
y (np.array) – 1-D array of y coordinates
sigmas (np.array) – 1-D array of errors on each pixel
- pyklip.covars.matern32(x, y, sigmas, corr_len)[source]
Generates a Matern (nu=3/2) covariance matrix that assumes x/y has the same correlation length
C_ij = sigma_i sigma_j (1 + sqrt(3) r_ij / l) exp(-sqrt(3) r_ij / l)
- Parameters:
x (np.array) – 1-D array of x coordinates
y (np.array) – 1-D array of y coordinates
sigmas (np.array) – 1-D array of errors on each pixel
corr_len (float) – correlation length of the Matern function
- Returns:
2-D covariance matrix parameterized by the Matern function
- Return type:
cov (np.array)
- pyklip.covars.sq_exp(x, y, sigmas, corr_len)[source]
Generates square exponential covariance matrix that assumes x/y has the same correlation length
C_ij = sigma_i sigma_j exp(-r_ij^2/[2 l^2])
- Parameters:
x (np.array) – 1-D array of x coordinates
y (np.array) – 1-D array of y coordinates
sigmas (np.array) – 1-D array of errors on each pixel
corr_len (float) – correlation length (i.e. standard deviation of Gaussian)
mode (string) – either “numpy”, “cython”, or None, specifying the implementation of the kernel.
- Returns:
2-D covariance matrix parameterized by the Matern function
- Return type:
cov (np.array)
pyklip.empca module
- pyklip.empca.np_calc_chisq(data, b, w, coef)[source]
Calculate chi squared
- Parameters:
im – nim x npix, single-precision numpy.ndarray. Data to be fit by the basis images
b – nvec x npts, double precision numpy.ndarray. The nvec basis images.
w – nim x npts, single-precision numpy.ndarray. Weights (inverse variances) of the data.
coef – nvec x npts, double precision numpy.ndarray. The coefficients of the basis image fits.
- Returns:
chisq, the total chi squared summed over all points and all images
- pyklip.empca.set_pixel_weights(imflat, rflat, ivar=None, mode='standard', inner_sup=17, outer_sup=66)[source]
- Parameters:
imflat – array of flattend images, shape (N, number of section indices)
rflat – radial component of the polar coordinates flattened to 1D, length = number of section indices
mode –
- ‘standard’: assume poission statistics to calculate variance as sqrt(photon count)
use inverse sqrt(variance) as pixel weights and multiply by a radial weighting
inner_sup – radius within which to supress weights
outer_sup – radius beyond which to supress weights
- Returns:
pixel weights for empca
- pyklip.empca.weighted_empca(data, weights=None, niter=25, nvec=5, randseed=1, maxcpus=1, silent=True)[source]
Perform iterative low-rank matrix approximation of data using weights.
Generated model vectors are not orthonormal and are not rotated/ranked by ability to model the data, but as a set they are good at describing the data.
- Parameters:
data – images to model
weights – weights for every pixel
niter – maximum number of iterations to perform
nvec – number of vectors to solve (rank of the approximation)
randseed – rand num generator seed; if None, don’t re-initialize
maxcpus – maximum cpus to use for parallel programming
silent – bool, whether to show chi_squared for each iteration
- Returns:
returns the best low-rank approximation to the data in a weighted least-squares sense (dot product of coefficients and basis vectors).
pyklip.fakes module
- pyklip.fakes.LSQ_gauss2d(planet_image, x_grid, y_grid, a, x_cen, y_cen, sig)[source]
Calculate the squared norm of the residuals of the model with the data. Helper function for least square fit. The model is a 2d symmetric gaussian.
- Parameters:
planet_image – stamp image (y,x) of the satellite spot.
x_grid – x samples grid as given by meshgrid.
y_grid – y samples grid as given by meshgrid.
a – amplitude of the 2d gaussian
x_cen – x center of the gaussian
y_cen – y center of the gaussian
sig – standard deviation of the gaussian
- Returns:
Squared norm of the residuals
- pyklip.fakes.PSFcubefit(frame, xguess, yguess, searchrad=10, psfs_func_list=None, wave_index=None, residuals=False, rmbackground=True, add_background2residual=False)[source]
Estimate satellite spot amplitude (peak value) by fitting a symmetric 2d gaussian. Fit parameters: x,y position, amplitude, standard deviation (same in x and y direction)
- Parameters:
frame – the data - Array of size (y,x)
xguess – x location to fit the 2d guassian to.
yguess – y location to fit the 2d guassian to.
searchrad – 1/2 the length of the box used for the fit
psfs_func_list – List of spline fit function for the PSF_cube.
wave_index – Index of the current wavelength. In [0,36] for GPI. Only used when psfs_func_list is not None.
residuals – If True (Default = False) then calculate the residuals of the sat spot fit (gaussian or PSF cube).
rmbackground – If true (default), remove any background slope to the data stamp.
add_background2residual – If True (default is false) and if rmbackground was true, it adds the background that was removed to the returned residuals.
- Returns:
- scalar, Estimation of the peak flux of the satellite spot.
ie Amplitude of the fitted gaussian.
- Return type:
returned_flux
- pyklip.fakes.airyfit2d(frame, xguess, yguess, searchrad=5, guessfwhm=3, guesspeak=1)[source]
Fits a 2d airy function to the data at point (xguess, yguess)
- Parameters:
frame – the data - Array of size (y,x)
xguess – location to fit the 2d guassian to (should be pretty accurate)
yguess – location to fit the 2d guassian to (should be pretty accurate)
searchrad – 1/2 the length of the box used for the fit
guessfwhm – approximate fwhm to fit to
- Returns:
the peakflux of the Airy function fwhm: diameter between the first minima along one axis xfit: x position yfit: y position
- Return type:
peakflux
- pyklip.fakes.convert_pa_to_image_polar(pa, astr_hdr)[source]
Given a position angle (angle to North through East), calculate what polar angle theta (angle from +X CCW towards +Y) it corresponds to
- Parameters:
pa – position angle in degrees
astr_hdr – wcs astrometry header (astropy.wcs)
- Returns:
polar angle in degrees
- Return type:
theta
- pyklip.fakes.convert_polar_to_image_pa(theta, astr_hdr)[source]
Reversed engineer from covert_pa_to_image_polar by JB. Actually JB doesn’t quite understand how it works…
- Parameters:
theta – parallactic angle in degrees
astr_hdr – wcs astrometry header (astropy.wcs)
- Returns:
polar angle in degrees
- Return type:
theta
- pyklip.fakes.gauss2d(x0, y0, peak, sigma)[source]
2d symmetric guassian function for guassfit2d
- Parameters:
x0 – center of gaussian
y0 – center of gaussian
peak – peak amplitude of guassian
sigma – stddev in both x and y directions
- pyklip.fakes.gaussfit2d(frame, xguess, yguess, searchrad=5, guessfwhm=3, guesspeak=1, refinefit=True)[source]
Fits a 2d gaussian to the data at point (xguess, yguess)
- Parameters:
frame – the data - Array of size (y,x)
xguess – location to fit the 2d guassian to (should be pretty accurate)
yguess – location to fit the 2d guassian to (should be pretty accurate)
searchrad – 1/2 the length of the box used for the fit
guessfwhm – approximate fwhm to fit to
guesspeak – approximate flux
refinefit – whether to refine the fit of the position of the guess
- Returns:
the peakflux of the gaussian fwhm: fwhm of the PFS in pixels xfit: x position (only chagned if refinefit is True) yfit: y position (only chagned if refinefit is True)
- Return type:
peakflux
- pyklip.fakes.gaussfit2dLSQ(frame, xguess, yguess, searchrad=5, fit_centroid=False, residuals=False)[source]
Estimate satellite spot amplitude (peak value) by fitting a symmetric 2d gaussian. Fit parameters: x,y position, amplitude, standard deviation (same in x and y direction)
- Parameters:
frame – the data - Array of size (y,x)
xguess – x location to fit the 2d guassian to.
yguess – y location to fit the 2d guassian to.
searchrad – 1/2 the length of the box used for the fit
fit_centroid – If False (default), disable the centroid fit and only fit the amplitude and the standard deviation
residuals – If True (Default = False) then calculate the residuals of the sat spot fit (gaussian or PSF cube).
- Returns:
- scalar, estimation of the peak flux of the satellite spot.
ie Amplitude of the fitted gaussian.
- Return type:
returned_flux
- pyklip.fakes.generate_dataset_with_fakes(dataset, fake_position_dict, fake_flux_dict, spectrum=None, PSF_cube=None, PSF_cube_wvs=None, star_type=None, mute=False, SpT_file_csv=None, real_planets_pos=None, sep_skip_real_pl=None, pa_skip_real_pl=None, dn_per_contrast=None, star_spec=None)[source]
Generate spectral datacubes with fake planets. It will do a copy of the cubes read in GPIData after having injected fake planets in them. This new set of cubes can then be reduced in the same manner as the campaign data.
Doesn’t work with remove slice: assumes that the dataset is made of a list of similar datacubes or images.
- Parameters:
dataset – An object of type GPIData. The fakes are injected directly into dataset so you should make a copy of dataset prior to running this function. In order for the function to query simbad for the spectral type of the star, the attribute object_name needs to be defined in dataset.
fake_position_dict –
- Dictionary defining the way the fake planets are positionned
- fake_position_dict[“mode”]=”sector”: Put a planet in each klip sector. Can actually generate several
datasets in which the planets will be shifted in separation and position angle with respect to one another. It can be usefull for fake based contrast curve calculation. Several parameters needs to be defined. - fake_position_dict[“annuli”]: Number of annulis in the image - fake_position_dict[“subsections”]: Number of angular sections in the image - fake_position_dict[“sep_shift”]: separation shift from the center of the sectors - fake_position_dict[“pa_shift”]: position angle shift from the center of the sectors
- fake_position_dict[“mode”]=”custom”: Put planets at given (separation, position angle).
The following parameter needs to be defined - fake_position_dict[“pa_sep_list”]: List of tuple [(r1,pa1),(r2,pa2),…] with each tuple giving
the separation and position angle of each planet to be injected.
fake_position_dict[“mode”]=”ROC”: Generate fake for ROC curves calculation. Use hard-coded parameters.
- fake_flux_dict:
Dictionary defining the way in which the flux of the fake is defined. - fake_flux_dict[“mode”]=”contrast”: Defines the contrast value of the fakes.
fake_flux_dict[“contrast”]: Contrast of the fake planets
- fake_flux_dict[“mode”]=”SNR”: Defines the brightness of the fakes relatively to the satellite spots.
fake_flux_dict[“SNR”]: SNR wished for the fake planets.
fake_flux_dict[“sep_arr”]: Separation sampling of the contrast curve in pixels.
fake_flux_dict[“contrast_arr”]: 5 sigma contrast curve (planet to star ratio).
PSF_cube – the psf of the image. A numpy array with shape (wv, y, x)
PSF_cube_wvs – the wavelegnths that correspond to the input psfs
spectrum –
spectrum name (string) or array - “host_star_spec”: The spectrum from the star or the satellite spots is directly used.
It is derived from the inverse of the calibrate_output() output.
”constant”: Use a constant spectrum np.ones(nl).
- other strings: name of the spectrum file in #pykliproot#/spectra/*/ with pykliproot the
directory in which pyklip is installed. It that case it should be a spectrum from Mark Marley or one following the same convention. Spectrum will be corrected for transmission.
ndarray: 1D array with a user defined spectrum. Spectrum will be corrected for transmission.
star_type – Spectral type of the current star. If None, Simbad is queried.
mute – If True prevent printed log outputs.
suffix – Suffix to be added at the end of the filenames.
SpT_file_csv – Filename of the table (.csv) contaning the spectral type of the stars.
real_planets_pos – list of position of real point sources in the dataset that should be avoided when injecting fakes. [(sep1,pa1),(sep2,pa2),…] with the separation in pixels and the position angle in degrees.
sep_skip_real_pl – Limit in seperation of how close a fake can be injected of a known GOI.
pa_skip_real_pl – Limit in position angle of how close a fake can be injected of a known GOI.
dn_per_contrast – array of the same size as spectrum giving the conversion between the peak flux of a planet in data number and its contrast.
star_spec – 1D array stellar spectrum sampling dataset.wvs. Otherwise uses a pickles spectrum in which the temperature in interpolated from the spectral type.
- pyklip.fakes.inject_disk(frames, centers, inputfluxes, astr_hdrs, pa, fwhm=3.5)[source]
Injects a fake disk into a dataset
- Parameters:
frames – array of (N,y,x) for N is the total number of frames
centers – array of size (N,2) of [x,y] coordiantes of the image center
intputfluxes –
array of size N of the peak flux of the fake disk in each frame OR array of 2-D models (North up East left) to inject into the data.
(Disk is assumed to be centered at center of image)
astr_hdrs – array of size N of the WCS headers
pa – position angles angle (in degrees) of disk plane
fwhm – if injecting a Gaussian disk (i.e inputfluxes is an array of floats), fwhm of Gaussian
- Returns:
saves result in input “frames” variable
- pyklip.fakes.inject_planet(frames, centers, inputflux, astr_hdrs, radius, pa, fwhm=3.5, thetas=None, stampsize=None, field_dependent_correction=None)[source]
Injects a fake planet into a dataset either using a Gaussian PSF or an input PSF
- Parameters:
frames – array of (N,y,x) for N is the total number of frames
centers – array of size (N,2) of [x,y] coordiantes of the image center
inputflux – EITHER array of size N of the peak flux of the fake planet in each frame (will inject a Gaussian PSF) OR array of size (N,psfy,psfx) of template PSFs. The brightnesses should be scaled and the PSFs should be centered at the center of each of the template images
astr_hdrs – array of size N of the WCS headers
radius – separation of the planet from the star
pa – position angle (in degrees) of planet
fwhm – fwhm (in pixels) of gaussian
thetas – ignore PA, supply own thetas (CCW angle from +x axis toward +y) array of size N
stampsize – in pixels, the width of the square stamp to inject the image into. Defaults to 3*fwhm if None
field_dependent_correction – a function of the form ``output_stamp = correction(input_stamp, input_dx, input_dy) where, input_stamp is a 2-D image of shape (y_stamp, x_stamp). input_dx and input_dy have the same shape and represent the offset of each pixel from the star (in pixel units). The function returns an output_stamp of the same shape, but corrected for any field dependent throughputs or distortions.
- Returns:
saves result in input “frames” variable
- pyklip.fakes.retrieve_planet(frames, centers, astr_hdrs, sep, pa, searchrad=7, guessfwhm=3.0, guesspeak=1, refinefit=True, thetas=None)[source]
Retrives the planet properties from a series of frames given a separation and PA
- Parameters:
frames – frames of data to retrieve planet. Can be a single 2-D image ([y,x]) for a series/cube ([N,y,x])
centers – coordiantes of the image center. Can be [2]-element lst or an array that matches array of frames [N,2]
astr_hdrs – astr_hdrs, can be a single one or an array of N of them
sep – radial distance in pixels
PA – parallactic angle in degrees
searchrad – 1/2 the length of the box used for the fit
guessfwhm – approximate fwhm to fit to
guesspeak – approximate flux
refinefit – whether or not to refine the positioning of the planet
thetas – ignore PA, supply own thetas (CCW angle from +x axis toward +y) single number or array of size N
- Returns:
(peakflux, x, y, fwhm). A single tuple if one frame passed in. Otherwise an array of tuples
- Return type:
measured
- pyklip.fakes.retrieve_planet_flux(frames, centers, astr_hdrs, sep, pa, searchrad=7, guessfwhm=3.0, guesspeak=1, refinefit=False, thetas=None)[source]
Retrives the planet flux from a series of frames given a separation and PA
- Parameters:
frames – frames of data to retrieve planet. Can be a single 2-D image ([y,x]) for a series/cube ([N,y,x])
centers – coordiantes of the image center. Can be [2]-element lst or an array that matches array of frames [N,2]
astr_hdrs – astr_hdrs, can be a single one or an array of N of them
sep – radial distance in pixels
PA – parallactic angle in degrees
searchrad – 1/2 the length of the box used for the fit
guessfwhm – approximate fwhm to fit to
guesspeak – approximate flux
refinefit – whether or not to refine the positioning of the planet
thetas – ignore PA, supply own thetas (CCW angle from +x axis toward +y) single number or array of size N
- Returns:
- either a single peak flux or an array depending on whether a single frame or multiple frames
where passed in
- Return type:
peakflux
pyklip.fitpsf module
- class pyklip.fitpsf.FMAstrometry(guess_sep, guess_pa, fitboxsize, method='mcmc')[source]
Bases:
FitPSF
Specifically for fitting astrometry of a directly imaged companion relative to its star. Extension of
pyklip.fitpsf.FitPSF
.- Parameters:
guess_sep – the guessed separation (pixels)
guess_pa – the guessed position angle (degrees)
fitboxsize – fitting box side length (pixels)
method (str) – either ‘mcmc’ or ‘maxl’ depending on framework you want. Defaults to ‘mcmc’.
- guess_sep
(initialization) guess separation for planet [pixels]
- Type:
float
- guess_pa
(initialization) guess PA for planet [degrees]
- Type:
float
- guess_RA_offset
(initialization) guess RA offset [pixels]
- Type:
float
- guess_Dec_offset
(initialization) guess Dec offset [pixels]
- Type:
float
- raw_RA_offset
(result) the raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_Dec_offset
(result) the raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_flux
(result) factor to scale the FM to match the flux of the data
- Type:
- covar_params
(result) hyperparameters for the Gaussian process
- Type:
list of
pyklip.fitpsf.ParamRange
- raw_sep
(result) the inferred raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_PA
(result) the inferred raw result from the MCMC fit for the planet’s location [degrees]
- Type:
- RA_offset
(result) the RA offset of the planet that includes all astrometric errors [pixels or mas]
- Type:
- Dec_offset
(result) the Dec offset of the planet that includes all astrometric errors [pixels or mas]
- Type:
- sep
(result) the separation of the planet that includes all astrometric errors [pixels or mas]
- Type:
- PA
(result) the PA of the planet that includes all astrometric errors [degrees]
- Type:
- fm_stamp
(fitting) The 2-D stamp of the forward model (centered at the nearest pixel to the guessed location)
- Type:
np.array
- data_stamp
(fitting) The 2-D stamp of the data (centered at the nearest pixel to the guessed location)
- Type:
np.array
- noise_map
(fitting) The 2-D stamp of the noise for each pixel the data computed assuming azimuthally similar noise
- Type:
np.array
- padding
amount of pixels on one side to pad the data/forward model stamp
- Type:
int
- sampler
an instance of the emcee EnsambleSampler. Only for Bayesian fit. See emcee docs for more details.
- Type:
emcee.EnsembleSampler
- fit_astrometry(nwalkers=100, nburn=200, nsteps=800, save_chain=True, chain_output='bka-chain.pkl', numthreads=None)[source]
Fits the PSF of the planet in either a frequentist or Bayesian way depending on initialization.
- Parameters:
nwalkers – number of walkers (mcmc-only)
nburn – numbe of samples of burn-in for each walker (mcmc-only)
nsteps – number of samples each walker takes (mcmc-only)
save_chain – if True, save the output in a pickled file (mcmc-only)
chain_output – filename to output the chain to (mcmc-only)
numthreads – number of threads to use (mcmc-only)
Returns:
- generate_data_stamp(data, data_center, data_wcs=None, noise_map=None, dr=4, exclusion_radius=10)[source]
Generate a stamp of the data_stamp ~centered on planet and also corresponding noise map
- Parameters:
data – the final collapsed data_stamp (2-D)
data_center – location of star in the data_stamp.
data_wcs – sky angles WCS object. To rotate the image properly [NOT YET IMPLMETNED] if None, data_stamp is already rotated North up East left
noise_map – if not None, noise map for each pixel in the data_stamp (2-D). if None, one will be generated assuming azimuthal noise using an annulus widthh of dr
dr – width of annulus in pixels from which the noise map will be generated
exclusion_radius – radius around the guess planet location which doens’t get factored into noise estimate
Returns:
- generate_fm_stamp(fm_image, fm_center, fm_wcs=None, extract=True, padding=5)[source]
Generates a stamp of the forward model and stores it in self.fm_stamp :param fm_image: full imgae containing the fm_stamp :param fm_center: [x,y] center of image (assuing fm_stamp is located at sep/PA) corresponding to guess_sep and guess_pa :param fm_wcs: if not None, specifies the sky angles in the image. If None, assume image is North up East left :param extract: if True, need to extract the forward model from the image. Otherwise, assume the fm_stamp is already
centered in the frame (fm_image.shape // 2)
- Parameters:
padding – number of pixels on each side in addition to the fitboxsize to extract to pad the fm_stamp (should be >= 1)
Returns:
- propogate_errs(star_center_err=None, platescale=None, platescale_err=None, pa_offset=None, pa_uncertainty=None)[source]
Propogate astrometric error. Stores results in its own fields
- Parameters:
star_center_err (float) – uncertainity of the star location (pixels)
platescale (float) – mas/pix conversion to angular coordinates
platescale_err (float) – mas/pix error on the platescale
pa_offset (float) – Offset, in the same direction as position angle, to set North up (degrees)
pa_uncertainity (float) – Error on position angle/true North calibration (Degrees)
- class pyklip.fitpsf.FitPSF(fitboxsize, method='mcmc')[source]
Bases:
object
Base class to perform astrometry on direct imaging data_stamp using GP regression. Can utilize a Bayesian framework with MCMC or a frequentist framework with least squares.
- Parameters:
fitboxsize – fitting box side length (pixels)
method (str) – either ‘mcmc’ or ‘maxl’ depending on framework you want. Defaults to ‘mcmc’.
fmt (str) – either ‘seppa’ or ‘xy’ depending on how you want to input the guess coordiantes
- guess_x
(initialization) guess x position [pixels]
- Type:
float
- guess_y
(initialization) guess y positon [pixels]
- Type:
float
- guess_flux
guess scale factor between model and data
- Type:
float
- fit_x
(result) the result from the MCMC fit for the planet’s location [pixels]
- Type:
- fit_y
(result) the result from the MCMC fit for the planet’s location [pixels]
- Type:
- fit_flux
(result) factor to scale the FM to match the flux of the data
- Type:
- covar_params
(result) hyperparameters for the Gaussian processa
- Type:
list of
pyklip.fitpsf.ParamRange
- fm_stamp
(fitting) The 2-D stamp of the forward model (centered at the nearest pixel to the guessed location)
- Type:
np.array
- data_stamp
(fitting) The stamp of the data (centered at the nearest pixel to the guessed location) (2-D unless there were NaNs in which 1-D)
- Type:
np.array
- noise_map
(fitting) The stamp of the noise for each pixel the data computed assuming azimuthally similar noise (same dim as data_stamp)
- Type:
np.array
- padding
amount of pixels on one side to pad the data/forward model stamp
- Type:
int
- sampler
an instance of the emcee EnsambleSampler. Only for Bayesian fit. See emcee docs for more details.
- Type:
emcee.EnsembleSampler
- best_fit_and_residuals(fig=None)[source]
Generate a plot of the best fit FM compared with the data_stamp and also the residuals :param fig: if not None, a matplotlib Figure object :type fig: matplotlib.Figure
- Returns:
the Figure object. If input fig is None, function will make a new one
- Return type:
fig (matplotlib.Figure)
- fit_psf(nwalkers=100, nburn=200, nsteps=800, save_chain=True, chain_output='bka-chain.pkl', numthreads=None)[source]
Fits the PSF to the data in either a frequentist or Bayesian way depending on initialization.
- Parameters:
nwalkers – number of walkers (mcmc-only)
nburn – numbe of samples of burn-in for each walker (mcmc-only)
nsteps – number of samples each walker takes (mcmc-only)
save_chain – if True, save the output in a pickled file (mcmc-only)
chain_output – filename to output the chain to (mcmc-only)
numthreads – number of threads to use (mcmc-only)
Returns:
- generate_data_stamp(data, guess_loc, noise_map, radial_noise_center=None, dr=4, exclusion_radius=10)[source]
Generate a stamp of the data_stamp ~centered on planet and also corresponding noise map :param data: the final collapsed data_stamp (2-D) :param guess_loc: guess location of where to fit the model in the data :param noise_map: if not None, noise map for each pixel (either same shape as input data, or shape of data stamp)
if None, one will be generated assuming azimuthal noise using an annulus widthh of dr. radial_noise_center MUST be defined.
- Parameters:
radial_noise_center – if we assume the noise is azimuthally symmetric and changes radially, this is the [x,y] center for it
dr – width of annulus in pixels from which the noise map will be generated
exclusion_radius – radius around the guess planet location which doens’t get factored into the radial noise estimate
Returns:
- generate_fm_stamp(fm_image, fm_pos=None, fm_wcs=None, extract=True, padding=5)[source]
Generates a stamp of the forward model and stores it in self.fm_stamp :param fm_image: full image containing the fm_stamp :param fm_pos: [x,y] location of the forwrd model in the fm_image :param fm_wcs: if not None, specifies the sky angles in the image. If None, assume image is North up East left :param extract: if True, need to extract the forward model from the image. Otherwise, assume the fm_stamp is already
centered in the frame (fm_image.shape // 2)
- Parameters:
padding – number of pixels on each side in addition to the fitboxsize to extract to pad the fm_stamp (should be >= 1)
Returns:
- property guess_flux
- make_corner_plot(fig=None)[source]
Generate a corner plot of the posteriors from the MCMC :param fig: if not None, a matplotlib Figure object
- Returns:
the Figure object. If input fig is None, function will make a new one
- Return type:
fig
- set_bounds(dx, dy, df, covar_param_bounds, read_noise_bounds=None)[source]
Set bounds on Bayesian priors. All paramters can be a 2 element tuple/list/array that specifies the lower and upper bounds x_min < x < x_max. Or a single value whose interpretation is specified below If you are passing in both lower and upper bounds, both should be in linear scale! :param dx: Distance from initial guess position in pixels. For a single value, this specifies the largest distance
form the initial guess (i.e. x_guess - dx < x < x_guess + dx)
- Parameters:
dy – Same as dx except with y
df – Flux range. If single value, specifies how many orders of 10 the flux factor can span in one direction (i.e. log_10(guess_flux) - df < log_10(guess_flux) < log_10(guess_flux) + df
covar_param_bounds – Params for covariance matrix. Like df, single value specifies how many orders of magnitude parameter can span. Otherwise, should be a list of 2-elem touples
read_noise_bounds – Param for read noise term. If single value, specifies how close to 0 it can go based on powers of 10 (i.e. log_10(-read_noise_bound) < read_noise < 1 )
Returns:
- set_kernel(covar, covar_param_guesses, covar_param_labels, include_readnoise=False, read_noise_fraction=0.01)[source]
Set the Gaussian process kernel used in our fit
- Parameters:
covar – Covariance kernel for GP regression. If string, can be “matern32” or “sqexp” or “diag” Can also be a function: cov = cov_function(x_indices, y_indices, sigmas, cov_params)
covar_param_guesses – a list of guesses on the hyperparmeteres (size of N_hyperparams). This can be an empty list for ‘diag’.
covar_param_labels – a list of strings labelling each covariance parameter
include_readnoise – if True, part of the noise is a purely diagonal term (i.e. read/photon noise)
read_noise_fraction – fraction of the total measured noise is read noise (between 0 and 1)
Returns:
- class pyklip.fitpsf.ParamRange(bestfit, err_range)[source]
Bases:
object
Stores the best fit value and uncertainities for a parameter in a neat fasion
- Parameters:
bestfit (float) – the bestfit value
err_range – either a float or a 2-element tuple (+val1, -val2) and gives the 1-sigma range
- bestfit
the bestfit value
- Type:
float
- error
the average 1-sigma error
- Type:
float
- error_2sided
[+error1, -error2] 2-element array with asymmetric errors
- Type:
np.array
- class pyklip.fitpsf.PlanetEvidence(guess_sep, guess_pa, fitboxsize, sampling_outputdir, l_only=False, fm_basename='Planet', null_basename='Null')[source]
Bases:
FMAstrometry
Specifically for nested sampling of the parameter space of a directly imaged companion relative to its star. Extension of
pyklip.fitpsf.FitPSF
.- Parameters:
guess_sep – the guessed separation (pixels)
guess_pa – the guessed position angle (degrees)
fitboxsize – fitting box side length (pixels)
fm_basename – Prefix of the foward model sampling files multinest saves in /chains/
null_basename – Prefix of the null hypothesis model sampling files multinest saves in /chains/
- guess_sep
(initialization) guess separation for planet [pixels]
- Type:
float
- guess_pa
(initialization) guess PA for planet [degrees]
- Type:
float
- guess_RA_offset
(initialization) guess RA offset [pixels]
- Type:
float
- guess_Dec_offset
(initialization) guess Dec offset [pixels]
- Type:
float
- raw_RA_offset
(result) the raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_Dec_offset
(result) the raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_flux
(result) factor to scale the FM to match the flux of the data
- Type:
- covar_params
(result) hyperparameters for the Gaussian process
- Type:
list of
pyklip.fitpsf.ParamRange
- raw_sep
(result) the inferred raw result from the MCMC fit for the planet’s location [pixels]
- Type:
- raw_PA
(result) the inferred raw result from the MCMC fit for the planet’s location [degrees]
- Type:
- RA_offset
(result) the RA offset of the planet that includes all astrometric errors [pixels or mas]
- Type:
- Dec_offset
(result) the Dec offset of the planet that includes all astrometric errors [pixels or mas]
- Type:
- sep
(result) the separation of the planet that includes all astrometric errors [pixels or mas]
- Type:
- PA
(result) the PA of the planet that includes all astrometric errors [degrees]
- Type:
- fm_stamp
(fitting) The 2-D stamp of the forward model (centered at the nearest pixel to the guessed location)
- Type:
np.array
- data_stamp
(fitting) The 2-D stamp of the data (centered at the nearest pixel to the guessed location)
- Type:
np.array
- noise_map
(fitting) The 2-D stamp of the noise for each pixel the data computed assuming azimuthally similar noise
- Type:
np.array
- padding
amount of pixels on one side to pad the data/forward model stamp
- Type:
int
- sampler
function that runs the pymultinest sampling for both hypotheses
- Type:
pymultinest.run
- pyklip.fitpsf.lnlike(fitparams, fma, cov_func, readnoise=False, negate=False)[source]
Likelihood function :param fitparams: array of params (size N). First three are [dRA,dDec,f]. Additional parameters are GP hyperparams
dRA,dDec: RA,Dec offsets from star. Also coordianaes in self.data_{RA,Dec}_offset f: flux scale factor to normalizae the flux of the data_stamp to the model
- Parameters:
fma (FMAstrometry) – a FMAstrometry object that has been fully set up to run
cov_func (function) – function that given an input [x,y] coordinate array returns the covariance matrix e.g. cov = cov_function(x_indices, y_indices, sigmas, cov_params)
readnoise (bool) – If True, the last fitparam fits for diagonal noise
negate (bool) – if True, negatives the probability (used for minimization algos)
- Returns:
log of likelihood function (minus a constant factor)
- Return type:
likeli
- pyklip.fitpsf.lnprior(fitparams, bounds, readnoise=False, negate=False)[source]
Bayesian prior
- Parameters:
fitparams – array of params (size N)
bounds – array of (N,2) with corresponding lower and upper bound of params bounds[i,0] <= fitparams[i] < bounds[i,1]
readnoise (bool) – If True, the last fitparam fits for diagonal noise
negate (bool) – if True, negatives the probability (used for minimization algos)
- Returns:
0 if inside bound ranges, -inf if outside
- Return type:
prior
- pyklip.fitpsf.lnprob(fitparams, fma, bounds, cov_func, readnoise=False, negate=False)[source]
Function to compute the relative posterior probabiltiy. Product of likelihood and prior :param fitparams: array of params (size N). First three are [dRA,dDec,f]. Additional parameters are GP hyperparams
dRA,dDec: RA,Dec offsets from star. Also coordianaes in self.data_{RA,Dec}_offset f: flux scale factor to normalizae the flux of the data_stamp to the model
- Parameters:
fma – a FMAstrometry object that has been fully set up to run
bounds – array of (N,2) with corresponding lower and upper bound of params bounds[i,0] <= fitparams[i] < bounds[i,1]
cov_func – function that given an input [x,y] coordinate array returns the covariance matrix e.g. cov = cov_function(x_indices, y_indices, sigmas, cov_params)
readnoise (bool) – If True, the last fitparam fits for diagonal noise
negate (bool) – if True, negatives the probability (used for minimization algos)
Returns:
- pyklip.fitpsf.quick_psf_fit(data, psf, x_guess, y_guess, fitboxsize)[source]
A wrapper for a quick maximum likelihood fit to a PSF to the data.
- Parameters:
data (np.array) – 2-D data frame
psf (np.array) – 2-D PSF template. This should be smaller than the size of data and larger than the fitboxsize
x_guess (float) – approximate x position of the location you are fitting the psf to
y_guess (float) – approximate y position of the location you are fitting the psf to
fitboxsize (int) – fitting region is a square. This is the lenght of one side of the square
- Returns:
x_fit, y_fit, flux_fit x_fit (float): x position y_fit (float): y position flux_fit (float): multiplicative scale factor for the psf to match the data
pyklip.fm module
- pyklip.fm.calculate_fm(delta_KL_nospec, original_KL, numbasis, sci, model_sci, inputflux=None)[source]
Calculate what the PSF looks up post-KLIP using knowledge of the input PSF, assumed spectrum of the science target, and the partially calculated KL modes (Delta Z_k^lambda in Laurent’s paper). If inputflux is None, the spectral dependence has already been folded into delta_KL_nospec (treat it as delta_KL).
Note: if inputflux is None and delta_KL_nospec has three dimensions (ie delta_KL_nospec was calculated using pertrurb_nospec() or perturb_nospec_modelsBased()) then only klipped_oversub and klipped_selfsub are returned. Besides they will have an extra first spectral dimension.
- Parameters:
delta_KL_nospec – perturbed KL modes but without the spectral info. delta_KL = spectrum x delta_Kl_nospec. Shape is (numKL, wv, pix). If inputflux is None, delta_KL_nospec = delta_KL
orignal_KL – unpertrubed KL modes (array of size [numbasis, numpix])
numbasis – array of KL mode cutoffs If numbasis is [None] the number of KL modes to be used is automatically picked based on the eigenvalues.
sci – array of size p representing the science data
model_sci – array of size p corresponding to the PSF of the science frame
input_spectrum – array of size wv with the assumed spectrum of the model
If delta_KL_nospec does NOT include a spectral dimension or if inputflux is not None: :returns:
- array of shape (b,p) showing the forward modelled PSF
Skipped if inputflux = None, and delta_KL_nospec has 3 dimensions.
klipped_oversub: array of shape (b, p) showing the effect of oversubtraction as a function of KL modes klipped_selfsub: array of shape (b, p) showing the effect of selfsubtraction as a function of KL modes Note: psf_FM = model_sci - klipped_oversub - klipped_selfsub to get the FM psf as a function of K Lmodes
(shape of b,p)
- Return type:
fm_psf
If inputflux = None and if delta_KL_nospec include a spectral dimension: :returns: Sum(<S|KL>KL) with klipped_oversub.shape = (size(numbasis),Npix)
klipped_selfsub: Sum(<N|DKL>KL) + Sum(<N|KL>DKL) with klipped_selfsub.shape = (size(numbasis),N_lambda or N_ref,N_pix)
- Return type:
klipped_oversub
- pyklip.fm.calculate_fm_singleNumbasis(delta_KL_nospec, original_KL, numbasis, sci, model_sci, inputflux=None)[source]
Same function as calculate_fm() but faster when numbasis has only one element. It doesn’t do the mutliplication with the triangular matrix.
Calculate what the PSF looks up post-KLIP using knowledge of the input PSF, assumed spectrum of the science target, and the partially calculated KL modes (Delta Z_k^lambda in Laurent’s paper). If inputflux is None, the spectral dependence has already been folded into delta_KL_nospec (treat it as delta_KL).
Note: if inputflux is None and delta_KL_nospec has three dimensions (ie delta_KL_nospec was calculated using pertrurb_nospec() or perturb_nospec_modelsBased()) then only klipped_oversub and klipped_selfsub are returned. Besides they will have an extra first spectral dimension.
- Parameters:
delta_KL_nospec – perturbed KL modes but without the spectral info. delta_KL = spectrum x delta_Kl_nospec. Shape is (numKL, wv, pix). If inputflux is None, delta_KL_nospec = delta_KL
orignal_KL – unpertrubed KL modes (array of size [numbasis, numpix])
numbasis – array of (ONE ELEMENT ONLY) KL mode cutoffs If numbasis is [None] the number of KL modes to be used is automatically picked based on the eigenvalues.
sci – array of size p representing the science data
model_sci – array of size p corresponding to the PSF of the science frame
input_spectrum – array of size wv with the assumed spectrum of the model
If delta_KL_nospec does NOT include a spectral dimension or if inputflux is not None: :returns:
- array of shape (b,p) showing the forward modelled PSF
Skipped if inputflux = None, and delta_KL_nospec has 3 dimensions.
klipped_oversub: array of shape (b, p) showing the effect of oversubtraction as a function of KL modes klipped_selfsub: array of shape (b, p) showing the effect of selfsubtraction as a function of KL modes Note: psf_FM = model_sci - klipped_oversub - klipped_selfsub to get the FM psf as a function of K Lmodes
(shape of b,p)
- Return type:
fm_psf
If inputflux = None and if delta_KL_nospec include a spectral dimension: :returns: Sum(<S|KL>KL) with klipped_oversub.shape = (size(numbasis),Npix)
klipped_selfsub: Sum(<N|DKL>KL) + Sum(<N|KL>DKL) with klipped_selfsub.shape = (size(numbasis),N_lambda or N_ref,N_pix)
- Return type:
klipped_oversub
- pyklip.fm.calculate_validity(covar_perturb, models_ref, numbasis, evals_orig, covar_orig, evecs_orig, KL_orig, delta_KL)[source]
- Calculate the validity of the perturbation based on the eigenvalues or the 2nd order term compared
to the 0th order term of the covariance matrix expansion
- Parameters:
evals_perturb – linear expansion of the perturbed covariance matrix (C_AS). Shape of N x N
models_ref – N x p array of the N models corresponding to reference images. Each model should contain spectral information
numbasis – array of KL mode cutoffs
evevs_orig – size of [N, maxKL]
- Returns:
perturbed KL modes. Shape is (numKL, wv, pix)
- Return type:
delta_KL_nospec
- pyklip.fm.find_id_nearest(array, value)[source]
Find index of the closest value in input array to input value :param array: 1D array :param value: scalar value
- Returns:
Index of the nearest value in array
- pyklip.fm.klip_dataset(dataset, fm_class, mode='ADI+SDI', outputdir='.', fileprefix='pyklipfm', annuli=5, subsections=4, OWA=None, N_pix_sector=None, movement=None, flux_overlap=0.1, PSF_FWHM=3.5, minrot=0, padding=0, numbasis=None, maxnumbasis=None, numthreads=None, corr_smooth=1, calibrate_flux=False, aligned_center=None, psf_library=None, spectrum=None, highpass=False, annuli_spacing='constant', save_klipped=True, mute_progression=False, time_collapse='mean')[source]
Run KLIP-FM on a dataset object
- Parameters:
dataset – an instance of Instrument.Data (see instruments/ subfolder)
fm_class – class that implements the the forward modelling functionality
mode – some combination of ADI, SDI, and RDI (e.g. “ADI+SDI”, “RDI”). Note that note all FM classes support RDI.
anuuli – Annuli to use for KLIP. Can be a number, or a list of 2-element tuples (a, b) specifying the pixel bondaries (a <= r < b) for each annulus
subsections – Sections to break each annuli into. Can be a number [integer], or a list of 2-element tuples (a, b) specifying the positon angle boundaries (a <= PA < b) for each section [radians]
OWA – if defined, the outer working angle for pyklip. Otherwise, it will pick it as the cloest distance to a nan in the first frame
N_pix_sector –
Rough number of pixels in a sector. Overwriting subsections and making it sepration dependent. The number of subsections is defined such that the number of pixel is just higher than N_pix_sector. I.e. subsections = floor(pi*(r_max^2-r_min^2)/N_pix_sector) Warning: There is a bug if N_pix_sector is too big for the first annulus. The annulus is defined from
0 to 2pi which create a bug later on. It is probably in the way pa_start and pa_end are defined in fm_from_eigen(). (I am taking about matched filter by the way)
movement – minimum amount of movement (in pixels) of an astrophysical source to consider using that image for a refernece PSF
flux_overlap – Maximum fraction of flux overlap between a slice and any reference frames included in the covariance matrix. Flux_overlap should be used instead of “movement” when a template spectrum is used. However if movement is not None then the old code is used and flux_overlap is ignored. The overlap is estimated for 1D gaussians with FWHM defined by PSF_FWHM. So note that the overlap is not exactly the overlap of two real 2D PSF for a given instrument but it will behave similarly.
PSF_FWHM – FWHM of the PSF used to calculate the overlap (cf flux_overlap). Default is FWHM = 3.5 corresponding to sigma ~ 1.5.
minrot – minimum PA rotation (in degrees) to be considered for use as a reference PSF (good for disks)
padding – for each sector, how many extra pixels of padding should we have around the sides.
numbasis – number of KL basis vectors to use (can be a scalar or list like). Length of b If numbasis is [None] the number of KL modes to be used is automatically picked based on the eigenvalues.
maxnumbasis – Number of KL modes to be calculated from whcih numbasis modes will be taken.
corr_smooth (float) – size of sigma of Gaussian smoothing kernel (in pixels) when computing most correlated PSFs. If 0, no smoothing
numthreads – number of threads to use. If none, defaults to using all the cores of the cpu
calibrate_flux – if true, flux calibrates the regular KLIP subtracted data. DOES NOT CALIBRATE THE FM
aligned_center – array of 2 elements [x,y] that all the KLIP subtracted images will be centered on for image registration
psf_library – a rdi.PSFLibrary object with a PSF Library for RDI
spectrum –
- (only applicable for SDI) if not None, optimizes the choice of the reference PSFs based on the
spectrum shape.
an array: of length N with the flux of the template spectrum at each wavelength.
a string: Currently only supports “methane” between 1 and 10 microns.
Uses minmove to determine the separation from the center of the segment to determine contamination and the size of the PSF (TODO: make PSF size another quantity) (e.g. minmove=3, checks how much contamination is within 3 pixels of the hypothetical source) if smaller than 10%, (hard coded quantity), then use it for reference PSF
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
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
save_klipped – if True, will save the regular klipped image. If false, it wil not and sub_imgs will return None
mute_progression – Mute the printing of the progression percentage. Indeed sometimes the overwriting feature doesn’t work and one ends up with thousands of printed lines. Therefore muting it can be a good idea.
time_collapse – how to collapse the data in time. Currently support: “mean”, “weighted-mean”
- pyklip.fm.klip_math(sci, refs, numbasis, covar_psfs=None, model_sci=None, models_ref=None, spec_included=False, spec_from_model=False)[source]
linear algebra of KLIP with linear perturbation disks and point source
- Parameters:
sci – array of length p containing the science data
refs – N x p array of the N reference images that characterizes the extended source with p pixels
numbasis – number of KLIP basis vectors to use (can be an int or an array of ints of length b) If numbasis is [None] the number of KL modes to be used is automatically picked based on the eigenvalues.
covar_psfs – covariance matrix of reference images (for large N, useful). Normalized following numpy normalization in np.cov documentation
in (# The following arguments must all be passed) –
work (or none of them for klip_math to) –
models_ref – N x p array of the N models corresponding to reference images. Each model should be normalized to unity (no flux information)
model_sci – array of size p corresponding to the PSF of the science frame
Sel_wv – wv x N array of the the corresponding wavelength for each reference PSF
input_spectrum – array of size wv with the assumed spectrum of the model
- Returns:
- array of shape (p,b) that is the PSF subtracted data for each of the b KLIP basis
cutoffs. If numbasis was an int, then sub_img_row_selected is just an array of length p
KL_basis: array of KL basis (shape of [numbasis, p]) If models_ref is passed in (not None):
delta_KL_nospec: array of shape (b, wv, p) that is the almost perturbed KL modes just missing spectral info
- Otherwise:
evals: array of eigenvalues (size of max number of KL basis requested aka nummaxKL) evecs: array of corresponding eigenvectors (shape of [p, nummaxKL])
- Return type:
sub_img_rows_selected
- pyklip.fm.klip_parallelized(imgs, centers, parangs, wvs, IWA, fm_class, OWA=None, mode='ADI+SDI', annuli=5, subsections=4, movement=None, flux_overlap=0.1, PSF_FWHM=3.5, numbasis=None, maxnumbasis=None, corr_smooth=1, aligned_center=None, numthreads=None, minrot=0, maxrot=360, spectrum=None, psf_library=None, psf_library_good=None, psf_library_corr=None, padding=0, save_klipped=True, flipx=True, N_pix_sector=None, mute_progression=False, annuli_spacing='constant', compute_noise_cube=False)[source]
multithreaded KLIP PSF Subtraction
- Parameters:
imgs – array of 2D images for ADI. Shape of array (N,y,x)
centers – N by 2 array of (x,y) coordinates of image centers
parangs – N length array detailing parallactic angle of each image
wvs – N length array of the wavelengths
IWA – inner working angle (in pixels)
fm_class – class that implements the the forward modelling functionality
OWA – if defined, the outer working angle for pyklip. Otherwise, it will pick it as the cloest distance to a nan in the first frame
mode – one of [‘ADI’, ‘SDI’, ‘ADI+SDI’] for ADI, SDI, or ADI+SDI
anuuli – Annuli to use for KLIP. Can be a number, or a list of 2-element tuples (a, b) specifying the pixel bondaries (a <= r < b) for each annulus
subsections – Sections to break each annuli into. Can be a number [integer], or a list of 2-element tuples (a, b) specifying the positon angle boundaries (a <= PA < b) for each section [radians]
N_pix_sector –
Rough number of pixels in a sector. Overwriting subsections and making it sepration dependent. The number of subsections is defined such that the number of pixel is just higher than N_pix_sector. I.e. subsections = floor(pi*(r_max^2-r_min^2)/N_pix_sector) Warning: There is a bug if N_pix_sector is too big for the first annulus. The annulus is defined from
0 to 2pi which create a bug later on. It is probably in the way pa_start and pa_end are defined in fm_from_eigen().
movement – minimum amount of movement (in pixels) of an astrophysical source to consider using that image for a refernece PSF
flux_overlap – Maximum fraction of flux overlap between a slice and any reference frames included in the covariance matrix. Flux_overlap should be used instead of “movement” when a template spectrum is used. However if movement is not None then the old code is used and flux_overlap is ignored. The overlap is estimated for 1D gaussians with FWHM defined by PSF_FWHM. So note that the overlap is not exactly the overlap of two real 2D PSF for a given instrument but it will behave similarly.
PSF_FWHM – FWHM of the PSF used to calculate the overlap (cf flux_overlap). Default is FWHM = 3.5 corresponding to sigma ~ 1.5.
numbasis – number of KL basis vectors to use (can be a scalar or list like). Length of b If numbasis is [None] the number of KL modes to be used is automatically picked based on the eigenvalues.
maxnumbasis – Number of KL modes to be calculated from whcih numbasis modes will be taken.
corr_smooth (float) – size of sigma of Gaussian smoothing kernel (in pixels) when computing most correlated PSFs. If 0, no smoothing
aligned_center – array of 2 elements [x,y] that all the KLIP subtracted images will be centered on for image registration
numthreads – number of threads to use. If none, defaults to using all the cores of the cpu
minrot – minimum PA rotation (in degrees) to be considered for use as a reference PSF (good for disks)
maxrot – maximum PA rotation (in degrees) to be considered for use as a reference PSF (temporal variability)
spectrum – if not None, a array of length N with the flux of the template spectrum at each wavelength. Uses minmove to determine the separation from the center of the segment to determine contamination and the size of the PSF (TODO: make PSF size another quanitity) (e.g. minmove=3, checks how much containmination is within 3 pixels of the hypothetical source) if smaller than 10%, (hard coded quantity), then use it for reference PSF
padding – for each sector, how many extra pixels of padding should we have around the sides.
save_klipped – if True, will save the regular klipped image. If false, it wil not and sub_imgs will return None
flipx – if True, flips x axis after rotation to get North up East left
mute_progression – Mute the printing of the progression percentage. Indeed sometimes the overwriting feature doesn’t work and one ends up with thousands of printed lines. Therefore muting it can be a good idea.
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
compute_noise_cube – if True, compute the noise in each pixel assuming azimuthally uniform noise
- Returns:
- array of [array of 2D images (PSF subtracted)] using different number of KL basis vectors as
specified by numbasis. Shape of (b,N,y,x).
Note: this will be None if save_klipped is False
fmout_np: output of forward modelling. perturbmag: output indicating the magnitude of the linear perturbation to assess validity of KLIP FM aligned_center: (x, y) location indicating the star center for all images and FM after PSF subtraction
- Return type:
sub_imgs
- pyklip.fm.pertrurb_nospec(evals, evecs, original_KL, refs, models_ref)[source]
Perturb the KL modes using a model of the PSF but with no assumption on the spectrum. Useful for planets.
By no assumption on the spectrum it means that the spectrum has been factored out of Delta_KL following equation (4) of Laurent Pueyo 2016 noted bold “Delta Z_k^lambda (x)”. In order to get the actual perturbed KL modes one needs to multpily it by a spectrum.
This function fits each cube’s spectrum independently. So the effective spectrum size is N_wavelengths * N_cubes.
- Parameters:
evals – array of eigenvalues of the reference PSF covariance matrix (array of size numbasis)
evecs – corresponding eigenvectors (array of size [p, numbasis])
orignal_KL – unpertrubed KL modes (array of size [numbasis, p])
Sel_wv – wv x N array of the the corresponding wavelength for each reference PSF
refs – N x p array of the N reference images that characterizes the extended source with p pixels
models_ref – N x p array of the N models corresponding to reference images. Each model should be normalized to unity (no flux information)
model_sci – array of size p corresponding to the PSF of the science frame
- Returns:
- perturbed KL modes but without the spectral info. delta_KL = spectrum x delta_Kl_nospec.
Shape is (numKL, wv, pix)
- Return type:
delta_KL_nospec
- pyklip.fm.perturb_nospec_modelsBased(evals, evecs, original_KL, refs, models_ref_list)[source]
Perturb the KL modes using a model of the PSF but with no assumption on the spectrum. Useful for planets.
By no assumption on the spectrum it means that the spectrum has been factored out of Delta_KL following equation (4) of Laurent Pueyo 2016 noted bold “Delta Z_k^lambda (x)”. In order to get the actual perturbed KL modes one needs to multpily it by a spectrum.
Effectively does the same thing as pertrurb_nospec() but in a different way. It injects models with dirac spectrum (all but one vanishing wavelength) and because of the linearity of the problem allow one de get reconstruct the perturbed KL mode for any spectrum. The difference however in the pertrurb_nospec() case is that the spectrum here is the asummed to be the same for all
cubes while pertrurb_nospec() fit each cube independently.
- Parameters:
evals –
evecs –
original_KL –
refs –
models_ref –
- Returns:
delta_KL_nospec
- pyklip.fm.perturb_specIncluded(evals, evecs, original_KL, refs, models_ref, return_perturb_covar=False)[source]
Perturb the KL modes using a model of the PSF but with the spectrum included in the model. Quicker than the others
- Parameters:
evals – array of eigenvalues of the reference PSF covariance matrix (array of size numbasis)
evecs – corresponding eigenvectors (array of size [p, numbasis])
orignal_KL – unpertrubed KL modes (array of size [numbasis, p])
refs – N x p array of the N reference images that characterizes the extended source with p pixels
models_ref – N x p array of the N models corresponding to reference images. Each model should contain spectral informatoin
model_sci – array of size p corresponding to the PSF of the science frame
- Returns:
perturbed KL modes. Shape is (numKL, wv, pix)
- Return type:
delta_KL_nospec
pyklip.klip module
- pyklip.klip.align_and_scale(img, new_center, old_center=None, scale_factor=1, dtype=<class 'float'>)[source]
Helper function that realigns and/or scales the image
- Parameters:
img – 2D image to perform manipulation on
new_center – 2 element tuple (xpos, ypos) of new image center
old_center – 2 element tuple (xpos, ypos) of old image center
scale_factor –
how much the stretch/contract the image. Will we scaled w.r.t the new_center (done after relaignment). We will adopt the convention
>1: stretch image (shorter to longer wavelengths) <1: contract the image (longer to shorter wvs) This means scale factor should be lambda_0/lambda where lambda_0 is the wavelength you want to scale to
- Returns:
shifted and/or scaled 2D image
- Return type:
resampled_img
- pyklip.klip.calc_scaling(sats, refwv=18)[source]
Helper function that calculates the wavelength scaling factor from the satellite spot locations. Uses the movement of spots diagonally across from each other, to calculate the scaling in a (hopefully? tbd.) centering-independent way. This method is definitely temporary and will be replaced by better scaling strategies as we come up with them. Scaling is calculated as the average of (1/2 * sqrt((x_1-x_2)**2+(y_1-y_2))), over the two pairs of spots.
- Parameters:
sats – [4 x Nlambda x 2] array of x and y positions for the 4 satellite spots
refwv – reference wavelength for scaling (optional, default = 20)
- Returns:
Nlambda array of scaling factors
- Return type:
scaling_factors
- pyklip.klip.collapse_data(data, pixel_weights=None, axis=1, collapse_method='mean')[source]
Function to collapse multi-dimensional data along axis using collapse_method
- Parameters:
data – (multi-dimension)arrays of 2D images or 3D cubes.
pixel_weights – ones if collapse method is not weighted collapse
axis – axis index along which to collapse
collapse_method – currently support ‘median’, ‘mean’, ‘weighted-mean’, ‘trimmed-mean’, ‘weighted-median’
- Returns:
Collapsed data
- pyklip.klip.define_annuli_bounds(annuli, IWA, OWA, annuli_spacing='constant')[source]
Defines the annuli boundaries radially
- Parameters:
annuli – number of annuli
IWA – inner working angle (pixels)
OWA – outer working anglue (pixels)
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
- Returns:
array of 2-element tuples that specify the beginning and end radius of that annulus
- Return type:
rad_bounds
- pyklip.klip.estimate_movement(radius, parang0=None, parangs=None, wavelength0=None, wavelengths=None, mode=None)[source]
Estimates the movement of a hypothetical astrophysical source in ADI and/or SDI at the given radius and given reference parallactic angle (parang0) and reference wavelegnth (wavelength0)
- Parameters:
radius – the radius from the star of the hypothetical astrophysical source
parang0 – the parallactic angle of the reference image (in degrees)
parangs – array of length N of the parallactic angle of all N images (in degrees)
wavelength0 – the wavelength of the reference image
wavelengths – array of length N of the wavelengths of all N images
NOTE – we expect parang0 and parangs to be either both defined or both None. Same with wavelength0 and wavelengths
mode – one of [‘ADI’, ‘SDI’, ‘ADI+SDI’] for ADI, SDI, or ADI+SDI
- Returns:
- array of length N of the distance an astrophysical source would have moved from the
reference image
- Return type:
moves
- pyklip.klip.high_pass_filter(img, filtersize=10)[source]
A FFT implmentation of high pass filter.
- Parameters:
img – a 2D image
filtersize – size in Fourier space of the size of the space. In image space, size=img_size/filtersize
- Returns:
the filtered image
- Return type:
filtered
- pyklip.klip.klip_math(sci, ref_psfs, numbasis, covar_psfs=None, return_basis=False, return_basis_and_eig=False)[source]
Helper function for KLIP that does the linear algebra
- Parameters:
sci – array of length p containing the science data
ref_psfs – N x p array of the N reference PSFs that characterizes the PSF of the p pixels
numbasis – number of KLIP basis vectors to use (can be an int or an array of ints of length b)
covar_psfs – covariance matrix of reference psfs passed in so you don’t have to calculate it here
return_basis – If true, return KL basis vectors (used when onesegment==True)
return_basis_and_eig – If true, return KL basis vectors as well as the eigenvalues and eigenvectors of the covariance matrix. Used for KLIP Forward Modelling of Laurent Pueyo.
- Returns:
- array of shape (p,b) that is the PSF subtracted data for each of the b KLIP basis
cutoffs. If numbasis was an int, then sub_img_row_selected is just an array of length p
KL_basis: array of shape (max(numbasis),p). Only if return_basis or return_basis_and_eig is True. evals: Eigenvalues of the covariance matrix. The covariance matrix is assumed NOT to be normalized by (p-1).
Only if return_basis_and_eig is True.
- evecs: Eigenvectors of the covariance matrix. The covariance matrix is assumed NOT to be normalized by (p-1).
Only if return_basis_and_eig is True.
- Return type:
sub_img_rows_selected
- pyklip.klip.make_polar_coordinates(x, y, center=[0, 0])[source]
- Parameters:
x – meshgrid of x coordinates
y – meshgrid of y coordinates
center – new location of origin
- Returns:
polar coordinates centered at the specified origin
- pyklip.klip.meas_contrast(dat, iwa, owa, resolution, center=None, low_pass_filter=True)[source]
Measures the contrast in the image. Image must already be in contrast units and should be corrected for algorithm thoughput.
- Parameters:
dat – 2D image - already flux calibrated
iwa – inner working angle
owa – outer working angle
resolution – size of noise resolution element in pixels (for speckle noise ~ FWHM or lambda/D) but it can be 1 pixel if limited by pixel-to-pixel noise.
center – location of star (x,y). If None, defaults the image size // 2.
low_pass_filter – if True, run a low pass filter. Can also be a float which specifices the width of the Gaussian filter (sigma). If False, no Gaussian filter is run
- Returns:
tuple of separations in pixels and corresponding 5 sigma FPF
- Return type:
(seps, contrast)
- pyklip.klip.nan_gaussian_filter(img, sigma, ivar=None)[source]
Gaussian low-pass filter that handles nans
- Parameters:
img – 2-D image
sigma – float specifiying width of Gaussian
ivar – inverse variance frame for the image, optional
- Returns:
2-D image that has been smoothed with a Gaussian
- Return type:
filtered
- pyklip.klip.nan_map_coordinates_2d(img, yp, xp, mc_kwargs=None)[source]
scipy.ndimage.map_coordinates() that handles nans for 2-D transformations. Only works in 2-D!
Do NaN detection by defining any pixel in the new coordiante system (xp, yp) as a nan If any one of the neighboring pixels in the original image is a nan (e.g. (xp, yp) = (120.1, 200.1) is nan if either (120, 200), (121, 200), (120, 201), (121, 201) is a nan)
- Parameters:
img (np.array) – 2-D image that is looking to be transformed
yp (np.array) – 2-D array of y-coordinates that the image is evaluated out
xp (np.array) – 2-D array of x-coordinates that the image is evaluated out
mc_kwargs (dict) – other parameters to pass into the map_coordinates function.
- Returns:
2-D transformed image. Each pixel is evaluated at the (yp, xp) specified by xp and yp.
- Return type:
transformed_img (np.array)
- pyklip.klip.rotate(img, angle, center, new_center=None, flipx=False, astr_hdr=None)[source]
Rotate an image by the given angle about the given center. Optional: can shift the image to a new image center after rotation. Also can reverse x axis for those left
handed astronomy coordinate systems
- Parameters:
img – a 2D image
angle – angle CCW to rotate by (degrees)
center – 2 element list [x,y] that defines the center to rotate the image to respect to
new_center – 2 element list [x,y] that defines the new image center after rotation
flipx – reverses x axis after rotation
astr_hdr – wcs astrometry header for the image
- Returns:
new 2D image
- Return type:
resampled_img
pyklip.nmf_imaging module
- pyklip.nmf_imaging.NMFbff(trg, model, fracs=None)[source]
BFF subtraction. :param trg: :param model: :param fracs: (if need to be).
- Returns:
best frac
- pyklip.nmf_imaging.NMFcomponents(ref, ref_err=None, n_components=None, maxiters=1000.0, oneByOne=False, ignore_mask=None, path_save=None, recalculate=False)[source]
Returns the NMF components, where the rows contain the information. :param ref and ref_err should be: :type ref and ref_err should be: N, p :param ignore_mask: array of shape (N, p). mask pixels in each image that you don’t want to use. :param path_save: string, path to save the NMF components (at: path_save + ‘_comp.fits’) and coeffieients (at: path_save + ‘_coef.fits’) :param recalculate: boolean, whether to recalculate when path_save is provided
- Returns:
NMf components (n_components * p).
- pyklip.nmf_imaging.NMFmodelling(trg, components, n_components=None, mask_components=None, mask_data_imputation=None, trg_err=None, maxiters=1000.0, cube=False, trgThresh=0)[source]
NMF modeling. :param trg: 1D array, p pixels :param components: N * p, calculated using NMFcomponents. :param n_components: how many components do you want to use. If None, all the components will be used. :param cube: whether output a cube or not (increasing the number of components). :param trgThresh: ignore the regions with low photon counts. Especially when they are ~10^-15 or smaller. I chose 0 in this case.
- Returns:
NMF model of the target.
- pyklip.nmf_imaging.NMFsubtraction(trg, model, frac=1)[source]
NMF subtraction with a correction factor, frac.
- pyklip.nmf_imaging.data_masked_only(data, mask=None)[source]
Return the data where the same regions are ignored in all the data :param data: (p, N) where N is the number of references, p is the number of pixels in each reference :param mask: 1d array containing only 1 and 0, with 0 will be ignored in the output
- Returns:
(p_focused, N) where p_focused is the number of 1’s in mask
- Return type:
data_focused
- pyklip.nmf_imaging.data_masked_only_revert(data, mask=None)[source]
Return the data where the same regions were ignored in all the data :param data: (p_focused, N) where N is the number of references, p_focused is the number of 1’s in mask :param mask: 1d array containing only 1 and 0, with 0 was previously ignored in the input
- Returns:
(p, N) where p is the number of pixels in each reference
- Return type:
data_focused_revert
- pyklip.nmf_imaging.nmf_math(sci, ref_psfs, sci_err=None, ref_psfs_err=None, componentNum=5, maxiters=100000.0, oneByOne=True, trg_type='disk', ignore_mask=None, path_save=None, recalculate=False, mask_data_imputation=None)[source]
Main NMF function for high contrast imaging. :param trg: target image, dimension: height * width. :type trg: 1D array :param refs: reference cube, dimension: referenceNumber * height * width. :type refs: 2D array :param trg_err: uncertainty for trg and refs, repectively. If None is given, the squareroot of the two arrays will be adopted. :param ref_err: uncertainty for trg and refs, repectively. If None is given, the squareroot of the two arrays will be adopted. :param componentNum: number of components to be used. Default: 5. Caution: choosing too many components will slow down the computation. :type componentNum: integer :param maxiters: number of iterations needed. Default: 10^5. :type maxiters: integer :param oneByOne: whether to construct the NMF components one by one. Default: True. :type oneByOne: boolean :param trg_type (string: “disk” or “d” for circumsetllar disks by Bin Ren, the user can use “planet” or “p” for planets): are we aiming at finding circumstellar disks or planets? :param default: “disk” or “d” for circumsetllar disks by Bin Ren, the user can use “planet” or “p” for planets): are we aiming at finding circumstellar disks or planets?
- Returns:
NMF modeling result. Only the final subtraction result is returned.
- Return type:
result (1D array)
pyklip.parallelized module
- pyklip.parallelized.generate_noise_maps(imgs, aligned_center, dr, IWA=None, OWA=None, numthreads=None, pool=None)[source]
Create a noise map for each image. The noise levels are computed using azimuthally averaged noise in the images
- Parameters:
imgs – array of shape (N,y,x) containing N images
aligned_center – [x,y] location of the center. All images should be aligned to common center
dr (float) – how mnay pixels wide the annulus to compute the noise should be
IWA (float) – inner working angle (how close to the center of the image to start). If none, this is 0
OWA (float) – outer working angle (if None, it is the entire image.)
numthreads – number of threads to be used
pool – multiprocessing thread pool (optional). To avoid repeatedly creating one when processing a list of images.
- Returns:
array of shape (N,y,x) containing N noise maps
- Return type:
noise_maps
- pyklip.parallelized.high_pass_filter_imgs(imgs, numthreads=None, filtersize=10, pool=None)[source]
filters a sequences of images using a FFT
- Parameters:
imgs – array of shape (N,y,x) containing N images
numthreads – number of threads to be used
filtersize – size in Fourier space of the size of the space. In image space, size=img_size/filtersize
pool – multiprocessing thread pool (optional). To avoid repeatedly creating one when processing a list of images.
- Returns:
array of shape (N,y,x) containing the filtered images
- Return type:
filtered
- pyklip.parallelized.klip_dataset(dataset, mode='ADI+SDI', outputdir='.', fileprefix='', annuli=5, subsections=4, movement=3, numbasis=None, numthreads=None, minrot=0, calibrate_flux=False, aligned_center=None, annuli_spacing='constant', maxnumbasis=None, corr_smooth=1, spectrum=None, psf_library=None, highpass=False, lite=False, save_aligned=False, restored_aligned=None, dtype=None, algo='klip', skip_derot=False, time_collapse='mean', wv_collapse='mean', verbose=True)[source]
run klip on a dataset class outputted by an implementation of Instrument.Data
- Parameters:
dataset – an instance of Instrument.Data (see instruments/ subfolder)
mode – some combination of ADI, SDI, and RDI (e.g. “ADI+SDI”, “RDI”)
outputdir – directory to save output files
fileprefix – filename prefix for saved files
anuuli – number of annuli to use for KLIP
subsections – number of sections to break each annuli into
movement – minimum amount of movement (in pixels) of an astrophysical source to consider using that image for a refernece PSF
numbasis – number of KL basis vectors to use (can be a scalar or list like). Length of b
numthreads – number of threads to use. If none, defaults to using all the cores of the cpu
minrot – minimum PA rotation (in degrees) to be considered for use as a reference PSF (good for disks)
calibrate_flux – if True calibrate flux of the dataset, otherwise leave it be
aligned_center – array of 2 elements [x,y] that all the KLIP subtracted images will be centered on for image registration
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
maxnumbasis – if not None, maximum number of KL basis/correlated PSFs to use for KLIP. Otherwise, use max(numbasis)
corr_smooth (float) – size of sigma of Gaussian smoothing kernel (in pixels) when computing most correlated PSFs. If 0, no smoothing
spectrum – (only applicable for SDI) if not None, optimizes the choice of the reference PSFs based on the spectrum shape. - an array: of length N with the flux of the template spectrum at each wavelength. - a string: Currently only supports “methane” between 1 and 10 microns. Uses minmove to determine the separation from the center of the segment to determine contamination and the size of the PSF (TODO: make PSF size another quantity) (e.g. minmove=3, checks how much contamination is within 3 pixels of the hypothetical source) if smaller than 10%, (hard coded quantity), then use it for reference PSF
psf_library – if not None, a rdi.PSFLibrary object with a PSF Library for RDI
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
lite – if True, run a low memory version of the alogirhtm
save_aligned – Save the aligned and scaled images (as well as various wcs information), True/False
restore_aligned – The aligned and scaled images from a previous run of klip_dataset (usually restored_aligned = dataset.aligned_and_scaled)
dtype – data type of the arrays. Should be either ctypes.c_float(default) or ctypes.c_double
algo (str) – algorithm to use (‘klip’, ‘nmf’, ‘empca’, ‘none’). None will run no PSF subtraction.
skip_derot – if True, skips derotating the images. Note, the saved time-collapsed cubes may not make sense
time_collapse – how to collapse the data in time. Currently support: “mean”, “weighted-mean”, ‘median’, “weighted-median”
wv_collapse – how to collapse the data in wavelength. Currently support: ‘median’, ‘mean’, ‘trimmed-mean’
verbose (bool) – if True, print warning messages during KLIP process.
- Returns
Saved files in the output directory Returns: nothing, but saves to dataset.output: (b, N, wv, y, x) 5D cube of KL cutoff modes (b), number of images
(N), wavelengths (wv), and spatial dimensions. Images are derotated. For ADI only, the wv is omitted so only 4D cube
- pyklip.parallelized.klip_parallelized(imgs, centers, parangs, wvs, filenums, IWA, OWA=None, mode='ADI+SDI', annuli=5, subsections=4, movement=3, numbasis=None, aligned_center=None, numthreads=None, minrot=0, maxrot=360, annuli_spacing='constant', maxnumbasis=None, corr_smooth=1, spectrum=None, psf_library=None, psf_library_good=None, psf_library_corr=None, save_aligned=False, restored_aligned=None, dtype=None, algo='klip', compute_noise_cube=False, verbose=True)[source]
Multitprocessed KLIP PSF Subtraction
- Parameters:
imgs – array of 2D images for ADI. Shape of array (N,y,x)
centers – N by 2 array of (x,y) coordinates of image centers
parangs – N length array detailing parallactic angle of each image
wvs – N length array of the wavelengths
filenums (np.array) – array of length N of the filenumber for each image
IWA – inner working angle (in pixels)
OWA – outer working angle (in pixels)
mode – one of [‘ADI’, ‘SDI’, ‘ADI+SDI’] for ADI, SDI, or ADI+SDI
anuuli – number of annuli to use for KLIP
subsections – number of sections to break each annuli into
movement – minimum amount of movement (in pixels) of an astrophysical source to consider using that image for a refernece PSF
numbasis – number of KL basis vectors to use (can be a scalar or list like). Length of b
aligned_center – array of 2 elements [x,y] that all the KLIP subtracted images will be centered on for image registration
numthreads – number of threads to use. If none, defaults to using all the cores of the cpu
minrot – minimum PA rotation (in degrees) to be considered for use as a reference PSF (good for disks)
maxrot – maximum PA rotation (in degrees) to be considered for use as a reference PSF (temporal variability)
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
maxnumbasis – if not None, maximum number of KL basis/correlated PSFs to use for KLIP. Otherwise, use max(numbasis)
corr_smooth (float) – size of sigma of Gaussian smoothing kernel (in pixels) when computing most correlated PSFs. If 0, no smoothing
spectrum – if not None, a array of length N with the flux of the template spectrum at each wavelength. Uses minmove to determine the separation from the center of the segment to determine contamination and the size of the PSF (TODO: make PSF size another quanitity) (e.g. minmove=3, checks how much containmination is within 3 pixels of the hypothetical source) if smaller than 10%, (hard coded quantity), then use it for reference PSF
psf_library – array of (N_lib, y, x) with N_lib PSF library PSFs
psf_library_good – array of size N_lib indicating which N_good are good are selected in the passed in corr matrix
psf_library_corr – matrix of size N_sci x N_good with correlation between the target franes and the good RDI PSFs
save_aligned – Save the aligned and scaled images (as well as various wcs information), True/False
restore_aligned – The aligned and scaled images from a previous run of klip_dataset (usually restored_aligned = dataset.aligned_and_scaled)
dtype – data type of the arrays. Should be either ctypes.c_float(default) or ctypes.c_double
algo (str) – algorithm to use (‘klip’, ‘nmf’, ‘empca’)
compute_noise_cube – if True, compute the noise in each pixel assuming azimuthally uniform noise
- Returns:
- array of [array of 2D images (PSF subtracted)] using different number of KL basis vectors as
specified by numbasis. Shape of (b,N,y,x).
aligned_center: (x,y) specifying the common center the output images are aligned to
- Return type:
sub_imgs
- pyklip.parallelized.klip_parallelized_lite(imgs, centers, parangs, wvs, filenums, IWA, OWA=None, mode='ADI+SDI', annuli=5, subsections=4, movement=3, numbasis=None, aligned_center=None, numthreads=None, minrot=0, maxrot=360, annuli_spacing='constant', maxnumbasis=None, corr_smooth=1, spectrum=None, dtype=None, algo='klip', compute_noise_cube=False, **kwargs)[source]
multithreaded KLIP PSF Subtraction, has a smaller memory foot print than the original
- Parameters:
imgs – array of 2D images for ADI. Shape of array (N,y,x)
centers – N by 2 array of (x,y) coordinates of image centers
parangs – N length array detailing parallactic angle of each image
wvs – N length array of the wavelengths
filenums (np.array) – array of length N of the filenumber for each image
IWA – inner working angle (in pixels)
OWA – outer working angle (in pixels)
mode – one of [‘ADI’, ‘SDI’, ‘ADI+SDI’] for ADI, SDI, or ADI+SDI
anuuli – number of annuli to use for KLIP
subsections – number of sections to break each annuli into
movement – minimum amount of movement (in pixels) of an astrophysical source to consider using that image for a refernece PSF
numbasis – number of KL basis vectors to use (can be a scalar or list like). Length of b
annuli_spacing – how to distribute the annuli radially. Currently three options. Constant (equally spaced), log (logarithmical expansion with r), and linear (linearly expansion with r)
maxnumbasis – if not None, maximum number of KL basis/correlated PSFs to use for KLIP. Otherwise, use max(numbasis)
corr_smooth (float) – size of sigma of Gaussian smoothing kernel (in pixels) when computing most correlated PSFs. If 0, no smoothing
aligned_center – array of 2 elements [x,y] that all the KLIP subtracted images will be centered on for image registration
numthreads – number of threads to use. If none, defaults to using all the cores of the cpu
minrot – minimum PA rotation (in degrees) to be considered for use as a reference PSF (good for disks)
maxrot – maximum PA rotation (in degrees) to be considered for use as a reference PSF (temporal variability)
spectrum – if not None, a array of length N with the flux of the template spectrum at each wavelength. Uses minmove to determine the separation from the center of the segment to determine contamination and the size of the PSF (TODO: make PSF size another quanitity) (e.g. minmove=3, checks how much containmination is within 3 pixels of the hypothetical source) if smaller than 10%, (hard coded quantity), then use it for reference PSF
kwargs – in case you pass it stuff that we don’t use in the lite version
dtype – data type of the arrays. Should be either ctypes.c_float (default) or ctypes.c_double
algo (str) – algorithm to use (‘klip’, ‘nmf’, ‘empca’)
compute_noise_cube – if True, compute the noise in each pixel assuming azimuthally uniform noise
- Returns:
- array of [array of 2D images (PSF subtracted)] using different number of KL basis vectors as
specified by numbasis. Shape of (b,N,y,x).
- Return type:
sub_imgs
- pyklip.parallelized.rotate_imgs(imgs, angles, centers, new_center=None, numthreads=None, flipx=False, hdrs=None, disable_wcs_rotation=False, pool=None)[source]
derotate a sequences of images by their respective angles
- Parameters:
imgs – array of shape (N,y,x) containing N images
angles – array of length N with the angle to rotate each frame. Each angle should be CCW in degrees.
centers – array of shape N,2 with the [x,y] center of each frame
new_centers – a 2-element array with the new center to register each frame. Default is middle of image
numthreads – number of threads to be used
flipx – flip the x axis after rotation if desired
hdrs – array of N wcs astrometry headers
- Returns:
array of shape (N,y,x) containing the derotated images
- Return type:
derotated
pyklip.rdi module
- class pyklip.rdi.PSFLibrary(data, aligned_center, filenames, correlation_matrix=None, wvs=None, compute_correlation=False)[source]
Bases:
object
This is an PSF Library to use for reference differential imaging
- master_library
aligned library of PSFs. 3-D cube of dim = [N, y, x]. Where N is ALL files
- Type:
np.ndarray
- aligned_center
a (x,y) coordinate specifying common center the library is aligned to
- Type:
array-like
- master_filenames
array of N filenames for each frame in the library. Should match with pyklip.instruments.Data.filenames for cross-matching
- Type:
np.ndarray
- master_correlation
N x N array of correlations between each 2 frames
- Type:
np.ndarray
- master_wvs
N wavelengths for each frame
- Type:
np.ndarray
- nfiles
the number of files in the PSF library
- Type:
int
- dataset
- correlation
N_data x M array of correlations between each 2 frames where M are the selected frames and N_data is the number of files in the dataset. Along the N_data dimension, files are ordered in the same way as the dataset object
- Type:
np.ndarray
- isgoodpsf
array of N indicating which M PSFs are good for this dataset
- Type:
np.ndarray
- add_new_dataset_to_library(dataset, collapse=False, verbose=False)[source]
Add all the files from a new dataset to the PSF library and add them to the correlation matrix. If a mask was used for the correlation matrix, use it here too.
NOTE: This routine already assumes that the data has been centered.
- Parameters:
dataset (pyklip.instruments.Instrument.Data) –
- prepare_library(dataset, badfiles=None)[source]
Prepare the PSF Library for an RDI reduction of a specific dataset by only taking the part of the library we need.
- Parameters:
dataset (pyklip.instruments.Instrument.Data) –
badfiles (np.ndarray) – a list of filenames corresponding to bad files we want to also exclude
Returns:
- save_correlation(filename, overwrite=False, clobber=None, format='fits')[source]
Saves self.correlation to a file specified by filename :param filename: filepath to store the correlation matrix :type filename: str :param overwrite: if true overwrite the previous correlation matrix :type overwrite: bool :param clobber: same as overwrite, but deprecated in astropy. :type clobber: bool :param format: type of file to store the correlation matrix as. Supports numpy?/fits?/pickle? (TBD) :type format: str
pyklip.spectra_management module
- pyklip.spectra_management.calibrate_star_spectrum(template_spectrum, template_wvs, filter_name, magnitude, wvs)[source]
Scale the Pickles stellar spectrum of a star to the observed apparent magnitude and returns the stellar spectrum in physical units sampled at the specified wavelengths. Currently only supports 2MASS filters and magnitudes. TODO: implement a way to take magnitude error input and propogate the error to the final spectrum
- Parameters:
template_spectrum – 1D array, model spectrum of the star with arbitrary units
template_wvs – 1D array, wavelengths at which the template_spectrum is smapled, in units of microns
filter_name – string, 2MASS filter, ‘J’, ‘H’, or ‘Ks’
magnitude – scalar, observed apparent magnitude of the star
mag_error – scalar or 1D array with 2 elements, error(s) of the magnitude, not yet implemented
wvs – 1D array, the wvs at which to sample the scaled spectrum, in units of angstroms
- Returns:
1D array, scaled stellar spectrum sampled at wvs
- Return type:
scaled_spectrum
- pyklip.spectra_management.extract_planet_spectrum(cube_para, position, PSF_cube_para, method=None, filter=None, mute=True)[source]
- pyklip.spectra_management.find_lower_nearest(array, value)[source]
Find the lower nearest element to value in array.
- Parameters:
array – Array of value
value – Value for which one wants the lower value.
- Returns:
(low_value, id) with low_value the closest lower value and id its index.
- pyklip.spectra_management.find_nearest(array, value)[source]
Find the nearest element to value in array.
- Parameters:
array – Array of value
value – Value for which one wants the closest value.
- Returns:
(closest_value, id) with closest_value the closest lower value and id its index.
- pyklip.spectra_management.find_upper_nearest(array, value)[source]
Find the upper nearest element to value in array.
- Parameters:
array – Array of value
value – Value for which one wants the upper value.
- Returns:
(up_value, id) with up_value the closest upper value and id its index.
- pyklip.spectra_management.get_planet_spectrum(spectrum, wavelength, ori_wvs=None)[source]
Get the normalized spectrum of a planet for a GPI spectral band or any wavelengths array. Spectra are extraced from .flx files from Mark Marley et al’s models.
- Parameters:
spectrum – Path of the .flx file containing the spectrum.
wavelength – array of wavelenths in microns (or string with GPI band ‘H’, ‘J’, ‘K1’, ‘K2’, ‘Y’). (When using GPI spectral band wavelength samples are linearly spaced between the first and the last wavelength of the band.)
- Returns:
is the gpi sampling of the considered band in micrometer. spectrum: is the spectrum of the planet for the given band or wavelength array and normalized to unit mean.
- Return type:
wavelengths
- pyklip.spectra_management.get_specType(object_name, SpT_file_csv=None)[source]
Return the spectral type for a target based on Simbad or on the table in SpT_file
- Parameters:
object_name – Name of the target: ie “c_Eri”
SpT_file – Filename (.csv) of the table containing the target names and their spectral type. Can be generated from quering Simbad. If None (default), the function directly tries to query Simbad.
- Returns:
Spectral type
- pyklip.spectra_management.get_star_spectrum(wvs_or_filter_name=None, star_type=None, temperature=None, mute=None)[source]
Get the spectrum of a star with given spectral type interpolating the pickles database. The spectrum is normalized to unit mean. It assumes type V star.
- Inputs:
- wvs_or_filter_name: array of wavelenths in microns (or string with GPI band ‘H’, ‘J’, ‘K1’, ‘K2’, ‘Y’).
(When using GPI spectral band wavelength samples are linearly spaced between the first and the last wavelength of the band.)
- star_type: ‘A5’,’F4’,… Is ignored if temperature is defined.
If star_type is longer than 2 characters it is truncated.
temperature: temperature of the star. Overwrite star_type if defined.
- Output:
- (wavelengths, spectrum) where
wavelengths: Sampling in mum. spectrum: is the spectrum of the star for the given band.