Caltech-IPAC / Montage

Image Mosaics for Astronomers
Other
111 stars 48 forks source link

MontagePy Reproject not flux conserving #64

Closed mattcarv closed 1 year ago

mattcarv commented 1 year ago

Hi, hope you are doing well! 🙂

I am Matheus, an astronomy M.Sc. student at Western University (Canada) under the supervision of Professor Pauline Barmby. I'm sending this email to ask for help regarding the Montage reprojection code.

I have to reproject H-alpha maps from SITELLE (tangential projection) into the same as my Radio data (IRAM - Sanson-Flamsteed projection). Although it seems like the coordinates match properly, the flux gets out of hand, making the flux drop around 100 times. I was expecting some decrease of flux due to the reprojecting, considering that the small rotation would create NaN regions but not at this amount. I am attaching the images and the region to be used for reprojection.

This is the original image. At the top you can see the integrated flux. image

This is after the reprojection. image

This is the simple code I am using to reproject.

from MontagePy.main import mProject, mViewer

rtn = mProject('original.fits', 
               'projected5000_2.fits', 
               'region.hdr', expand=True)
print(rtn)

And this is the region I am trying to reproject into:

SIMPLE = T BITPIX = -64 NAXIS = 2 NAXIS1 = 67 NAXIS2 = 67 CTYPE1 = 'RA---SFL' CTYPE2 = 'DEC--SFL' CRVAL1 = 58.26047500001 CRVAL2 = 35.58937944444 CRPIX1 = 33.97131766181 CRPIX2 = 34.00082117819 CROTA2 = 0.000000 CDELT1 = -0.001472222283541 CDELT2 = 0.001472222283541 RADESYS = 'FK5' EQUINOX = 2000.0 END

If you have any tips, please let me know! Thank you very much!

Cheers, Matheus

JohnGood commented 1 year ago

Is there some reason you output uses much larger pixels than the input? It makes the two images rather hard to compare.

Montage is very strictly flux conserving and NaN values are handled correctly. I can't tell from what you sent: Is the original data total energy or flux density? Montage can handle either but you have to tell it. mProject also produces an "area" image to let you know how much of an output pixel is actually covered by the input. With your large pixels and NaNs that could be a significant effect.

Is there some specific reason you are using SFL projection? It can introduce significant spatial distortions and there are other projections that don't and are equal-area.

mattcarv commented 1 year ago

Hi, the image I am using as the reprojection has a (very) different resolution, so maybe that's why it binned like that. Below is the said image.

image

And I also attached its header, for reference. wcs_header.txt

And the original data (i.e. the integrated flux at the top) is the sum of all pixel values. I did not perform any sort of gaussian fitting to obtain the flux density as for my case this is not necessary. This H-alpha map is obtained from the LUCI code https://github.com/crhea93/LUCI and one of the outputs is a specific map for the H-alpha line, where each pixel value is already the flux density.

As for the SFL projection, unfortunately I had no choice in this case :(. My collaborator provided me the radio data from the IRAM 30m telescope so I don't know specifically why this projection is used.

Thanks! Matheus

JohnGood commented 1 year ago

The Montage command you used assumes the input is flux density and outputs flux density. So far so good. But to sum the pixels to get total energy you have to compensate for the pixel sizes, especially if the pixels vary in size due to projection. That shouldn't be a problem here since you are covering a small area and aren't too close to the pole (or the -180 / 180 wraparound longitude).

However if, as you say, there are multiple empty (NaN) input pixels, then there can be a lot of variability in the fractional pixel coverage in the output and to properly sum the values you have to convert them from flux density to total energy. For that you use the output "_area.fits" file, which contains the areas in steradians of the coverage of each pixel. Divide each flux density (.fits) pixel by the corresponding area pixel and sum the resultant "total energy" value.

Your "Integrated Flux" is in units of total energy, not flux density, so you have to have included pixel area in some form. Basically what I am saying is that you need to do the conversion on a pixel by pixel basis, taking into account partial pixel coverage.

mattcarv commented 1 year ago

I see, thanks for the input!

I appreciate the comments :)

Have a great day.