TRIPPy examples

Introduction: SExtractor and emcee

To perform photometry and source subtraction, in addition to having a good PSF (which trippy will generate) one needs three very important parameters: x, y, and m, or source position and amplitude.

When one has the PSF and TSF already generated, one can run a fitting routine to solve for these. For this purpose, we use emcee. emcee is an MCMC routine which allows for good estimates of (x,y,m) and their uncertainties. We use a likelihood definition as the natural log likelihood of the exponential flux, basically exactly what you'd expect. If you are uncertain of what this means, or care for more detail, please go read the emcee documentation.

If the PSF or TSF is not yet known, to get a centroid (x,y), we need to use some other software. We haven't included this inside trippy because there is no point in reinventing a wheel that has already been nearly perfected. For this purpose, we use the venerable SExtractor. All jokes on its name aside, sextractor does exactly what we need, as well as we would ever need it to be done.

Trippy includes a module trippy.scamp with functions defined in scamp.py and makeParFiles.py that mearly provide convenient wrappers to call sextractor. This has been done in a couple other packages, but not in a way that satisfies me. Hence my own implementation. A couple details to note: makeParFiles creates all the parameter files in the working directory (eg. makeParFiles.writeConv()), and scamp is responsible for sextractor execution and catalog reading (scamp.runSex() and scamp.getCatalog). Catalogs are stored in FITS_LDAC format. This choice was done to facilitate execution of the sextractor sister program scamp, though we won't need to know what that means for full use of trippy. If you are unfamiliar with sextractor and its use, don't adopt trippy as a blackbox. RTFM!

With that out of the way, on to actual business.

The trippy tutorial

The first thing to do is import all the necessary packages. Note that this notebook assumes you have the optional packages installed, as well as SExtractor available on your command line.

NOTE: proper use of psfStarChooser requires plot interaction. So for this tutorial you'd best comment out the first line, %matplotlib inline. But for my web presentation, I leave inline.

The function trim catalog is a convenience function to simply return only those sources that are well enough isolated for PSF generation. It rejects any sources within 30 pixels of another source, any sources with peak pixel above 70,000, and any sources that sextractor has flagged for what ever reason. We may fold this into psfStarChooser in the future.

Get the image this tutorial assumes you have. If wget fails then you are likely on a mac, and should just download it manually

First load the fits image and get out the header, data, and exposure time.

Next run sextractor on the images, and use trimCatalog to create a trimmed down list of isolated sources.

makeParFiles handles the creation of all the sextractor files, including the .sex file which we call example.sex, the default.conv, the param file which is saved as def.param.

.runSex creates example.cat which is read by .getCatalog. getCatalog takes as input the catalog name and the parameter file "def.param".

The parameters that are actually used by psfStarChooser and psf.genLookupTable are XWIN_IMAGE, YWIN_IMAGE, FLUX_AUTO, and FLUXERR_AUTO, which are the x,y coordinates, the flux, and the flux uncertainty estimate respectively. The latter two are used in the SNR cut that psfStarChooser makes.

Finally, find the source closest to 811, 4005 which is the bright asteroid, 2006 Polonskaya. Also, set the rate and angle of motion. These were found from JPL horizons. The 1 degree increase is to account for the slight rotation of the image.

Note: in this image, the asteroid is near (4005,811) and we apply a distance sort to the catalog to find correct catalog entry, and the source centroid, which we store in (xt,yt).

Setting the important asteroid parameters. xt,yt contain the location of the asteroid itself (near 811,4005), rate and angle are the rate and angle of traililng, in "/hr and degrees. We find the actual centroid as the location closest to that point.

Now use psfStarChooser to select the PSF stars. The first and second parameters to starChooser are the fitting box width in pixels, and the SNR minimum required for a star to be considered as a potential PSF star.

Optional but important inputs are autoTrim and noVisualSelection. The former, when True, uses bgFinder.fraserMode to attempt to determine what FWHM corresponds to actual stars, and rejects all sources with FWHM outside +-0.5 pixels of the modal value. noVisualSelection determines if manual input is required. When set to false, all stars are considered. Until you know the software, I suggest you use noVisualSelection=True for manual selection, and autoTrim=False to see all sources in the plot window.

