reproject_exact() is not conserving flux

[acikota]

Dear all,

I have an issue with the reproject_exact() function, which doesn't properly conserve the flux.

To demonstrate the issue, I prepared an example which you can download from my Dropbox: https://www.dropbox.com/s/x1a2upsxft5iru0/example.zip?dl=1

I have a highly sampled image (input.fits), which contains 4 pixels of values 10, 20, 30 and 40. When I reproject the image onto a template (templates/template.fits), the sum of the reprojected pixels is smaller by a factor of ~700.

I thought that the issue is produced by the difference in the area between the input image and the reprojected output image, i.e. that the flux is conserved per area, however, the pixel area ratio between template and input image is only ~226. Furthermore, there is also an inconsistency in the pixel values, e.g. the pixels with the values 40 and 30 are equally bright in the reprojected image.

I would very appreciate your help or advice! Thank You!

Best regards, Aleksandar

[astrofrog]

@acikota - thanks for the bug report! I spent a little while investigating this, and it looks like this is a precision issue due to the very small pixels. I did some simple tests and it looks like things start to not work properly below a few times 0.01" in pixel sizes:

For the record this plot can be made with:

# Check how good the reprojection is as a function of resolution

import numpy as np
from astropy.wcs import WCS
from astropy.io import fits
from reproject import reproject_exact

resolutions = np.logspace(-7, 0, 1000)
ratios = []

for res in resolutions:

    wcs1 = WCS()
    wcs1.wcs.ctype = 'RA---TAN', 'DEC--TAN'
    wcs1.wcs.crval = 13, 20
    wcs1.wcs.crpix = 10., 10.
    wcs1.wcs.cdelt = res, res

    wcs2 = WCS()
    wcs2.wcs.ctype = 'RA---TAN', 'DEC--TAN'
    wcs2.wcs.crval = 13, 20
    wcs2.wcs.crpix = 3, 3
    wcs2.wcs.cdelt = 3 * res, 3 * res

    array = np.zeros((19, 19))
    array[9, 9] = 1

    result, _ = reproject_exact((array, wcs1), wcs2, shape_out=(5, 5))

    ratios.append(9 * np.nansum(result) / np.nansum(array))

import matplotlib.pyplot as plt

plt.semilogx(resolutions * 3600, ratios)
plt.xlabel('resolution (arcsec)')
plt.ylabel('actual / expected')
plt.ylim(0, 5)
plt.savefig('ratios.png')

Of course we should try and fix this and determine where precision is being lost, but in the mean time I think I should definitely include a warning in the documentation as well as when running reproject_exact if the resolution of the target or final image is too small, so I'll do that in a PR now.

[acikota]

@astrofrog Thank you a lot for investigating the bug!

What would you suggest I should do in order to reproject my image onto larger pixels? Do you think modifying the TOLERANCE in overlapArea.c would be useful? (I am trying to find out how to re-compile the module after modifying the overlapArea.c in the site-packages/reproject/spherical_intersect/)

Thanks!

[astrofrog]

@acikota - I think someone needs to investigate what could be going wrong with the precision in the module, as I'm not sure if changing the tolerance is going to be sufficient. If you modify overlapArea.c, the easiest way to recompile is to run:

pip install -e .

at the top level of the repository.

[astrofrog]

I will try and find time this month to investigate the precision issues, but let me know if you manage to get things working in the mean time!

[mtakahiro]

@astrofrog I would like to second acikota's first point --- the functionality seems to return an image that conserves surface flux density (e.g., count/pixel) but not total flux per area in the sky. So if users change a pixel scale by a factor of x, then the new array needs to be divided by x^2.

What I was doing was opposite (low-res image reprojected to high-res), and conversion scale x is rather large (x=7).

[thomaswilliamsastro]

@astrofrog did this ever get fixed? For JWST images, the NIRCam short pixel scale is 0.03", exactly the regime where this might be problematic

[astrofrog]

@thomaswilliamsastro - it looks like this is still an issue. I haven't had a chance to investigate but I would welcome any help!

[astrofrog]

I have been assuming in the past that the issue is in the C overlap code but it could very well be in the Python code setting up the pixel vertices - I will look into this more.

[astrofrog]

This does appear to be an issue in the C overlap code which comes from Montage (and Montage shows the same precision issues).