Investigate differences between user and pipeline match algorithm in skymatch

[stscijgbot-jp]

Issue JP-3843 was created on JIRA by Ned Molter:

Internal pipeline user Nicolas Flagey ran the skymatch step on the data mentioned above using the "match" algorithm, and compared it to their custom implementation of the same algorithm, which is in the attached Jupyter notebook.  Although the two are supposedly doing the same thing, the Jupyter notebook version seemed to outperform the pipeline version, with the pipeline version showing much larger mismatches in the sky level between exposures than the notebook version.  This remained the case even after persistence and 1/F noise were handled prior to skymatch; the cleaned version of the data is what is in the box folder.  This also remained the case playing with all the various input parameters to skymatch.

It should be investigated what is going on here.

• Are they indeed the same algorithm?
• Can this dataset be processed successfully using the existing pipeline version by setting parameters in a way that Nicolas didn't already try?
• Is there a bug in the pipeline version?
• If not a bug, should improvements be made to the pipeline version?

[emolter]

@mcara What is your take on this?

[stscijgbot-jp]

Comment by Nicolas Flagey on JIRA:

Mosaic obtained with pipeline without our own background matching: https://stsci.box.com/s/h486zcatdcrrq4i5bh3t5eew0lxqgfbq

Mosaic obtained with pipeline after our own background matching: https://stsci.box.com/s/3sillvgnnblobxqik7hrmals06g7n8ez

These two mosaics only use the NIRCam B3 and B4 cals but if needed, I can also send you the links for the entire field.

[stscijgbot-jp]

Comment by David Law on JIRA:

Looks like I can't access any of those box links?

[stscijgbot-jp]

Comment by Nicolas Flagey on JIRA:

Can you access this folder: https://stsci.box.com/s/zj8te9xb1447gjuqovgabqmpj01gagt7 ?

And then check for the mosaics with a filename that contain B3B4?

[stscijgbot-jp]

Comment by Nicolas Flagey on JIRA:

I send direct links to the mosaics to Ned Molter, David Law and Mihai Cara via email

[stscijgbot-jp]

Comment by David Law on JIRA:

I can't access the folder, but I can now see the individual mosaics linked above.

[stscijgbot-jp]

Comment by David Law on JIRA:

Looking at the data I assume it's getting thrown off by the bright diffraction spike from the 4th mag star nearby, but couldn't say why offhand.  Will be useful to hear Mihai Cara 's opinion on this case.

[stscijgbot-jp]

Comment by Nicolas Flagey on JIRA:

David, that might be true, but why is it that the code we run off the pipeline works then? My understanding is that the pipeline is doing the exact same thing but fails. The devil might be in the details here, about how the overlap regions between CALs is computed and what kind of filtering, if any, is done ...

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

It would be very helpful if I could somehow get access to the data and possibly a notebook/script indicating how stage 3 was run (skymatch, tweakreg, resample) so that I could reproduce the original ({}"Mosaic obtained with pipeline without our own background matching"{}) resampled image. 

[stscijgbot-jp]

Comment by Ned Molter on JIRA:

Hi Mihai, these should all be provided already.  See the "data location" field at the top of this Jira ticket for the data in a Box folder, and see the attachment to this Jira ticket for a notebook showing how Nicolas ran the custom match algorithm

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

Nicolas Flagey Please see attached screenshot.

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

Maybe you could put these files on central store with appropriate permissions for me to get access to them?

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

To reproduce Nicolas Flagey results do the following:

Assign a unique "group_id" to {}each image{}.

run skymatch step with the following parameters:

'dqbits': '~0',  # every pixel is good
'nclip': 0,  # disable clipping
'skymethod': 'match',
'skystat': 'median',```
Poor performance of skymatch using default grouping could be caused by [hereafter it is just speculations]: poor calibration of detector's sensitivities or some random variation in these sensitivities (hardware issues???) or, maybe, [Nicolas Flagey](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nflagey) images have been manipulated in a way that resulted in unequal subtraction from image data (for example, maybe that 1/f noise removal procedure causes this, etc), internal reflections, ghosting, etc.

The default grouping assumes that all "chips" in an exposure (nrcb1, nrcb2, nrcb3, nrcb4 + the A side) have been calibrated well such that a common background level is applicable to all chips/detectors. Unfortunately this does not appear to be true for [Nicolas Flagey](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nflagey) images that have been provided to me. Please see the attached image showing the comparison jw02143001001_02101_00001-2.png shows that nrcb* detectors in jw02143001001_02101_00001 have very different levels compared to the same detectors in jw02143001001_02101_00002 (right side) for almost the same "scene". That is, on the right side nrcb detectors are almost uniform while on the left side nrcb2 and 3 appear darker while nrcb1 and 4 appear brighter than the same detectors on right. Similar kind of variations can be seen in jw02143001001_02101_00005-6.png

By "ungrouping" all these detectors, skymatch is allowed to compute a separate sky value for each individual detector (which should not be needed).

Otherwise, at the core, both methods attempt to minimize overall sky differences. The differences between the methods are:
 # skymatch uses spherical_geometry to compute polygon intersections on the sky. [Nicolas Flagey](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nflagey) method resamples one image onto another one and then fills non-overlaps with NaN with the only purpose to find the overlaps between two images. Most likely polygon intersections should be significantly faster than resampling every pixel in every input image for a total of about N**2 / 2 resamples (where N is the total number of images).
 # skymatch solves a system of linear equations that one gets once analytical minimization is performed while [Nicolas Flagey](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nflagey)  uses non-linear optimization. Probably the former is numerically more stable and accurate (accuracy differences should be negligible though).
 # The method in this notebook is median and there is no clipping. Skymatch offers other methods + clipping.

While the differences listed above are mostly affecting performance and "versatility", the method shown here is not incorrect but it also does not add anything new.

I will add a link to the fits file created with skymatch & resample from the pipeline using non-grouped detectors for you to compare with images from [Nicolas Flagey](https://jira.stsci.edu/secure/ViewProfile.jspa?name=nflagey) .

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

Link to the mosaic obtained with skymatch using independent treatment of all detectors: https://stsci.box.com/s/rfyhcrb2k4efi5cvpkpah5p7rxfxdzmx

[bhilbert4]

This is really useful information. I didn't realize that skymatch acts like some of the other stage 3 steps and assumes that the sky signal is the same between detectors in a single exposure. I wasn't able to find this in the skymatch documentation. I think it would be very valuable to add it, along with a description of how to force skymatch to treat each detector independently. I've seen multiple help desk tickets where this strategy probably would have helped.

[stscijgbot-jp]

Comment by Ned Molter on JIRA:

Thanks so much Mihai, this is indeed a really helpful analysis!  I'm glad it sounds like algorithmic changes are not necessary in the skymatch step.  Hopefully Nicolas can confirm that.

I like Bryan's suggestion to update the docs with a note about this.  I'm happy to take that on, unless you would prefer to add it yourself Mihai.  Let's keep this ticket open to cover that work.

As far as why, for a nearly identical scene, the same detector has a different background level, it would be nice to figure out but it's really a separate issue.  I would suggest that we check whether that was still the case before the custom 1/f handling, and if it was, then open a separate ticket.

Does that sound good to everyone? David Law Mihai Cara Nicolas Flagey Bryan Hilbert 

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

This is from actual docs:

.. note::
    Group ID (``group_id``) is used by both ``tweakreg`` and ``skymatch`` steps and so modifying it for one step will affect the results in another step. If it is desirable to apply different grouping strategies to the ``tweakreg`` and ``skymatch`` steps, one may need to run each step individually and provide a different ASN as input to each step. ```
This is How skymatch works for HST too (WFC/WF).

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

https://jwst-pipeline.readthedocs.io/en/latest/jwst/tweakreg/README.html#grouping

[bhilbert4]

Ah ok. It's on the tweakreg page. But for someone coming to the skymatch page, I don't see any mention of it. And I'm not sure we should rely on users knowing to then go look at the tweakreg page. I may need more coffee this morning, but that didn't occur to me.

[emolter]

The note is also not written in a way that would encourage a user to try completely un-grouping all the files and running skymatch that way. The suggestion to "modify the grouping" if it's "desirable to apply different grouping strategies" doesn't give any suggestions of when a modified grouping might be desirable nor what kind of grouping to perhaps try.

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

