Various spectroscopy routines

Contents:

Suite to reduce spectroscopic data.

subfunctions:
calibrate
setheaders – exptime, gain, readnoise, etc. makeflat – make median flat and noisy pixel map makedark – make median dark, and estimate noise in each pixel.

clean – clean and replace bad pixels extract

trace – trace spectral orders makeprofile – compute mean spectral PSF (a spline) for an order fitprofile – fit given spline-PSF to a spectral cross-section
Utilities:
pickloc fitPSF
spec.add_partial_pixel(x0, y0, z0, z)[source]
INPUTS:
x0 : int or sequence of ints

first index of z at which data will be added

y0 : int or sequence of ints

second index of z at which data will be added

z0 : scalar or sequence

values which will be added to z

z : 2D NumPy array

initial data, to which partial-pixels will be added

spec.atmosdisp(wave, wave_0, za, pressure, temp, water=2.0, fco2=0.0004, obsalt=0.0)[source]
NAME:

atmosdisp

PURPOSE:

Compute the atmosperic dispersion relative to lambda_0.

CATEGORY:

Spectroscopy

CALLING SEQUENCE:
result = atmosdisp(wave,wave_0,za,pressure,temp,[water],[obsalt],$

CANCEL=cancel)

INPUTS:

wave - wavelength in microns wave_0 - reference wavelength in microns za - zenith angle of object [in degrees] pressure - atmospheric pressure in mm of Hg temp - atmospheric temperature in degrees C

OPTIONAL INPUTS:

water - water vapor pressure in mm of Hg. fco2 - relative concentration of CO2 (by pressure) obsalt - The observatory altitude in km.

KEYWORD PARAMETERS:

CANCEL - Set on return if there is a problem

OUTPUTS:

Returns the atmospheric disperion in arcseconds.

PROCEDURE:

Computes the difference between the dispersion at two wavelengths. The dispersion for each wavelength is derived from Section 4.3 of Green’s “Spherical Astronomy” (1985).

EXAMPLE:

MODIFICATION HISTORY:

2000-04-05 - written by M. Cushing, Institute for Astronomy, UH 2002-07-26 - cleaned up a bit. 2003-10-20 - modified formula - WDV 2011-10-07 15:51 IJMC: Converted to Python, with some unit conversions

class spec.baseObject[source]

Empty object container.

spec.calibrate_stared_mosfire_spectra(scifn, outfn, skycorrect, pixcorrect, subreg_corners, **kw)[source]
Correct non-dithered WIDE-slit MOSFIRE spectral frames for:

pixel-to-pixel nonuniformities (i.e., traditional flat-fielding)

detector nonlinearities

tilted spectral lines

non-uniform slit widths (which cause non-smooth backgrounds)

Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frame.

INPUTS:
scifn : str

Filename of raw, uncalibrated science frame (in ADUs, not electrons)

outfn : str

Name into which final, calibrated file should be written.

skycorrect : str or 2D NumPy array

Slitloss correction map (i.e., for slits of nonuniform width), such as generated by make_spectral_flats().

pixcorrect : str or 2D NumPy array

Pixel-by-pixel sensitivity correction map (i.e., flat field), such as generated by make_spectral_flats().

subreg_corners : sequence of 2- or 4-sequences

Indices (or merely starting- and ending-rows) for each subregion of interest in ‘sky’ and ‘domeflat’ frames. Syntax should be that of tools.extractSubregion(), or such that subreg=sky[subreg_corners[0]:subreg_corners[1]]

linearize : bool

Whether to linearity-correct the data.

If linearizing: linearity correction is computed & applied AFTER applying the pixel-by-pixel (flatfield) correction, but BEFORE the slitloss (sky) correction.

bkg_radii : 2-sequence

Inner and outer radius for background computation and removal; measured in pixels from the center of the profile. The values [11,52] seems to work well for MOSFIRE K-band spectra.

locs : None, or sequence
Row-index in each subregion of the location of the spectral-trace-of-interest. Only used for rectifying of sky frame; leaving this at None will have at most a mild effect. If not None, should be the same length as subreg_corners. If subreg_corners[0] = [800, 950] then locs[0] might be set to, e.g., 75 if the trace lies in the middle of the subregion.
gain: scalar
Detector gain, in electrons per ADU
readnoise: scalar
Detector readnoise, in electrons
polycoef : str, None, or sequence
Polynomial coefficients for detector linearization: see ir.linearity_correct() and ir.linearity_mosfire().
unnormalized_flat : str or 2D NumPy array
If not ‘None’, this is used to define the subregion header keywords (‘subreg0’, ‘subreg1’, etc.) for each slit’s two-dimensional spectrum. These keywords are required by much of my extraction machinery!
EXAMPLE:
try:
    from astropy.io import fits as pyfits
except:
    import pyfits

import spec
import ir

obs = ir.initobs('20121010')
allinds = [[7, 294], [310, 518], [532, 694], [710, 960], [976, 1360], [1378, 1673], [1689, 2022]]
locs = [221, 77, 53, 88, 201, 96, 194]
unnorm_flat_fn = obs._raw + 'mudflat.fits'
if False:
   sky = pyfits.getdata(obs._raw + 'masktersky.fits')
   unnorm_domeflat = pyfits.getdata(unnorm_flat_fn)
   skycorrect, pixcorrect = spec.make_spectral_flats(sky, unnorm_domeflat, allinds, badpixelmask=obs.badpixelmask, locs=locs)
else:
   skycorrect = obs._raw + 'skycorrect.fits'
   pixcorrect = obs._raw + 'pixcorrect.fits'

linearcoef='/Users/ianc/proj/transit/data/mosfire_linearity/linearity/mosfire_linearity_cnl_coefficients_bic-optimized.fits'
rawfn = obs.rawsci
outfn = obs.procsci
spec.calibrate_stared_mosfire_spectra(rawfn, outfn, skycorrect, pixcorrect, allinds, linearize=True, badpixelmask=obs.badpixelmask, locs=locs, polycoef=linearcoef, verbose=True, clobber=True, unnormalized_flat=unnorm_flat_fn)
NOTES:

This routine is slow, mostly because of the call to defringe_sinusoid(). Consider running multiple processes in parallel, to speed things up!

SEE_ALSO:

ir.mosfire_speccal(), defringe_sinusoid()

spec.defringe_sinusoid(subframe, badpixelmask=None, nmed=5, dispaxis=0, spatial_index=None, period_limits=[20, 100], retall=False, gain=None, readnoise=None, bictest=False, sinonly=False)[source]

Use simple fitting to subtract fringes and sky background.

