PSF looks wild

[laldoroty]

After updating the PSF retrieval and rotation to take the WCS into account, the coadded PSF looks wild.

Originally posted by @laldoroty in https://github.com/Roman-Supernova-PIT/diff-img/issues/13#issuecomment-2237682761

[laldoroty]

Copying and pasting @wmwv's comment from the other issue:

Does Swarp know how to coadd PSFs? I'm sensing the answer is no, it doesn't know anything special about them so instead you have to feed it a thing that is an image exactly aligned. You don't really want a full WCS, but rather just the rotation matrix.

When I look at the WCS for the PSF files, I see things like:

GS_WCS  = 'AffineTransform'    / GalSim WCS name                                
CTYPE1  = 'LINEAR  '           / name of the world coordinate axis              
CTYPE2  = 'LINEAR  '           / name of the world coordinate axis              
CRVAL1  =  -347.00793890391725 / world coordinate at reference pixel = u0       
CRVAL2  =    243.3111579925501 / world coordinate at reference pixel = v0       
CRPIX1  =     3167.59201485829 / image coordinate of reference pixel = x0       
CRPIX2  =   2247.2583100556526 / image coordinate of reference pixel = y0       
CD1_1   = -0.10799018359621737 / CD1_1 = dudx                                   
CD1_2   = -0.00219783174873922 / CD1_2 = dudy                                   
CD2_1   = 0.001247801388696990 / CD2_1 = dvdx                                   
CD2_2   =  0.10651140156275885 / CD2_2 = dvdy

and

GS_WCS  = 'AffineTransform'    / GalSim WCS name                                
CTYPE1  = 'LINEAR  '           / name of the world coordinate axis              
CTYPE2  = 'LINEAR  '           / name of the world coordinate axis              
CRVAL1  =  -189.75909248410068 / world coordinate at reference pixel = u0       
CRVAL2  =   -275.8411189309953 / world coordinate at reference pixel = v0       
CRPIX1  =   3100.0282045535396 / image coordinate of reference pixel = x0       
CRPIX2  =    240.1685487576857 / image coordinate of reference pixel = y0       
CD1_1   = -0.05421338482127558 / CD1_1 = dudx                                   
CD1_2   = -0.09033685129072228 / CD1_2 = dudy                                   
CD2_1   = -0.09335584816625218 / CD2_1 = dvdx                                   
CD2_2   =  0.05647968276683095 / CD2_2 = dvdy

But I think you want CRPIX1, CRPIX2, CRVAL1, CRVAL2 to have the same values for each PSF image, such that CRPIX1, CRPIX2 point to the middle of the array of the PSF. And then the CD matrix contains the rotation about this point.

[laldoroty]

Okay, this feels illegal, but I'm going to try manually setting all CRPIXs to NAXIS/2, and all CRVAL to 0.

[wmwv]

Yes, something like that is correct. Figure out a small test to make sure you get the pixel, half-pixel alignment correct for both even and odd. E.g., for even NAXIS=40 you want in between the middle pixels, and for NAXIS=40 you want the middle pixel. But check if you're 0-, 0.5-, or 1-indexed when reading the FITS file (I don't remember off the top of my head).

[wmwv]

"feels illegal" Just to explicitly reassure you. Yes, the issue here is that you want to make a relative coordinate expression but are forced to do so in an absolute coordinate system context. So the "feels illegal" part is exactly making the absolute coordinates agree, which is conceptual making it a relative coordinate system between the PSF stamps.

[laldoroty]

Ok, updated code and tested get_imsim_psf() with PSF size = 31 and 30. Then, applied rotate_psf() to each. Results (pre-rotate, after rotate):

So, looks like rotate_psf() is working correctly but get_imsim_psf() only works for NAXIS%2 != 0.

As for indexing, astropy says "It should be noted that astropy uses zero-based indexing when referring to HDUs and header cards, though the FITS standard (which was designed with Fortran in mind) uses one-based indexing." What I'm not sure about is where the location of the integer is, if that makes sense--i.e., is it "1" before the pixel, at the center of the pixel, or after the pixel?

[laldoroty]

Here's what our SWarp coadded PSF looks like with stamp size = 201:

Probably should have done it with a smaller stamp size so it would be more obvious if there was some weird offset. I'll try 30 next.

[laldoroty]

Here's what we got for 30x30:

And 31x31:

[laldoroty]

They do look slightly different around the edges, so I definitely need to write a test, but I think i need to understand the problem better first. I am going to go ahead and let it run all the way through with 31x31 and just see what happens.

[laldoroty]

hmmm yeah difference image wild. Crashes on that decorrelation kernel-generating step with "Exception: The kernel can't be normalized, because its sum is close to zero. The sum of the given kernel is < 0.01. For a zero-sum kernel, set normalize_kernel=False or pass a custom normalization function to normalize_kernel."

[laldoroty]

The reference and science images do appear to be aligned to each other:

[laldoroty]

However... PSF and image axes are flipped.

The original science PSF has a different orientation:

