How to keep FITS header when stacking?

[Luca]

Hello everyone, I am using Astroimagej for the photometry utilities.

The photometry tool uses parameters written inside the FITS header to find error bars and to correctly set the time at which each picture has been taken. Sometimes though, we need to STACK multiple images to reduce the signal to noise ratio, and this can be done with the "Grouped Z Project" tool under Image\Stacks\Tools.

Problem is, once the stacking is done (example: reduce a 700 image stack into a 350 stack by summing 2 subsequent images) there's nothing in the FITS header anymore. This makes sense because the stacking might change some of the parameters written inside, but it makes impossible to do photometry with the new stack. And there's no way, as far as I know, to make a new FITS header with the informations we need (by the way, the Fits_header_standalone_editor plugin doesn't work). I guess the solution would be to make a macro that automatically adjusts the time of exposure by summing the one found in each image but there's also no way to make that happen since there's no command as far as I know that enables to automatically edit the FITS header.

Any solution?

[karenacollins]

Hi Luca,

When images are combined using underlying ImageJ functionality (i.e. the "Grouped Z Project" feature), the FITS header information is not copied forward or adjusted as you have discovered. Maybe someday I will find a way to work around that, but for now we are stuck with the current functionality.

Since you are doing photometry, another approach is to run photometry on the full (unstacked) image set. Then in Multi_Plot, use the "bin" feature to combine the photometric data to reduce the light curve scatter. The "bin" feature not only combines the photometry, but also combines the time of the individual photometric points.

You can also use the fits header editor:

to add back a hand-calculated EXPTIME and DATE-OBS keyword and associated values, but that would likely not be reasonable here if you are processing 10's or 100's of images. If you do this, AIJ will add all of the other appropriate file description related keywords when the file is saved (i.e. # 1-8 below).

Karen

[Luca]

Thanks for the reply!

The bin funtionality as far as I know changes the plot only, its not a true sum of the images so the error is not really reduced. That being said, it would be a great improvement for photometry to get that fits header upgraded with z grouped project.

[karenacollins]

Hi Luca,

The error data is indeed recalculated as part of the binning and would be reduced when more points are combined (you can see it when you plot the error bars). You would need to save the binned data back into the measurements table to see the binned time, flux, and error numerically. Use the "New Col" button to save the binned data as new columns in the measurements table when the data are plotted and binned as you want them:

Karen

[Luca]

Alright then I guess I was wrong after all!

Since you're here, do you mind me asking why the scintillation noise hasn't been added to the error formula? Most of the time the exposure time has to be low to avoid things like saturation and sky background error accumulating, so it ends up being the most important form of error there. I've been taking photometry images with a 50cm (effective) diameter telescope, and the error bars calculated by astroimagej are less than half the size of the scatter. Scintillation is also "easy" to estimate since astroimagej already has all the parameters required: aperture, exposure time, altitude and star coordinates.

[karenacollins]

Hi Luca,

I've seen formulae for estimating scintillation noise for single aperture photometry, but I'm not sure how to reliably
do it for differential photometry. My thinking is that if comparison stars are near the target star, the scintillation
noise will be less of an effect, and if the comparison stars are far from the target stars, the effect will be greater.
Do you have a reference on how to calculate the effect of scintillation noise in differential photometry? Beyond that,
ground based light curves almost always have correlated (red) noise of some sort too, which can't really be calculated.
When doing modeling using light curves and RV data, MCMC algorithms typically scale errors such that the reduced
Chi-square of the data fit to the model is 1.0, so as long as the relative values of the calculated errors within a
dataset are reasonable correct, the scaled errors will be reasonably good representations of the actual errors
(including scintillation and red noise).

If you are taking short exposures of focused images, inter-pixel sensitivity variations may be the dominate factor in
the scatter, unless you have perfect guiding. If you can defocus slightly and extend the exposure time a bit to 30-60
seconds, the light curve scatter will likely be relatively close to the formal calculated error. If you are observing a
crowded field, then obviously this technique should be avoided.

All of that being said, I have considered adding the single aperture scintillation noise calculations to AIJ and then
just propagating those through the differential photometry error calculations, but I just haven't had time, and it would
likely over-estimate the errors for high-quality differential photometry. Maybe someday I can get this done when I have
spare time after my "day job".

Karen

[Luca]

Thank you for answering! I have a couple of questions.

The reduced chi-square for my observation is 13, way higher than 1.. so it's obvious that the error bars are underestimated but it's also true that the MCMC algorithm did not scale them. What you're saying is that the relative values of the calculated errors are not correct? Wouldn't an option to manually scale them be very useful?

I'm also curious about why you think that the scintillation noise would change based on the proximity to the target star? I thought it was supposed to depend on airmass, telescope effective diameter and exposure time only. By the way: https://arxiv.org/pdf/0903.2139.pdf (eq. 13) they treat the single aperture scintillation noise as the other noises and don't propagate because the error is similar when using many comparison stars (or so they say, at least in this measure).

Last question, does Astroimagej use this algorithm: https://tinyurl.com/y8dge5gn ?

P.S. Nice way of estimating the red noise: https://arxiv.org/pdf/astro-ph/0606395.pdf (eq 5 and 6).

[karenacollins]

Hi Luca,

Your chi2r does seem quite high. Do you have the CCD noise settings specified properly in Aperture Settings for your CCD?

(these numbers are for one of my CCDs and shouldn't be interpreted as generally correct for all CCDs and gain settings)

Have you calculated by hand a representative scintillation noise to see if it is on the right order of magnitude for the difference in your expected uncertainty? Of course the shorter the exposure time, the more the scintillation component dominates. Following the Southworth paper you referenced, you will have significant gains in photometric precision if your field is not crowded and you can perform defocused photometry.

All that said, I'd be happy to analyze your dataset to understand why the chi2r is so high (and I might learn something along the way). We could start by you forwarding me your measurements table directly (in confidence) to kcollins atat cfa dotdot harvard dotdot edu, if you are interested. If that leaves open questions, we might need to investigate a way for you to transfer your calibrated images.

Regarding MCMC, I am referring to full system modeling tools such as EXOFAST that are typically used by professional astronomers to model the full system dataset (multiple light curves, multiple RV datasets, Doppler tomography, and measured stellar parameters). In those cases, AIJ is used to formally calculate the errors based on the CCD equations specified in the AIJ paper, and then the MCMC in the global fitting tools (not AIJ) should scale the errors as needed. AIJ itself does not implement MCMC at this time.

For now, you can always quickly load the AIJ measurements table into a spreadsheet program as a tab delimited file and scale the uncertainties as needed, if you are not using a separate tool such as EXOFAST to scale the errors.

I am not arguing that AIJ ideally calculates uncertainties, but it is adequate for most transiting exoplanet work.  However, I think it is a good idea to provide an option within AIJ to allow error scaling. I'll do my best to add that feature when I get time. I will also probably provide an option to allow AIJ to calculate errors to include contributions from scintillation. However, it will likely be in the December time frame when I have some time off from my daily job demands to get it finished. To do this properly, I'll need to make a separate panel for the user to enter/verify noise calculation data (aperture size, airmass keyword name in fits header, telescope altitude [which can default to the one shown in Coordinate Converter], and exposure time fits keyword name) {notes to self when implementing}.

To answer your question "I'm also curious about why you think that the scintillation noise would change based on the proximity to the target star?", the point of differential photometry is that the affects of the atmosphere are partially cancelled out, and the closer the target and comparison stars on the sky, the better the cancellation. I don't have data to prove it, but I expect that same principle applies to scintillation noise as well.

Regarding the Southworth reference (https://arxiv.org/pdf/0903.2139.pdf (eq. 13)), my quick-look interpretation is that they are only calculating the photometric error for the target star, because they claim that calculating the error for the target star is good enough and that the contribution from the comp stars does not need to be included:

(after equation 15)

This is indeed not always the case in my experience, and is highly dependent on the field and availability of good comp stars. It is probably as incorrect as AIJ not calculating scintillation noise in many cases.

Regarding the reference to "https://tinyurl.com/y8dge5gn": No, AIJ does not implement that  detrending algorithm. I believe it works best with very long datasets, rather than short single night observations. We use something similar in the KELT data reduction pipeline (i.e. ~5,000-10,000 exposures over many years). You could approximate this in AIJ by detrending using one or more comp star light curves (i.e. rel_flux_Cxx, where xx = 02, 03, ...) as detrend parameters (in addition to airmass, etc).

Karen

[Luca]

You're amazing! Thanks. I'm sending you the measurements file right now, but if you want I can also just upload the aligned fits files tomorrow on dropbox once I'm home.

Yes CCD noise settings are correct, or at least they're taken from the camera sensor manual. We're going to test every feature in September to see if they're somehow different from the average:
CCD gain: 1.5
CCD readout noise: 15
CCD dark current per sec: 0.1

I've actually made a huge Excel spreadsheet which calculates all forms of errors for a single star and comparison stars (with some input), and then finds the best exposure time and defocussing in order to obtain the least amount of error per unit time. In this case each image had 10 sec exposure time, and by my calcs the scintillation noise was dominant. Only considering that, the error bars should be 3x the size that AIJ finds. Then I also made an estimation of the correlated noise (1.3 mmag or so), and that should increase the overall error by quite a bit as well.

P.S. I find that most of the times, with the setup we have, it's better to keep the exposure time below 60 seconds (even more so because more images you get for the transit, less error you have on the calculated parameters), so the scintillation noise is a big issue that doesn't go away with less bright targets.
This is the errors (sigma mag) I find for HD189733b (7.7 visible mag!) using 3 comparison stars, over exposure times:

[karenacollins]

Hi Luca,

That's actually quite a nice light curve. If I bin the exposures by 6 to simulate a minute exposure, your RMS is right
at 1.0 mmag, which is on the better end of ground-based results.

I understand now why your error bars are so underestimated. The photon noise contribution is tiny because the target
star (HD189733) aperture has nearly 10^7 counts. In this extreme case of very high photon count, and short exposures,
the AIJ errors will be way underestimated (as you have discovered) because of scintillation (3x as you have pointed
out), and there are other sources of noise in that regime that will probably be hard to pinpoint.

There are also systematics (the red noise you mention I assume) that are causing excursions from the model, which could
be caused by changes in the sky background, changes in FWHM, or movement of the field on the detector, etc. It's hard
for me to assess the movement on the detector since apparently the images were aligned before running photometry. It
would be best to remove any known sources of systematics using detrending.

If I detrend using airmass, time, sky_background, and FWHM, the fit cleans up a good bit. I wonder if being able to
detrend on x- and y- centroid would resolve some more. Also, I noticed that some of your exposures are very near
saturation. Even if you are not saturating, you are likely into the higher end of the non-linear regime. If your
detector is like most I have worked with, I'd try to keep the peak pixel counts to 40-45K or less (or even 30K if the
star is not so bright).

The bottom line is that I agree that having AIJ handle the scintillation errors would dramatically help in this special
situation. I will indeed add this feature when I can find time to do so. By the way, I understand that making the
calculations should be an easy addition, but the biggest part of the work is tying it into the user interface and
testing to make sure the results are correct and that I haven't broken other stuff in the process.

Karen

[Luca]

I just followed your advice to manually change the error bars through Excel and basically multiplied them by 4 (3 times for the scintillation noise plus 1 for the red noise as calculated through my calc sheet) and found an almost too perfect "1" value for the CHI^2/dof. Now I can definitely say I am happy with the result. Obviously when I tell AstroimageJ to bin the data, it doesn't know that the red noise exists (meaning it doesn't scale the error bars correctly) so the CHI^2/dof value rises, but not that much. Adding an option for red noise (either for calculating it or just adding it manually) would be incredible, just saying.. ahah .

This time I used non-aligned images as you implied and found that using the WIDTH_T1 parameter to detrend improves the BIC value (and the CHI^2/dof) by quite a lot... almost as much as AIRMASS actually**. A bit of improvement also comes from SKY/PIXEL_T1 and ANGLE_T1. The other parameters don't really help (FWHM_mult now is a constant value). Could you point me to a source where the detrending calculations done here are explained? I have a vague idea of what the software is doing but I would like to be sure about that.

Since you were wondering about what filter we used, I am linking it here, combining this with the CCD bandwidth we're basically catching light from 610 to 1000nm or so. By the way, the fact that the depth of the light curve is a couple of mmag low compared to the depths in the literature is due to HD189733 companion, a ~10 v_mag star that is definitely inside the aperture and as you probably know, adds a constant mag value that causes this effect. Doing the math (subtracting the light of the companion) I find that the depth is consistent with the other ones.

Thank you for your incredible and continuous support. I would gladly help to update the software if I could but you have done an amazing job already.

P.S. Yes we had some "saturation alarms" at the end but consider that our CCD is linear up to 60k counts, so it doesn't matter that much. But yes it was a first attempt and we didn't consider that during the night the focusing was going to change (increase in this case) and therefore each pixel got more counts than expected towards the end. **Thinking about it now, since the focusing changed (~sinusoidally) during the observation, might this be the cause of the high detrending effect of the WIDTH_T1 value?

P.P.S. Do you generally apply the detrending parameters to check stars as well?

[karenacollins]

Hi Luca,

Detrending is described in the AIJ paper in section 4.4 "Light Curve Fitting and Detrending". Let me know if you have specific questions after you take a look at that section. If you are fitting and detrending an exoplanet transit light curve using the () fit mode, the model in the chi2 equation (equation 5) is the transit model. If you are detrending a comp star using the () fit mode, a constant model is assumed in the chi2 equation.

Ah, I recall the neighbor star, and indeed that makes sense if it was in your aperture.

I appreciate your attention to detail!

Karen

On 8/15/2017 12:54 PM, Luca [via AstroImageJ] wrote:

I just followed your advice to manually change the error bars through Excel and basically multiplied them by 4 (3 times for the scintillation noise plus 1 for the red noise as calculated through my calc sheet) and found an almost too perfect "1" value for the CHI^2/dof. Now I can definitely say I am happy with the result. Obviously when I tell AstroimageJ to bin the data, it doesn't know that the red noise exists (meaning it doesn't scale the error bars correctly) so the CHI^2/dof value rises, but not that much.

Oh by the way, this time I used non-aligned images as you implied and found that using the WIDTH_T1 parameter to detrend, it improves the BIC value (and the CHI^2/dof) by a lot... almost as much as Airmass actually. A bit of improvement also comes from SKY/PIXEL_T1 and ANGLE_T1. The other parameters don't really help (FWHM_mult now is a constant value). Could you point me to a source where the detrending calculations done here are explained? I have a vague idea of what the software is doing but I would like to be sure about that.

Since you were wondering about the filter we used, I am linking it here, combining this with the CCD bandwidth we're basically catching light from 610 to 1000nm or so. By the way, the fact that the depth of the light curve is a couple of mmag low compared to the depths in the literature is due to HD189733 companion, a ~10 v_mag star that is definitely inside the aperture and as you probably know the companion adds a constant mag value that causes this effect. Doing the math (subtracting the light of the companion) the depth is consistent with the other ones.

Thank you for your incredible and continuous support. I would gladly help to update the software if I could but you have done an amazing job already.

P.S. Yes we had some "saturation alarms" at the end but consider that our CCD is linear up to 60k counts, so it doesn't matter that much. But yes it was a first attempt and we didn't consider that during the night the focusing was going to change (increase) and therefore each pixel got more pixels than expected going to the end.

---

If you reply to this email, your message will be added to the discussion below:

http://astroimagej.1065399.n5.nabble.com/How-to-keep-FITS-header-when-stacking-tp725p749.html

To start a new topic under AstroImageJ, email [hidden email]
To unsubscribe from AstroImageJ, click here.
NAML

[Luca]

Alright, thanks!



This is how the focus changed during the night, so basically what WIDTH_T1 detrending does is linearly correct the change in FWHM size.
If so, the linear detrend did wonders even with that messy focus trend, implying that this might be a big source of error in our data.