PSF Subtraction Step Function

[ell-bogat]

Describe your changes

• New step function to perform PSF subtraction via pyKLIP
• New subclass of pyKLIP_data in data.py
• Tests for do_psf_subtraction() and PyKLIP_data in tests/test_psfsub.py

Results using OS11 simulations:

Testing using gaussian stellar & planetary PSFs:
Note: I'm only comparing to the "analytical" result outside of the IWA, as there is increased deviation as you approach the mask center. Also, there is offset between the corgidrp and analytical result equal to the median of the corgidrp output, so if you subtract the median (on the order of -1e-4) then they match much more precisely (order 1e-6 or less).

Output header keywords:

• We keep by default all the header keywords from the first science frame in the dataset
• We then add all the output pyklip header keywords to the ext_hdr.

Still to do:

• Handle DQ arrays properly (may need to update what value we assign to the nans pixels after pyklip runs)
• Double check that roll angle direction convention is correct!

Future PRs:

• Calibrate for algorithm throughput in next PR (use OS11 mock calibration for now, eventually use a calibration input file that contains off-axis psfs)
  • Account for core throughput transmission for point sources near the detection limit
• Incorporate WCS header when ready
• How to track mask position after PSF subtraction? If star-mask alignment is possibly different in each input image?
• Handle error arrays properly

Type of change

• New feature (non-breaking change which adds functionality)