subframe : NumPy array
data subframe containing sky data to be subtracted (and, perhaps, an object’s spectral trace). Assumed to be in ADU, not electrons (see gain and readnoise)
nmed : int
Size of window for 2D median filter (to reject bad pixels, etc.)
dispaxis : int
set dispersion axis: 0 = horizontal and 1 = vertical
spatial_index : None, or 1D NumPy array of type bool
Which spatial rows (if dispaxis=0) to use when fitting the tilt of sky lines across the spectrum. If you want to use all, set to None. If you want to ignore some (e.g., because there’s a bright object’s spectrum there) then set those rows’ elements of spatial_index to ‘False’.
period_limits : 2-sequence
Minimum and maximum periods (in pixels) of fringe signals to accept as ‘valid.’ Resolution elements with best-fit periods outside this range will only be fit by a linear trend.
gain : scalar
Gain of detector, in electrons per ADU (where ‘subframe’ is in units of ADUs).
readnoise : scalar
Readnoise of detector, in electrons (where ‘subframe’ is in units of ADUs).
bictest : bool
If True, use ‘gain’ and ‘readnoise’ to compute the Bayesian Information Criterion (BIC) for each fit; a sinusoid is only removed if BIC(sinusoid fit) is lower than BIC(constant fit).
sinonly : bool
If True, the output “model 2D spectrum” will only contain the sinusoidal component. Otherwise, it will contain DC and linear-trend terms.
NOTES:

Note that (in my experience!) this approach works better when you set all weights to unity, rather than using the suggested (photon + read noise) variances.

REQUIREMENTS:

Numerical analysis routines (for analysis.fmin())

Planetary phase curve routines (for phasecurves.errfunc())

lomb (for Lomb-Scargle periodograms)

SciPy ‘signal’ module (for median-filtering)

Returns the best-fit sky frame.

spec.extractSpectralProfiles(args, **kw)[source]

Extract spectrum

INPUTS:
args : tuple or list

either a tuple of (splineProfile, profileStack, errorStack, profileMask), or a list of such tuples (from makeprofile).

OPTIONS:
bkg_radii : 2-sequence

inner and outer radii to use in computing background

extract_radius : int

radius to use for both flux normalization and extraction

RETURNS:
3-tuple:

[0] – spectrum flux (in electrons)

[1] – background flux

[2] – best-fit pixel shift

EXAMPLE:
::

out = spec.traceorders(‘aug16s0399.fits’,nord=7)

out2 = spec.makeprofile(‘aug16s0399.fits’,out,retall=True)

out3 = spec.extractSpectralProfiles(out2)

SEE_ALSO:

optimalExtract()

NOTES:

Note that this is non-optimal with highly tilted or curved spectra, for the reasons described by Marsh (1989) and Mukai (1990).

spec.fit2Gaussian(vec, err=None, verbose=False, guess=None, holdfixed=None)[source]

Fit two Gaussians simultaneously to an input data vector.

INPUTS:
vec : sequence

1D array or list of values to fit to

err : sequence

uncertainties on vec

guess : sequence

guess parameters: [area1, sigma1, cen1, area2, sig2, cen2, constant]. Note that parameters are in pixel units.

holdfixed : sequence

parameters to hold fixed in analysis, _IF_ guess is passed in.

SEE ALSO: analysis.gaussian(), fitGaussian()

spec.fitGaussian(vec, err=None, verbose=False, guess=None)[source]

Fit a Gaussian function to an input data vector.

Return the fit, and uncertainty estimates on that fit.

SEE ALSO: analysis.gaussian()

spec.fitGaussiann(vec, err=None, verbose=False, guess=None, holdfixed=None)[source]

Fit a Gaussian function to an input data vector.

Return the fit, and uncertainty estimates on that fit.

SEE ALSO: analysis.gaussian()

spec.fitPSF(ec, guessLoc, fitwidth=20, verbose=False, sigma=5, medwidth=6, err_ec=None)[source]

Helper function to fit 1D PSF near a given region. Assumes spectrum runs horizontally across the frame!

ec : 2D numpy array
echellogram array, with horizontal dispersion direction
guessLoc : 2-tuple
A slight misnomer for this (x,y) tuple: y is a guess and will
be fit, but x is the coordinate at which the fitting takes place.
fitwidth : int
width of cross-dispersion direction to use in fitting
medwidth : int
number of columns to average over when fitting a profile
verbose : bool
verbosity/debugging printout flag
sigma : scalar
sigma scale for clipping bad values
spec.fitProfileSlices(splineProfile, profileStack, stdProfile, goodPixels, verbose=False, bkg_radii=None, extract_radius=None)[source]

Fit a given spatial profile to a spectrum

Helper function for extractSpectralProfiles()

spec.fitTophat(vec, err=None, verbose=False, guess=None)[source]

Fit a 1D tophat function to an input data vector.

Return the fit, and uncertainty estimates on that fit.

SEE ALSO: analysis.gaussian()

spec.fitprofile(spc, profile, w=None, verbose=False, guess=None, retall=False)[source]

Fit a spline-class spatial profile to a spectrum cross-section

spec.gaussint(x)[source]
PURPOSE:

Compute the integral from -inf to x of the normalized Gaussian

INPUTS:
x : scalar

upper limit of integration

NOTES:

Designed to copy the IDL function of the same name.

spec.humidpressure(RH, T)[source]
  • NAME:

    humidpressure

    PURPOSE:

    To convert relative humidity into a H2O vapor partial pressure

    CATEGORY:

    Spectroscopy

    CALLING SEQUENCE:

    humidpressure(RH, 273.15)

    INPUTS:

    RH - relative humidity, in percent T - temperature, in Kelvin

    OUTPUTS:

    h2o_pp : water vapor partial pressure, in Pascals

    PROCEDURE:

    As outlined in Butler (1998): “Precipitable Water at KP”, MMA Memo 238 (which refers in turn to Liebe 1989, “MPM - An Atmospheric Millimeter-Wave Propagation Model”). Liebe claims that this relation has an error of <0.2% from -40 C to +40 C.

    EXAMPLE:

    None

    MODIFICATION HISTORY:

    2011-10-08 17:08 IJMC: Created.

spec.lightloss(objfile, wguide, seeing, press=None, water=None, temp=None, fco2=None, wobj=None, dx=0, dy=0, retall=False)[source]
  • NAME:

    lightloss

    PURPOSE:

    To determine the slit losses from a spectrum.

    CATEGORY:

    Spectroscopy

    CALLING SEQUENCE:

    ### TBD lightloss, obj, std, wguide, seeing, out, CANCEL=cancel

    INPUTS:

    obj - FITS file of the object spectrum wguide - wavelength at which guiding was done seeing - seeing FWHM at the guiding wavelength

    OPTIONAL INPUTS:

    press - mm Hg typical value (615 for IRTF, unless set) water - mm Hg typical value (2 for IRTF, unless set) temp - deg C typical value (0 for IRTF, unless set) fco2 - relative concentration of CO2 (0.004, unless set) wobj - wavelength scale for data dx - horizontal offset of star from slit center dy - vertical offset of star from slit center retall- whether to return much diagnostic info, or just lightloss.

    NOTES:

    ‘seeing’, ‘dx’, and ‘dy’ should all be in the same units, and also the same units used to define the slit dimensions in the obj FITS file header

    KEYWORD PARAMETERS:

    CANCEL - Set on return if there is a problem

    OUTPUTS:

    array : fractional slit loss at each wavelength value

    OR:

    tuple of arrays: (slitloss, disp_obj, diff, fwhm, dx_obj, dy_obj)

    PROCEDURE:

    Reads a Spextool FITS file.

    EXAMPLE:

    None

    REQUIREMENTS:

    Aperture and simple PSF-fitting photometrymem pyfits

    MODIFICATION HISTORY:

    2003-10-21 - Written by W D Vacca 2011-10-07 20:19 IJMC: Converted to Python, adapted for single objects 2011-10-14 14:01 IJMC: Added check for Prism mode (has

    different slit dimension keywords) and different pyfits header read mode.

    2011-11-07 15:53 IJMC: Added ‘retall’ keyword

