FITS概要

1各种图像格式

一般的天体照片等的“二维”图像数据，通常，以规定的格式作为文件被存储(因为多数在存储的时候压缩，所以也包含压缩方式)。这种格式有以下几种。

• jpeg:用于照片图像等。以每个像素24位的数据存储。
• bmp:Windows标准的格式。每个像素24位的数据。
• gif:用于交换图像数据的格式。每个像素8位。(也有gif动画等简易动画的扩展)
• png:作为代替gif的免费授权格式被开发。每个像素最多48位。(也有mng等用于动画的扩展)

2什么是FITS格式?

相对于上述的一般图像格式，FITS格式原本也是作为天文领域的图像交换用格式而被开发的。这在FITS (Flexible ImageTransport System)上也有所体现。因为最初被提出是在1979年，所以反而比一般的图像格式成立得早。但是，就像Flexible的名字一样，目标是不限于图像，也能用于数据交换的格式。也就是说，不仅限于图像。目前压缩还没有包含在规格中。例如，用编辑器打开FITS文件(下左)和bmp文件(下右)，就会变成下面这样。

3FITS文件的结构

FITS文件的一般结构如下。

以下，同样的扩展FITS元素(optional FITS structure)也可以接着。

3.1报头记录和数据记录

在报头中，用文本描述与FITS文件中的数据有关的信息。
关键词=值/注释
例如：

BITPIX的关键字是16，每像素数据的比特数是16,NAXIS的关键字是2，是二维数据。跟随报头的数据记录中有遵从报头描述的二进制数据。

3.2主要信息的例子

各关键字中记述的值按照逻辑值、整数值、实数值、字符串等的类别决定写法。举一个代表性的例子。

SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                   16 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                 1024
NAXIS2  =                 1024
EXTEND  =                    T
BIAS    =                 1000
FOCALLEN=                    0
APTAREA =                    0
APTDIA  =                    0
OBSERVER= ''
DATE-OBS= '2020-07-08T23:59:33.352574+00:00'
TIME-OBS= '23:59:33.352574'
SWCREATE= 'SATSIM 0.9.1'
SET-TEMP=                    0
COLORCCD=                    0
DISPCOLR=                    1
IMAGETYP= 'Light Frame'
CCDSFPT =                    1
XORGSUBF=                    0
YORGSUBF=                    0
CCDSUBFL=                    0
CCDSUBFT=                    0
XBINNING=                    2
CCDXBIN =                    2
YBINNING=                    2
CCDYBIN =                    2
CCD-TEMP=                    0
TEMPERAT=                    0
OBJECT  = ''
OBJCTRA = '21 40 49.6872 '
OBJCTDEC= '+61 13 38.9387 '
TELTKRA =   -934.7532981349218
TELTKDEC=    167.8345302059228
CENTAZ  =    350.9034367524328
CENTALT =    4.307733550377353
TELHA   = ''
LST     = ''
AIRMASS =                    0
SITELAT = '20.746111 N'
SITELONG= '156.431667 W'
ORCHCOMM= 'rate#[3:32]@[sim]%[filter wheel disabled by user]'
INSTRUME= ''
FILTER  = ''
EXPTIME =                  0.4
EXPOSURE=                  0.4
CBLACK  =                  958
CWHITE  =                 4000
SH-UT001= '2020-07-08T23:59:33.352574+00:00'
SH-UT002= '2020-07-08T23:59:33.752574+00:00'
CTYPE1  = 'RA---TAN' / 第一个坐标轴的类型
CTYPE2  = 'DEC--TAN' / 第二个坐标轴的类型
CRVAL1  =    340.2070298626763 / 第一个坐标轴的参考值
CRVAL2  =    71.22748298103676 / 第二个坐标轴的参考值
CRPIX1  =                512.0 / 第一个坐标轴的参考像素坐标
CRPIX2  =                512.0 / 第二个坐标轴的参考像素坐标
CDELT1  =              0.00034 / 第一个坐标轴的增量
CDELT2  =              0.00034 / 第二个坐标轴的增量
CD1_1   =              0.00034
CD1_2   =                  0.0
CD2_1   =                  0.0
CD2_2   =              0.00034
CROTA1  =                  0.0
CROTA2  =                  0.0
EQUINOX =               2000.0
BZERO   =                32768
BSCALE  =                    1
END

