JP-3035: Fix computation of the photometric keywords PIXAR_* and scale resampled intensities according to actual pixel scale ratio.

[mcara]

Resolves JP-3035

Closes #7668
Closes #7390

This PR re-designed how resample computes photometric keywords PIXAR_* as well as corrects image data so that the flux of the resampled image (when image intensity is multiplied by PIXAR_SR) matches the flux of the input image. In my brief testing on NIRCAM images, flux error could have been up to 1.6% using previous computations and now it is down to about 0.005% (difference between the computed flux in the input and output images).

This required undoing some of the changes made in #5389.

I expect a large number of regression (and some unit) tests to fail due to this PR. I will will deal with failing unit tests in next commits.

This PR also allows some additional info to be passed with custom user-provided output_wcs as well it modified how output resampled image size is computed from such a WCS.

CC: @nden @hbushouse

Checklist for maintainers

• added entry in CHANGES.rst within the relevant release section
• updated or added relevant tests
• updated relevant documentation
• added relevant milestone
• added relevant label(s)
• ran regression tests, post a link to the Jenkins job below.
  How to run regression tests on a PR
• Make sure the JIRA ticket is resolved properly

Regression tests:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/931
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/934
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/944
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/945
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/953
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/974

[codecov]

Codecov Report

Attention: 26 lines in your changes are missing coverage. Please review.

Files	Coverage Δ
jwst/resample/gwcs_drizzle.py	80.00% <ø> (ø)
jwst/resample/resample_spec.py	78.46% <100.00%> (+0.09%)	⬆️
jwst/resample/resample_step.py	89.30% <88.88%> (-0.06%)	⬇️
jwst/model_blender/blendmeta.py	79.33% <33.33%> (-1.34%)	⬇️
jwst/resample/resample.py	90.00% <81.74%> (-6.09%)	⬇️

📢 Thoughts on this report? Let us know!.

[hbushouse]

I think this looks OK. I worry a bit, though, about rescaling the intensity of the output images based on a scale factor derived from one of the PIXAR_xx keywords in the input images. Since those keyword values are just some kind of average or nominal pixel area, and are only as accurate as the person who computed them and put them into a ref file chose to be, it seems like it could be risky.

[mcara]

I think these are some very very good questions!

