overzealous outlier detection for NIRSpec IFU

[stscijgbot-jp]

Issue JP-2050 was created on JIRA by James Muzerolle:

In the course of processing simulated IFU exposures of a dithered point source, we're encountering problems with the outlier detection step.  A majority of usable pixels are being flagged as outliers, despite the source being bright and no cosmic rays were included.  As a result, the output combined cube is mostly zero.  For the full set of 4 dithered exposures, cube_build actually crashes because there is no unflagged data.  We used only the default parameters for outlier_detection; it's not clear if this is just a matter of tuning them, or if there's a bug in the code.

[jdavies-st]

The median is not standard output, just optional, right?

So probably not super important, but it should probably follow the same pattern as for general spectral products in level 3. There's a product name, and an association name that get folded into the final name.

So the following NIRCam imaging _cal files

jw42424001001_01101_0000?_nrc?5_cal.fits

get combined into a final

jw42424-o002_t001_nircam_clear-f444w_i2d.fits

which is a product name generated by the association.

    "products": [{
        "name": "jw42424-o002_t001_nircam_clear-f444w",
        "members": [{
                "expname": "jw42424001001_01101_00001_nrca5_cal.fits",
                "exptype": "science"
            },
            {
                "expname": "jw42424001001_01101_00001_nrcb5_cal.fits",
                "exptype": "science"
            },
            {
                "expname": "jw42424001001_01101_00002_nrca5_cal.fits",
                "exptype": "science"
            },
            {
                "expname": "jw42424001001_01101_00002_nrcb5_cal.fits",
                "exptype": "science"
            },
            {
                "expname": "jw42424001001_01101_00003_nrca5_cal.fits",
                "exptype": "science"
            },
            {
                "expname": "jw42424001001_01101_00003_nrcb5_cal.fits",
                "exptype": "science"
            }
        ]
    }]
}

And ideally the median output name would also be

jw42424-o002_t001_nircam_clear-f444w_median.fits

indicating it is a median made from combined images in that same association. Though maybe the things that get recorded by default in the product name are probably different for MIRI MRS than for NIRCam imaging? For NIRCam imaging it is the pupil and filter (clear and f444w).

[stscijgbot-jp]

Comment by James Davies [X] on JIRA:

This looks like simulated data? Are the error arrays in the rate file realistic?

Would be useful to have the _rate and _cal files in the above data dump as well.

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

James Davies [X] cal files are now in the directory, along with one of the original count rate inputs.  I don't know the details of what exactly went into the simulations, but it looks like the input err arrays contain some estimate of the photon noise and possibly read noise.  If anything, they are probably underestimated, mostly ~0.1%, so that doesn't jive with the outlier flagging.  Another possibility might be that the pointing is slightly off, so that the PSF pattern doesn't exactly overlap from dither to dither; however, I did a run with 2 dithers skipping outlier rejection, and the combined product qualitatively looks reasonable.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

Looking at the values in the SCI and ERR arrays of the cal files indicates an SNR ~1000, which seems unrealistically high for data with relatively low signal levels. I found an old IFU exposure with actual lamp signal in it (from CV3?) and reprocessed from scratch, so that all of the variance and error arrays would be recomputed properly. This exposure, with LOTS of signal from the lamp, yields SCI/ERR ratios in the cal file ~100. So if high levels of signal give you SNR~100, this simulation of low to moderate signal with SNR~1000 is clearly unrealistic. So I'm guessing the overzealous CR flagging is due to the fact that the use of the low ERR values is making everything look like an outlier. I'm rerunning outlier_detection using much higher values of params like "snr" and "scale", to see if this helps to make up for the unrealistically low ERR values.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

James Muzerolle could you put the original rate (cts) files into the data directory listed above, so that I can rerun them through calwebb_spec2 and try to generate more realistic error estimates? I'm sure that this problem with outlier_detection is simply due to the fact that the existing ERR values in the cal files are way too low. They imply a SNR~1000, whereas when I use the 1 cts file that had been provided and manually calculate VAR_RNOISE and VAR_POISSON arrays (using the same formulae as in the ramp_fit step) and then run through calwebb_spec2, which also then computes VAR_FLAT, I get cal files with implied SNR more in the 30-40 range.