spec.lightloss2(wobj, slitwd, slitht, slitPA, targetPA, zenith_angle, wguide, seeing, press=615.0, water=2.0, temp=0.0, fco2=0.004, obsalt=4.16807, teldiam=3.0, dx=0, dy=0, retall=False, ydisp=None, xdisp=None, fwhm=None)[source]
  • NAME:

    lightloss2

    PURPOSE:

    To determine the slit losses from an observation (no FITS file involved)

    CATEGORY:

    Spectroscopy

    CALLING SEQUENCE:

    ### TBD lightloss, obj, std, wguide, seeing, out, CANCEL=cancel

    INPUTS:

    wobj - wavelength scale for data slitwd - width of slit, in arcsec slitht - height of slit, in arcsec slitPA - slit Position Angle, in radians targetPA - Parallactic Angle at target, in radians zenith_angle - Zenith Angle, in radians wguide - wavelength at which guiding was done seeing - seeing FWHM at the guiding wavelength

    OPTIONAL INPUTS:

    press - mm Hg typical value (615, unless set) water - mm Hg typical value (2 , unless set) temp - deg C typical value (0 , unless set) fco2 - relative concentration of CO2 (0.004, unless set) obsalt- observatory altitude, in km teldiam- observatory limiting aperture diameter, in m dx - horizontal offset of star from slit center dy - vertical offset of star from slit center retall- whether to return much diagnostic info, or just lightloss.

    ydisp - The position of the spectrum in the slit at all

    values of wobj. This should be an array of the same size as wobj, with zero corresponding to the vertical middle of the slit and positive values tending toward zenith. In this case xdisp will be computed as XXXX rather than from the calculated atmospheric dispersion; dx and dy will also be ignored.

    fwhm - Full-width at half-maximum of the spectral trace, at

    all values of wobj. This should be an array of the same size as wobj, measured in arc seconds.

    NOTES:

    ‘slitwidth’, ‘slitheight’, ‘seeing’, ‘dx’, ‘dy’, and ‘fwhm’ (if used) should all be in the same units: arc seconds.

    OUTPUTS:

    array : fractional slit loss at each wavelength value

    OR:

    tuple of arrays: (slitloss, disp_obj, diff, fwhm, dx_obj, dy_obj)

    PROCEDURE:

    All input-driven. For the SpeXTool-version analogue, see lightloss()

    EXAMPLE:

    import numpy as np import spec w = np.linspace(.5, 2.5, 100) # Wavelength, in microns d2r = np.deg2rad(1.) #targetPA, za = spec.parangle(1.827, 29.67*d2r, lat=20.*d2r) targetPA, za = 105.3, 27.4 slitPA = 90. * d2r

    spec.lightloss2(w, 3., 15., slitPA, targetPA*d2r, za*d2r, 2.2, 1.0)

    REQUIREMENTS:

    Aperture and simple PSF-fitting photometrymem

    MODIFICATION HISTORY:

    2003-10-21 - Written by W D Vacca 2011-10-07 20:19 IJMC: Converted to Python, adapted for single objects 2011-10-14 14:01 IJMC: Added check for Prism mode (has

    different slit dimension keywords) and different pyfits header read mode.

    2011-11-07 15:53 IJMC: Added ‘retall’ keyword 2011-11-07 21:17 IJMC: Cannibalized from SpeXTool version 2011-11-25 15:06 IJMC: Added ydisp and fwhm options.

spec.makeSpexSlitlessSky(skyfns, scatcen=[980, 150], scatdim=[60, 300])[source]

Generate a normalized Sky frame from SpeX slitless spectroscopy data.

INPUTS:
skyfns : list of strs

Filenames of slitless-spectroscopy sky frames

scatcen : 2-sequence

center of region to use in median-normalizing the frame.

scatdim : 2-sequence

full width of region to use in median-normalizing the frame.

OUTPUTS:

(sky, skyHeader)

spec.make_spectral_flats(sky, domeflat, subreg_corners, badpixelmask=None, xord_pix=[15, 2], xord_sky=[2, 1], yord=2, minsnr=5, minfrac_pix=0.7, minfrac_sky=0.5, locs=None, nsigma=3)[source]

Construct appropriate corrective frames for multi-object spectrograph data. Specifically: corrections for irregular slit widths, and pixel-by-pixel detector sensitivity variations.

sky : 2D NumPy array

Master spectral sky frame, e.g. from median-stacking many sky frames or masking-and-stacking dithered science spectra frames. This frame is used to construct a map to correct science frames (taken with the identical slit mask!) for irregular sky backgrounds resulting from non-uniform slit widths (e.g., Keck/MOSFIRE).

Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frames.

domeflat : 2D NumPy array

Master dome spectral flat, e.g. from median-stacking many dome spectra. This need not be normalized in the dispersion or spatial directions. This frame is used to construct a flat map of the pixel-to-pixel variations in detector sensitivity.

Note that the dispersion direction should be ‘horizontal’ (i.e., parallel to rows) in this frames.

subreg_corners : sequence of 2- or 4-sequences
Indices (or merely starting- and ending-rows) for each subregion of interest in ‘sky’ and ‘domeflat’ frames. Syntax should be that of tools.extractSubregion(), or such that subreg=sky[subreg_corners[0]:subreg_corners[1]]
badpixelmask : 2D NumPy array, or str
Nonzero for any bad pixels; these will be interpolated over using nsdata.bfixpix().
xord_pix : sequence
Polynomial orders for normalization in dispersion direction of pixel-based flat (dome flat) on successive iterations; see makexflat().
xord_sky : sequence
Polynomial orders for normalization in dispersion direction of slit-based correction (sky frame) on successive iterations; see makexflat().
yord : scalar
Polynomial order for normalization of pixel-based (dome) flat in spatial direction.
locs : None, or sequence
Row-index in each subregion of the location of the spectral-trace-of-interest. Only used for rectifying of sky frame; leaving this at None will have at most a mild effect. If not None, should be the same length as subreg_corners. If subreg_corners[0] = [800, 950] then locs[0] might be set to, e.g., 75 if the trace lies in the middle of the subregion.
RETURNS:

wideslit_skyflat, narrowslit_domeflat

EXAMPLE:
try:
    from astropy.io import fits as pyfits
except:
    import pyfits

import spec
import ir

obs = ir.initobs('20121010')
sky = pyfits.getdata(obs._raw + 'masktersky.fits')
domeflat = pyfits.getdata(obs._raw + 'mudflat.fits')
allinds = [[7, 294], [310, 518], [532, 694], [710, 960], [976, 1360], [1378, 1673], [1689, 2022]]
locs = [221, 77, 53, 88, 201, 96, 194]

skycorrect, pixcorrect = spec.make_spectral_flats(sky, domeflat, allinds, obs.badpixelmask, locs=locs)
spec.makemodel(params, xvec, specvec, wtemplate)[source]

Helper function for wavelengthMatch(): generate a scaled, interpolative model of the template.

spec.makeprofile(filename, trace, **kw)[source]

Make a spatial profile from a spectrum, given its traced location. We interpolate the PSF at each pixel to a common reference frame, and then average them.

filename : str _OR_ 2D numy array
2D echellogram
trace : 1D numpy array
set of polynomial coeficients of order (P-1)
dispaxis : int
set dispersion axis: 0 = horizontal and = vertical
fitwidth : int
Total width of extracted spectral sub-block. WIll always be increased to at least twice the largest value of bkg_radii.
neg : bool scalar
set True for a negative spectral trace
nsigma : scalar
Sigma-clipping cut for bad pixels (beyond read+photon noise). Set it rather high, and feel free to experiment with this parameter!
xylims : 4-sequence
Extract the given subset of the data array: [xmin, xmax, ymin, ymax]
retall : bool
Set True to output several additional parameters (see below)
rn : scalar
Read noise (electrons)
g : scalar
Detector gain (electrons per data unit)
bkg_radii : 2-sequence
Inner and outer radius for background computation and removal; measured in pixels from the center of the profile.
bkg_order : int > 0
Polynomial order of background trend computed in master spectral profile
interp : bool
Whether to (bi-linearly) interpolate each slice to produce a precisely-centered spectral profile (according to the input ‘trace’). If False, slices will only be aligned to the nearest pixel.
OUTPUT:
if retall:

a spline-function that interpolates pixel locations onto the mean profile

a stack of data slices

estimates of the uncertainties

good pixel flag

list of splines

else:
a spline-function that interpolates pixel locations onto the mean profile
spec.makexflat(subreg, xord, nsigma=3, minsnr=10, minfrac=0.5, niter=1)[source]

Helper function for XXXX.

INPUTS:
subreg : 2D NumPy array

spectral subregion, containing spectral background, sky, and/or target flux measurements.

xord : scalar or sequence

Order of polynomial by which each ROW will be normalized. If niter>0, xord can be a sequence of length (niter+1). A good approach for, e.g., spectral dome flats is to set niter=1 and xord=[15,2].

nsigma : scalar

Sigma-clipping level for calculation of column-by-column S/N

minsnr : scalar

Minimum S/N value to use when selecting ‘good’ columns for normalization.

minfrac : scalar, 0 < minfrac < 1

Fraction of columns to use, selected by highest S/N, when selecting ‘good’ columns for normalization.

niter : int

Number of iterations. If set to zero, do not iterate (i.e., run precisely once through.)

NOTES:

Helper function for functions XXXX

spec.message(text)[source]

Display a message; for now, with text.

spec.modelSpectralTrace(param, shape=None, nscat=None, npw=None, npy=None, noy=None, now=None, ndist=None, x=None, y=None, transpose=False)[source]

Model a raw spectral trace!

INPUTS:

NOTE that most inputs should be in the _rectified_ frame.

Trace background pedestal level : 1D array

Width of background pedestarl level : scalar (for now)

Center of trace : 1D array

Offset of object spectrum, relative to center : scalar

width of 1D PSF : scalar

Area of 1D psf : 1D array

Distortion (x and y, somehow???)

Scattered light background : scalar (for now)

spec.model_resel(param, x)[source]

Model a spectral resolution element.

INPUTS:

param : sequence

param[0, 1, 2] - amplitude, sigma, and central location of Gaussian line profile (cf. analysis.gaussian()).

param[3, 4, 5] - amplitude, width, and central location of top-hat-like background (cf. tophat()).

param[6::] - additional (constant or polynomial) background components, for evaluation with numpy.polyval()

x : sequence

Values at which to evaluate model function (i.e., pixels). Typically 1D.

OUTPUTS:
line : NumPy array

model of the resolution element, of same shape as x.

DESCRIPTION:

Model a spectral resolution element along the spatial direction. This consists of a (presumably Gaussian) line profile superimposed on the spectral trace’s top-hat-like background, with an additional constant (or polynomial) out-of-echelle-order background component.

spec.modelline(param, prof2, dv)[source]

Generate a rotational profile, convolve it with a second input profile, normalize it (simply), and return.

INPUTS:
param : 1D sequence

param[0:3] – see rotationalProfile() param[4] – multiplicative scaling factor

dv : velocity grid

prof2 : second input profile

spec.normalizeSpecFlat(flatdat, nspec=1, minsep=50, median_width=51, readnoise=40, badpixelmask=None, traces=None)[source]

Trace and normalize a spectroscopic flat field frame.

INPUTS:
flatdat : 2D NumPy array

Master, unnormalized flat frame: assumed to be measured in photoelectrons (for computing uncertainties).

nspec : int

Number of spectral orders to find and normalize

minsep : int

Minimum separation, in pixels, between spectral orders that will be found.

median_width : int

Width of median-filter kernel used to compute the low-

readnoise : scalar

Detector read noise, in electrons. For computing uncertainties.

badpixelmask : 2D NumPy array

bad pixel mask: 1 at bad pixel locations, 0 elsewhere.

traces : 2D NumPy Array
(nord, pord) shaped numpy array representing the polynomial

coefficients for each order (suitable for use with np.polyval), as produced by traceorders()

spec.optimalExtract(*args, **kw)[source]

Extract spectrum, following Horne 1986.

INPUTS:
data : 2D Numpy array

Appropriately calibrated frame from which to extract spectrum. Should be in units of ADU, not electrons!

variance : 2D Numpy array

Variances of pixel values in ‘data’.

gain : scalar

Detector gain, in electrons per ADU

readnoise : scalar

Detector readnoise, in electrons.

OPTIONS:
goodpixelmask : 2D numpy array

Equals 0 for bad pixels, 1 for good pixels

bkg_radii : 2- or 4-sequence

If length 2: inner and outer radii to use in computing background. Note that for this to be effective, the spectral trace should be positions in the center of ‘data.’

If length 4: start and end indices of both apertures for background fitting, of the form [b1_start, b1_end, b2_start, b2_end] where b1 and b2 are the two background apertures, and the elements are arranged in strictly ascending order.