For each star provided to psfStarChooser, it will print a line to screen of x,y and best fit alpha, beta, and FWHM of the moffat profile fit.

Then psfStarChooser will pop-up a multipanel window. Top left: histogram of fit chi values. Top right: chi vs. FWHM for each fitted source. Middle right: histogram of FWHM. Bottom right: image display of the currently selected source. Bottom left: Radial profiles of all sources displayed in the top right scatter plot.

The point of this window is to select only good stars for PSF generation, done by zooming to the good sources, and rejecting those that are bad.

Use the zoom tool to select the region containing the stars. In this image, that's a cluser at FWHM~3.5 pixels.

Left and right clicks will select a source, now surrounded by a diamond, displaying the radial profile bottom left, and the actual image bottom right.

Right click will oscillate between accepted source and rejected source (blue and red respectively).

Keyboard funcitonality is now also implemented. Use the left/right arrow keys (or a/d) to cycle through each source, and the up/down keys (or w/d) to mark a source as rejected (red) or accepted (blue). This is probably the fastest way to cycle through sources. Note that for some mac python installs, key presses won't be recognized inside a pylab window. To solve this, invoke your trippy script with pythonw instead of python.

When the window is closed, only those sources shown as blue points, and within the zoom of the top right plot will be used to generate the PSF.

The array goodFits is returned for convenience and contains the moffat fit details of each accepted source. Each entry is [FWHM, chi, alpha, beta, x, y, local background value].

The array goodMeds is just the median of goodFits, and provides the median moffat alpha and beta of the selected stars.

Note on a couple starChooser options:

--bgRadius is the radius outside of which the image background level is sampled. The fitting is relatively insensitive to this value, however, if you happen to know what the FWHM is approximately, then the best fitting results can be had with bgRadius>~3xFWHM in pixels.

--ftol is the least squares fitting tolerance parameter passed to the scipy least sqaures fitter. Increasing this number can result in dramatic performance improvements. Default is 1.4e-8 to provide an extremely accurate fit. Good enough fits can be had with 1.e-7 or even 1.e-6 if one has a need for speed.

--repFact defaults to 5. If you want to run faster but still preserve most accuracy in the fitting procedure, use repFact = 3

--quickFit = True will provide the fastest moffat fitting. The speed improvement over quickFit = False is dramatic, but results in slightly less accurate moffat fit parameters. For the majority of use cases, where the number of good psf stars are more than a few, the degredation in PSF accuracy will not be appreciable because of the fact that a lookup table is used. But the user should confirm this be comparing PSFs generated in both circumstances.

--printStarInfo = True will display an inset in the starchooser plot that shows the parameters of the selected source, such as alpha, beta, and FWHM, among others.

Generate the PSF. We want a 61 pixel wide PSF, adopt a repFactor of 10, and use the mean star fits chosen above.

always use odd values for the dimensions. Even values (eg. 60 instead of 61) result in off centered lookup tables.

Repfactors of 5 and 10 have been tested thoroughly. Larger is pointless, smaller is inaccurate. 5 is faster than 10, 10 is more accurate than 5.

The PSF has to be wide/tall enough to handle the trailing length and the seeing disk. For Polonskaya, the larger is trailing, at ~19"/hr*480s/3600/0.185"/pix = 14 pixels. Choose something a few times larger. Also, stick with odd width PSFs, as the even ones have some funny centroid stuff that I haven't fully sorted out.

The full PSF is created with instantiation, and running both genLookupTable and genPSF.

Now generate the TSF, which we call the line/long PSF interchangeably through the code...

Rate is in units of length/time and pixScale is in units of length/pixel, time and length are in units of your choice. Sanity suggests arcseconds and hours. Then rate in "/hr and pixScale in "/pix. Angle is in degrees counter clockwise from horizontal between +-90 degrees.

This can be rerun to create a TSF with different rate/angle of motion, though keep in mind that the psf class only contains one longPSF (one rate/angle) at any given time.