Given the existing (very low) ERR values, I have to set the outlier_detection snr parameter to ~1000 in order to prevent it from flagging everything in sight. The detection threshold uses the ERR values as the noise estimate, so if those aren't realistic, you don't get reliable outlier detection.

P.S. the outlier_detection code seems to work fine on simulated MIRI MRS data, because it has more realistic ERR estimates.

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Howard Bushouse done.

Do you happen to have a script for the manual noise term calculations?  We'll need to do the same thing for our other simulations, since all of them start at level 2a.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

Howard Bushouse can you remind me how outlier detection is supposed to work when we cover more than 1 band of data. This NIRSpec data covers NRS1 and NRS2. I thought it worked on each band individually. So take all the NRS1 data - build single type cubes so they cover the same wavelength range and then build a median cube from these singles cubes  again covering the same wavelength range. When the median cube is blotted back it covers only the range of wavelengths on NRS1. But when I run calspec3 on this data it creates single type cubes covering the wavelength range of NRS1 and NRS2. I am not sure how the median cube would be made from this data because 1/2 of single type cubes have 0's in the non covered wavelength. Then I am not sure how the median cube covering the wavelength range of NRS1 and NRS2 should be blotted back to just one detector. I am a little confused how this is supposed to work.  Is the way calspec3 is working on this set set up correctly - or should it split outlier detection into at least 2 parts - one for each detector ?

The same is true about MIRI data that cover all 12 bands. I think I have some Miri data covering all 12 bands I will see how it does outlier detection. 

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

Not sure how it was supposed to be doing things for NIRSpec IFU, but from what you describe it does in fact sound like it should be working on only 1 band of data at a time, as it does for MRS. There's no point in having the median constructed from the entire wavelength range, since only half the data contribute, which will affect the median. So I'm guessing the NIRSpec side needs some "hooks" like the MRS side, so that outlier_detection works on data from just one detector (wavelength range) at a time. The associations are already setup to restrict the data to just one grating/filter, so you shouldn't need to worry about that domain.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

Looking at outlier_detection_ifu.py the hooks for MIRI are setup on channel. I am running a test case of 12 bands of data and it works on the channels separately - building a median for each channel. For NIRSpec the hook is on grating. But that is not fine enough it should be grating and detector because the grating is the G140H so it covers both detectors

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

David Law

This ticket is about a NIRSpec outlier detection problem but I am not sure that outlier detection is working correctly for MIRI.  You have more experience than I do with testing MIRI MRS and outlier detection. So maybe you can set me straight. I have an association that covers the full wavelength of MIRI. I have 3 dithers of each data - so I have a total of 18 files. To be very clear 3 files for: IFULONG_LONG, IFULONG_MEDIUM, IFULONG_SHORT, IFUSHORT_LONG, IFUSHORT_MEDIUM, IFUSHORT_SHORT.

The first step is to create single type ifu cubes. These have a size of the entire 18 exposures on the sky - but are for each band. So with 18 files I get 36 single IFU cubes. Next a median cube is created that will be blotted back to the detector. MY CONFUSSION is that currently there is a median cube created for each channel. CHANNEL ? So for channel 1 it is a median of channel 1 Short, Channel 1 Medium and Channel 1 Long. Should it not be a median of each band - so a median of channel 1 short, another one for channel 1 medium and a third for channel 1 long. Currently I just have 4 median fits files. 

The output names from using --steps.outlier_detection.save_intermiedate_results = TRUE are just plain confusing.  I really think this needs improvement. For example I get a file 

det_image_seq1_MIRIFUSHORT_12MEDIUMexp1_cal_ch1-mediumshortlong-_single_a3001__outlier_s3d.fits.fits

 

I find it confusing that mediumshortlong is in the name. I think that is not correct. Single IFU cube have just 1 band, so either I am very confused or those extra bands in the name really should not be there. 

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

David Law I tracked down 1 error for MIRI MRS outlier detection. The set of input_images needs to be consistent between outlier_detection and blot_cube. Blot_cube was redefining the list of images to be only those contained on the channel that is being blotted. This results in incorrect data arrays being added together in outlier_detection. I want to fix some of the output name strangeness and also FIX NIRSPEC and then I will put in a pr for you to look at.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