Sure, docs can be added to skymatch as well. A note for programmers, in the code, SkyImage and SkyGroup classes exist for this reason. I'll add a note to the docs.

[emolter]

@penaguerrero I think the helpdesk ticket could probably be closed too

[stscijgbot-jp]

Comment by David Law on JIRA:

Ned Molter That sounds good to me.  I'll take a quick look at the individual files and see if there's an obvious reason for the sky level offsets coming out of the pipeline.  We can treat that as a separate issue.

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

So this was a help desk issue?

[emolter]

originally yes, but we upgraded it to a JIRA ticket because of the possibility it could be a bug

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

A bug is always a possibility

[emolter]

yes, glad you figured out the difference in this case though, and it doesn't appear to be related to the code itself. Or have I misunderstood that?

[mcara]

yes, glad you figured out the difference in this case though, and it doesn't appear to be related to the code itself. Or have I misunderstood that?

Is this a question for me? If yes, "glad you figured out the difference" what difference are we talking about? I may have missed this in previous messages.

[stscijgbot-jp]

Comment by David Law on JIRA:

Looking at the rate files I don't think this is an issue with 1/f correction; there are large offsets even in the rate files on MAST.  1/f correction unrelatedly looks like it's doing a good job of correcting the striping, and doesn't impact the pedestal difference much.

example.png attached shows the 4 nrca MAST rate images for jw02143001001_02101_00006

nrca4 typical counts in a source-free region are about 0.020 DN/s while nrca2 is about 0.013 DN/s median.  My guess would be stray light coming off the very bright star nearby, but I defer to Bryan Hilbert on whether this is normal for NIRCam or if a new ticket needs to be opened to investigate further.

[emolter]

@mcara yes it's a question for you. The "difference" I'm talking about was the original ticket itself, there was an apparent difference between the results of the pipeline skymatch step and the user skymatch implementation. But it sounds like the "difference" was related to the way the files were being input and grouped, and not a problem with the pipeline skymatch itself. Can you please confirm whether this is a true statement?

[mcara]

@emolter OK. I will confirm that this is not a bug in the code.

[mcara]

Actually, I should reword my answer: while I cannot guarantee lack of bugs in the code, at this moment I believe the observed difference between standard pipeline processing and the alternative user code is due to mainly issues with data and due to default grouping by the pipeline.

[stscijgbot-jp]

Comment by Nicolas Flagey on JIRA:

Thank you all for investigating the images!

I will look into more details at the parameters that Mihai Cara changed in the pipeline. I think there is a path for my team to stick to the pipeline and simply ungroup the detectors but still use the sigma-clipping and other advantages that Mihai Cara listed in their initial comparison.

A couple of things that I wanted to highlight though:

• the CAL images that I shared with you have been through 2 steps: a 1/f correction and a persistence correction. For the former, we used Chris Willott's code. For the persistence correction, it is important to understand our field is quite faint and diffuse, and observed through a narrow filter. Right before our observation, a very bright field was pointed at. We managed to remove most of the persistence effect but it is likely that it affected the "relative background level" of each detector.
• however, we saw the mismatch in background between detectors even using the CAL images from MAST, without 1/f and persistence correction.

Again, in the end, it looks like the best path forward (for my science team) is to go back to the pipeline and assign unique group_id for each image while keeping all other parameters nominal.

Would that be a fair statement?

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

??Would that be a fair statement???

I can't really answer this question as everyone should decide for themselves what's best for their science projects. Fundamentally, there is not much difference between your method or the ungrouped pipeline method.

For the median or mode statistics, clipping has either no effect or minimal effect but again it depends on the scene. For example, in a field with mostly a constant background (+ noise) and some stars, background pixels will dominate in the median and it will be affected minimally by clipping. The same goes for mode since the peak of the mode will be centered around the background level and that peak is not affected by clipping the wings of the histogram. In a complicated scene, there could be multiple modes, or a poorly defined peak, etc. - clipping may matter but most likely mode would not be a good choice. Again, up to you to decide what is a good statistics.

Some additional features of the pipeline code are: limiting min/max range of "good" data for background, use of DQ flags. If you don't need these - your code is perfectly fine.

[stscijgbot-jp]

Comment by Mihai Cara on JIRA:

Documentation updated in #9080