extract_radius : int or 2-sequence

radius to use for both flux normalization and extraction. If a sequence, the first and last indices of the array to use for spectral normalization and extraction.

dispaxis : bool

0 for horizontal spectrum, 1 for vertical spectrum

bord : int >= 0

Degree of polynomial background fit.

bsigma : int >= 0

Sigma-clipping threshold for computing background.

pord : int >= 0

Degree of polynomial fit to construct profile.

psigma : int >= 0

Sigma-clipping threshold for computing profile.

csigma : int >= 0

Sigma-clipping threshold for cleaning & cosmic-ray rejection.

finite : bool

If true, mask all non-finite values as bad pixels.

nreject : int > 0

Number of pixels to reject in each iteration.

RETURNS:
3-tuple:

[0] – spectrum flux (in electrons)

[1] – uncertainty on spectrum flux

[1] – background flux

EXAMPLE:
SEE_ALSO:

superExtract().

NOTES:

Horne’s classic optimal extraction algorithm is optimal only so long as the spectral traces are very nearly aligned with detector rows or columns. It is not well-suited for extracting substantially tilted or curved traces, for the reasons described by Marsh 1989, Mukai 1990. For extracting such spectra, see superExtract().

spec.optspecextr_idl(frame, gain, readnoise, x1, x2, idlexec, clobber=True, tempframefn='tempframe.fits', specfn='tempspec.fits', scriptfn='temp_specextract.pro', IDLoptions="adjfunc='adjgauss', adjoptions={center:1, centerfit:1, centerdeg:3}, bgdeg=3", inmask=None)[source]

Run optimal spectral extraction in IDL; pass results to Python.

INPUTS:
frame : str

filename, or 2D Numpy Array, or list of filenames containing frames from which spectra will be extracted. This should be in units of ADU (not electrons) for the noise properties to come out properly.

Also, the spectral trace must run vertically across the frame.

gain : scalar

Detector gain, in electrons / ADU

readnoise : scalar

Detector read noise, in electrons

x1, x2 : ints, or lists of ints

Start and stop indices of the spectral trace across the frame. If multiple frames are input and a single x1/x2 is input, the same value will be used for each frame. Note however that multiple x1/x2 can also be input (one for each frame).

idlexec : str

Path to the IDL executable. OPTSPECEXTR.PRO and its associated files must be in your IDL path. If set to None, then it will be set to: os.popen(‘which idl’).read().strip()

OPTIONS:
clobber : bool

Whether to overwrite files when writing input data to TEMPFRAMFN.

tempframefn : str

If input ‘frame’ is an array, it will be written to this filename in order to pass it to IDL.

specfn : str

IDL will write the spectral data to this filename in order to pass it back to Python.

scriptfn : str

Filename in which the short IDL script will be written.

IDLoptions : str

Options to pass to OPTSPECEXTR.PRO. For example: “adjfunc=’adjgauss’, adjoptions={center:1,centerfit:1,centerdeg:3}, bgdeg=3”

Note that this Python code will break if you _don’t_ trace the spectrum (adjoptions, etc.); this is an area for future work if I ever use a spectrograph with straight traces.

inmask : None or str

Name of the good pixel mask for OPTSPECEXTR.PRO. Equal to 1 for good pixels, and 0 for bad pixels.

OUTPUTS:
For each input frame, a list of four items:

[0] – Extracted spectrum, ADU per pixel [1] – Uncertainty (1 sigma) of extracted spectrum [2] – Location of trace (in pixels) across the frame [3] – Width of trace across the frame

NOTES:

Note that this more closely follows Horne et al. than does optimalExtract(), and is faster than both that function and (especially!) extractSpectralProfiles(). The only downside (if it is one) is that this function requires IDL.

TO-DO:

Add options for user input of a variance frame, or of sky variance.

Allow more flexibility (tracing, input/output options, etc.)

REQUIREMENTS:

IDL

OPTSPECEXTR

spec.parangle(HA, DEC, lat)[source]
  • NAME:

    parangle

    PURPOSE:

    To compute the parallactic angle at a given position on the sky.

    CATEGORY:

    Spectroscopy

    CALLING SEQUENCE:

    eta, za = parangle(HA, DEC, lat)

    INPUTS:

    HA - Hour angle of the object, in decimal hours (0,24) DEC - Declination of the object, in degrees lat - The latitude of the observer, in degrees

    KEYWORD PARAMETERS:

    CANCEL - Set on return if there is a problem

    OUTPUTS:

    eta - The parallactic angle za - The zenith angle

    PROCEDURE:

    Given an objects HA and DEC and the observers latitude, the zenith angle and azimuth are computed. The law of cosines then gives the parallactic angle.

    EXAMPLE:

    NA

    MODIFICATION HISTORY:

    2000-04-05 - written by M. Cushing, Institute for Astronomy,UH 2002-08-15 - cleaned up a bit. 2003-10-21 - changed to pro; outputs zenith angle as well - WDV 2011-10-07 17:58 IJMC: Converted to Python

spec.pickloc(ax=None, zoom=10)[source]
INPUTS:

ax : (axes instance) – axes in which to pick a location

zoom : int – zoom radius for target confirmation

: 2-tuple – (x,y) radii for zoom confirmation.

spec.profileError(param, spc, profile, w, xaxis=None, retresidual=False)[source]

Compute the chi-squared error on a spectrum vs. a profile

spec.reconstitute_gmos_roi(infile, outfile, **kw)[source]

Convert GMOS frames taken with custom ROIs into standard FITS frames.

INPUTS:
in : string or sequence of strings

Input filename or filenames.

out : string or sequence of strings

Output filename or filenames.

OPTIONS:
clobber : bool

Passed to PyFITS; whether to overwrite existing files.

spec.resamplespec(w1, w0, spec0, oversamp=100)[source]

Resample a spectrum while conserving flux density.

INPUTS:
w1 : sequence

new wavelength grid (i.e., center wavelength of each pixel)

w0 : sequence

old wavelength grid (i.e., center wavelength of each pixel)

spec0 : sequence

old spectrum (e.g., flux density or photon counts)

oversamp : int

factor by which to oversample input spectrum prior to rebinning. The worst fractional precision you achieve is roughly 1./oversamp.

NOTE:

Format is the same as numpy.interp()

REQUIREMENTS:

Plotting and analysis tools

spec.rotationalProfile(delta_epsilon, delta_lam)[source]

Compute the rotational profile of a star, assuming solid-body rotation and linear limb darkening.

This uses Eq. 18.14 of Gray’s Photospheres, 2005, 3rd Edition.

INPUTS:

delta_epsilon : 2-sequence

[0] : delta_Lambda_L = lambda * V * sin(i)/c; the rotational

displacement at the stellar limb.

[1] : epsilon, the linear limb darkening coefficient, used in

the relation I(theta) = I0 + epsilon * (cos(theta) - 1).

[2] : OPTIONAL! The central location of the profile (otherwise

assumed to be located at delta_lam=0).