Reference any relevant issues (don't forget the #)

Issue #8 : Produce Calibrated PSF Subtracted Images with pyKLIP

Checklist before requesting a review

• I have linted my code
• I have verified that all unit tests pass in a clean environment and added new unit tests, as needed
• I have verified that all docstrings are properly formatted and added new documentation, as needed

[ell-bogat]

I'm trying to confirm the convention for which direction is a positive roll angle. I was assuming the names of the OS11 files referred to the roll of the telescope, but when I do a straight subtraction of the negative direction roll from the positive direction roll, I get this:

If the positive roll direction is from North towards East (left), then the planet should appear to move in the opposite direction, aka clockwise, right? So is it possible the filenames are really referring to the planet position angle, or possibly mislabeled? Or is the convention different? If I don't swap the sign of the rolls when I input them into pyKLIP, I get a negative planet, and I want to make sure the convention is correct.

[semaphoreP]

Hi Ell, I frequently have a negative sign uncertainty on the roll angles. I think it's ok if, for OS11, we need to apply a negative sign to get the roll angles right. Hopefully we can get this right after we figure out the ROLL keywords in real CGI headers, but for now, I think having to apply a negative sign for OS11 is fine.

[ell-bogat]

Thanks @semaphoreP !
The main issue I'm running into now is in the WCS info. My current understanding is there are two ways to store the coordinate transformation, either using a CD-matrix or a PC-matrix + a CDELT matrix (still figuring out exactly what those are). And the fits headers default to storing the PC version, while pyKLIP expects the CD version. For some reason it was working anyway before I updated to the same software versions that the github testing environment uses, but now I'm getting a "no cd present" angry message from pyKLIP causing the tests to fail. Also if I just overwrite the WCS header to None so that pyKLIP ignores it, we don't get the correct psf subtraction result anymore. Is there a simple way to convert from PC to CD? Do we need to tackle generating the correct WCS header info first?

[semaphoreP]

Can you clarify the statement "fits headers default to storing the PC version"? How are the fits headers being created? Pyklip can internally make wcs headers from roll angles and centers using the wcsgen module it has, in case that's useful.

…

________________________________ From: Ell Bogat ***@***.***> Sent: Friday, January 31, 2025 10:48 AM To: roman-corgi/corgidrp ***@***.***> Cc: Jason Wang ***@***.***>; Mention ***@***.***> Subject: Re: [roman-corgi/corgidrp] PSF Subtraction Step Function (PR #280) Thanks @semaphoreP<https://urldefense.com/v3/__https://github.com/semaphoreP__;!!Dq0X2DkFhyF93HkjWTBQKhk!Wgj0KvmVAgbRCM4nmi31-Kct5FbAprhkAjlqpDX73mSSU8iecnPZKn-m3n69R8GFVbG9RgB7gc9wiGoJ66KfYHzRN6_2stgU$> ! The main issue I'm running into now is in the WCS info. My current understanding is there are two ways to store the coordinate transformation, either using a CD-matrix or a PC-matrix + a CDELT matrix (still figuring out exactly what those are). And the fits headers default to storing the PC version, while pyKLIP expects the CD version. For some reason it was working anyway before I updated to the same software versions that the github testing environment uses, but now I'm getting a "no cd present" angry message from pyKLIP causing the tests to fail. Also if I just overwrite the WCS header to None so that pyKLIP ignores it, we don't get the correct psf subtraction result anymore. Is there a simple way to convert from PC to CD? Do we need to tackle generating the correct WCS header info first? — Reply to this email directly, view it on GitHub<https://urldefense.com/v3/__https://github.com/roman-corgi/corgidrp/pull/280*issuecomment-2627791200__;Iw!!Dq0X2DkFhyF93HkjWTBQKhk!Wgj0KvmVAgbRCM4nmi31-Kct5FbAprhkAjlqpDX73mSSU8iecnPZKn-m3n69R8GFVbG9RgB7gc9wiGoJ66KfYHzRN1yDgTNH$>, or unsubscribe<https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AAREI3DKAIMMVUIQPGD5H4D2NOSOTAVCNFSM6AAAAABUV6NUBCVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMMRXG44TCMRQGA__;!!Dq0X2DkFhyF93HkjWTBQKhk!Wgj0KvmVAgbRCM4nmi31-Kct5FbAprhkAjlqpDX73mSSU8iecnPZKn-m3n69R8GFVbG9RgB7gc9wiGoJ66KfYHzRNx9Yw4EB$>. You are receiving this because you were mentioned.Message ID: ***@***.***>

[ell-bogat]

Basically any time you use the method to convert a WCS object to a fits header, it is saved using the PC formalism.
e.g. even if I generate an object using the pyKLIP generate_wcs function which uses the CD matrix:

center = (30,30)
roll = 13.
wcs_generated = generate_wcs(roll,center) # WCS object
display(wcs_generated)

> WCS Keywords
> 
> Number of WCS axes: 2
> CTYPE : 'RA---TAN' 'DEC--TAN' 
> CRVAL : 0.0 0.0 
> CRPIX : 30.0 30.0 
> CD1_1 CD1_2  : -2.7065835132923e-06 6.24864039844069e-07 
> CD2_1 CD2_2  : 6.24864039844069e-07 2.70658351329232e-06 
> NAXIS : 0  0

If I then call the wcs.to_header() function to store it in the extension header it produces this:

wcs_header = wcs_generated.to_header()
display(wcs_header)

> WCSAXES =                    2 [/](https://file+.vscode-resource.vscode-cdn.net/) Number of coordinate axes                      
> CRPIX1  =                 30.0 [/](https://file+.vscode-resource.vscode-cdn.net/) Pixel coordinate of reference point            
> CRPIX2  =                 30.0 [/](https://file+.vscode-resource.vscode-cdn.net/) Pixel coordinate of reference point            
> PC1_1   = -2.7065835132923E-06 [/](https://file+.vscode-resource.vscode-cdn.net/) Coordinate transformation matrix element       
> PC1_2   =  6.2486403984407E-07 [/](https://file+.vscode-resource.vscode-cdn.net/) Coordinate transformation matrix element       
> PC2_1   =  6.2486403984407E-07 [/](https://file+.vscode-resource.vscode-cdn.net/) Coordinate transformation matrix element       
> PC2_2   =  2.7065835132923E-06 [/](https://file+.vscode-resource.vscode-cdn.net/) Coordinate transformation matrix element       
> CDELT1  =                  1.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Coordinate increment at reference point  
> CDELT2  =                  1.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Coordinate increment at reference point  
> CUNIT1  = 'deg'                [/](https://file+.vscode-resource.vscode-cdn.net/) Units of coordinate increment and value        
> CUNIT2  = 'deg'                [/](https://file+.vscode-resource.vscode-cdn.net/) Units of coordinate increment and value        
> CTYPE1  = 'RA---TAN'           [/](https://file+.vscode-resource.vscode-cdn.net/) Right ascension, gnomonic projection           
> CTYPE2  = 'DEC--TAN'           [/](https://file+.vscode-resource.vscode-cdn.net/) Declination, gnomonic projection               
> CRVAL1  =                  0.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Coordinate value at reference point      
> CRVAL2  =                  0.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Coordinate value at reference point      
> LONPOLE =                180.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Native longitude of celestial pole       
> LATPOLE =                  0.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [deg] Native latitude of celestial pole        
> MJDREF  =                  0.0 [/](https://file+.vscode-resource.vscode-cdn.net/) [d] MJD of fiducial time                       
> RADESYS = 'ICRS'               [/](https://file+.vscode-resource.vscode-cdn.net/) Equatorial coordinate system

[maxwellmb]

Hi all,
This is looking overall nice to me. Please see the comments attached here and let me know if anything is unclear or if you want to discuss any of them further. Always happy to iterate.
Cheers,
Max

[ell-bogat]

Ready for re-review! Right now here are the tasks that I think should be follow-up issues:

• Calibrate for algorithm throughput in next PR (use OS11 mock calibration for now, eventually use a calibration input file that contains off-axis psfs)
  • Account for core throughput transmission for point sources near the detection limit
• Incorporate WCS header when ready
• Track mask position in each frame after PSF subtraction?
• Handle error arrays properly