This is my first question here, so I'll try to explain it well.
I am working with some all-sky fits astronomical images, which are in galactic coordinates (in the AITOF projection). I would like to transform these fits images into RA,Dec coordinates, so that the image is rotated in the proper way (i.e. the plane of the galaxy will not be horizontal at the center of the map anymore, but kind of twisted).
Does anyone have any idea of how to do it? I was trying to do it in a not-so-elegant way, that I explain here:
I open the fits file, and the hdu extension where the image is. I read the header's appropriate keywords I create an array with the same shape, but every element will be now a couple of values, which will be the coordinates of each pixel. Then I transform each element into the new coordinates that I want (using the astropy.coordinates package) and finally I have to move each element in order to sort them in the new coordinates (first they were sorted in galactic coordinates, and now they must be sorted in celestial ones).
The piece of code that I had until now is very slow and I am sure there must be a better way to do what I intend to do:
import numpy as np
import pyfits as pf
from astropy.coordinates import ICRS, Galactic
from astropy import units as u
file = pf.open('IC443.fits')
crpixx = file[0].header["CRPIX1"] # X-axis reference pixel number
crpixy = file[0].header["CRPIX2"] # Y-axis reference pixel number
crvalx = file[0].header["CRVAL1"] # X-axis reference pixel value in coordinates (galactic longitude in this case)
crvaly = file[0].header["CRVAL2"] # Y-axis reference pixel value in coordinates (galactic latitude in this case)
stepx = file[0].header["CDELT1"] # X-axis increment per pixel
stepy = file[0].header["CDELT2"] # Y-axis increment per pixel
numpixelsx = file[0].header["NAXIS1"] # number of pixels in X axis
numpixelsy = file[0].header["NAXIS2"] # number of pixels in Y axis
tab = np.zeros([numpixelsx,numpixelsy,2]) # array of zeros with the same
# number of elements than the image, where each element is now a copuple
# of two values
def assign_coords(i,j):
# function to calculate coordinate correspondent to any pixel
coord_i = (-crpix1+i)*step1+crval1
coord_j = (-crpix2+j)*step2 + crval2
return coord_i, coord_j
for k, z in product(np.arange(numpixelsx), np.arange(numpixelsy)):
tab[k][z][0], tab[k][z][1] = assign_coords(k,z)
tab_temp_x = Galactic(tab[k][z][0], tab[k][z][1], unit=(u.degree, u.degree)).fk5.ra.value
tab_temp_y = Galactic(tab[k][z][0], tab[k][z][1], unit=(u.degree, u.degree)).fk5.dec.value
tab[k][z][0] = tab_temp_x # X-coord value in the new coordinates
tab[k][z][1] = tab_temp_y # Y-coord value in the new coordinates
After that, what I do not know is how to do the re-sort of the array in the proper way...
Thanks a lot guys!
