I would like if someone could have a look at how I do ensemble photometry using astropy and astropy's afiliated package photutils. I'd like to know if I'm doing anything wrong or if it could be improved.
I have put together a jupyter notebook as an example of photometry of KIC08462852. The link is here > https://nbviewer.jupyter.org/gist/dokeeffe/416e214f134d39db696c7fdac261…
Thanks in advance for any comments.
I assume that you are doing DSLR photometry, and that you do preprocessing to de-Beyer, apply bias, dark and flat field corrections off line, using some other software. If these assumptions are not correct, some of the following comments will be rendered inappropriate.
My first concern is about the comparison star selection process. It appears that by going directly to Simbad you are bypassing the AAVSO's standard comparison star sets. (No big deal.) On the other hand, the selection process does not appear to consider spectral match between the target and comparison stars and check star, and (of lesser concern) the lack of consideration of estimated errors in the photometry in choosing them.
Another concern is the fixed choice of a 6 pixel radius for the measurement aperture, and the apparent lack of a background annulus. My own practice is to look at the profile of the target signature (and, at least initially, of the other stars) to select an aperture that maximizes SNR while avoiding including any extraneous background stars, and ensuring that none of the stars are in saturation. My experience is that a single aperture simply does not work over the range of seeing and defocus conditions encountered in the real world. So, I suggest you implement a visualization routine to allow display of the ADU profiles and SNRs to facilitate proper choice of measurement aperture. For that you will need to implement the background annulus, also variable, to provide a noise estimate.
I notice that you do not correct for, or even calculate, air mass. IMO, that is a serious deficiency.
There is some question in my mind as to whether the residual on the linear regression is the correct measure of error. Perhaps someone else would like to clear this up? (In my own practice, the magnitude is estimated based on a zero residual (i.e., the target is placed on the curve) and the error is given as the standard deviation of the multi-image measurement process.) Speaking of which, it appears that your script is for a single image. Of course, it would not take much to loop over an ensemble of images and use the result to estimate magnitude.
Finally, the report needs to comport to AAVSO requirements. You need to indicate air mass and filter. Results should be reported to three decimal places only to avoid the mis-impression of unrealistic precision.
SFS, Thanks for taking the time to review it. I really appreeciate it.
Yes you are correct, the fits files are already calibrated as per AAVSO guidelines. Its not DSLR though, Its a CCD. I used the AAVSO documentation and ccd course material to calibrate the data as best I can.
Actually I'm only using Simbad to get the target star RA and DEC as I only know the name of the star at the start of the notebook. Once I have the RA and DEC I download the comparison star data from AAVSO and dont use Simbad anymore.
Thanks for the comments on the aperture radius. I need to look into this more and see how I can set a background annulus in photutils. Looking at the profile is best and I'll see how I can include that.
Air mass, yes actually the capture software I use has recently added AIRMASS as a fits header so I can include this in future observations. Although how to correct for it I'm not sure yet.
Regarding looping over several images and using the std dev as the error. This is a great idea! I'll try this as I generally take 5 images of a single source anyway. Well it will work for stars that dont change too quickly I suppose.
Oh yes I manually enter the results into AAVSO then and use 3decimal places.
Regarding filiters, I currently only have a V filter so my data is not transformed. I have ordered a B filter though and hope to begin the process of being able to transform my data too.
Thanks again for your time and comments.
I quite enjoyed what you have done :-) Very nice! I'll add also couple of ideas/suggestions that could be useful for you:
- just for sake of simplicity and reliability, I'd use circular annulus for sky estimation, mode or median value (if you can choose) would be good
- to ensure the best sky measurements, creating frame mask for stars and using that mask for sky annulus measurements should give you the best estimate of local sky - photutils has a function for that
- determine/estimate FWHM automatically (fitting?) and define your aperture radiae based on that (e.g. stellar aperture 1.5xFWHM, sky annulus from 3-5 FWHM, your may like different numbers) - that helps somewhat with varying seeing/trailing
- always estimate uncertainties as well as you can (i.e. include your ensemble members, target(s) etc ), IMHO there was a function for that. It is often said in exact sciences that without uncertainty estimate there is no measurement...
- you may like to convert all magnitudes to unit exposure time (e.g. 1 second): mag+2.5*math.log10(exposure) - then it is easier to compare e.g. results when you have taken frames using multiple exposure times, also to create airmass-instrumental magnitude graphs for multiple nights etc etc
- if your target and ensemble stars do not have dramatically different colours, you could just combine ensemble stars into "one" comparison star, most probably a correct way would be adding up aperture sums (in your tables) and deriving magnitude from that.
I also wanted to ask if I may use that notebook for a student's astronomy course to explain basic photometry?
Thanks for taking the time to read my notebook and comment.
Yes of course you can use it in a course! Glad to help other people learn. Although do take into consideration the other comments in this thread since it can be improved.
Looks like I need to work on my aperture selection instead of just using 6 pixels. Using the FWHM is a good idea and I could do it dynamically then too.
Regarding the colours of the ensemble stars, It looks like this is an area that I could improve here. I currently only have a V filter but soon will also have a B filter. Once I have this and collect data to enable transformation of results I think this will correct for the colours of comparison/ensemble stars. At least that is my understanding.... I could be totally wrong :-)
P.S. If you need the actual fits data file to run the notebook locally let me know and I can send it....
For the past couple of years, I've been using a Python Jupyter notebook that I wrote to perform aperture photometry on CCD time series. Like you, I also use the WCS for each image to determine the positions needed for the photometric apertures. I have a few initial comments, in no particular order.
- I think that the magnitude error is underestimated. The linear regression gives equal weight to the stars, regardless of their SNR. The uncertainty will increase as the SNR decreases, which becomes a problem if your ensemble and target have a wide range of magnitudes. Also, this method of error estimation assumes that the spectral characteristics of the comp stars are similar to those of the target star. For example, if your comp stars are in the range 0.4 < B-V < 0.8 and the target star has B-V = 1.0, the relative redness of the target would probably result in a larger error than the comp-star residuals would indicate.
- I use PythonPhot to do the actual photometry. I never really tried the aperture photometry routines within photutils. One of the nice things about PythonPhot's aperture photometry routine is that it allows you to specify a background-measurement annulus so that your star flux is background-subtracted.It also gives you reliable uncertainties on your star magnitudes.
- If you were to modify your code to analyze a time series, I've found that it is computationally inefficient to call wcs = WCS(fits_header). Instead, I wrote a custom function that creates an empty WCS object and then populates it with the relevant WCS keys from the fits header. This is shown in the Astropy WCS documentation. If I remember correctly, this sped up my for loop by a factor of 5.
- If you already know the equatorial coordinates of your target and comp stars, you don't necessarily need to do the source extraction unless you're worried about proper motion, inaccurate catalog coordinates, or a poor WCS solution. In practice, I rarely encounter these issues. As a safety check, I simply plot the photometric apertures over a sample CCD image to verify that they are positioned correctly.
- I think that a better way of matching sources to target/comp stars would be to eliminate the nested for loop, find the Pythagorean distance to the nearest source, and to use np.argmin() to return the index of the matched source from your sources table. This should be faster than the current method, which I believe could also produce a false match if there are multiple sources within 4 pixels.
Thanks for your suggestions Colin,
Regarding the error, yes I think I'll try some other way and compare. I was thinking of estimating all comparison stars' mags with the ensemble, compare with AAVSO published mag and use the std dev as the error. This may be better for a single frame. Otherwise the stddev of multiple images taken around the same time is another option.
Yes I have done time series with much the same workflow and it is.... slow. Very slow. So thanks for the performance improvement tips.
On the source extraction, the reason I extract the sources and then do the pythagorean match is that I have a feeling its more accurate. The source extraction gives sub-pixel accuracy and would cater for field curvature etc. I think using WCS on its own to get the aperture centroid may be a tiny bit off. Although I have not tested this.. Another thing to check.
Looks like the spectral characteristics of my comparison stars is something I really need to work on as a top priority.
Derek, I am trying to learn photometry with Photutils and have your notebook, but trying to run it, I could not locate
FITS_FILE = '/home/dokeeffe/Pictures/CalibratedLight/2017-10-02/KIC08462852-2017-10-02-20-56-22Light_005.fits'
Did you remove it? could you put it back or send it to me?
Thanks, Richard Post
I've been using your excellent Jupyter notebook as a tutorial on photometry using astropy and photutils. Its very clearly organized and it has been a great learning tool!
In order to get it to work locally, I've been using DSS fits images. I was wonder if I could have a copy of the image you used in order to duplicate the results exactly: