def flatten(f, channel=0, freqaxis=0):
""" Flatten a fits file so that it becomes a 2D image. Return new header and data """
from astropy import wcs
naxis=f[0].header['NAXIS']
if naxis<2:
raise RadioError('Can\'t make map from this')
if naxis==2:
return f[0].header,f[0].data
w = wcs.WCS(f[0].header)
wn=wcs.WCS(naxis=2)
wn.wcs.crpix[0]=w.wcs.crpix[0]
wn.wcs.crpix[1]=w.wcs.crpix[1]
wn.wcs.cdelt=w.wcs.cdelt[0:2]
wn.wcs.crval=w.wcs.crval[0:2]
wn.wcs.ctype[0]=w.wcs.ctype[0]
wn.wcs.ctype[1]=w.wcs.ctype[1]
header = wn.to_header()
header["NAXIS"]=2
copy=('EQUINOX','EPOCH')
for k in copy:
r=f[0].header.get(k)
if r:
header[k]=r
slice=[]
for i in range(naxis,0,-1):
if i<=2:
slice.append(np.s_[:],)
elif i==freqaxis:
slice.append(channel)
else:
slice.append(0)
# slice=(0,)*(naxis-2)+(np.s_[:],)*2
return header,f[0].data[slice]
评论列表
文章目录