...and when that PSF goes through rotate_psf(), it loses the N, E axes entirely:

[laldoroty]

I think that if I'm using SWarp to rotate and coadd, I should be using SWarp to rotate the science PSF as well, and not the "custom" rotate_psf() function. I will rewrite this function to do that instead.

[laldoroty]

Underneath rotate_psf() is sfft.utils.ImageZoomRotate.IZR(). This uses SWarp to rotate. I will leave it, but the issue is that IZR only returns the data array, not the data array AND the header. Added the header back in (in sfft) to see if it fixes the dropped yellow axis in the above comment. It does, but the rotated and un-rotated PSFs have different orientations again. Un-rotated: Rotated:

The coadded and original get_imsim_psf() images also show "LINEAR-LINEAR" for WCS instead of FK5 (everything else, including the rotated get_imsim_psf(), has FK5).

Still getting the kernel can't be normalized error, with the same bad-looking difference image.

I also changed CRVALs to SN RA, Dec instead of 0., 0.

[laldoroty]

The coadded and original get_imsim_psf() images also show "LINEAR-LINEAR" for WCS instead of FK5 (everything else, including the rotated get_imsim_psf(), has FK5).

I think this is because of galsim.affine. https://galsim-developers.github.io/GalSim/_build/html/wcs.html

This preserves the transformation near the location of image_pos, but it is linear, so the transformed values may not agree as you get farther from the given point.

[laldoroty]

ah-ha... opened everything as a mosaic in DS9. it still thinks the PSFs belong at 0, 0 in the corner of the images.

[laldoroty]

Whenever I try to open any of the PSFs as a mosaic segment on top of an image, DS9 throws the error "!! AST: Error in routine wcsSkyFrame at line 37 in file frame/wcsast.C. ! astSetAttrib(Frame): Invalid System description "FK5".".

[laldoroty]

Okay, this runs when I take the PSF coadd out and just use a single PSF as the template PSF (but leaving the template coadd in). Difference image still has some stuff in it, but at least it actually generates the decorrelation kernel.

[laldoroty]

I think I should rewrite to use Masao's suggested approach, which is: do single-epoch DIA on each image N times, where N is the number of template images. Then, coadd the difference images at the end. So, I would not coadd anything until the very end. It is more computationally expensive, but it builds on code that already works, and I think I'm wasting a lot of time trying to hunt this down. I could also simultaneously rewrite so that I rotate the reference to match the science instead of the other way around (like I'm currently doing).

To make the science image where I get the ZPT (i.e., the science image that's un-subtracted but has been cross-convolved with the reference and the decorrelation kernel applied), I could just make N of these, do photometry on the stars in all N images, and take the median exactly as I was doing before. This also may allow for SFFT on smaller stamp sizes (1K)--the issue with making the stamps smaller was there weren't enough stars in the stamps to get a ZPT. This would multiply the number of stars by N.

Thoughts? @wmwv

[dscolnic]

Im not sure that this math works for getting the zpt. I think you can do this path, but there's a different math

[wmwv]

@laldoroty If you're going to go done the N-combinations path, I might offer:

"The NN2 Flux Difference Method for Constructing Variable Object Light Curves" https://ui.adsabs.harvard.edu/abs/2005AJ....130.2272B/abstract Barris, Tonry, Novicki, Wood-Vasey """ We present a new method for optimally extracting point-source time variability information from a series of images. Differential photometry is generally best accomplished by subtracting two images separated in time, since this removes all constant objects in the field. By removing background sources such as the host galaxies of supernovae, such subtractions make possible the measurement of the proper flux of point-source objects superposed on extended sources. In traditional difference photometry, a single image is designated as the ``template'' image and is subtracted from all other observations. This procedure does not take all the available information into account and for suboptimal template images may produce poor results. Given N total observations of an object, we show how to obtain an estimate of the vector of fluxes from the individual images using the antisymmetric matrix of flux differences formed from the N(N-1)/2 distinct possible subtractions and provide a prescription for estimating the associated uncertainties. We then demonstrate how this method improves results over the standard procedure of designating one image as a template and differencing against only that image. """

I think going back to single-image subtractions make sense. Coadding the difference images doesn't seem necessary. I would photometry and then combine in catalog space.

[wmwv]

I agree that the coaddition and construction of templates is its own project.

[dscolnic]

Yeah, I think can do that as @wmwv said - the co-addition in catalog space. So with three templates and one science image, you can three zeropointed fluxes. Scale the fluxes once more to a common zpt, average them. Im not sure though what the proper uncertainty from this process is? Cause you cant do scatter/sqrt(images) cause there is separate noise coming from template and science image. Thoughts?

[laldoroty]

Okay, coaddition in catalog space makes more sense to me.

The uncertainty is in eqns 9-13 in the paper!

[laldoroty]

Reorganized the code structure plan according to this: https://github.com/Roman-Supernova-PIT/diff-img/discussions/14

I am going to open a new issue for this.

[laldoroty]

Ok to close this issue since the problem is no longer relevant?