CRVAL代表坐标系中的参考值（coordinate reference value），通常是一个数值，用于确定图像中每个像素的实际物理坐标。CRVAL通常是指定天体坐标系的起始点，如赤道坐标系中的赤经或赤纬，或者是银道坐标系中的银道经或银道纬度。它通常以度为单位给出，并且可以是一个浮点数或整数
计算CRVAL

物理坐标 = CD * (像素坐标 - CRPIX) + CRVAL
CD是一个2×2矩阵，它描述了像素坐标和物理坐标之间的线性变换关系
CRVAL = CRPIX * CD + CRVAL_p

CD1_1代表在坐标系中X方向上相邻两个像素之间的物理距离。CD1_1的单位是度/像素
CDELT1是坐标轴上相邻两个像素之间的距离（以像素为单位），CROT1是图像中坐标轴的旋转角度（以度为单位），cos表示余弦函数。
CRVAL1和CRVAL2：代表参考点（即中心像素）在天球坐标系下的赤经和赤纬，单位为度。
CRPIX1和CRPIX2：代表参考点在像素坐标系下的坐标，即中心像素的x和y坐标。
CD1_1、CD1_2、CD2_1和CD2_2：代表像素坐标系和天球坐标系之间的关系，即坐标变换矩阵，单位为度/像素。通常情况下，这些值是对应的CD矩阵元素值的近似值。述像素坐标系旋转和缩放的矩阵元素

3.3关于FITS的EXTENSION(扩展FITS元素)

FITS最初的提案是以单纯的二维图像数据为对象，但是为了适应之后观测数据等的多样化，能够存储各种各样的数据，不断地扩展了规格。这就是所谓的EXTENSION。现在被标准化并正式成为FITS格式一部分的EXTENSION有以下这些。

• ASCII表:目录等表格形式的数据的存储
• 图像EXTENSION:同时存储图像数据和辅助信息
• 二进制表:表格形式包含二进制数据并存储
  除此之外，虽然不是EXTENSION，但是作为为了使天球坐标和数据排列对应的规约，WCS(World Coordinate System)达成了协议，制定了相关的关键字。

4fits文件图像

4.1一维FITS文件

一维FITS文件，即在黑器中描述为NAXIS =1的文件的代表性示例是频谱。通过数据分析软件等将分光器的数据记录为每个波长的强度的一维光谱数据时，比起图像，绘制成图表更容易理解其含义。例如，在下左的图像中，故意使存储了这种光谱数据的一维FITS文件具有宽度来显示，而在下右的图像中，将其作为图表来显示。

4.2二维FITS文件

二维FITS文件，即在报头中描述为NAXIS = 2的文件的代表性例子是二维图像，也是最常见的类型。例如，下左是星系的摄像数据的例子，下右是直方图分布。

4.3其他FITS文件

作为到此为止举出的代表性的FITS文件以外的FITS文件的例子有如下的东西。

• 作为NAXIS = 3的数据的示例，包括收集了二维图像数据的三维数据(右图)。
  
• 在NAXIS = 4的示例中，将(RA, DEC, Freq.， Stokes)的电波数据集中在一个FITS文件中。
• 3. NAXIS =0的例子不包含图像有目录等ASCI表。
• 4.即使NAXIS = 2也不是单纯的图像的例子是X射线卫星数据等，用二维二进制表存储事件的到达时间和脉冲高度等。

5fits文件