delta_lam : scalar or sequence

Wavelength minus offset: Lambda minus lambda_0. Grid upon which computations will be done.

EXAMPLE:
import pylab as py
import spec

dlam = py.np.linspace(-2, 2, 200) # Create wavelength grid
profile = spec.rotationalProfile([1, 0.6], dlam)

py.figure()
py.plot(dlam, profile)
spec.runlblrtm(lamrange, pwv=2.0, zang=0.0, alt=4.2, co2=390, res=200, dotrans=True, dorad=False, pwv_offset=4.0, verbose=False, _save='/Users/ianc/temp/testatmo.mat', _wd='/Users/ianc/proj/atmo/aerlbl_v12.0_package/caltech_wrapper_v02/', scriptname='runtelluric.m', command='/Applications/Octave.app/Contents/Resources/bin/octave')[source]

Run LBLRTM to compute atmospheric transmittance and/or radiance.

INPUTS:
lamrange : 2-sequence

approximate minimum and maximum wavelengths

OPTIONS:
pwv : float

mm of Precipitable Water Vapor above observation site. If negative, then the value abs(pwv_offset-pwv) will be used instead.

pwv_offset : float

Only used if (pwv < 0); see above for description.

zang : float

observation angle from zenith, in degrees

alt : float

Observation elevation, in km.

co2 : float

CO2 concentration in ppm by volume. Concentration is assumed to be uniform throughout the atmosphere.

res : float

approximate spectral resolution desired

dotrans : bool

compute atmospheric transmittance

dorad : bool

compute atmospheric radiance. NOT CURRENTLY WORKING

_save : str

path where temporary MAT save file will be stored

_wd : str

path where MATLAB wrapper scripts for LBLRTM are located

scriptname : str

filename for temporary MATLAB/OCTAVE script (saved after exit)

command : str

path to MATLAB/OCTAVE executable

OUTPUTS:

A 2- or 3-tuple: First element is wavelength in microns, second element is transmittance (if requested). Radiance will (if requested) always be the last element, and in f_nu units: W/cm2/sr/(cm^-1)

REQUIREMENTS:

SciPy

OCTAVE or MATLAB

LBLRTM

D. Feldman’s set of MATLAB wrapper scripts

spec.scaleSpectralSky_cor(subframe, badpixelmask=None, maxshift=20, fitwidth=2, pord=1, nmed=3, dispaxis=0, spatial_index=None, refpix=None, tord=2)[source]

Use cross-correlation to subtract tilted sky backgrounds.

subframe : NumPy array
data subframe containing sky data to be subtracted (and, perhaps, an object’s spectral trace).
badpixelmask : None or NumPy array
A boolean array, equal to zero for good pixels and unity for bad pixels. If this is set, the first step will be a call to nsdata.bfixpix() to interpolate over these values.
nmed : int
size of 2D median filter for pre-smoothing.
pord : int
degree of spectral tilt. Keep this number low!
maxshift : int
Maximum acceptable shift. NOT YET IMPLEMENTED!
fitwidth : int
Maximum radius (in pixels) for fitting to the peak of the cross-correlation.
nmed : int
Size of window for 2D median filter (to reject bad pixels, etc.)
dispaxis : int
set dispersion axis: 0 = horizontal and 1 = vertical
spatial_index : None, or 1D NumPy array of type bool
Which spatial rows (if dispaxis=0) to use when fitting the tilt of sky lines across the spectrum. If you want to use all, set to None. If you want to ignore some (e.g., because there’s a bright object’s spectrum there) then set those rows’ elements of spatial_index to ‘False’.
refpix : scalar
Pixel along spatial axis to which spectral fits should be aligned; if a spectral trace is present, one should set “refpix” to the location of the trace.
tord : int
Order of polynomial fits along spatial direction in aligned 2D-spectral frame, to account for misalignments or irregularities of tilt direction.
RETURNS:a model of the sky background, of the same shape as ‘subframe.’
spec.scaleSpectralSky_dsa(subframe, variance=None, badpixelmask=None, nk=31, pord=1, nmed=3, gain=3.3, readnoise=30, dispaxis=0, spatial_index=None)[source]

Use difference-imaging techniques to subtract moderately tilted sky background. Doesn’t work so well!

subframe : NumPy array
data subframe containing sky data to be subtracted (and, perhaps, an object’s spectral trace). Assumed to be in ADU, not electrons (see gain and readnoise)
variance : None or NumPy array
variance of each pixel in ‘subframe’
nmed : int
size of 2D median filter
pord : int
degree of spectral tilt. Keep this number low!
nk : int
Number of kernel pixels in dia.dsa()
nmed : int
Size of window for 2D median filter (to reject bad pixels, etc.)
dispaxis : int
set dispersion axis: 0 = horizontal and 1 = vertical
gain, readnoise : ints
If ‘variance’ is None, these are used to estimate the uncertainties.
spatial_index : None, or 1D NumPy array of type bool
Which spatial rows (if dispaxis=0) to use when fitting the tilt of sky lines across the spectrum. If you want to use all, set to None. If you want to ignore some (e.g., because there’s a bright object’s spectrum there) then set those rows’ elements of spatial_index to ‘False’.
NOTES:Note that (in my experience!) this approach works better when you set all weights to unity, rather than using the suggested (photon + read noise) variances.

Returns the best-fit sky frame.

spec.scaleSpectralSky_pca(frame, skyframes, variance=None, mask=None, badpixelmask=None, npca=3, gain=3.3, readnoise=30)[source]

Use PCA and blank sky frames to subtract

frame : str or NumPy array
data frame to subtract sky from. Assumed to be in ADU, not electrons (see gain and readnoise)
npca : int
number of PCA components to remove

f0 = pyfits.getdata(odome.procsci[0]) mask = pyfits.getdata(odome._proc + ‘skyframes_samplemask.fits’).astype(bool) badpixelmask = pyfits.getdata( odome.badpixelmask).astype(bool)

The simplest way to fit sky to a ‘frame’ containing bright spectral is to include the spectral-trace regions in ‘mask’ but set the ‘variance’ of those regions extremely high (to de-weight them in the least-squares fit).

To use for multi-object data, consider running multiple times (once per order)

Returns the best-fit sky frame as determined from the first ‘npca’ PCA components.

spec.slittrans(*varargin)[source]
NAME:

slittrans

PURPOSE:

Compute flux passing through a slit assuming a gaussian PSF.

CATEGORY:

Spectroscopy

CALLING SEQUENCE:
 

result = slittrans(width,height,fwhm,xoffset,yoffset,CANCEL=cancel)

INPUTS:

width - Width of slit. height - Height of slit. fwhm - Full-width at half-maximum of the gaussian image. xoffset - Offset in x of the image from the center of the slit. yoffset - Offset in y of the image from the center of the slit.

Note: the units are arbitrary but they must be the same across

all of the input quantities.