Howard Bushouse  I want to update how the names for some intermediate files  are formed in outlier_detection. The  name of median of the single IFU cubes is  formed by the first name in the association. For MIRI a median image is formed for each channel. If the association contains multiple channels then the output median name is confusing. For example - I have an association for containing all four channels. The output medians names are: 

det_image_seq1_MIRIFULONG_34MEDIUMexp1_a3001_band1_median.fits - this is for channel 1

det_image_seq1_MIRIFULONG_34MEDIUMexp1_a3001_band2_median.fits- this is for channel 2

det_image_seq1_MIRIFULONG_34MEDIUMexp1_a3001_band3_median.fits - this is for channel 3

det_image_seq1_MIRIFULONG_34MEDIUMexp1_a3001_band4_median.fits - this is for channel 4

Is there some requirement that the median contain the initial part of the first file. Or can I change to just a3001_ch1_median.fits .... a3001_ch4_median.fits - where a3001 is in the association.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

I agree with James Davies [X] that the file names of the intermediate products (single s3d cubes, median cube, blotted 2d images) should at least contain some of the original programmatic info from the input cal files. Remember that the names being used in your test ("det_image_seq1...") are totally bogus relative what actually gets used in normal ops processing. Normal input cal file names for MIRI MRS are going to be something like "jw00623036001_02105_00001_mirifushort_cal.fits", where the "02105" and "00001" fields are going to vary across the multiple inputs that go into a given median cube. Meanwhile, the output product root name defined in the spec3 ASN file will be a succinct value like "jw00623-o036_t008_miri" (no additional detector, filter, pupil, channel, band info, because we know the final spec3 product will cover multiple settings of all of those). So you could start with the ASN product root name and append something like "ch[n]_median.fits". Will a given median cube also be restricted to one band or could it cover multiple bands for the given channel? If restricted to one band, you could indicate both the channel and band in the appended suffix.

It probably wouldn't hurt to also include one of our typical product type indicators in the suffix, like "s3d" for the median cubes, because they are 3D spectral products, so "ch[n]_median_s3d.fits" for the full suffix.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

David Law

General question. Outlier detection is currently organizing data by channel. It gathers all the exposures according to the channel number and first creates the "single" mode ifu cubes. Let me explain this using an example. If we have data for channel 1 and channel 2 that covers sub channels short, medium and long. The outlier detection would first gather all the data from channel 1 - the size of the single mode IFU cube would be that covered by channel 1, sub-channel short, medium, long. But when building 'single mode' ifus for outlier detection each exposure of channel 1 would create a separate single IFU file. Then all the channel 1 single IFU cubes are combined together to give us the median image that will be blotted back to each exposure and used for flagging outliers.  At the edges between short- medium- long the wavelengths overlap a little. Do you think this is beneficial - when creating the medium images if some wavelengths are a median of 2 bands (will it help flag outliers) or since the resolution is just slightly different - will it create a problem of flagging too many outliers ?

[stscijgbot-jp]

Comment by David Law on JIRA:

Hm, tricky.  I'm inclined to say that outlier detection should work only on a per-band basis.  I.e., the only information used to determined whether to flag things as an outlier in 1A should be other 1A data.  That also avoid the complication of what to do when using crossed-band observations, in which case you really don't want to be combining bands to identify things to reject.  Requiring the same band/channel for outlier detection seems like a clean solution to avoid potentially introducing problems.

Things like residual background matching should certainly look across bands to ensure continuity, but that's a different issue.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

James Muzerolle  In working on this ticket I have found an oversight in both cube_build and outlier detection. IFU cubes are built  in the spec3 pipeline by combining data from the same band - which for NIRSpec we defined as the same grating and same filter. However we did not select on detector.  The data that was for this ticket contains NRS1, NRS2 exposures for grating G140H and filter F100LP.  The default IFU cube for spec3 is one IFU cube containing NRS1 and NRS2.  Is that what is desired ? Or should it be two IFU cubes covering the wavelength region covered by each detector ? 