处理FITS文件需要使用对应的软件。天文用的软件大多能读取FITS文件进行简易分析，一般用的图像处理软件也有一部分能显示FITS文件，但大多只支持1、2维的FITS文件。看。另外，面向研究人员的分析软件(+对应这些的浏览器)也能显示和分析FITS文件。随着PC UNIX的普及，专业的分析软件在一般的PC上也能运行，门槛降低了。
FITS是一种灵活的图像传输格式，主要用于存储和传输天文图像。天文研究人员不能将图像以FITS的格式传输出去，需要将图像格式进行转换（PNG，JPEG）,转换后的图像将丢失重要的观测参数、WCS坐标参数等信息。需要将图像文件头中的信息提取出来，用AVM标准化的元数据标签进行转换，并添加到PNG或JPEG的文件头中。最终获得的PNG或者JPEG格式图像，不仅包含像素信息，而且包含了坐标参数信息。

FITS图像格式标准

一个FITS文件中包含整数条逻辑记录，一个FITS的逻辑记录长度为23040bits，每个记录以标题记录开始，标题和图像数据都在一个新的逻辑记录里开始。一个FITS文件包含三部分：主头文件和主数据单元（HDU）、扩展（Exten Extention） 和特殊记录（Special Record）,后两个是可选的，体现fits格式的灵活和可拓展性。
FITS标题用ASICII码编辑成每一个开始用8字符的关键字定义的类型（80字符）。标题记录语法为：keyword-value/comment。每一个头文件必须按照顺序含有SIMPLE、BITPIX、NAXIS、NAXISn和END五个关键词作为基本关键字，它指明了图像的基本信息：
SIMPLE：逻辑值，只能在主头文件头中，T表示该数据是标准FITS文件，F为否。
BITPIX：整数值，表示一个像素占的位数，8，16，32，-32，-64，表示图像数据的位宽，8表示字符串或无符号二进制整数，16，32分别表示16位和32位补码二进制整数，-32，-64分别表示IEEE单精度和双精度浮点数。
NAXIS:整数值，表示数据矩阵的纬度，0表示头文件后没有数据区。
NAXISn:表示数据矩阵大小即n轴的像素个数，若NAXI为0，则没有这一个关键词。
END:头文件结束标记。

FITS格式转换分析

FITS图像中每个像素的灰度范围通常很大，操作系统通常只能显示256级的灰度范围，首先要把灰度的变化范围压缩为0-255，压缩方式多种：线性变换、对数变换】

AVM标准的PNG(JPEG图像)的WCS信息添加

AVM天文可视化元数据标准是VAMP(Virtual Astronomy Multimedia Project)项目的底层架构。该标准以XMP（Extensible Metadata Platform）标准为基础开发而来，

5.1处理FITS的软件示例(Windows)

Qfits: JAHOU(Japanese Hands On Universe)用开发，IE插件也可以
Stellalmage:AstroArts公司的商用软件
AvisFV免费
Aladin

Astropy

可以使用Astropy库中的fits模块来读取Fits图像。下面是一个简单的示例：

from astropy.io import fits
# Open the FITS file
hdulist = fits.open('image.fits')
# Extract the image data
image_data = hdulist[0].data
# Close the FITS file
hdulist.close()

这样就可以将图像数据读取到image_data中，可以使用numpy库来操作图像。

如果你想要更高效的读取，可以使用fits_tools库中的fits_image_reader来读取，fits_image_reader能够支持多线程读取，能够大大提高读取速度。

!pip install fits_tools
from fits_tools import fits_image_reader

image_data = fits_image_reader('image.fits')

Astropy库中的WCS模块来计算FITS图像中的时角坐标。下面是一个简单的示例

from astropy.io import fits
from astropy.wcs import WCS
import numpy as np

# Open the FITS file
hdulist = fits.open('image.fits')

# Extract the WCS information
w = WCS(hdulist[0].header)

# Create a set of pixel coordinates
x, y = np.meshgrid(np.arange(hdulist[0].header['NAXIS1']),
                   np.arange(hdulist[0].header['NAXIS2']))

# Convert the pixel coordinates to sky coordinates
ra, dec = w.all_pix2world(x, y, 0)

# Close the FITS file
hdulist.close()

创建一个FITS文件保存图像和header信息

from astropy.io import fits
import numpy as np

# Generate some dummy data for the image
data = np.random.rand(512, 512)

