python处理卫星云图获取亮温值

在处理卫星云图hdf文件时，matlab非常好用，但随着python的流行，还是想学习一下python处理卫星云图数据获取亮温值。

我处理数据的前提是：

1. 根据最佳路径数据集，获取经纬度和时间。
2. 根据时间找对应的hdf文件，获取其中的经纬度，然后根据hdf文件中的经度确定是哪颗卫星的。
3. 在全圆盘标对称文件经纬度对照表dat文件获取经纬度时加上相应的经度。
4. 最后根据最佳路径的经纬度获取dat文件中最近格点的经纬度索引，再通过hdf文件获得亮温值。

1.处理hdf文件

import h5py
hdfFile = h5py.File('F:/新建文件夹/Satellite_Imagery/IR_data/韦森特/FY2D_FDI_ALL_NOM_20120723_1945.hdf', 'r')
db1 = hdfFile['/CALChannelIR1']
hw1 = hdfFile['/NOMChannelIR1']
db2 = hdfFile['/CALChannelIR2']
hw2 = hdfFile['/NOMChannelIR2']
db3 = hdfFile['/CALChannelIR3']
hw3 = hdfFile['/NOMChannelIR3']
db4 = hdfFile['/CALChannelIR4']
hw4 = hdfFile['/NOMChannelIR4']
infoh = hdfFile['/NomFileInfo']

a = hdfFile['NomFileInfo'][()]
print(a[0][3])  # 获取纬度
print(a[0][4])  # 获取经度

对于风云卫星的数据格式，可以查阅相应的资料，获取静止气象卫星的数据格式。这里不再详细描述。

亮温值的获取分为四个IR通道和一个可见光通道，不同通道的用处不同，根据自己的需求选取数据。

对于获取卫星经纬度等数据，需要，根据NomFileInfo的数据集分析。

a = hdfFile['NomFileInfo'][()]
print(a)

[(b'FY2D', b'NOM', b'Nomalized Projection Full Disc Image', 0., 86.5, 35785864., b'fit', 2012, 7, 10, 12, 34, 0, 2012, 7, 23, 19, 45, 1, 2012, 7, 23, 19, 58, 58, 2012, 7, 23, 20, 0, 27, 6378137., 3.5e-05, 0.00014, 298.25722356, b'NOM Fit HDF5')]

对应信息见下图数据集。

获取亮温值：

获取通道1的亮温值，代码如下：

hw = hw1[()]  # 通道IR1，定标表
db = db1[()]  # 获取定标表的值
tb = np.zeros(shape=(2288, 2288))  # 2288*2288的图像每个具体的亮温值
for i in range(2288):
    for j in range(2288):
        if hw[i][j] == 65535 or hw[i][j] == 65534:  # 判断缺省值
            tb[i][j] = 0
        else:
            a = hw[i][j]
            tb[i][j] = db[0][a]
print(tb)
tb = tb.T  # 注意需要转置

2. 处理dat文件

对于FY2气象卫星而言，经纬度查找表文件（NOM_ITG_2288_2288(0E0N)_LE.zip）可以从网上下载。数据解压之后，里面有三个文件，分别是：
在数据说明文件中，详细说明了经纬度查找表数据的使用方法，同时针对FY2系列卫星的经度进行了说明，由于上述经纬度数据是以中心点为0度经度生成的，所以对于FY2的数据，需要在经度数据上加上卫星所在的经度。

5	FY-2G	104.5°E
6	FY-2H	79°E

对于不同的星可以根据hdf文件中的经度，判断对应的查找表的经度。

下面以H星为例，读取dat文件，获取网格经纬度

lonlatfile = 'F:/Satellite_Imagery/Code/NOM_ITG_2288_2288(0E0N)_LE.dat'
    with open(lonlatfile, 'rb') as f:
        #lon_fy = np.fromfile(f, count=2288 * 2288, dtype='float32') + 79  # 先存经度，根据卫星的不同加上对应的经度值
        #lat_fy = np.fromfile(f, count=2288 * 2288, dtype='float32')  # 再存纬度
        data = np.fromfile(f, dtype='float32')
        data = data.reshape([2288, 2288, 2], order='F')
    #lon = lon_fy.reshape([2288, 2288], order='F')
    #lat = lat_fy.reshape([2288, 2288], order='F')
    lon = data[:, :, 0] + 104.5
    lat = data[:, :, 1]

这里需要注意的是，需要先取出来经度，再取出纬度，这是经纬度在dat文件中的存储顺序，然后加上卫星的经度，这里是79，故经度加79，获取经纬度。

# 经度表
[[379. 379. 379. ... 379. 379. 379.]
 [379. 379. 379. ... 379. 379. 379.]
 [379. 379. 379. ... 379. 379. 379.]
 ...
 [379. 379. 379. ... 379. 379. 379.]
 [379. 379. 379. ... 379. 379. 379.]
 [379. 379. 379. ... 379. 379. 379.]]
# 纬度表 
 [[300. 300. 300. ... 300. 300. 300.]
 [300. 300. 300. ... 300. 300. 300.]
 [300. 300. 300. ... 300. 300. 300.]
 ...
 [300. 300. 300. ... 300. 300. 300.]
 [300. 300. 300. ... 300. 300. 300.]
 [300. 300. 300. ... 300. 300. 300.]]

这样就获取了dat文件中的网格经纬度，即2288*2288的对照表，这样根据最佳路径数据集中的经纬度，找到网格中最近的经纬度。

注：因为网格存储的经纬度不是等间距的，所以不能直接近似找最近的点，需要根据角度计算距离，找到网格上离所求经纬度最近的网格点，找到对应的索引值，再根据索引值找到亮温值。

以上就是我通过python读取hdf，dat文件，根据经纬度获取卫星云图的亮温值的过程，欢迎大家讨论。

绘制登陆时的卫星云图亮温值请参考python绘制登陆时的卫星云图（TBB）。