Flux conservation for spectral resampling

[stscijgbot-jp]

Issue JP-3547 was created on JIRA by Melanie Clarke:

I am looking into a couple of potential issues with flux conservation in spectral resampling for NIRSpec.

One, using the pixel_scale or pixel_scale_ratio options for the resample_spec step does not seem to conserve flux in the output.  Testing with MOS data designated as both point source and extended source, the output spectra have noticeably different flux when resample is run with pixel_scale_ratio = 1.0 (default) and with  pixel_scale_ratio = 1.2 (output pixels have larger spatial dimension than input).  Hacking the iscale at line 320 in  jwst/resample/resample.py to use 1 / np.sqrt(self.pscale_ratio) instead of 1.0 produces much better agreement for both source types.

Two, investigating some MOS calibration data taken in long slit mode, we found that the flux in a point source spectrum after resampling was significantly different when we modified the MSA metafile to use different slitlet sizes: wider slitlets produced less flux in the extracted spectrum.  A 5-shutter slitlet has significantly less flux than a 3-shutter slitlet, and a 30-shutter slitlet has a little less than the 5-shutter slitlet.  I suspect this is because the effective spatial pixel size on the sky is different in these different cases.

It looks like JP-3035 implemented some improvements for pixel area changes and flux conservation for imaging modes, but the changes are explicitly not implemented for spectral modes.  Can/should this work be extended to cover the spectral case?

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Test data for both cases is here, using calibration file jw01128019001_03102_00002_nrs1_rate.fits, reduced with pipeline v1.13.4, and modified MSA files:

██████████████████████████████████████████████████████████████

See the README in each subdirectory for the pipeline command to run.

 

[stscijgbot-jp]

Comment by David Law on JIRA:

Now that you've filed this ticket it's tickling my memory; we did run into a similar issue with NIRSpec IFU flux conservation before that I'd forgotten about.  This was a discussion via slack/email, and I think the only final ticket coming out of that was https://jira.stsci.edu/browse/JP-3163  As I recall, the issue was that all other modes assume that the *cal.fits files are in surface brightness units and process them accordingly, but since NIRSpec point source data was in flux units it didn't conserve flux when going to anything other than the default output spaxel scale.  That ticket fixed it for IFU mode, but not any other NIRSpec mode.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks David Law - I remembered that one, and initially thought this might be a surface brightness vs flux units issue too.  But testing with the same data propagated as a point source and as an extended source shows the same problem in both modes.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Double checking that the issue is the resampling and not the extraction aperture, I added a quick notebook to my test area to sum the flux over the whole slit in both the s2d and cal files.  The total flux shows the same effect in the s2d files: the total flux is different in the 4 different cases.  Total flux in the cal files matches very well.

Attaching a couple plots from the notebook to demonstrate.

[stscijgbot-jp]

Comment by David Law on JIRA:

I've been taking a look into this, and there seem to be a few issues going on.  I've tried running the example NIRSpec data and get similar results, and have also tried some LRS slit data as well.

First, a comment on how I believe the resample step is working; line 130 of resample_step.py calls resample_spec to figure out how to structure the WCS of the output grid, then line 131 of resample_step.py calls resample.py to actual do the drizzle resampling.  WCS restructuring is done with one set of routines for NIRSpec data, and another set of routines for non-NIRSpec data.

• For whatever reason the NIRSpec s2d resampled files have weird BUNIT keywords on the SCI array.  In Melanie's example data, the BUNIT keyword on the s2d files is DN/s, even though the BUNIT on the cal data was MJy.  When I reprocess the NIRSpec data I don't get a BUNIT keyword at all on these files.  LRS slit s2d files have the expected BUNIT of MJy/sr.

• It looks like the pixel_scale and pixel_scale_ratio keywords have no effect for LRS slit data, even though no output_wcs is provided at the user input.  Instead, initialization of the ResampleSpecData class sets pscale_ratio=1.0 and user inputs have no effect on this.

• If I force a different pscale_ratio for LRS slit data by changing the hardcoded value, I see weird behaviour similar to what Melanie saw for the NIRSpec data.  Somewhat surprisingly, it seems that the pscale_ratio keyword is scaling both the spatial and the spectral pixel size for the output grid for LRS data (though not for NIRSpec data).  It's not obvious to me that this is a desirable thing; changing the PSF sampling in the spatial domain is very different than changing the spectral sampling.  Hacking out this spectral resampling doesn't change the extracted 1d spectral mismatch though.