KEYWORD PARAMETERS:
CANCEL - Set on return if there is a problem
OUTPUTS:
Returned is the fraction of the total gaussian included in the slit.
EXAMPLE:

result = slittrans(0.3,15,0.6,0,0)

Computes the fraction of the flux transmitted through a slit 0.3x15 arcseconds with a PSF of 0.6 arcseconds FWHM. The PSF is centered on the slit.

MODIFICATION HISTORY:

Based on M Buie program, 1991 Mar., Marc W. Buie, Lowell Observatory Modified 2000 Apr., M. Cushing to include y offsets. 2011-10-07 15:45 IJMC: Converted to Python 2011-11-14 16:29 IJMC: Rewrote to use erf() rather than

spec.spexsxd_scatter_fix(fn1, fn2, **kw)[source]

Fix scattered light in SpeX/SXD K-band and write a new file.

INPUTS:
fn1 : str

file to be fixed

fn2 : str

new filename of fixed file.

OPTIONS:
clobber : bool

whether to overwrite existing FITS files

Other options will be passed to spexsxd_scatter_model()

OUTPUTS:
status : int

0 if a problem, 1 if everything is Okay

spec.spexsxd_scatter_model(dat, halfwid=48, xlims=[470, 1024], ylims=[800, 1024], full_output=False, itime=None)[source]

Model the scattered light seen in SpeX/SXD K-band frames.

INPUTS:
dat : str or numpy array

filename of raw SXD frame to be corrected, or a Numpy array containing its data.

OPTIONS:
halfwid : int

half-width of the spectral orders. Experience shows this is approximately 48 pixels. This value is not fit!

xlims : list of length 2

minimum and maximum x-pixel values to use in the fitting

ylims : list of length 2

minimum and maximum y-pixel values to use in the fitting

full_output : bool

whether to output only model, or the tuple (model, fits, chisq, nbad)

itime : float

integration time, in seconds, with which to scale the initial guesses

OUTPUT:
scatter_model : numpy array

Model of the scattered light component, for subtraction or saving.

OR:

scatter_model, fits, chis, nbad

REQUIREMENTS:

pyfits, numpy, fit_atmo, Numerical analysis routines, Planetary phase curve routines

TO_DO_LIST:

I could stand to be more clever in modeling the scattered light components – perhaps fitting for the width, or at least allowing the width to be non-integer.

spec.spextractor(frame, gain, readnoise, framevar=None, badpixelmask=None, mode='superExtract', trace=None, options=None, trace_options=None, verbose=False)[source]

Extract a spectrum from a frame using one of several methods.

INPUTS:
frame : 2D Numpy array or filename

Contains a single spectral trace.

gain : None or scalar

Gain of data contained in ‘frame;’ i.e., number of collected photoelectrons equals frame * gain.

readnoise : None or scalar

Read noise of detector, in electrons.

framevar : None, 2D Numpy array, or filename

Variance of values in ‘frame.’

If and only if framevar is None, use gain and readnoise to compute variance.

badpixelmask : None, 2D Numpy array, or filename

Mask of bad pixels in ‘frame.’ Bad pixels are set to 1, good pixels are set to 0.

mode : str

Which spectral extraction mode to use. Options are:

superExtract – see superExtract()

optimalExtract – see optimalExtract()

spline – see extractSpectralProfiles()

Must also input trace_options.

trace : None, or 1D Numpy Array

Spectral trace location: fractional pixel index along the entire spectral trace. If None, traceorders() will be called using the options in ‘trace_options.’

options : None or dict

Keyword options to be passed to the appropriate spectral extraction algorithm. Note that you should be able to pass the same sets of parameters to superExtract() and optimalExtract() (the necessary parameter sets overlap, but are not identical).

trace_options : None or dict

Keyword options to be passed to traceorders() (if no trace is input, or if mode=’spline’)

RETURNS:

spectrum, error or variance of spectrum, sky background, ...

NOTES:

When ‘optimalextract’ is used: if len(bkg_radii)==2 then the background apertures will be reset based on the median location of the trace. If extract_radius is a singleton, it will be similarly redefined.

EXAMPLE:
frame = pyfits.getdata('spectral_filename.fits')
gain, readnoise = 3.3, 30
output = spec.spextractor(frame, gain, readnoise, mode='superExtract',                options=dict(bord=2, bkg_radii=[20, 30], extract_radius=15,                polyspacing=1./3, pord=5, verbose=True, trace=trace,                qmode='slow-linear'))

output2 = spec.spextractor(frame, gain, readnoise, mode='optimalExtract',                 options=dict(bkg_radii=[20,30], extract_radius=15, bord=2,                 bsigma=3, pord=3, psigma=8, csigma=5, verbose=1))
spec.superExtract(*args, **kw)[source]

Optimally extract curved spectra, following Marsh 1989.

INPUTS:
data : 2D Numpy array

Appropriately calibrated frame from which to extract spectrum. Should be in units of ADU, not electrons!

variance : 2D Numpy array

Variances of pixel values in ‘data’.

gain : scalar

Detector gain, in electrons per ADU

readnoise : scalar

Detector readnoise, in electrons.

OPTIONS:
trace : 1D numpy array

location of spectral trace. If None, traceorders() is invoked.

goodpixelmask : 2D numpy array

Equals 0 for bad pixels, 1 for good pixels

npoly : int

Number of profile polynomials to evaluate (Marsh’s “K”). Ideally you should not need to set this – instead, play with ‘polyspacing’ and ‘extract_radius.’ For symmetry, this should be odd.

polyspacing : scalar

Spacing between profile polynomials, in pixels. (Marsh’s “S”). A few cursory tests suggests that the extraction precision (in the high S/N case) scales as S^-2 – but the code slows down as S^2.

pord : int

Order of profile polynomials; 1 = linear, etc.

bkg_radii : 2-sequence

inner and outer radii to use in computing background

extract_radius : int

radius to use for both flux normalization and extraction

dispaxis : bool

0 for horizontal spectrum, 1 for vertical spectrum

bord : int >= 0

Degree of polynomial background fit.

bsigma : int >= 0

Sigma-clipping threshold for computing background.

tord : int >= 0

Degree of spectral-trace polynomial (for trace across frame – not used if ‘trace’ is input)

csigma : int >= 0

Sigma-clipping threshold for cleaning & cosmic-ray rejection.

finite : bool

If true, mask all non-finite values as bad pixels.

qmode : str (‘fast’ or ‘slow’)

How to compute Marsh’s Q-matrix. Valid inputs are ‘fast-linear’, ‘slow-linear’, ‘fast-nearest,’ ‘slow-nearest,’ and ‘brute’. These select between various methods of integrating the nearest-neighbor or linear interpolation schemes as described by Marsh; the ‘linear’ methods are preferred for accuracy. Use ‘slow’ if you are running out of memory when using the ‘fast’ array-based methods. ‘Brute’ is both slow and inaccurate, and should not be used.

nreject : int

