Revise how resample handles VAR arrays

[stscijgbot-jp]

Issue JP-3569 was created on JIRA by David Law:

At the moment, the resample step (for both imaging and spectroscopy) is using inverse variance weighting to combine the variance arrays (RNOISE, POISSON, FLAT).  It's doing so by first drizzling each variance from a single exposure to the output frame using uniform weighting, and then doing an inverse-variance weighted combination across all of the exposures to be combined together.  Crucially though, these arrays use themselves as the inverse variance weights.  I.e., for VAR_FLAT, 1/VAR_FLAT is being used as the weighting to combine together the individual exposures to create the final VAR_FLAT.  Likewise, VAR_POISSON is being used as weights for creating the final VAR_POISSON.

In general, for the weighted combination of two images

f3 = (w1 * f1 + w2 * x2)/(w1 + w2)

and the variance V3 of the output image is

V3 = (w1w1V1 + w2w2V2)/((w1+w2)*(w1+w2))

It's only in the special case when weighting by the variance of these inputs that this reduces to

V3 = 1/( (1/V1) + (1/V2))

Generally we don't want to do that as the weights should always be derived from the read noise variance, not whichever variance is being propagated.

We should therefore modify resample.py in two ways:

1. When using ivm weighting, resample_variance_array should be sure to use the read noise variance to weight the combination of the drizzled variances (Poisson, RN, and Flat) from different exposures.

2. When using exposure time weighting, resample_variance_array should actually use the exposure times (e.g., output_model.meta.exposure.duration) for weighting rather doing an inverse variance weight.

This should address the outstanding issue of zero-valued flatfield variances, and allow such values in reference files to propagate as expected to output products.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

David Law - while we are revising error propagation in the resampling step, should we also consider tracking covariance for the propagating the flat errors?  Generally, the same flat is used for all input exposures, so flat variance is not statistically independent for each exposure.

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke It's a thought, but I'm not sure how helpful it would be or the best way to implement it.  If the flatfield errors are stochastic and on scales small compared to the dither offsets then they should average down like normal.  I'd imagine this would be a bigger issue if working with undithered data, or if there are systematic flatfield errors on large scales.  Do we have an estimate of that which you're aware of, or any NIRSpec programs that don't dither around enough to sample different locations?

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

NIRSpec flat field errors include systematic flux calibration errors, which are not small compared to the statistical errors and should not average down. Charles Proffitt might have a better estimate of the scale; we have discussed this recently in the context of flux calibration from Cycle 1 data.

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke Ok, let's separate these issues though.  This ticket should be a fairly straight forward fix, and if you can please file a new ticket on flatfield covariance we can continue the discussion there.  I'll be keen to see what Charles is finding.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks, I filed ticket JP-3575 for further discussion.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

We also recently got a Help Desk inquiry about this issue of resampling the variances, which I'm including here just as additional input. From INC0199428

I have been working on calculating accurate uncertainties for data in JWST 2391 and have run
into an issue with the method of resampling variances in the JWST pipeline in the resample
step with drizzle. I am using JWST version 1.13.4 (the most up to date version) and looking at 
line 436 and 458 in resample.py. Here, the variances are resampled by the drizzle algorithm 
with no weights, similar to how the science data is resampled if you provide --
weight_type='None' to the pipeline. However, if a user provides weights when resampling the 
science data (as it is by default, using IVM weighting), shouldn't the variances be resampled 
the same way to properly propagate the uncertainties? I'm concerned that the variance terms 
provided in the final level 3 data product are therefore inaccurate.

Investigating this further, I also found that the drizzle paper (Frechter & Hook 2002) suggests 
that the variances should be resampled by the square of the weights, including squaring both 
the drizzle weighting kernel and the input weights (such as IVM or exposure time). See 
equation 7 here: https://iopscience.iop.org/article/10.1086/338393#df6. I don't believe the 
JWST pipeline is resampling the variances this way, which adds to my unease of using the 
variances as an accurate estimate of the uncertainties in the data.

Could the help desk let me know if I am mistaken in the logic of resampling variances? I am 
curious as to why the JWST pipeline resamples variances in this way. Further, if the help desk 
has any advice on how to calculate proper uncertainties for my data or future improvements 
to the pipeline that would be greatly appreciated.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Question for David Law -