• Adopting Melanie Clarke 's line 320 revision of iscale = 1 / np.sqrt(self.pscale_ratio) fixes the s2d file for the NIRSpec point source data, but not for the MIRI LRS data.  To fix the MIRI LRS s2d file I need to use iscale = np.sqrt(self.pscale_ratio) instead.  This may be related to the MIRI data being in surface brightness units while the NIRSpec point source data is in flux units.  I'd be curious to know which works best for the NIRSpec extended source data.

• Even with a corrected MIRI LRS s2d file the extracted 1d spectrum is still wrong, which may be telling us something about problems with the LRS spectral extraction step as well.

Tagging Greg Sloan, Katherine Murray and Sarah Kendrew as well since LRS data is being brought into the test now.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

David Law - I think the weird BUNIT for the s2d file is related to this issue: JP-3088

[stscijgbot-jp]

Comment by David Law on JIRA:

Thanks Melanie Clarke ; I agree JP-3088 looks probable for that issue.

Regarding my final point above, the reason the LRS extracted 1d spectrum seems wrong in any resampled output seems to be because 1d extraction is always happening between pixels 27 and 34.  There's no accounting for dithers, and no accounting for whether any rescaling would mean that the trace was now at a different pixel number, which explains why I wasn't seeing a simple multiplicative factor difference.  Other tickets forthcoming on the LRS issues so that they don't clog up the generic issues at hand in this ticket.

[stscijgbot-jp]

Comment by Sarah Kendrew on JIRA:

I can elaborate on the extraction issue David Law but it's basically related to my email of 14 Feb on which you were cc'd.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

David Law - is this ticket in the plan for Build 11? It is high priority for NIRSpec.

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke I think it should be possible, though we probably need to run down some of the discrepancies above first.  Specifically, you found iscale = 1 / np.sqrt(self.pscale_ratio) for point source data above, but what scale do you need when working with extended source data?

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

That value was just a local hack, for testing - I'm not recommending implementing it like that! 

I think what this actually needs is to extend the handling that already exists for the imaging resampling, from JP-3035, to account for spatial pixel size changes regardless of whether they come from input arguments or from the WCS rectification.

[stscijgbot-jp]

Comment by David Law on JIRA:

Yes, though if you find a different hack necessary for pt vs extended source data that tells us something about what unit issues are involved in the problem.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

I can check for extended source data.  The relevant difference for point source and extended data should only be in the value of the area keywords, used to convert to Jy in the extraction.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Okay, I updated my test data in ██████████████████████████████████████████████████████████████ to include equivalent reductions of the same source, reduced with a 3-shutter slitlet, identified as extended instead of a point source, resampling with and without pixel_ratio=1.2.  I modified the original reductions as well to turn off any steps that behave differently for point and extended sources.

The upshot is that with the current resampling implementation (1.13.4) the point source and extended source extractions look identical, and flux is not conserved in an identical way.  These extractions look the same:

  3_shutter_output/jw01128019001_03102_00002_nrs1_x1d.fits

  3_shutter_output_extended/jw01128019001_03102_00002_nrs1_x1d.fits

and these extractions look the same:

  3_shutter_output_pixel_scale_ratio_1.2/jw01128019001_03102_00002_nrs1_x1d.fits

  3_shutter_output_extended_pixel_scale_ratio_1.2/jw01128019001_03102_00002_nrs1_x1d.fits 

Hacking the iscale to 1 / np.sqrt(self.pscale_ratio) works in both cases to bring the spectra with a non-standard pixel_scale_ratio in line with the default reduction.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

As David said, and as a slight understatement, I think there are a few different things going on here.

First, related to explicitly specifying a non-default pixel_scale_ratio or pixel_scale:

• There is currently no correction implemented for changes between input and output pixel size, when the user specifies a different scale. If the units are in flux density, this correction needs to be applied to the flux. If the units are in surface brightness, it needs to be applied to the pixel area keywords (PIXAR_SR, PIXAR_AR).
• The sense of the application for pixel_scale_ratio is opposite for NIRSpec and MIRI. For NIRSpec, the input pixel scale is divided by pixel_scale_ratio to get the output pixel scale.  For MIRI, it is multiplied.  I think the MIRI implementation needs to change to match the documented intent (pixel_scale_ratio is "Ratio of input to output pixel scale").
• When this ticket was originally filed, pixel_scale_ratio was ignored if specified for MIRI, but with recent changes to remove the drizpars file, it is now passed along. When specified with a value other than 1.0, the step sometimes crashes because some size assumptions no longer apply.  When it does not crash, the spatial sampling still doesn't look like it's correctly implemented.
• MIRI attempts to use the pixel_scale_ratio to change the wavelength array; NIRSpec does not. I think the MIRI implementation needs to change here: changing the pixel scale in the spectral dimension is a very different thing than changing it in the spatial one.
• Directly specifying pixel_scale is not implemented for either instrument. It should either be implemented or documented as unavailable.

Next, related to implicitly changing the pixel areas during resampling with default values:

• For MIRI LRS, looking at some sample data, the input and output pixel sky areas, according to the spatial WCS, are consistent, within a few thousandths of a percent.
• For the NIRSpec sample data from this ticket, they are not. For the 3-shutter case, the average area per pixel is about 10% different in the output than in the input.  For the 5-shutter case, it's about 7% different.  These implicit changes to pixel area are not accounted for in either flux scaling or area keywords.
• I think the difference between NIRSpec and MIRI is that NIRSpec output is linearly sampled in slit units rather than in arcsec along the slit. I think this implementation probably needs to change to handle all NIRSpec slits consistently and to match user expectations that by default, resampled pixels are about the same size as input pixels.

Finally, related to differences between the nominal and actual pixel areas:

• For imaging, the difference between nominal area (meta.photometry.pixelarea_steradians) and average sky area according to the WCS for each image is accounted for with a scaling factor for the flux (as of JP-3035). This accounts for both an offset between the nominal value and the average value in the image, as well as any differences between input images.
• For spectroscopy, these differences are not checked or accounted for. For my sample MIRI and NIRSpec data, it looks like the nominal value is off from the computed value by about a factor of 2 for MIRI and a factor of 4 for NIRSpec.
• I'm not sure if it's desirable to account for these offsets from nominal in spectroscopic resampling, at least for NIRSpec, since the data is calibrated to MJy/pixel prior to the photom step, and the nominal area value is used to convert the data to and from surface brightness in photom and extraction steps. I would like to know, however, if the nominal or computed values are more correct – either the reference files in photom or the computed WCS may need updating. For NIRSpec, JP-1341 may be relevant.

I will start work on addressing all of these things, starting by looking at how the spatial sampling is constructed for the output WCS for NIRSpec.  Please let me know if there are any objections or additional suggestions, especially if any of this sounds urgent and small enough to split out into a separate ticket.  I think this work is complex enough that it won't all be ready in time for Build 11.0.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

I have a draft PR started here:

#8596

This needs a lot more testing, especially for edge cases, but the fix as currently implemented fixes all my sample use cases from this ticket.  NIRSpec spectra extracted with all these input conditions:

• 3 shutters, point and extended
• 3 shutters with pixel_scale_ratio or pixel_scale specified
• 5 shutters

are now directly comparable to each other and to spectra extracted directly from the input cal file.  The issues with pixel_scale_ratio are also fixed, for both NIRSpec and MIRI.

I am currently using the nominal pixel area as set by photom, rather than the spatial area computed directly from the spectral WCS, on the assumption that's the more accurate value for flux calibration purposes.

Christian Hayes James Muzerolle - I think this change will impact NIRSpec flux calibration, since the last F-flats were derived from resampled spectra.  For my sample data, which was part of the calibration data set, the old resampled 3 shutter spectrum is higher than a direct extraction from a cal file by about 5%.  The new one is the same, to within ~.03%.

[stscijgbot-jp]

Comment by Christian Hayes on JIRA:

Thanks Melanie Clarke, I will take a look and can run some additional tests for the NIRSpec side.  

I can also pass on that this will likely require updated F-flats.  We may want to discuss if we need to coordinate delivering new reference files with merging these changes.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

I'm still working on testing for this issue and pushing minor fixes as I come across them, but I want to note a couple things.

First, thinking some more about the nominal pixel area, I think there is no way around assuming its value, rather than trying to compute it from the spatial WCS.  The spatial WCS does not encode the slit width, and that value does not seem to be recorded anywhere else in the data model either.  The solution I settled on is to just use the nominal area as specified, but modify it in the output for any non-default pixel scale ratio specified by the user. If the user specifies a direct pixel scale in arcsec instead, it is translated into a pixel scale ratio, using an approximate spatial sampling in the cross-dispersion direction, for consistent flux and area scaling everywhere in the code.