Now calculate aperture corrections for the PSF and TSF. Store for values of r=1.4*FWHM.

Note that the precision of the aperture correction depends lightly on the sampling from the compute functions. 10 is generally enough to preserve 1% precision in the .roundAperCorr() and lineAperCorr() functions which use linear interpolation to get the value one actually desires.

NOTE: Set useLookupTable=False if one wants to calculate from the moffat profile alone. Generally, not accuarate for small apertures however.

Store the PSF. In TRIPPy v1.0 we introduced a new psf save format which decreases the storage requirements by roughly half, at the cost of increase CPU time when restoring the stored PSF. The difference is that the moffat component of the PSF was originally saved in the fits file's first extension. This is no longer saved, as it's pretty quick to calculate.

Default behaviour is the old PSF format, but the new format can be flagged with psfV2=True as shown below.

If we've already done the above once, we could doing it again by restoring the previously constructed PSF by the following commented out code.

And we could generate a new line psf by recalling .line with a new rate and angle

Now let's do some pill aperture photometry. Instantiate the class, then call the object you created to get photometry of Polonskaya. Again assume repFact=10.

pillPhot takes as input the same coordinates as outputted by sextractor.

First example is of a round star which I have manually taken the coordinates from above. Second example is for the asteroid itself.

New feature! The input radii can either be singletons like in the example below, or a numpy array of radii. If photometry of the same source using multiple radii are needed, the numpy array is much much faster than passing individual singletons.

enableBGselection=True will cause a popup display of the source, in which one can zoom to a section with no background source.

The detault background selection technique is "smart". See bgFinder documentation for what that means. If you want to change this away from 'fraserMode', take a look at the options in bgFinder.

display=True to see the image subsection

r is the radius of the pill, l is the length, a is the angle. Sky radius is the radius of a larger pill aperture. The pixels in this larger aperture, but outside the smaller aperture are ignored. Anything outside the larger pill, but inside +-width is used for background estimation.

Trimbghighpix is mostly made not important if mode=smart. But if you want to use a mean or median for some reason, then this value is used to reject pixels with values trimBGhighPix standard deviations above the mean of the cutout.

The SNR function calculates the SNR of the aperture,as well as provide an estiamte of the magnitude/flux uncertainties. Select useBGstd=True if you wish to use the background noise level instead of sqrt of the background level in your uncertainty estimate. Note: currently, this uncertainty estimate is approximate, good to a few percent. Future improvements will be made to get this a bit more accurate.

If the photometry radius was an array, then so are the products created using the SNR function.

verbose=True puts some nice terminal output in your face. These values can be accessed with their internal names.

Let's get aperture corrections measured directly from a star.

Finally, let's do some PSF source subtraction. This is only possible with emcee and sextractor installed.

First get the cutout. This makes everything faster later. Also, remove the background, just because.

This also provides an example of how to use zscale now built into trippy and astropy.visualization to display an astronomy image using the zscale scaling.

Now instantiate the MCMCfitter class, and then perform the fit. Verbose=False will not put anything to terminal. Setting to true will dump the result of each step. Only good idea if you insist on seeing what's happening. Do you trust black boxes?

Set useLinePSF to True if you are fitting a trailed source, False if a point source.

Set useErrorMap to True if you care to use an estimate of the poisson noise in each pixel during your fit. This produces honest confidence ranges.

I personally like nWalkers=nBurn=nStep=40. To get a reasonable fit however, that's overkill. But to get the best... your mileage will vary.

This will take a while on a computer. ~1 minute on a modern i5 processor, much longer if you computer is a few years old. You can reduce the number of walkers, nBurn and nStep to ~10 each if you are impatient. This will drop the run time by ~4x

Now get the fits results, including best fit and confidence region using the input value. 0.67 for 1-sigma is shown

Finally, lets produce the model best fit image, and perform a subtraction. Plant will plant a fake source with the given input x,y,amplitude into the input data. If returnModel=True, then no source is planted, but the model image that would have been planted is returned.

remove will do the opposite of plant given input data (it actually just calls plant).

Now show the image and the image with model removed for comparison.