Reading through the comment Howard Bushouse passed on, and the referenced drizzle paper, it looks to me like the only way to get truly accurate variances would be to modify the drizzle code to accept squared input.  As the user says, the input weights for the variance should match those used for science, and all three of the input weight, pixel scale, and the kernel overlap should be squared for each pixel before adding its contribution.  We could provide squared weights and pixel scale, but the kernel overlap is computed internally in the drizzle code.

The current process is to resample variance components for each input image with uniform weights and a squared scale factor, then to add the resampled variances with a weighted sum.  The proposal here is to modify the weighting used in the final sum: either use the read noise variance or the exposure time weights for each image, as appropriate for the science. 

To get closer to accurate variance propagation, should we also modify the resampling for the variance components? If we resampled them as errors instead of variances, with weights matching the science resampling, the kernel overlap, weight, and scale would not need to be squared. We could then square and sum them as proposed.  I think this would effectively be an overestimate of the variance for each pixel, but the relative pixel variations might be more accurately accounted for.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

I have both options (resampling as variance images with uniform weight, resampling as error images with weights matching science) coded up in a branch on my fork, along with the update to the weights in the sum over images. Draft PR is here:

#8437

Let me know if you want to test one or both options - it's currently set to resample variance arrays in the same way it was previously.

I spot checked both options for NIRSpec FS regression test data, and the resampling-as-error option shows consistently lower overall error values than resampling-as-variance, at the few percent level.  This is true whether using uniform weights for the resampled error image or weights that match science, for ivm and exptime weight_types. Using ivm weights for resampling the error image additionally makes a difference to the error on low-weight pixels: it can be significantly lower for isolated values.

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke Thanks for digging into this!

Repeating back to you my understanding (just to check that we're on the same page): as in the general equation for V3 above, the proper thing to do is use squared weights when combining variances.  Our variance propagation has two steps though; drizzling each image to the output frame (which uses normal non-square weights from the fractional pixel overlaps), and combining the different drizzled images together (which you've now proper weighting for).  Actually doing squared variance weights from the fractional pixel overlaps is hard, since that would be buried deep within the drizzle code, but we might be able to do a bit better than weight-combining the variances in drizzle and instead weight-combining the errors.

Offhand I don't have a good feeling for which would be a better approximation, but it might be possible to determine empirically.

Looks like I should be able to switch between the two methods by replacing resample_one_variance_array with resample_one_variance_array_as_err; let me get back to you.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks David Law - that is an accurate summary. 

I did work out explicitly what the variance would be for an output pixel with two input pixels in two input images in the exact case, the case of drizzling error images with user weights matching science, and the case of drizzling variance images with uniform weights. It's too much to type out here, but I'm happy to discuss offline if you like.

Yes, for testing, you can just replace resample_one_variance_array with resample_one_variance_array_as_err.  Let me know if you run into difficulties!

[stscijgbot-jp]

Comment by David Law on JIRA:

Adding some notes here for posterity summarizing an offline discussion with Melanie Clarke.

First the straight forward: the PR above successfully fixes the issue with zero-variance propagation described here and at https://jira.stsci.edu/browse/JP-3250

I agree that Method 2 above (sending ERR values into drizzle instead of VAR values) gives better performance than Method 1 (sending VAR values into drizzle where they will be weighted by the pixel area overlap instead of the proper square of pixel area overlap).  As a minimum example, consider the case of only two pixel inputs to a given output pixel.  As above, the format proper answer should be

f3 = (w1 * f1 + w2 * x2)/(w1 + w2)

V3 = (w1w1V1 + w2w2V2)/((w1+w2)*(w1+w2))

If f1=1, f2=7, v1=0.5, v2=50, w1=0.99, and w2=0.01 (i.e., pixel 2 is much more uncertain, but the output is almost entirely dominated by pixel 1 area of overlap), then

f3=1.06, v3=0.495

Method1 (sending variances into drizzle) is effectively doing

v3_method1 = (w1V1 + w2V2)/(w1+w2) = 0.995

Method2 (sending ERR into drizzle and squaring them) is effectively doing

v3_method2 = (w1sqrt(V1) + w2sqrt(V2))/(w1+w2))**2 = 0.59

I.e., method2 is high by about 20%, but method1 is high by about 100%.

Given the difficulties of reworking the core drizzle infrastructure, Method2 is clearly the best option.  We should be sure to update readthedocs to reflect this decision so that we can remember what we did in future years.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

David Law - I added some documentation updates to my PR.  Let me know if it's unclear or you think it needs more.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

Fixed by #8437

[stscijgbot-jp]

Comment by David Law on JIRA:

Now working as intended, closing.