Second, the resample_spec step borrows its parameter spec from the resample step, but the resample step for imaging allows more customizations for the output WCS than are sensible for spectral resampling.  Currently, if crval, crpix, or rotation are specified for spectral resampling, they are ignored, but it is possible to specify output_shape.  If it is specified, the output array is simply truncated or expanded with padding.  I think the output shape handling is not doing anything useful for spectral data, since the only other customization that is usefully applied is the spatial scaling.

For simplicity, I'd like to remove output_shape from the resample_spec handling.  I think it also might make more sense to give the ResampleSpecStep class its own spec, instead of directly inheriting the one from ResampleStep, to make it clearer which parameters are actually used in spectral resampling. 

David Law Nadia Dencheva - do you have thoughts on that?  It would change the parameters listed for the resample_spec step, but none of the ones currently specified in the resampling parameter files would be impacted.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Since it's not necessary for the flux conservation work, I filed a separate ticket to discuss separating the resample_spec interface from resample: JP-3679.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Christian Hayes James Muzerolle - does NIRSpec want to hold off on merging this until new reference files can be delivered, or is it okay to merge to the dev pipeline first?  We've just started work on the new build (11.1), so there is plenty of time before these changes would go to production.

[stscijgbot-jp]

Comment by Christian Hayes on JIRA:

Melanie Clarke we had briefly discussed how to handle the need for new reference files and synchronizing with this merge.  We were leaning toward wanting this merged into dev first so we could use the dev pipeline to produce new reference files.  We may have to hold the new reference files from going to CRDS ops until the Build 11.1 candidate is released to avoid having Build 11 process data (or users reprocess data) with the new reference files.  

Let me confirm that we are okay with merging this in when it is ready and I can get back to you.  

[stscijgbot-jp]

Comment by Christian Hayes on JIRA:

Melanie Clarke I confirmed with the NIRSpec team that we are okay with merging this into the dev pipeline when it is ready and will work on the new reference files once that happens.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Great, thanks Chris!

[stscijgbot-jp]

Comment by Tyler Pauly on JIRA:

Fix included in #8596

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke I'm finally getting to test this a little for LRS and seeing some unexpected behavior.

I'm working with jw01536027001_03102 (two dithered LRS slit exposures), and by default the output s2d file measures 63x388.  If I use pixel_scale_ratio=3 it looks like it's oversampling the science area about as expected (roughly 185 pixels wide or so), but the total array width is 3x this again (561 pixels).  If I use pixel_scale_ratio=0.5 the science area looks like it's scaling to the right width, but the array is being truncated at half the expected size.  See lrs_resample_drl.png attached (cases in order left to right).  Note that Y axis is now not scaling, which is good and as intended with these changes, even though the stretch in the image makes this unclear.

Any thoughts?

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks for testing, David Law - that sounds like a bug in the MIRI array sizing.  I'll take a look.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks for reporting that!  The issue was exactly what it looked like: the pixel_scale_ratio was mistakenly applied twice when computing the output array extent.  I tested pixel_scale_ratio for MIRI with spec2, but the array size is computed differently for multiple inputs in spec3 for MIRI, and I neglected to test that case.  I have a fix in this PR, along with a new unit test for the case I overlooked:

#8727

[stscijgbot-jp]

Comment by David Law on JIRA:

Melanie Clarke Fix looks good to me, thanks!  Empirically, using absurd scalings (e.g., pixel_scale_ratio=10) also looks like it is behaving as expected in terms of keeping information in surface brightness units.

There will be an issue with LRS spectral extraction from s2d files with such resampled pixel scales, since the extraction relies on specific pixel positions.  However, that's a different issue and shouldn't hold up this work.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Thanks for testing, David. 

It's now noted in the documentation that extraction apertures need to change if pixel scales change, but in the long run, it might be worth considering if we should specify extraction apertures in arcsec instead of pixels.

[stscijgbot-jp]

Comment by David Law on JIRA:

Agreed- extractions in pixel units is fragile if applied to resampled data where the pixel size can change.

[stscijgbot-jp]

Comment by Melanie Clarke on JIRA:

Additional fix in #8727

[stscijgbot-jp]

Comment by David Law on JIRA:

Tested again after both PRs merged and still looks ok to me, closing.