[AstroPy] Question regarding WCS coordinates
Eric L. N. Jensen
ejensen1 at swarthmore.edu
Wed Jul 3 16:53:14 EDT 2013
Hi John,
I had to do something like this recently, converting images from a simulation into FITS files with WCS in the header so I could compare directly with ALMA data.
Below is the code I used - it's a little more complicated than your case since you aren't creating a FITS header from scratch, and maybe you don't need a third (velocity/frequency) dimension in your image, but hopefully it will give you an idea of how to proceed.
To give credit where it's due, this code derives largely from an IDL routine supplied by Kees Dullemond with his radmc3d code - I translated it to Python and modified it somewhat.
If you need more information about WCS keywords in the FITS standard, this paper is a good reference:
Pence et al. 2010 A&A 524, 40, http://www.aanda.org/articles/aa/full_html/2010/16/aa15362-10/aa15362-10.html
Hope this helps,
Eric
P.S. In the code below, the 'read image' routine returns a dictionary with various bits of information about the image, as well as the numpy array of the image itself - adapt as needed.
def radmcimage_to_fits(imagename,
fitsname,
dpc,
arcsec=False,
mas=False,
restfreq=None,
vel_offset=None,
radeg=None,
decdeg=None,):
import numpy as np
from astropy.io import fits
cc = 2.9979245800000e10 # Light speed [cm/s]
pc = 3.08572e18 # Parsec [cm]
im = readimage(imagename)
pixdeg_x = 180.0*(im['sizepix_x']/(dpc*pc))/np.pi
pixdeg_y = 180.0*(im['sizepix_y']/(dpc*pc))/np.pi
freqhz = 1e4*cc/im['lambda']
image = im['image']
n_freq = len(freqhz)
# Because of the way that Python stores arrays, we need to swap around
# our data array before writing it, so that it is indexed in the order
# (freq, y, x), rather than (x, y, freq) - the last index in the Numpy
# array will become AXIS1 in the FITS image, and so on.
if im['stokes']:
(nx, ny, nstokes, nfreq) = np.shape(image)
newimage = np.transpose(image, (3, 2, 1, 0))
else:
(nx, ny, nfreq) = np.shape(image)
newimage = np.transpose(image, (2, 1, 0))
#
# Compute the conversion factor from erg/cm^2/s/Hz/ster to Jy/pixel
#
pixsurf_ster = pixdeg_x*pixdeg_y * (np.pi/180.)**2
factor = 1e+23 * pixsurf_ster
# And scale the image array accordingly:
image_in_jypix = factor * newimage
#
# Make FITS header information:
#
header = fits.Header()
header['BTYPE'] = 'Intensity'
header['BSCALE'] = 1
header['BZERO'] = 0
header['BUNIT'] = 'JY/PIXEL'
if radeg is not None:
header['EPOCH'] = 2000.
header['LONPOLE'] = 180.
header['CTYPE1'] = 'RA---SIN'
header['CRVAL1'] = radeg
# Figure out what units we are using per pixel:
if arcsec:
unit = 'arcsec'
multiplier = 3600.
elif mas:
unit = 'mas'
multiplier = 3.6e6
else:
unit = 'deg'
multiplier = 1
header['CDELT1'] = multiplier*pixdeg_x
header['CUNIT1'] = unit
#
# ...Zero point of coordinate system
#
header['CRPIX1'] = 1.0*((nx+1)/2)
if decdeg is not None:
header['CTYPE2'] = 'DEC--SIN'
header['CRVAL2'] = decdeg
header['CDELT2'] = multiplier*pixdeg_y
header['CUNIT2'] = unit
#
# ...Zero point of coordinate system
#
header['CRPIX2'] = 1.0* ((ny+1)/2)
# If the keyword is set for rest frequency, add that to the header as
# it gives a reference point for, e.g., the velocity of the frequency
# steps relative to some line rest frequency. Note that even though
# the keyword used here is 'restfreq', the FITS standard for this
# field is RESTFRQ, i.e. with no second 'E'.
if restfreq is not None:
header['RESTFRQ'] = (float(restfreq), 'Rest frequency of this transition')
#
# ...Frequency
#
# If they have specified a velocity offset, we shift all of the
# specified frequencies by that amount. This allows, e.g., mapping a
# transition of a particular CO line in radmc3d, but then accounting
# for the fact that a real object might have a center-of-mass radial
# velocity offset. Note that we do *not* shift the input 'restfreq'
# keyword, as that is presumed to be the real frequency of that
# transition.
if vel_offset is not None:
# Assume input velocity is in km/s, so
# specify speed of light in those
# units:
c = 2.9979245800000e5 # Light speed [km/s]
freqhz *= (1. - vel_offset/c)
header['COMMENT'] = "Velocity offset of %0.2f km/s applied to model." % vel_offset
if n_freq > 1:
# multiple frequencies - set up the header keywords to define the
# third axis as frequency
header['CTYPE3'] = 'FREQ'
header['CUNIT3'] = 'Hz'
header['CRPIX3'] = 1.0
header['CRVAL3'] = freqhz[0]
# Calculate the frequency step, assuming equal steps between all:
delta_freq = freqhz[1] - freqhz[0]
header['CDELT3'] = delta_freq
else: # only one frequency
header['RESTFREQ'] = freqhz
# Make a FITS file!
#
fits.writeto(fitsname, image_in_jypix, header, output_verify='fix')
On Jul 3, 2013, at 2:36 PM, John ZuHone wrote:
> Hi all,
>
> This question is probably more about WCS than it is about Astropy, but I intend to use Astropy to do it so I thought I would ask.
>
> I have some projected maps of gas temperature of a galaxy cluster that I made from FLASH simulation data using yt that I have put into FITS files. I also have some simulated Chandra event files made from the MARX Chandra simulator, of the same system.
>
> I would like to match the WCS coordinates of the two files. The simulated events file already has WCS coordinates in RA and DEC, and this is the coordinate system I would like to use for the temperature map. The information I have for the temperature map is the cell spacing in kpc, and the redshift of the simulated cluster, so with the chosen cosmology I also have the distance.
>
> How can I go about this? Sorry, I realize this is a bit involved (at least it seems that way to me).
>
> Best,
>
> John ZuHone
> _______________________________________________
> AstroPy mailing list
> AstroPy at scipy.org
> http://mail.scipy.org/mailman/listinfo/astropy
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/astropy/attachments/20130703/5798c809/attachment.html>
More information about the AstroPy
mailing list