Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unable to get WCS information from fits files #46

Open
lsmercadante opened this issue Feb 2, 2024 · 6 comments
Open

Unable to get WCS information from fits files #46

lsmercadante opened this issue Feb 2, 2024 · 6 comments

Comments

@lsmercadante
Copy link

I have downloaded the difference images and read in the fits files to a jupyter notebook using astropy. The problem I am having is I cannot convert the pixels to WCS coordinates. Whenever I attempt this I get an error. Here is how the image is read in:

image

And here is the error I encounter:

image image

I'm not sure how to fix this problem.

@genghisken
Copy link
Collaborator

Hi Lilah

Thanks for raising this. It helps us keep track of any solutions and time spent on them.

As you may have seen there is quite strong feeling about astropy standards compliance! I'm not a WCS expert, but one solution I've found that kind-of works is below. I'm not sure what effect deleting the corner cards have, so I'd need to double check with local experts. But try that and see what happens. The solution will still whine about SIP, but should only produce a warning, not an error this time.

from astropy.io import fits
from astropy.wcs import WCS
h = fits.open('01a57801o0580o_diff.fits')
header = h[0].header

w = WCS(header)
# whines about SIP but then will also barf out with an error.

# Delete the corner info

del header['CNPIX1']
del header['CNPIX2']

# Now try again - still whines about SIP but should't crash this time.

w = WCS(header)

@lsmercadante
Copy link
Author

This seems to have solved the problem, thank you! As a follow up question, I noticed that when the pixel position changes from say (200,200) to (199,200), instead of just the RA value changing, both the RA and DEC values change.

This seems unusual to me and I was wondering if there is a reason for this?

@genghisken
Copy link
Collaborator

Nothing unusual I don't think - let us know if you see something particularly unusual. Is the dec change particularly large? Although it's an equatorial mount there will be slight variations in the orientation of the camera, plus the source image is a very large field of view (29 sq degrees), so astrometric corrections will have been applied to account for distortion, etc.

@djbrout
Copy link

djbrout commented Feb 8, 2024

Excellent. Thank you Ken!

@rsiverd
Copy link

rsiverd commented Mar 26, 2024

In addition to using a few keywords in ways we shouldn't, the reduced ATLAS images actually contain two different WCS solutions. A first pass is performed using astrometry.net that installs SIP terms. A second solve adjusts the solution and saves its result using the PV distortion convention (TPV CTYPE). Unless you modify the headers, you will be getting the TPV solution.

I use the following code to clean things up on-the-fly and take the SIP solution:

import astropy.io.fits as pf
import astropy.wcs as awcs

## Load image header:
ihdr = pf.getheader(context.input_file, extname='image')

## Modify header for compatibility:
whdr = ihdr.copy(strip=True)                            # make a working copy
toss = ['RP_SCHIN', 'CDANAM', 'CDSKEW', 'RADECSYS']     # these prompt warnings
toss += ['CNPIX1', 'CNPIX2']
toss += [x for x in whdr if x.startswith('PV')]         # these are PV terms
for kk in toss:
    whdr.remove(kk, ignore_missing=True, remove_all=True)

## Update CTYPEx to prompt SIP interpretation:
whdr['CTYPE1'] = 'RA---TAN-SIP'                         # overwrite RA---TPV
whdr['CTYPE2'] = 'DEC--TAN-SIP'                         # overwrite DEC--TPV

## This should work correctly and without tons of warnings:
twcs = awcs.WCS(whdr)
```

@genghisken
Copy link
Collaborator

Rob's solution above is the definitive one. Please use it, rather than my hack!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

4 participants