Resample rescales image intensity on wrong assumptions

[stscijgbot-jp]

Issue JP-3907 was created on JIRA by Mihai Cara:

Please refer to #7894 for details and links to other issues that that PR was attempting to address.

In summary, that PR was trying to fix two issues: correct computation of PIXAR_* keywords for the resampled images and also minimize the difference between computed fluxes in input cal image(s) and output resampled image. Unfortunately, at the foundation of the math used to justify the use of iscale to rescale resampled image intensity to get a match between computed fluxes in CAL and I2D images - see #7894 (comment) - is the assumption that flux in the CAL image is obtained by multiplying the surface brightness by PIXAR_SR. I now believe this is wrong.

I do not know for sure why I thought (at that time) that this was intended usage for PIXAR_SR (it is possible it was a combination of factors: private communication, example notebooks that compared CAL fluxes with resampled fluxes, ...) In particular, I think this comment in the section "Calibrated image units and conversions" in https://jwst-docs.stsci.edu/jwst-near-infrared-camera/nircam-performance/nircam-absolute-flux-calibration-and-zeropoints is incorrect:

??JWST pipeline-calibrated NIRCam images are in surface brightness units, specifically MJy/sr (that is, MegaJanskys / steradian) {}per average pixel{}, as noted in the FITS "SCI" extension header keyword BUNIT in stage 2 pipeline products: "cal.fits" and "i2d.fits" images. It is appropriate to use average pixels since flat fielding corrects for variations in responsivity and areal coverage between pixels.??

??To measure photometry in MJy, convert the image to units of MJy/pixel by multiplying by the average area of a pixel in sr, a constant provided in the FITS header keyword PIXAR_SR. Note this value varies slightly (by a few percent) among the 8 short wavelength detectors and between the 2 long wavelength detectors. They correspond approximately to the nominal pixel sizes: 0.031" and 0.063", for the short and long wavelength channels, respectively.??

If one were to follow this recommendation, one would get a reasonably good agreement between aperture photometry in CAL and I2D images only at locations where CAL pixels have area close to PIXAR_SR and only for relatively small apertures. In most other locations, the discrepancy between CAL and resampled aperture photometry would be larger (probably up to 1.5% for NIRCAM) due to distortions in the CAL images.

With regard to the above cited paragraphs: "JWST pipeline-calibrated NIRCam images are in surface brightness units, specifically MJy/sr (that is, MegaJanskys / steradian) per average pixel" directly contradicts the notion of flat-fielding mentioned in the next sentence because flat-fielding specifically accounts for differences in the image intensity due to variable area on sky corresponding to different pixels.

In order to obtain correct photometry in flat-fielded CAL images it is necessary to multiply those images by pixel area maps (PAM) in order to convert to flux. The same is true for resampled images although, due to instrument distortions being removed (but not WCS projection effects!), one can assume that PAM is almost constant, especially for mosaics that are not very large.

If PAM is used in CAL images to convert surface brightness to flux instead of PIXAR_SR and those fluxes are compared to fluxes computed in resampled images, then there is no need to use iscale (that is, we can set iscale=1). Not only that, but this approach insures a very good agreement between CAL and I2D fluxes at any location in the image.

I believe the part related to iscale of #7894 must be undone and documentation should be updated to provide correct advice on how to convert surf. brightness to flux.

I also want to add that data models have .area arrays that seem they might be used as PAM. Unfortunately, in my experiments for a NIRCAM image, (area * pixelarea_steradians) != true PAM, probably because pixelarea_steradians is not defined where the area array was normalized to 1.