Affecter des coordonnées WCS à une image FITS

Question

J'ai cherché la documentation comme un fou et je ne trouve pas de réponse pour celle-ci.

Je génère des images FITS en python et j'ai besoin d'assigner des coordonnées WCS à l'image. Je sais qu'il existe de nombreuses façons de faire cela en faisant correspondre des sources ponctuelles avec un catalogue connu, mais dans ce cas, je génère une carte de poussière, donc la correspondance de sources ponctuelles ne fonctionnera pas (pour autant que je sache). L'image est donc un tableau numérique Numpy de forme (240, 240). Il est écrit comme si (x et y coordonnées affectations sont un peu bizarre, il fonctionne en quelque sorte):

H, xedges, yedges = np.histogram2d(glat, glon, bins=[ybins, xbins], weights=Av) 
count, x, y = np.histogram2d(glat, glon, bins=[ybins, xbins]) 
H/=count 
hdu = pyfits.PrimaryHDU(H) 
hdu.writeto(filename) 

>>> print H.shape 
(240,240) 

Que tout fonctionne très bien lui-même. Pour l'attribution de coordonnées galactiques semble que tout ce que vous aurez besoin de faire est quelque chose comme:

glon_coords = np.linspace(np.amin(glon), np.amax(glon), 240) 
glat_coords = np.linspace(np.amin(glat), np.amax(glat), 240) 

Mais je ne comprends pas comment l'image FITS stocke ces coordonnées, donc je ne sais pas comment les écrire. J'ai essayé de les assigner à SAO DS9, mais pas de chance. J'ai juste besoin d'un moyen simple d'assigner ces coordonnées à l'image.

Merci pour toute aide que vous pouvez fournir.

Answer 1

Je vous conseille de commencer à utiliser astropy. Pour les besoins de votre projet, le package astropy.wcs peut vous aider à écrire un en-tête WCS FITS, et l'API astropy.io.fits est fondamentalement identique aux pyfits que vous utilisez actuellement. De plus, les pages d'aide sont excellentes, et tout ce que je suis sur le point de faire est de traduire leur page de construction WCS pour correspondre à votre exemple.

À votre question: FITS ne "étiquette" pas chaque pixel avec une coordonnée. Je suppose qu'il est possible de créer une table de recherche de pixels ou quelque chose comme ça, mais le actual WCS est une traduction algorithmic de X, Y pixels aux coordonnées astrométriques (dans votre cas "galactique"). A nice page is here.

L'exemple que je voudrais vous indiquer est ici:

http://docs.astropy.org/en/latest/wcs/index.html#building-a-wcs-structure-programmatically

Et voici mon non testé pseudocode pour votre projet:

# untested code 

from __future__ import division # confidence high 

# astropy 
from astropy.io import fits as pyfits 
from astropy import wcs 

# your code 
H, xedges, yedges = np.histogram2d(glat, glon, bins=[ybins, xbins], weights=Av) 
count, x, y = np.histogram2d(glat, glon, bins=[ybins, xbins]) 
H/=count 

# characterize your data in terms of a linear translation from XY pixels to 
# Galactic longitude, latitude. 

# lambda function given min, max, n_pixels, return spacing, middle value. 
linwcs = lambda x, y, n: ((x-y)/n, (x+y)/2) 

cdeltaX, crvalX = linwcs(np.amin(glon), np.amax(glon), len(glon)) 
cdeltaY, crvalY = linwcs(np.amin(glat), np.amax(glat), len(glat)) 

# wcs code ripped from 
# http://docs.astropy.org/en/latest/wcs/index.html 

w = wcs.WCS(naxis=2) 

# what is the center pixel of the XY grid. 
w.wcs.crpix = [len(glon)/2, len(glat)/2] 

# what is the galactic coordinate of that pixel. 
w.wcs.crval = [crvalX, crvalY] 

# what is the pixel scale in lon, lat. 
w.wcs.cdelt = numpy.array([cdeltX, cdeltY]) 

# you would have to determine if this is in fact a tangential projection. 
w.wcs.ctype = ["GLON-TAN", "GLAT-TAN"] 

# write the HDU object WITH THE HEADER 
header = w.to_header() 
hdu = pyfits.PrimaryHDU(H, header=header) 
hdu.writeto(filename)