# Create a new HDU (Header Data Unit)
hdu = fits.PrimaryHDU(data)

# Update the header with some dummy information
hdu.header['OBJECT'] = 'Test Object'
hdu.header['TELESCOP'] = 'Test Telescope'

# Write the HDU to a FITS file
hdu.writeto('image.fits', overwrite=True)

FITS文件的header信息中写入注释

hdu.header['OBJECT'].comment = 'The name of the observed object'

FITS图像中DATE-OBS
FITS图像的DATE-OBS关键字代表的是观测的开始时间。可以使用Astropy库中的fits模块读取FITS文件的header信息，然后获取DATE-OBS的值。
DATE-OBS的格式通常是YYYY-MM-DDTHH:MM:SS，其中T是分隔符。可以使用Python的内置库datetime来生成这样的格式时间。
下面是一个简单的示例：

import datetime

# Get the current date and time
now = datetime.datetime.now()

# Format the date and time as a string in the required format
date_obs = now.strftime('%Y-%m-%dT%H:%M:%S')

FITS文件中的DATE-OBS通常是由观测仪器（如望远镜）的控制软件自动生成的。该信息表示观测数据的开始时间，是由观测仪器内部的时钟生成的，通常是准确的。

如果生成的FITS文件不包含DATE-OBS信息，则可以使用上面我提供的代码手动生成。这些信息是可以在后期更新的，因此即使生成时不准确也可以在后续修正。

=================================
**FITS文件中的不需要手动写入END，Astropy库会自动在FITS文件的末尾写入END
**
卫星目标的时角坐标通常是通过以下方法计算的：

1. 计算卫星目标在天文坐标系（如赤道坐标系或银河坐标系）中的位置。
2. 计算观测站的位置，包括纬度、经度和高度。
3. 计算卫星目标在观测站观测时的方位角，通常使用空间三维坐标。
4. 计算时角，即太阳系经过某一点的方向与卫星目标的方向之间的夹角。
5. 将时角转换为卫星目标的时角坐标。
   这是一个简化版的计算过程，实际上，在计算过程中还需要考虑到地球自转、太阳偏转等问题。因此，实际计算中，需要使用相关的天文学模型和算法，并在计算过程中使用相关的软件工具，如Astropy库中的skycoord模块。

from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time

# 卫星目标在天文坐标系中的位置
coord = SkyCoord(ra=10.0*u.degree, dec=45.0*u.degree, frame='icrs')

# 观测站位置
observing_location = EarthLocation(lat=35.0*u.deg, lon=-105.0*u.deg, height=2.0*u.km)

# 观测时间
observing_time = Time('2022-12-22T00:00:00')

# 计算卫星目标的AltAz坐标
altaz = coord.transform_to(AltAz(obstime=observing_time, location=observing_location))

# 计算时角
hour_angle = altaz.ra - observing_location.lon

# 时角转换为度
hour_angle = hour_angle.to(u.degree)

在上面的代码中，我们首先定义了卫星目标在天文坐标系中的位置，并定义了观测站的位置和观测时间。然后，我们计算了卫星目标的AltAz坐标，并使用AltAz坐标计算出了时角，最后，将时角转换为度。

要计算在FITS图像中检测到的卫星位置的时角坐标，需要先获得该图像的WCS信息，并使用WCS信息将卫星的位置从像素坐标转换为天球坐标。然后，可以将天球坐标转换为时角坐标。

from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u

# 读取FITS图像数据
hdul = fits.open('your_image.fits')
data = hdul[0].data

# 读取FITS图像的WCS信息
wcs = WCS(hdul[0].header)

# 假设卫星位置在像素坐标(x, y)
x = ...
y = ...

# 将像素坐标转换为天球坐标
sky_coord = wcs.pixel_to_world(x, y)

# 获取时角坐标
ra = sky_coord.ra.to(u.hourangle)
dec = sky_coord.dec

print(f"卫星位置的时角坐标: {ra}, 赤纬坐标: {dec}")

计算图像中某一个像素位置的天球坐标