1. Why doesn't the drizzling process take care of this?

   It actually does but not in a way that you expect. You probably expect to see resampled values to be 1/4 of input pixel values if output pixel scale is 1/2 of the input. That is how drizzlepac works. But it works like that because C code is supplied with the pixel scale ratio. But if you do not provide pixel scale ratio (i.e., leave it at the default value of 1 as JWST pipeline does after JP-1400: Add pixel_scale_ratio parameter to resample #5389) then drizzle still takes all of this into account because it will compute output pixel value as a weighted mean of a several input pixels with weights being proportional to pixel area of the intersection of input pixels with the output pixels (areas of overlaps). Since pixel values in JWST are in units of MJy/sr, pixel intensity should stay the same as that of the input image regardless of the pixel scale ratio. That is in theory... I will explain this in a little bit.

2. why the output image flux level has to be rescaled manually by the ratio of pixel areas.

   Normally, you would not need to do it manually if the desired output is in surface brightness units. However, for JWST there is a weird quantity meta.photometry.pixelarea_steradians that does not represent any real input pixel area (or maybe it does but only at one location in the input image which I could not figure out where - I asked this question on JIRA). But even if it were at a known, well defined location such as CRPIX, it would still be the pixel scale at one single location. When you work with surface brightness, drizzle will preserve pixel intensity (so, if all input pixels were 1 then all output pixels will be 1 if scale parameter to crizzle is set to 1). But ... ... measured flux will be different and it is different in the current implementation of all these scales and re-scalings. The reason for that is that WCS transformations will map input pixels to a different area in the output image that is not proportional to the ratio (meta.photometry.pixelarea_steradians / output_scale)**2 because meta.photometry.pixelarea_steradians in the best case is valid at one single location but the WCS will take into account distortions and pixel scales at every pixel location. The way this is currently done, total flux of the resampled images is 1.5-2% off from the total flux in the input image. IMO, this is not correct and this is a pretty large systematic error.

Simple Math

For simplicity, let's imagine that all input pixels have value of 1 or some other constant value (MJy/sr). I will use "prime" symbol to denote output (resampled) quantities. Based on my understanding of the documentation, total flux captured in the input image would be:

F = N * I * P

where N is the number of pixels in the input image and I is pixel intensity assumed to be constant across the image and also on the sky (thanks to flat-fielding) and P is the meta.photometry.pixelarea_steradians computed "somehow". Let's also assume that output image size is large enough to include all pixels from the input image (no flux is lost).

The flux in the resampled image would have a similar formula:

F' = N' * I * P'

(I is the same as in input because it is a surface brightness independent on pixel area).

Ideally flux should be preserved and so F'=F. Let's introduce average pixel scales <A> and <A'>. Then the relationship between N and N' is:

N' * <A'> = N * <A>

Therefore we can write:

F = N * I * P = (N * <A>) * I * (P / <A>) = (N' * <A'>) * I * (P / <A>)
  = N' * [I * (P / <A>)] * <A'> = N' * I' * <A'> = F' (if P' is defined as <A'>)

(just some tautology). Please observe how I (aptly I may say 😄 ) introduced the notation I' = I * (P / <A>). That is the pixel intensity that we must have in the output image in order for the flux to be conserved and also for the meta.photometry.pixelarea_steradians of the output image to actually be the average pixel scale of the output image (and since resampled image pixels have almost the same pixel scale due to lack of distortions - with exception of usually negligible projection effects - it will be the actual pixel area of every output pixel).

A more complicated math would not assume a constant intensity but would deal with the number of terms in summation (N * would become sum_{i}^{N}), individual pixel areas, etc., but I think the basic result would be the same.

Previously the code was computing something like this:

F' = N' * I * P * R

where R was the pixel ratio computed as the ratio of the output pixel scale (same for all pixels so <A'>) to the input pixel scale at CRPIX (let's call it A0) of the first of the input images. So:

F' = N' * I * P * (<A'> / A0) = (N * <A> / <A'>) * I * P * (<A'> / A0)
   = (N * I * P) * <A> / A0 = F * <A> / A0  <>  F

In addition, P' = P * R which was used for the "nominal pixel area" did not correspond to the actual pixel area of the resampled images. You can estimate the error in the flux obtained using the previous code as:

|<A> / A0 - 1| * 100%

Array A can be taken as Pixel Area Map of any instrument.

NOTE

The following links: https://jwst-pipeline.readthedocs.io/en/latest/jwst/photom/reference_files.html, https://jwst-pipeline.readthedocs.io/en/latest/jwst/photom/main.html (see "Pixel Area Data" section) all say something like this:

The step also populates the keywords PIXAR_SR and PIXAR_A2 in the science data product, which give the average pixel area in units of steradians and square arcseconds, respectively.

This is incorrect. No matter what method I used to compute average pixel area, it was alway quite different from PIXAR_SR. Documentation should be corrected.

[mcara]

Another benefit of the new code and scaling of individual input image intensities is that it allows to correctly resample input images that have different distortion models (i.e., if the distortion model has changed between visits, or photometric calibration results in different values in meta.photometry.pixelarea_steradians keyword in different input images). The old code was assuming that all input images had exactly the same distortion models, PIXAR_* values, etc.

[mcara]

Many of the regression tests are erroring with these errors:

ValueError: Unable to compute input pixel area from WCS of input image ...

These arise from _compute_image_pixel_area(wcs) returning None and it happens for spectral WCS. I will need to investigate and learn specifics of these instruments.

[mcara]

Restarted regression tests after this latest commit:
https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/945

[hbushouse]

Latest regtest run: https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/951/

Lots of unrelated failures. PR branch needs a rebase to pick up other changes from master.

[hbushouse]

@mcara What's the status of this? Are you done with your attempts to handle spectroscopic modes or still investigating?

[mcara]

Yes, I am done. Spectroscopic data have never had PIXAR_* modified, so I chose to not handle spectroscopic data, at least for now. The current algorithm will not work as it is tricky to derive pixel scale (or pixel area from one axis) from a spectroscopic WCS (and there are many types with different specifics).

[jdavies-st]

I know this is already implemented and working and seems correct. I just want to make a point here that in the _cal files that have been flat fielded, the variation in pixel area across the detector that affects the amount of flux a pixel receives has already been corrected to be the same.

Imagine that our QE of the detector is 100% for every pixel, and I take an image of a completely uniform white light source with exactly the same flux at all points. My detector will not register the same value at each pixel as some pixels have larger or smaller area than others due to distortions. So uniform QE and uniformly-illuminated astronomical scene still yields different per-pixel flux because of distortions that change the size of the pixels.

But that is only in the raw data. _cal files have been flat fielded, and with our assumption that QE is 1 at all pixels, our counts from a flat lamp will also vary with pixel area, even though the sensitivity of pixels is exactly 1 everywhere, in exactly the same way as above. So when I divide by the normalized flat, then I will divide out the differences in pixel area.

So once data has been flattened, any flux difference due to pixel area has been removed. And this is the input into the drizzle code for both HST and JWST.

Are we all on the same page?

It would be nice if the drizzle code did the correct thing when passed such flat-fielded data, i.e. that corrections such as this iscale business happens automatically within whatever drizzle API we are using, instead of being tucked into jwst.resample. drizzle itself should automatically know how to preserve surface brightness without it needing to be pre-computed here. It points to an oversight over in the drizzle package.

[hbushouse]

@jdavies-st I agree (as we have for years) that once images have been flat-fielded they have the correct surface brightness across all pixels, which I believe is the form that drizzle is expecting.

On a related note, when PR #8187 was merged recently, which shifted the photom step from retrieving pixel area keywords from the PHOTOM ref files to the AREA ref files, we noticed some differences in resampled results. Some of the AREA ref files had slightly different values for the pixel area keywords (PIXAR_SR, PIXAR_A2) than what had been in the PHOTOM ref files and this seems to have affected the pixel values in resampled images. I was surprised by this, because I assumed that those keyword values were not used at all in resampling (just information only). I assumed all computations from input to output WCS grids relied on the information contained in the gwcs object only. But apparently I was wrong. It appears that the pixel area keyword value is in fact being used somewhere within resampling. Can any drizzle aficionados explain to me where and why that's used?

[mcara]

@jdavies-st

I just want to make a point here that in the _cal files that have been flat fielded, the variation in pixel area across the detector that affects the amount of flux a pixel receives has already been corrected to be the same.
..............................
It would be nice if the drizzle code did the correct thing when passed such flat-fielded data, i.e. that corrections such as this iscale business happens automatically within whatever drizzle API we are using, instead of being tucked into jwst.resample. drizzle itself should automatically know how to preserve surface brightness without it needing to be pre-computed here. It points to an oversight over in the drizzle package.

drizzle actually works fine (i.e., does exactly what you expect it to do on flat-fielded data) modulo bugs, of course. This PR was not prompted by inability of drizzle to do the right thing with flat-fielded images but rather inconsistencies between where (position-wise) flat field is 1 and where the "nominal" pixel scale represented by PIXAR_SR (PIXAR_A2) were computed AND where, in the input image, when building the output WCS, the code is choosing reference pixel for computing the pixel scale of the output WCS which is neither the nominal pixel nor the pixel where flat field was 1... And all this needs to be made to work in such a way that when multiplying output image by PIXAR_SR (of the output image) one would get the correct flux. To do this, we need to compute a correction factor that gets multiplied to image data. Basically, we can just take resampled image and use numpy to multiply it by a constant factor. However, because drizzle already uses scale (or iscale) as a factor that it multiplies to image data (by default all data are multiplied to 1), I thought I could let drizzle do this and save an extra array multiplication simply for improved performance.

[mcara]

And this correction that I described in the previous message is specific to JWST pipeline: how it defines PIXAR_SR (at "nominal" pixel), how is flat-field normalization performed, and how it constructs the output WCS (which input pixel is chosen as reference for pixel scale computation). drizzle would not be able to guess what users do with their calibrations or how pixmap is constructed, etc.

[mcara]

@hbushouse

On a related note, when PR #8187 was merged recently, which shifted the photom step from retrieving pixel area keywords from the PHOTOM ref files to the AREA ref files, we noticed some differences in resampled results. Some of the AREA ref files had slightly different values for the pixel area keywords (PIXAR_SR, PIXAR_A2) than what had been in the PHOTOM ref files and this seems to have affected the pixel values in resampled images. I was surprised by this, because I assumed that those keyword values were not used at all in resampling (just information only). I assumed all computations from input to output WCS grids relied on the information contained in the gwcs object only. But apparently I was wrong. It appears that the pixel area keyword value is in fact being used somewhere within resampling. Can any drizzle aficionados explain to me where and why that's used?

It is not surprising at all that the result has changed. The whole purpose of this PR was to compute a scale factor (iscale) by which image data are multiplied. The value of this iscale factor depends on pixel area keywords, see:

https://github.com/mcara/jwst/blob/a6cac587e555865ea03f2b2a0cde4a91eeec24cf/jwst/resample/resample.py#L234

Basically,

iscale = np.sqrt(img.meta.photometry.pixelarea_steradians / input_pixel_area_computed_from_WCS)

(ratio of reported "nominal" area to the area computed from WCS).

[mcara]

Xref: #9246