I’ve tried to programatically set up WCS in python using data aquired from astrometry (i know i have ready wcs files there, but i want to understand this topic fundamentally). Unfortunetly my best try has 29" error in determination of Neptunes coordinates but when using wcs from astrometry i have only 2.12" error. I think i may got something wrong with rotation becouse there are multpile definitions and i don’t know which one to use (i’m still an astronomy newbie). Below i show screenshots of astrometry measurments and my python code. Funny thing is when i input 55 deg into my rotation matrix the grid lines are off by many degrees. Can someone please point my errors? I want to be able to create good WCS objects in the future:)
Astrometry calibrated image ^
My code to create WCS object:
import numpy as np
from astropy import wcs
from astropy.io import fits
# Otwórz plik FITS
hdulist = fits.open('_V_30.00s _10.10_2023-11-04_0000.fits')
# Utwórz nowy obiekt WCS
w = wcs.WCS(naxis=2)
rotation_angle = 145.0 # rotation_angle w stopniach
rotation_angle_rad = np.radians(rotation_angle) # Convert to radians
# Oblicz środek obrazu
image_width = hdulist[0].header['NAXIS1']
image_height = hdulist[0].header['NAXIS2']
image_center = [image_width / 2, image_height / 2]
# Ustawienia obiektu WCS
w.wcs.crpix = image_center
w.wcs.cdelt = np.array([-0.0000933333, -0.0000933333]) # pixel scale 0.336"/px
# Set the pc matrix for rotation
w.wcs.pc = np.array([[np.cos(rotation_angle_rad), -np.sin(rotation_angle_rad)],
[np.sin(rotation_angle_rad), np.cos(rotation_angle_rad)]])
w.wcs.crval = [355.776, -3.193] # RA, DEC środka obrazu
w.wcs.ctype = ["RA---AIR", "DEC--AIR"]
w.wcs.set_pv([(2, 1, 45.0)])
# Przekształcenie obiektu WCS na nagłówek FITS
header = w.to_header()
# Zaktualizuj nagłówek głównego HDU
hdulist[0].header.update(header)
hdulist.writeto('Neptune_with_coordinates3.fits', overwrite=True)
First off, my personal preference: rather than setting the attributes after the WCS object is instantiated, I almost always instead first create a dictionary of header keyword/values, and then instantiate the WCS object from that header. The attribute setters are convenient for parsing computed output (e.g., when the PCij matrix is the output of a matrix expression), but in your script, you are essentially setting each individual matrix element separately, so I think a dictionary header would be clearer and facilitate comparisons.
Now, it’s difficult to be certain about the errors because you haven’t provided much information/output to go on beyond the Python code. But, here are my guesses:
A definite error: Your values for CRPIXi are almost certainly off by half a pixel each (add 0.5 to each value), because I think you intend to be specifying the center of your data array. (For an even number of pixels, the center of the data array falls on a pixel edge, so should be on a half pixel, not a whole pixel.) Given your plate scale, this error should be miniscule, but it’s an error nonetheless.
A probable error: Your CDELTi values raise a yellow flag because they are both negative. It’s common enough for CDELT1 (RA) to be negative, but it’s not common for CDELT2 (declination) to be negative. Now, it’s totally possible that your image is flipped, so this may not be an error. But, combined with the fact that you are rotating your image by 145 degrees rather than 55 degrees, this might indicate that you are incorrectly compensating for a CDELT with the wrong sign by adding an extra 90 degrees of rotation. You ought to double-check your image with this WCS object and confirm that RA and declination increase in the correct directions for the image.
A maybe error: I’m surprised that you’re not only using the Airy projection (AIR), but using it with a theta_b value of exactly 45 degrees. This doesn’t feel like a projection that would naturally match a telescope image, either ideally or as determined from calibration. What motivates this choice of projection? Perhaps the gnomonic projection (TAN) would be a more appropriate choice? I haven’t checked, but perhaps the difference between these projections is insignificant over your FOV.
Right, so as i said some of decision made were just from example code that’s froma AIRy’s projection came from (I also checked TAN but it made no difference). The 145 degree rotation is becouse it matched the rotation of the grid in astrometry (and it also made a little bit of sense becouse 90+55=145) but i think it’s not the way to go. In cdelt the signs were added so they were minimizing the error. But image is almost certainly fliped in RA. I have no idea how rotation in WCS is defined and how rotation in astrometry is defined so that’s why these weird values show up. This is the code for displaying the image:
import warnings
from matplotlib import pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS, FITSFixedWarning
from astropy.utils.data import get_pkg_data_filename
from matplotlib.colors import LogNorm
from astropy.coordinates import SkyCoord
from astropy import units as u
filename = get_pkg_data_filename('Neptune_WCS.fits')
hdu = fits.open(filename)[0]
with warnings.catch_warnings():
warnings.filterwarnings('ignore', message="'datfix' made the change",
category=FITSFixedWarning)
wcs = WCS(hdu.header)
# Wyodrębnij datę i godzinę obserwacji
obs_date = hdu.header.get('DATE-OBS', 'Nieznana data')
fig = plt.figure()
ax = fig.add_subplot(111, projection=wcs)
plt.imshow(hdu.data, origin='lower', cmap=plt.cm.gray, norm=LogNorm())
ax.coords.grid(color='white', alpha=0.5, linestyle='solid')
# Dodaj datę i godzinę obserwacji do tytułu
plt.title(f'Data obserwacji: {obs_date} UTC')
# Oblicz środek obrazu
image_width = hdu.header['NAXIS1']
image_height = hdu.header['NAXIS2']
image_center = [image_width / 2, image_height / 2]
# Funkcja onclick jest wywoływana, gdy użytkownik kliknie na obraz.
def onclick(event):
if event.inaxes is None:
return
ra, dec = wcs.all_pix2world(event.xdata, event.ydata, 0)
coord = SkyCoord(ra, dec, unit='deg')
print(f'RA: {coord.ra.to_string(unit=u.hour, sep=":")}, Dec: {coord.dec.to_string(unit=u.degree, sep=":")}')
# Calculate the center of the image in world coordinates
ra_center, dec_center = wcs.all_pix2world(image_center[0], image_center[1], 0)
center_coord = SkyCoord(ra_center, dec_center, unit='deg')
# Calculate the separation in arcseconds
separation = coord.separation(center_coord).arcsec
print(f'Distance from center: {separation} arcsec')
plt.plot(event.xdata, event.ydata, 'rx')
plt.draw()
# Połącz zdarzenie kliknięcia z funkcją onclick
cid = fig.canvas.mpl_connect('button_press_event', onclick)
plt.xlabel('RA')
plt.ylabel('Dec')
plt.show()
And the output with marked Neptune:
RA: 23:43:03.67426488, Dec: -3:13:32.68333391
Distance from center: 124.22115118636745 arcsec
Made the changes.
Now the output is:
RA error grown to 10s
but DEC is almost non present so that’s great.
Wondering what can i improve.
I’ll try to investigate deeper into WCS header made by astrometry:)
Few things i noticed.
CRVAL - is much more precise
CTYPE - im not sure what is it
CRPIX - that’s not the middle at all
CD(n) - i don’t know how they calculated those but they made it instead rotation matrix (How is grid rotated then)
This means an RA/dec coordinate system, using the gnomonic projection (TAN), plus Simple Imaging Polynomial (SIP) distortion coefficients. Presumably astrometry.net attempts to fit the actual distortion in the image if there is enough information in the star field.
The CRVAL values do not match what is shown in the screenshot – it’s not merely more precision – which at least is consistent with the fact that CRPIX is not in the center of data array. CRPIX defines the reference point for the WCS projection, which for the gnomonic projection and an ideal telescope would be the optical axis of the telescope. It’s of course rare for the optical axis of the telescope to be perfectly aligned with the center of the detector, so presumably astrometry.net attempts to determine the optimal reference point.
The CDij matrix is simply another way to encode rotation/scaling. For square pixels (CDELT1 = CDELT2), it’s equivalent to CDELT * PCij.