For outlier detection to we make 'single' type IFU cubes. These are IFU cubes for a single band. Each exposure for the band is mapped to the entire sky covered by the data but the single IFU only contains the data from a single exposure. All the single IFUs are combined together to create a median IFU cube. This median IFU cube is mapped back (blotted)  to the detector space to flag outliers. The problem is that since we only look at the grating and the filter the 3D region the single IFUs cover include both NRS1 and NRS2. This results in a problem. When we take each exposure and mapped it to this single 3d IFU region the NRS1 exposures have data for the ~first 1/2 of the wavelengths while the NRS2 exposures have data for the ~second 1/2 of the wavelengths. Combining all the single IFU cubes into a median image is not correct. So we need to create single IFU cube for NIRSPec also looking at the detector. As I am making this change I wanted to see if we need to also look at the detector when making the default spec3 IFU cubes . I think we do - but I wanted to confirm with you.

[stscijgbot-jp]

Comment by David Law on JIRA:

Jumping in with an opinion quickly, though of course I defer to James on NIRSpec issues.  It seems overly complicated for cube building to have to worry about splitting data coming from NRS1 and NRS2, especially since the wavelength transition between the detectors changes as a function of location in the FOV.  Adding such a split would introduce quite a bit of extra complexity in the files that needed to be handled in many different places.

Rather, I would think that the median combination of the single IFU cubes should be able to handle cases like the one described above.  It shouldn't matter if NRS1 exposures only have data for the first 1/2 of the wavelength range, as these spaxels should be flagged with NODATA (or the equivalent) and ignored by the median routine.  The median combination shouldn't just bring in a zero when there is no data (or bad data), it should be respecting the DQ cube in determining what goes into the median for each spaxel.

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Jane Morrison we do want a single cube including data from both detectors, when applicable - note that only this is only relevant for the high-resolution gratings (except for the specific combination G140H+F070LP).  In those cases, there is some overlap in the wavelength coverage on both detectors among some of the slices because of the wavelength tilt, so from a resampling standpoint, it's preferable to include input pixels from both detectors. It's also preferable from a science standpoint, so that there is always a single level 3 cube regardless of which disperser was used.

I'm not sure I understand the problem with outlier detection; is it that you can't determine the mapping back to detector space when both detectors are included? I wouldn't have expected any ambiguity as long as you can tell the IFU slice correspondence for a given output spaxel.

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Just saw David's last comment, I agree with the reasoning in his first paragraph.

[stscijgbot-jp]

Comment by David Law on JIRA:

Interesting- it looks like the DQ array is entirely zero for intermediate-stage outlier detection products.  That would need to be fixed before we could use such information of course.

[stscijgbot-jp]

Comment by James Davies [X] on JIRA:

When making a median cube, the input median cubes need not have fill values of zero. They can have fill values of np.nan where there is no input data maps onto the cube and then np.nanmedian can be used to create the median.

Spaxels themselves cannot actually have DQ values, right? There's no way to map a partial DQ or distribute input DQ values amongst adjacent spaxels.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

James ok great news on the default IFU cubes for NRS1 and NRS2 - so no change to cube build. 

 

On outlier detection - one solution is to enhance the median combination routine. I will look into this. For single mode IFU cubes the DQ array was set to zero - mainly because it was not being used for these special type of cubes and since it was not being used I did not want to waste the time to set up the DQ plane. It takes a rather long time to set the DQ plane for NIRSpec. We could also use the weight map instead. I think it might be a better plane to use to create the median. The DQ plane setting for NIRSpec really needs more testing. I will need to explain sometime in detail what I did so you can critique the algorithm. One MORE THING - if the some of the wavelength slices overlap between the NRS1 and NRS2 - it could confuse the outlier detection algorithm. The blotted image that mapped back to the detector will be a combination of both detections for these wavelengths.  Will this be a problem ? 

[stscijgbot-jp]

Comment by David Law on JIRA:

Edit: Just realized that this comment crossed with Jane's, in which she answers why the extension is zero.

The data cube files have a DQ extension that's 3 dimensional, and contain information about the quality of individual spaxels.  Generally you're right, and it's non-trivial to figure out how to map detector-based DQ values to DQ cubes, but the cubes do at least flag spaxels outside the footprint of the input data with a 513 value {'DO_NOT_USE', 'NON_SCIENCE'}

At least, the cubes that I'm used to looking at do- for some reason the cubes that outlier detection is creating do not contain anything except zeros in the DQ extension.

[stscijgbot-jp]

Comment by David Law on JIRA:

Jane Morrison Using the weight map to select what goes into the median could work too since this array is already being populated.  Anything with non-zero WMAP would be a good starting point.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

James Muzerolle  The issue about having the median cube being made from both NRS1 and NRS2 is not a problem when it comes to mapping (blotting) back to the two detectors. The only issue I think is for those wavelengths that overlap between NRS1 and NRS2. The median cube will be a combination of data from both detectors (which could be a good think since I think the resolution is the same since it is the same grating).  These overlap regions will be backed back to both detectors and outlier detection will then be performed. Are there any differences in the data in these overlap regions  from one detector to another. Are the noise properties between the detectors that same ? If not then it I will just improve the way the median image is created and not worry about separating the single mode IFU by detector. 

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Jane Morrison are you able to back out the IFU slice (i.e., in the slicer plane, not cube wavelength slice) to which a spaxel corresponds? That would break the wavelength degeneracy with detector since the overlap occurs between different slices.

[stscijgbot-jp]

Comment by Jane Morrison on JIRA:

James Muzerolle Ok what am I concerned about is the MEDIAN image for the the wavelengths that are found on both NRS1 and NRS2. This median in this wavelength region  is will be derived using data from both detectors. This median image will be mapped back to every exposure that is was created from (this results in blot images).

There is no problem blotting it to NRS 1 or NRS2 that works fine. The blotting images are used to flag outliers. I just want to make sure that having a median image created using data from 2 different detectors (in the small overlap region)  is not going cause the outlier detection any problems.  I think if the noise properties of one detector was much different than the other one - it might. But It sounds like this may not be the case.  - so I will just go with that and we can test it and see what happens. 

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Jane Morrison ok, I think I understand now. The median cube spaxels with wavelengths in the overlap region were calculated using pixels from both detectors, so there may be some inconsistencies in flagging outliers compared to the rest of the pixels. I suspect this shouldn't be a problem since the detectors have similar total noise (NRS2 is ~15-25% higher, depending on the readout pattern). Also, the number of spaxels involved is relatively small, conservatively maybe 0.1%, and again this is only relevant for the high-resolution gratings.

[stscijgbot-jp]

Comment by Howard Bushouse on JIRA:

James Muzerolle After fixing the blotting algorithm used in outlier_detection for NIRSpec IFU data and computing estimated variances/errors from the original countrate images that I hope is closer to what would be coming out of the ramp_fit step for real exposures, I think we're very close to reasonable results. When using the default outlier_detection snr and scale param values, there are still quite a few pixels flagged as outliers along the centers of the spectral traces in many slices, but I believe this is due to 2 things. 1) the estimated errors are still not quite realistic and 2) possible slight mismatches between the simulated data in each dither position and the assigned pointing information. If I raise the snr and scale param values, the number of flagged pixels decreases accordingly. Also, when doing a close inspection of the blotted image for each exposure against the original sci image, I can often see what appear to be slight offsets in the spectral trace in the areas where a lot of outlier flags occur (the offset causes the difference between the original and blotted data to exceed the detection threshold). This could be due to slight inconsistencies between the pointing information for each exposure and the actual location of the spectral data in the images.

I've attached some images showing the blotted, original sci, and resulting DQ map for part of one exposure. You probably can't see it from these images, but blinking the blotted and sci images shows a small, but distinct, offset in the centroid of the spectral trace in some slices, right where the most pixels are flagged in the DQ array. [Note that in the DQ image the "outlier" flag values show up as grey; the brighter pixels are very high DQ values, such as those pixels affected by stuck open shutters.]

So overall, I think we've done as much as we can do for now to the actual code, based on testing with the simulated data. Also, if the NGROUPS=80 value used in these simulated images is representative of the way real observations will be taken, I don't think you'll need to be relying on the outlier_detection step in level 3 at all to catch remaining CR's. With such a large value of NGROUPS, the jump step in calwebb_detector1 should have no problem at all catching pretty much everything. !blot_image.png|thumbnail! !sci_image.png|thumbnail! !crf_dq_image.png|thumbnail!

[stscijgbot-jp]

Comment by James Muzerolle on JIRA:

Verified using the simulations and high numbers for the snr and scale parameters. We will need to do more analysis using real data to determine the optimum values for these parameters, especially for point sources.