Number of outlier-pixels to reject at each iteration.

retall : bool

If true, also return the 2D profile, background, variance map, and bad pixel mask.

RETURNS:
object with fields for:

spectrum

varSpectrum

trace

EXAMPLE:
import spec
import numpy as np
import pylab as py

def gaussian(p, x):
   if len(p)==3:
       p = concatenate((p, [0]))
   return (p[3] + p[0]/(p[1]*sqrt(2*pi)) * exp(-(x-p[2])**2 / (2*p[1]**2)))

# Model some strongly tilted spectral data:
nx, nlam = 80, 500
x0 = np.arange(nx)
gain, readnoise = 3.0, 30.
background = np.ones(nlam)*10
flux =  np.ones(nlam)*1e4
center = nx/2. + np.linspace(0,10,nlam)
FWHM = 3.
model = np.array([gaussian([flux[ii]/gain, FWHM/2.35, center[ii], background[ii]], x0) for ii in range(nlam)])
varmodel = np.abs(model) / gain + (readnoise/gain)**2
observation = np.random.normal(model, np.sqrt(varmodel))
fitwidth = 60
xr = 15

output_spec = spec.superExtract(observation, varmodel, gain, readnoise, polyspacing=0.5, pord=1, bkg_radii=[10,30], extract_radius=5, dispaxis=1)

py.figure()
py.plot(output_spec.spectrum.squeeze() / flux)
py.ylabel('(Measured flux) / (True flux)')
py.xlabel('Photoelectrons')
TO_DO:

Iterate background fitting and reject outliers; maybe first time would be unweighted for robustness.

Introduce even more array-based, rather than loop-based, calculations. For large spectra computing the C-matrix takes the most time; this should be optimized somehow.

SEE_ALSO:
spec.tophat2(param, x)[source]

Grey-pixel tophat function with set width param: [cen_pix, amplitude, background] newparam: [amplitude, full width, cen_pix, background] x : must be array of ints, arange(0, size-1) returns the model.

spec.tophat_alt(param, x)[source]

Standard tophat function (alternative version).

INPUTS:
p : sequence

p[0] – Amplitude p[1] – full width dispersion p[2] – central offset (mean location) p[3] – vertical offset (OPTIONAL)

x : scalar or sequence

values at which to evaluate function

OUTPUTS:
y : scalar or sequence

1.0 where |x| < 0.5, 0.5 where |x| = 0.5, 0.0 otherwise.

spec.traceorders(filename, pord=5, dispaxis=0, nord=1, verbose=False, ordlocs=None, stepsize=20, fitwidth=20, plotalot=False, medwidth=6, xylims=None, uncertainties=None, g=5, rn=25, badpixelmask=None, retsnr=False, retfits=False)[source]

Trace spectral orders for a specified filename.

filename : str OR 2D array
full path and filename to a 2D echelleogram FITS file, _OR_ a 2D numpy array representing such a file.

OPTIONAL INPUTS: pord : int

polynomial order of spectral order fit
dispaxis : int
set dispersion axis: 0 = horizontal and = vertical
nord : int
number of spectral orders to trace.
ordlocs : (nord x 2) numpy array
Location to “auto-click” and
verbose: int
0,1,2; whether (and how much) to print various verbose debugging text
stepsize : int
number of pixels to step along the spectrum while tracing
fitwidth : int
number of pixels to use when fitting a spectrum’s cross-section
medwidth : int
number of columns to average when fitting profiles to echelleograms
plotalot : bool
Show pretty plot? If running in batch mode (using ordlocs) default is False; if running in interactive mode (ordlocs is None) default is True.
xylims : 4-sequence
Extract the given subset of the data array: [xmin, xmax, ymin, ymax]
uncertainties : str OR 2D array

full path and filename to a 2D uncertainties FITS file, _OR_ a 2D numpy array representing such a file.

If this is set, ‘g’ and ‘rn’ below are ignored. This is useful if, e.g., you are analyzing data which have already been sky-subtracted, nodded on slit, or otherwise altered. But note that this must be the same size as the input data!

g : scalar > 0
Detector gain, electrons per ADU (for setting uncertainties)
rn : scalar > 0
Detector read noise, electrons (for setting uncertainties)
retsnr : bool
If true, also return the computed S/N of the position fit at each stepped location.
retfits : bool
If true, also return the X,Y positions at each stepped location.
RETURNS:
(nord, pord) shaped numpy array representing the polynomial

coefficients for each order (suitable for use with np.polyval)

NOTES:

If tracing fails, a common reason can be that fitwidth is too small. Try increasing it!

spec.wavelengthMatch(spectrum, wtemplate, template, etemplate, guess=None, nthread=1)[source]

Determine dispersion solution for a spectrum, from a template.

INPUTS:
spectrum : 1D sequence

Spectrum for which a wavelength solution is desired.

wtemplate : 1D sequence

Known wavelength grid of a template spectrum.

template : 1D sequence

Flux (e.g.) levels of the template spectrum with known wavelength solution.

etemplate : 1D sequence

Uncertainties on the template values. This can be important in a weighted fit!

OPTIONS:
guess : sequence

Initial guess for the wavelength solution. This is very helpful, if you have it! The guess should be a sequence containing the set of Chebychev polynomial coefficients, followed by a scale factor and DC offset (to help in scaling the template).

If guess is None, attempt to fit a simple linear dispersion relation.

order : int > 0

NOT YET IMPLEMENTED! BUT EVENTUALLY: if guess is None, this sets the polynomial order of the wavelength solution.

FOR THE MOMENT: if guess is None, return a simple linear solution. This is likely to fail entirely for strongly nonlinear dispersion solutions or poorly mismatched template and spectrum.

nthread : int > 0

Number of processors to use for MCMC searching.

RETURNS:

(wavelength, wavelength_polynomial_coefficients, full_parameter_set)

NOTES:

This implementation uses a rather crude MCMC sampling approach to sample parameters space and ‘home in’ on better solutions. There is probably a way to do this that is both faster and more optimal...

Note that if ‘spectrum’ and ‘template’ are of different lengths, the longer one will be trimmed at the end to make the lengths match.

REQUIREMENTS:

emcee

spec.zenith(ra, dec, ha, obs)[source]

Compute zenith angle (in degrees) for an observation.

INPUTS:
ra : str

Right Ascension of target, in format: HH:MM:SS.SS

dec : str

Declination of target, in format: +ddd:mm:ss

ha : str

Hour Angle of target, in format: +HH:MM:SS.SS

obs : str

Name of observatory site (keck, irtf, lick, lapalma, ctio, andersonmesa, mtgraham, kpno) or a 3-tuple containing (longitude_string, latitude_string, elevation_in_meters)

OUTPUTS:

Zenith angle, in degrees, for the specified observation

REQUIREMENTS:

Aperture and simple PSF-fitting photometrymem

Numpy

Previous topic

Observing Tools

Next topic

NIRSPEC data reduction

This Page