Python: 如何批量预处理FY4A L1 DISK和REGC产品?(辐射定标/裁剪/GLT校正/HDF5转TIFF文件等)

目录

01 前言

1.1 想要说

1.2 Requirements

1.3 程序适用数据集

02 函数说明

2.1 读取HDF5文件某一数据集

2.2 读取HDF5文件数据集属性

2.3 对FY4A数据集进行辐射定标

2.4 基于官方地理对照表获取经纬度数据(仅适用DISK)

2.5 依据行列号计算经纬度数据(仅适用DISK)

2.6 基于*GEO*.HDF获取行列号数据集以计算经纬度数据(仅适用REGC)

2.7 裁剪

2.8 GLT校正

2.9 输出TIFF文件

03 完整程序

3.1 基于官方地理对照表的DISK(全圆盘)FY4A L1产品的批量预处理

3.2 基于GEO文件的REGC(中国区域)的FY4A L1产品的批量预处理

3.3 基于行列号计算经纬度数据集(DISK/全圆盘)的FY4A L1产品的单个文件预处理

04 结果展示 

4.1 ArcGIS显示

4.2 ENVI 显示

---

01 前言

1.1 想要说

这篇博客,其实我想了很久要不要写，因为关于如何去预处理FY4A产品我已经花了很久的时间去收集和整理资料，花费的时间和精力其实不少。这方面数据的处理教程本来就少(包括我自己去搜集资料都是困难重重)，我自己去处理产品也有受益于其它博客的影响。还有很多想说，暂且憋着吧。

另外，由于博客参考了部分其它博客和资料，时间跨度有点长，若有遗漏很抱歉。

1.2 Requirements

很抱歉我并没有放到Github上，因为时间和精力，所以需要大家自行安装以下模块，版本仅供参考，若你进行处理出现问题，请查看是否是因为版本更新问题所致.

numpy 1.23.5

h5py 3.9.0

gdal 3.3.3

scipy 1.8.1

猜测gdal比较难以安装, 因为它因为一些原因被放到osgeo模块，如果实在难以安装，可以尝试直接下载.whl文件进行安装，这里不再解释，自行搜索。若无法安装，也可以使用rasterio等模块，人不能在一棵树上吊死，但是你需要付出代价 ==> 修改部分代码。

1.3 程序适用数据集

你的数据集或许也是FY4A L1产品，但是需要注意，FY4A有很多产品，你需要注意自己的产品是否可以进行当前处理，如果不能，请查阅其它方法，本博客仅供参考。

DISK(全圆盘)文件示例

FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200601000000_20200601001459_4000M_V0001.HDF

形如AGRI*DISK*FDI*MULT*NOM*4000M*.HDF

解释如下：

 其中，全圆盘很重要，如果你的文件非DISK，那么你需要考虑当前程序是否适用这些方法；另外是MULT多通道，对于我的文件，一共有14个通道，形如下：

 若你的文件不存在多个波段，你需要了解你的文件属性，或许它并非一级产品；

最后，当前程序仅仅适用于4000M分辨率的产品，若你的产品非4000M产品，你可参考本程序，但是在关键代码中我将提醒你替换关键参数（参数从官方文档中查找）。

REGC（中国区域）文件示例

FY4A-_AGRI--_N_REGC_1047E_L1-_FDI-_MULT_NOM_20200601003000_20200601003417_4000M_V0001.HDF

形如：*AGRI*REGC*FDI*MULT*NOM*4000M*.HDF

对于DISK的GLT校正，需要使用到官方地理对照表；对于REGC的GLT校正，需要使用到官方的*GEO*.HDF定位配准文件.（后续会详细说明，请勿用担心）

关于FY4A数据集的下载请参考官方文档： 风云四号气象卫星数据存档文件名命名规范_20170509.pdf

如打不开，请参照如下步骤：
 

02 函数说明

在这些自定义函数上我会参杂一些官方文档等供大家理解；函数内可能会调用其它函数，为了确保程序正常运行，请将理解所有函数尽管有时候它很简单。另外，我并没有把函数进行提炼写成类，因为我认为这会增加阅读的困难，不方便大家直观的理解其中的原理。

2.1 读取HDF5文件某一数据集

1. # 读取HDF5文件数据集
2. def h5_data_get(hdf_path, dataset_name):
3. """
4. 读取HDF5文件中的数据集
5. :param hdf_path: HDF5文件的路径
6. :param dataset_name: HDF5文件中数据集的名称
7. :return: 数据集
8. """
10. with h5py.File(hdf_path, 'r') as f:
11. dataset = f[dataset_name][:]
13. return dataset

对于数据集名称，你需要通过HDFExp，Panoply等软件查看会更为方便，如下：

 当然，你也可以通过h5py模块进行查看，这里不再详细说明。

2.2 读取HDF5文件数据集属性

1. # 读取HDF5文件数据集属性
2. def h5_attr_get(hdf_path, dataset_name, attr_name):
3. """
4. 获取HDF5文件数据集的属性
5. :param hdf_path: HDF5文件的路径
6. :param dataset_name: HDF5文件中数据集的名称
7. :param attr_name: HDF5文件中数据集属性的名称
8. :return: 数据集属性
9. """
11. with h5py.File(hdf_path, 'r') as f:
12. attr = f[dataset_name].attrs[attr_name]
14. return attr

2.3 对FY4A数据集进行辐射定标

1. # 对FY4A数据集进行辐射定标
2. def fy4a_calibration(hdf_path, nom_channel_name, cal_channel_name):
3. """
4. 获取FY4A数据集并对其进行辐射定标
5. :param hdf_path: HDF5文件的路径
6. :param nom_channel_name: 待定标的数据集名称
7. :param cal_channel_name: 用于定标的数据集名称
8. :return: 辐射定标后的数据集
9. """
11. # 读取数据集
12. nom_channel = h5_data_get(hdf_path, nom_channel_name)
13. cal_channel = h5_data_get(hdf_path, cal_channel_name)
15. # 读取数据集属性
16. nom_min, nom_max = h5_attr_get(hdf_path, nom_channel_name, 'valid_range')
17. cal_min, cal_max = h5_attr_get(hdf_path, cal_channel_name, 'valid_range')
18. nom_fill_value = h5_attr_get(hdf_path, nom_channel_name, 'FillValue')[0]
19. cal_fill_value = h5_attr_get(hdf_path, cal_channel_name, 'FillValue')[0]
21. # 数据集掩码和填充值预准备
22. nom_mask = (nom_channel > nom_min) & (nom_channel < nom_max) & (nom_channel != nom_fill_value)
23. cal_mask = (cal_channel > cal_min) | (cal_channel < cal_max)
24. cal_channel[cal_channel == cal_fill_value] = int(cal_min - 10)
26. # 辐射定标
27. target_channel = np.zeros_like(nom_channel, dtype=np.float32)
28. target_channel[nom_mask] = cal_channel[cal_mask][nom_channel[nom_mask]]
30. # 无效值处理(包括不在范围及填充值)
31. target_channel[~nom_mask] = np.nan
32. target_channel[target_channel == int(cal_min - 10)] = np.nan
34. return target_channel

我认为有必要对这一段代码进行解释方便大家的理解。

首先请让我说明常规的辐射定标：

 对于Landsat部分系列有如下公式：

 那么对于FY4A L1数据集，它是如何进行辐射定标的呢？

首先，他主要给定两种数据集，其一是波段影像：

 以及对应各个波段的辐射定标文件（每一个文件均为一维，类似列表）：

 他们的关系是：由于NOMChannel中DN值为整数，因此将NOMChannel中DN值作为CALChannel的索引值，通过索引值在CALChannel可查找到该DN值对应的辐射定标之后的值。

因此实际上辐射定标后的值已经在CALChannel已经计算出来，需要自行查找。

当然，我的程序中进行了很多填充值和有效值范围的限定，这看起来并不清晰，对于为什么最后进行np.nan的赋值是因为原始的DN值为整数，而nan为浮点数因此需要预先创建mask掩码。

当然，如果你的文件存在辐射定标的参数数据集，如下：

 第0列为前文的Gain增益系数，第1列为偏置Offset；不同的行代表不同的波段，第0行通常代表NOMChannel01数据集。

为什么在一个HDF5文件中存在两种辐射定标的方法呢？于是我尝试通过两种方法（CALchannel和CALIBRATION_COEF(SCALE+OFFSET)）去对同一数据集进行辐射定标，对于辐射定标后的数据集进行相关系数的计算以及目视查看，发现相关系数达到0.9999987，目视查看发现二者完全一致。因此认为两种方法没有任何区别

2023-11-11修正：上述想法错误，实际上CALChannel在前面的可见光波段和一些近红外波段与CaliBration_COEF(SCALE+OFFSET)是完全一致的；但是对于后面的热红外波段就不一样，CALChannel计算得到的是辐射亮温（我们熟知的地表温度），而CaliBration_COEF(SCALE+OFFSET)计算得到的是辐射亮度。

由于基于CALIBRATION_COEF(SCALE+OFFSET）的辐射定标较为简单，代码我也没有保存，因此这里就不贴了，需要注意辐射定标时顺便考虑填充值和有效值范围。

2.4 基于官方地理对照表获取经纬度数据(仅适用DISK)

1. # 读取地理对照表(.raw, FY4A)
2. def read_glt(raw_path, shape):
3. """
4. 该函数用于读取地理对照表
5. :param raw_path: 地理查找表的路径
6. :param shape: 地理查找表的形状
7. :return: 二维数组(lon, lat)
8. """
9. # 读取二进制数据，将其转换为两个2D数组（经度和纬度）
10. with open(raw_path, 'rb') as f:
11. raw_data = np.fromfile(f, dtype=np.float64) # 读取二进制数据, '
12. raw_lat = raw_data[::2].reshape(shape)
13. raw_lon = raw_data[1::2].reshape(shape)
14. raw_lat[(raw_lat < -90.0) | (raw_lat > 90.0)] = np.nan
15. raw_lon[(raw_lon < -180.0) | (raw_lon > 180.0)] = np.nan
17. """
18. 官方文档中的说明：文件从北向南按行（从西到东）填写，每个网格对应 16 字节，前8字节为经度值，后 8 字节为纬度值，double 类型，高字节在前
19. -----------------------------------------------------------------------------------------------------------------
20. | 通过查看输出的lon, lat发现, lon的范围为[-90, 90], lat的范围为[-180, 180];
21. | 通过读取二进制流时，如果以高字节在前(Big-Endian/小端模式)的方式(dtype=‘
22. | 即默认是小端模式)读取数据时，会出现错误的结果;
23. -----------------------------------------------------------------------------------------------------------------
24. 综上所述, 官方文档的说明存在错误: 应该为前8字节为纬度值，后 8 字节为经度值，double 类型,低字节在前;
25. 参考: https://blog.csdn.net/modabao/article/details/83834580
26. """
28. return raw_lon, raw_lat

在了解上方的代码前，我们需要了解地理查找表，在早期的博客中，我曾经写过：
ENVI: 如何创建GLT文件并基于GLT对图像进行几何校正？

但是其实那时候我还是不太能理解地理查找表（地理对照表，或者它就是GLT文件因为地理位置查找表文件 (Geographic Lookup Table, GLT) ）的含义。
我所下载的地理对照表文件(.raw后缀)为：

FullMask_Grid_4000.raw

现在我们来查看官方文档对于地理对照表的原理诠释：

但是需要注意是，官方文档存在勘误（代码中也有详细说明）。这是在经过反复确认得到的结论，当然也参考加以确认：
Python3.FY-4A标称（NOM）数据提取

正确说明应当是：每一个像元值(网格值)共16个字节，前8字节是纬度值，后8字节是经度值，double类型，低字节在前(np默认即是小端模式即低字节在前)；

但是4KM影像的行列数是没有问题的，为2748×2748.

 至此，我基本理解地理对照表中含有什么信息了，我之前一直认为地理对照表含有经纬度、投影等诸多信息，现在看来仅仅只有经纬度信息。上述代码就是将地理对照表转换为两个2748×748的数据集包括经度和纬度。

关于地理对照表的下载：FullMask_Grid_4000.raw

2.5 依据行列号计算经纬度数据(仅适用DISK)

需要说明的是，依据行列号计算经纬度数据与前文基于官方地理对照表获取经纬度数据 功能都是一致的，都是要获取经纬度数据集。

这里仅仅是从更底层去探讨经纬度的获取（当然如果你的数据集已经有经纬度数据集那么这些处理对于你来说都是无需的）。

参考：FY-4A建立中国区域图像行列号转经纬度的经纬度查找表进行几何校正

但是上文是关于REGC（中国区域） 而且是基于*GEO*HDF获取中国区域的行列号数据集。并不适用DISK，因此我尝试以类似的思路去进行DISK的经纬度数据集计算。经过和前文的基于官方地理对照表获取经纬度数据相比，二者的相关性系数在0.95上下波动，所以还是不错的。

1. def get_lon_lat(dataset):
2. """
3. FY4A数据集经纬度计算(基于FY4A数据格式说明书的公式)
4. :param dataset: FY4A数据集
5. :return: 经纬度
6. """
7. # 行列号 ==> 经纬度, 公式由FY4A数据格式说明书给出
9. # 获取数据集的行列号
10. rows, cols = dataset.shape
11. # 生成行列号矩阵
12. col_mesh, row_mesh = np.meshgrid(np.arange(cols), np.arange(rows))
14. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
15. ea = 6378.137 # 地球长半轴, 单位km
16. eb = 6356.7523 # 地球短半轴, 单位km
17. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
18. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
19. coff = 1373.5 # 列偏移
20. cfac = 10233137 # 列比例因子
21. loff = 1373.5 # 行偏移
22. lfac = 10233137 # 行比例因子
24. # step1. 求x, y
25. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
26. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
28. # step2. 求sd, sn, s1, s2, s4, sxy
29. sd = sqrt(
30. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
31. )
32. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
33. s1 = h - sn * cos(x) * cos(y)
34. s2 = sn * sin(x) * cos(y)
35. s3 = -sn * sin(y)
36. sxy = sqrt(power(s1, 2) + power(s2, 2))
38. # step3. 求lon, lat
39. lon = lon0 + arctan(s2 / s1) * 180 / pi
40. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
42. return lon, lat

(若上述参数与官方文档存在不符 ，请一定告知我)

行列号数据集我并非如REGC数据集获取*GEO*HDF中的行列号数据集，而是直接由波段NOMChannel的行列数生成行列号数据集，这是由于官方文档说明行列转换经纬度说明书的附录说明4KM范围的行列号范围为0~2747，这与DISK数据集的行列数2748完全匹配，因此我尝试如此操作具有一定的依据，事实上两种方法的区别差异比较小，也有一定的说服力。这是二者生成的图像对比查看：

上述代码仅适用于4KM分辨率。关于参数的选择和计算公式(由行列计算经纬度)请查阅官方文档：
这贴出当前部分参数说明：

 官方文档的下载查阅：FY4A成像仪标称上行列号和经纬度的互相转换方法.V2.pdf

 

 2.6 基于*GEO*.HDF获取行列号数据集以计算经纬度数据(仅适用REGC)

关于*GEO*.HDF的示例：
FY4A-_AGRI--_N_REGC_1047E_L1-_GEO-_MULT_NOM_20200606194500_20200606194917_4000M_V0001.HDF
其中GEO表示定位配准产品：

 适用Panoply打开如下：

 或许你已经发现，其中包含行列号数据集(二维)，而且*GEO*.HDF还包含具体的时间，说明应当是对应不同*REGC*.HDF文件，有不同的*GEO*.HDF文件。但是我通过多幅影像查看，发现似乎所有的REGC的行列号数据集是没有区别的，因此在后续的script中我对于所有REGC均适用同一GEO文件的行列号数据集。至于存在不同时间的GEO文件，大概存在于其它数据集例如其中不仅包含行列号数据集还包含NOMSatelliteAzimuth卫星方位角、NOMSunGlintAngle太阳高度角等，这些应该与大气校正有关。所以很明显，因为能力有限，我并没有对数据集进行大气校正，这是非常遗憾。

1. def read_geo(geo_path, row_name, col_name):
2. """
3. 基于FY4A的GEO文件(定位配准产品中的对地图像定位产品)获取像元级的行列号矩阵, 并通过公式计算获取对应的经纬度数据集
4. :param geo_path: FY4A的GEO文件路径
5. :param row_name: 行号数据集名称
6. :param col_name: 列号数据集名称
7. :return: 经纬度数据集, 形式(lon, lat)
8. """
9. # 读取数据集及其属性
10. original_row_mesh = h5_data_get(geo_path, row_name)
11. original_col_mesh = h5_data_get(geo_path, col_name)
12. row_fill_value = h5_attr_get(geo_path, col_name, 'FillValue')[0]
13. col_fill_value = h5_attr_get(geo_path, row_name, 'FillValue')[0]
14. # 有效值掩码
15. mask = (original_row_mesh != row_fill_value) & (original_col_mesh != col_fill_value)
16. # 有效值
17. row_mesh = original_row_mesh[mask]
18. col_mesh = original_col_mesh[mask]
20. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
21. ea = 6378.137 # 地球长半轴, 单位km
22. eb = 6356.7523 # 地球短半轴, 单位km
23. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
24. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
25. coff = 1373.5 # 列偏移
26. cfac = 10233137 # 列比例因子
27. loff = 1373.5 # 行偏移
28. lfac = 10233137 # 行比例因子
30. # step1. 求x, y
31. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
32. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
34. # step2. 求sd, sn, s1, s2, s4, sxy
35. sd = sqrt(
36. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
37. )
38. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
39. s1 = h - sn * cos(x) * cos(y)
40. s2 = sn * sin(x) * cos(y)
41. s3 = -sn * sin(y)
42. sxy = sqrt(power(s1, 2) + power(s2, 2))
44. # step3. 求lon, lat
45. lon = lon0 + arctan(s2 / s1) * 180 / pi
46. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
48. # 输出
49. out_lon = np.zeros_like(original_row_mesh, dtype=np.float32)
50. out_lat = np.zeros_like(original_row_mesh, dtype=np.float32)
51. out_lon[mask] = lon
52. out_lat[mask] = lat
53. out_lon[~mask] = np.nan
54. out_lat[~mask] = np.nan
56. return out_lon, out_lat

此处的代码处理与此前DISK的依据行列号数据集计算经纬度数据集一致，所用文档也是一致的，仅仅在于行列号数据集的获取不同，DISK是通过波段的行列数自行计算的每像元行列号数据集，而此处是REGC没法如此处理（因为REGC为DISK范围的子集），因此我们需要通过GEO文件获取行列号数据集。

2.7 裁剪

这一步的裁剪很关键，在大多数博客的GLT校正中，或多或少都是通过代码或者ENVI生成GLT地理查找表，然后适用ENVI IDL二次打开接口或者ENVI现有工具实现GLT校正。大多数博客都提及GLT校正中在创建GLT文件前需要经纬度数据集进行裁剪，因为DISK全圆盘的缘故，如下圆外均为NoData值。

(源: FY-4A 静止卫星圆盘数据几何校正)

如果直接进行GLT校正会导致GLT校正失败，因此我们需要对经纬度数据集进行选定。

例如博客：FY-4A建立中国区域图像行列号转经纬度的经纬度查找表进行几何校正

1. def clip(original_band, original_lon, original_lat, lon_lat_range=None):
2. """
3. 裁剪, 非常规意义的矩形裁剪, 而是基于经纬度范围的裁剪并进行重组
4. :param original_band: 原始数据集
5. :param original_lon: 对应的经度矩阵
6. :param original_lat: 对应的纬度矩阵
7. :param lon_lat_range: 经纬度范围, 默认为中国区域(形式:[lon_min, lon_max, lat_min, lat_max])
8. :return: 裁剪并重组好的数据集, 包括(band, lon, lat)
9. """
10. if lon_lat_range is None:
11. lon_lat_range = [73, 136, 18, 54] # 中国区域(为兼容REGC, 仅南至海南省主岛屿)
12. lon_min, lon_max, lat_min, lat_max = lon_lat_range
14. # 获取经纬度范围内的掩码
15. mask = (original_lon >= lon_min) & (original_lon <= lon_max) & (original_lat >= lat_min) & (original_lat <= lat_max)
16. valid_lon = original_lon[mask]
17. valid_lat = original_lat[mask]
18. valid_band = original_band[mask]
19. valid_num = np.sum(mask) # True为1, False为0
20. valid_num_sqrt = int(sqrt(valid_num))
21. # 重组数据
22. reform_lon = valid_lon[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
23. reform_lat = valid_lat[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
24. reform_band = valid_band[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
26. return reform_band, reform_lon, reform_lat

此处的裁剪并不与其它博客裁剪一致，而是指定裁剪的经纬度范围（默认为中国区域），而非指定矩形范围（即类似行列号的裁剪） 。另外，为方便进行GLT校正，我将裁剪后的经纬度数据集和波段数据集进行重组，这部分理解较为困难，与后面的校正函数相关联，这里不予解释。

2.8 GLT校正

1. def glt_warp(original_dataset, original_lon, original_lat, out_res, method='nearest'):
2. """
3. 基于地理查找表的校正
4. :param original_dataset: 待校正的数据
5. :param original_lon: 原始数据的经度
6. :param original_lat: 原始数据的纬度
7. :param out_res: 输出数据的分辨率
8. :param method: 插值方法, 默认为最近邻插值('linear', 'nearest', 'cubic')
9. :return: 校正后的数据
10. """
12. # 校正数据的最大最小经纬度
13. lon_min, lon_max = np.nanmin(original_lon), np.nanmax(original_lon)
14. lat_min, lat_max = np.nanmin(original_lat), np.nanmax(original_lat)
16. # 生成目标数据的经纬度矩阵
17. grid_lon, grid_lat = np.meshgrid(
18. np.arange(lon_min, lon_max, out_res),
19. np.arange(lat_max, lat_min, -out_res),
20. )
22. # 插值
23. interp_dataset = griddata(
24. (original_lon.ravel(), original_lat.ravel()), # 原始数据的经纬度, ravel()将多维数组转换为一维数组
25. original_dataset.ravel(), # 原始数据
26. (grid_lon, grid_lat), # 目标数据的经纬度
27. method=method, # 插值方法
28. fill_value=np.nan, # 用于填充超出边界的点
29. )
31. return interp_dataset

此GLT校正思路参考某IDL思路，由于IDL代码并非我写且并未询问公开，这里不展示。（IDL生成的影像与Python实现的影像的相关性系数在0.96，不同时间点有差异，差异可能在于插值算法的具体实现有关）。

其思路为通过传入的经纬度数据集获取最大最小的经纬度值，并通过传入的输出分辨率进行输出经纬度数据集的创建（例如最大最小经度为100， 90，而输出分辨率为1°，那么每一行的纬度均为90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100；对于纬度亦是如此。其实类似于前文的行列号数据集的创建）。创建好输出经纬度数据集，接下来就是依据原始经纬度数据集(传入的)以及对应的波段数据集将波段中每一个像元的值填入输出经纬度数据集对应的位置，对于缺失位置的像元通过插值进行填充。

为了方便理解，这里贴出未进行插值但是已经将输入的波段数据集填入输出经纬度数据集对应位置的区域的输出影像：

 2.9 输出TIFF文件

1. # 输出TIF
2. def write_tiff(out_path, dataset, transform, nodata=np.nan):
3. """
4. 输出TIFF文件
5. :param out_path: 输出文件的路径
6. :param dataset: 待输出的数据
7. :param transform: 坐标转换信息(形式:[左上角经度, 经度分辨率, 旋转角度, 左上角纬度, 旋转角度, 纬度分辨率])
8. :param nodata: 无效值
9. :return: None
10. """
12. # 创建文件
13. driver = gdal.GetDriverByName('GTiff')
14. out_ds = driver.Create(out_path, dataset[0].shape[1], dataset[0].shape[0], len(dataset), gdal.GDT_Float32)
16. # 设置基本信息
17. out_ds.SetGeoTransform(transform)
18. out_ds.SetProjection('WGS84')
20. # 写入数据
21. for i in range(len(dataset)):
22. out_ds.GetRasterBand(i + 1).WriteArray(dataset[i]) # GetRasterBand()传入的索引从1开始, 而非0
23. out_ds.GetRasterBand(i + 1).SetNoDataValue(nodata)
24. out_ds.FlushCache()

由于此处的输出TIFF函数进行了改写方便多波段影像的输出，因此如果你是单波段的输出，应该传入[band]而非单波段band数据集，因为我们的形式为[band1, band2, band3,···]。详细信息查看后文主程序。

另外需要说明的是，由于输出的数据集中所有像元值的范围普遍小于1，因此当你使用ArcGIS打开进行显示可能会出现全灰色的情况：

 拉伸方法如下：

 如果在ENVI中打开，则可以正常显示，因为其会默认进行2%线性拉伸。

03 完整程序

为避免错误，对于不同的处理我均将函数搬运在同一个脚本中，为了规范和后续代码的维护，我建议你将函数专门放置于专门的位置，例如：

 然后再行调用。

3.1 基于官方地理对照表的DISK(全圆盘)FY4A L1产品的批量预处理

1. # @炒茄子 2023-07-04
3. """
4. 当前脚本用于处理FY4A数据，包括辐射定标、基于地理查找表(.raw文件)的校正、GeoTIFF文件的创建等, 最后输出为GeoTIFF文件(WGS84)
5. --批量处理
6. """
8. # 函数
9. # @炒茄子 2023-07-09
11. import os
12. import glob
13. import h5py
14. import numpy as np
15. from numpy import sqrt, sin, cos, power, arctan, pi
16. from osgeo import gdal
17. from scipy.interpolate import griddata
20. # 读取HDF5文件数据集
21. def h5_data_get(hdf_path, dataset_name):
22. """
23. 读取HDF5文件中的数据集
24. :param hdf_path: HDF5文件的路径
25. :param dataset_name: HDF5文件中数据集的名称
26. :return: 数据集
27. """
29. with h5py.File(hdf_path, 'r') as f:
30. dataset = f[dataset_name][:]
32. return dataset
35. # 读取HDF5文件数据集属性
36. def h5_attr_get(hdf_path, dataset_name, attr_name):
37. """
38. 获取HDF5文件数据集的属性
39. :param hdf_path: HDF5文件的路径
40. :param dataset_name: HDF5文件中数据集的名称
41. :param attr_name: HDF5文件中数据集属性的名称
42. :return: 数据集属性
43. """
45. with h5py.File(hdf_path, 'r') as f:
46. attr = f[dataset_name].attrs[attr_name]
48. return attr
51. # 对FY4A数据集进行辐射定标
52. def fy4a_calibration(hdf_path, nom_channel_name, cal_channel_name):
53. """
54. 获取FY4A数据集并对其进行辐射定标
55. :param hdf_path: HDF5文件的路径
56. :param nom_channel_name: 待定标的数据集名称
57. :param cal_channel_name: 用于定标的数据集名称
58. :return: 辐射定标后的数据集
59. """
61. # 读取数据集
62. nom_channel = h5_data_get(hdf_path, nom_channel_name)
63. cal_channel = h5_data_get(hdf_path, cal_channel_name)
65. # 读取数据集属性
66. nom_min, nom_max = h5_attr_get(hdf_path, nom_channel_name, 'valid_range')
67. cal_min, cal_max = h5_attr_get(hdf_path, cal_channel_name, 'valid_range')
68. nom_fill_value = h5_attr_get(hdf_path, nom_channel_name, 'FillValue')[0]
69. cal_fill_value = h5_attr_get(hdf_path, cal_channel_name, 'FillValue')[0]
71. # 数据集掩码和填充值预准备
72. nom_mask = (nom_channel > nom_min) & (nom_channel < nom_max) & (nom_channel != nom_fill_value)
73. cal_mask = (cal_channel > cal_min) | (cal_channel < cal_max)
74. cal_channel[cal_channel == cal_fill_value] = int(cal_min - 10)
76. # 辐射定标
77. target_channel = np.zeros_like(nom_channel, dtype=np.float32)
78. target_channel[nom_mask] = cal_channel[cal_mask][nom_channel[nom_mask]]
80. # 无效值处理(包括不在范围及填充值)
81. target_channel[~nom_mask] = np.nan
82. target_channel[target_channel == int(cal_min - 10)] = np.nan
84. return target_channel
87. def get_lon_lat(dataset):
88. """
89. FY4A数据集经纬度计算(基于FY4A数据格式说明书的公式)
90. :param dataset: FY4A数据集
91. :return: 经纬度
92. """
93. # 行列号 ==> 经纬度, 公式由FY4A数据格式说明书给出
95. # 获取数据集的行列号
96. rows, cols = dataset.shape
97. # 生成行列号矩阵
98. col_mesh, row_mesh = np.meshgrid(np.arange(cols), np.arange(rows))
100. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
101. ea = 6378.137 # 地球长半轴, 单位km
102. eb = 6356.7523 # 地球短半轴, 单位km
103. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
104. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
105. coff = 1373.5 # 列偏移
106. cfac = 10233137 # 列比例因子
107. loff = 1373.5 # 行偏移
108. lfac = 10233137 # 行比例因子
110. # step1. 求x, y
111. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
112. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
114. # step2. 求sd, sn, s1, s2, s4, sxy
115. sd = sqrt(
116. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
117. )
118. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
119. s1 = h - sn * cos(x) * cos(y)
120. s2 = sn * sin(x) * cos(y)
121. s3 = -sn * sin(y)
122. sxy = sqrt(power(s1, 2) + power(s2, 2))
124. # step3. 求lon, lat
125. lon = lon0 + arctan(s2 / s1) * 180 / pi
126. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
128. return lon, lat
131. def read_geo(geo_path, row_name, col_name):
132. """
133. 基于FY4A的GEO文件(定位配准产品中的对地图像定位产品)获取像元级的行列号矩阵, 并通过公式计算获取对应的经纬度数据集
134. :param geo_path: FY4A的GEO文件路径
135. :param row_name: 行号数据集名称
136. :param col_name: 列号数据集名称
137. :return: 经纬度数据集, 形式(lon, lat)
138. """
139. # 读取数据集及其属性
140. original_row_mesh = h5_data_get(geo_path, row_name)
141. original_col_mesh = h5_data_get(geo_path, col_name)
142. row_fill_value = h5_attr_get(geo_path, col_name, 'FillValue')[0]
143. col_fill_value = h5_attr_get(geo_path, row_name, 'FillValue')[0]
144. # 有效值掩码
145. mask = (original_row_mesh != row_fill_value) & (original_col_mesh != col_fill_value)
146. # 有效值
147. row_mesh = original_row_mesh[mask]
148. col_mesh = original_col_mesh[mask]
150. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
151. ea = 6378.137 # 地球长半轴, 单位km
152. eb = 6356.7523 # 地球短半轴, 单位km
153. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
154. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
155. coff = 1373.5 # 列偏移
156. cfac = 10233137 # 列比例因子
157. loff = 1373.5 # 行偏移
158. lfac = 10233137 # 行比例因子
160. # step1. 求x, y
161. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
162. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
164. # step2. 求sd, sn, s1, s2, s4, sxy
165. sd = sqrt(
166. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
167. )
168. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
169. s1 = h - sn * cos(x) * cos(y)
170. s2 = sn * sin(x) * cos(y)
171. s3 = -sn * sin(y)
172. sxy = sqrt(power(s1, 2) + power(s2, 2))
174. # step3. 求lon, lat
175. lon = lon0 + arctan(s2 / s1) * 180 / pi
176. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
178. # 输出
179. out_lon = np.zeros_like(original_row_mesh, dtype=np.float32)
180. out_lat = np.zeros_like(original_row_mesh, dtype=np.float32)
181. out_lon[mask] = lon
182. out_lat[mask] = lat
183. out_lon[~mask] = np.nan
184. out_lat[~mask] = np.nan
186. return out_lon, out_lat
189. # 读取地理对照表(.raw, FY4A)
190. def read_glt(raw_path, shape):
191. """
192. 该函数用于读取地理对照表
193. :param raw_path: 地理查找表的路径
194. :param shape: 地理查找表的形状
195. :return: 二维数组(lon, lat)
196. """
197. # 读取二进制数据，将其转换为两个2D数组（经度和纬度）
198. with open(raw_path, 'rb') as f:
199. raw_data = np.fromfile(f, dtype=np.float64) # 读取二进制数据, '
200. raw_lat = raw_data[::2].reshape(shape)
201. raw_lon = raw_data[1::2].reshape(shape)
202. raw_lat[(raw_lat < -90.0) | (raw_lat > 90.0)] = np.nan
203. raw_lon[(raw_lon < -180.0) | (raw_lon > 180.0)] = np.nan
205. """
206. 官方文档中的说明：文件从北向南按行（从西到东）填写，每个网格对应 16 字节，前8字节为经度值，后 8 字节为纬度值，double 类型，高字节在前
207. -----------------------------------------------------------------------------------------------------------------
208. | 通过查看输出的lon, lat发现, lon的范围为[-90, 90], lat的范围为[-180, 180];
209. | 通过读取二进制流时，如果以高字节在前(Big-Endian/小端模式)的方式(dtype=‘
210. | 即默认是小端模式)读取数据时，会出现错误的结果;
211. -----------------------------------------------------------------------------------------------------------------
212. 综上所述, 官方文档的说明存在错误: 应该为前8字节为纬度值，后 8 字节为经度值，double 类型,低字节在前;
213. 参考: https://blog.csdn.net/modabao/article/details/83834580
214. """
216. return raw_lon, raw_lat
219. def clip(original_band, original_lon, original_lat, lon_lat_range=None):
220. """
221. 裁剪, 非常规意义的矩形裁剪, 而是基于经纬度范围的裁剪并进行重组
222. :param original_band: 原始数据集
223. :param original_lon: 对应的经度矩阵
224. :param original_lat: 对应的纬度矩阵
225. :param lon_lat_range: 经纬度范围, 默认为中国区域(形式:[lon_min, lon_max, lat_min, lat_max])
226. :return: 裁剪并重组好的数据集, 包括(band, lon, lat)
227. """
228. if lon_lat_range is None:
229. lon_lat_range = [73, 136, 18, 54] # 中国区域(为兼容REGC, 仅南至海南省主岛屿)
230. lon_min, lon_max, lat_min, lat_max = lon_lat_range
232. # 获取经纬度范围内的掩码
233. mask = (original_lon >= lon_min) & (original_lon <= lon_max) & (original_lat >= lat_min) & (original_lat <= lat_max)
234. valid_lon = original_lon[mask]
235. valid_lat = original_lat[mask]
236. valid_band = original_band[mask]
237. valid_num = np.sum(mask) # True为1, False为0
238. valid_num_sqrt = int(sqrt(valid_num))
239. # 重组数据
240. reform_lon = valid_lon[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
241. reform_lat = valid_lat[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
242. reform_band = valid_band[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
244. return reform_band, reform_lon, reform_lat
247. def glt_warp(original_dataset, original_lon, original_lat, out_res, method='nearest'):
248. """
249. 基于地理查找表的校正
250. :param original_dataset: 待校正的数据
251. :param original_lon: 原始数据的经度
252. :param original_lat: 原始数据的纬度
253. :param out_res: 输出数据的分辨率
254. :param method: 插值方法, 默认为最近邻插值('linear', 'nearest', 'cubic')
255. :return: 校正后的数据
256. """
258. # 校正数据的最大最小经纬度
259. lon_min, lon_max = np.nanmin(original_lon), np.nanmax(original_lon)
260. lat_min, lat_max = np.nanmin(original_lat), np.nanmax(original_lat)
262. # 生成目标数据的经纬度矩阵
263. grid_lon, grid_lat = np.meshgrid(
264. np.arange(lon_min, lon_max, out_res),
265. np.arange(lat_max, lat_min, -out_res),
266. )
268. # 插值
269. interp_dataset = griddata(
270. (original_lon.ravel(), original_lat.ravel()), # 原始数据的经纬度, ravel()将多维数组转换为一维数组
271. original_dataset.ravel(), # 原始数据
272. (grid_lon, grid_lat), # 目标数据的经纬度
273. method=method, # 插值方法
274. fill_value=np.nan, # 用于填充超出边界的点
275. )
277. return interp_dataset
280. # 输出TIF
281. def write_tiff(out_path, dataset, transform, nodata=np.nan):
282. """
283. 输出TIFF文件
284. :param out_path: 输出文件的路径
285. :param dataset: 待输出的数据
286. :param transform: 坐标转换信息(形式:[左上角经度, 经度分辨率, 旋转角度, 左上角纬度, 旋转角度, 纬度分辨率])
287. :param nodata: 无效值
288. :return: None
289. """
291. # 创建文件
292. driver = gdal.GetDriverByName('GTiff')
293. out_ds = driver.Create(out_path, dataset[0].shape[1], dataset[0].shape[0], len(dataset), gdal.GDT_Float32)
295. # 设置基本信息
296. out_ds.SetGeoTransform(transform)
297. out_ds.SetProjection('WGS84')
299. # 写入数据
300. for i in range(len(dataset)):
301. out_ds.GetRasterBand(i + 1).WriteArray(dataset[i]) # GetRasterBand()传入的索引从1开始, 而非0
302. out_ds.GetRasterBand(i + 1).SetNoDataValue(nodata)
303. out_ds.FlushCache()
306. # 准备工作
307. in_path = r'===>\DISK' # FY4A L1(DISK全圆盘)产品所在文件夹
308. raw_path = r'===>\FullMask_Grid_4000.raw' # 地理对照表(地理查找表)文件路径
309. out_path = r'===>\output' # 输出文件夹
310. out_res = 0.036 # 0.036为输出的分辨率(°/度), 约为4km
312. # 检索HDF5文件
313. paths = glob.iglob(in_path + r'\*DISK*FDI*.HDF')
314. # 读取地理查找表
315. lon, lat = read_glt(raw_path, (2748, 2748)) # 地理查找表的shape与band的shape相同,或可从官方文档中获取((2748, 2748))
316. # 循环处理每一个HDF5文件
317. for path in paths:
318. # HDF5文件的文件名
319. hdf_name = os.path.basename(path).split('.')[0]
320. # 读取HDF5文件
321. hdf = h5py.File(path, 'r')
322. nom_channels_name = [key for key in hdf.keys() if 'NOMChannel' in key] # 读取各个波段的通道名
323. cal_channels_name = [key for key in hdf.keys() if 'CALChannel' in key] # 读取对应各个波段的辐射定标通道名
324. channel_names = zip(nom_channels_name, cal_channels_name)
325. # 循环处理每一个通道
326. bands = []
327. transform = None
328. for nom_channel_name, cal_channel_name in channel_names:
329. # 辐射定标
330. band = fy4a_calibration(path, nom_channel_name, cal_channel_name)
331. # 裁剪
332. clip_band, clip_lon, clip_lat = clip(band, lon, lat)
333. # GLT校正
334. processed_band = glt_warp(clip_band, clip_lon, clip_lat, out_res, method='linear')
335. # 仿射变换参数
336. transform = (
337. np.nanmin(clip_lon), # 左上角经度
338. out_res, # x方向分辨率
339. 0, # 旋转角度
340. np.nanmax(clip_lat), # 左上角纬度
341. 0, # 旋转角度
342. -out_res, # y方向分辨率, 由于自左上角开始(纬度往下逐渐减小), 因此为负;
343. )
344. bands.append(processed_band)
345. # 将处理后的数据和地理查找表的信息写入GeoTIFF文件
346. hdf_out_path = out_path + '\\' + hdf_name + '.tif'
347. write_tiff(hdf_out_path, bands, transform)

3.2 基于GEO文件的REGC(中国区域)的FY4A L1产品的批量预处理

1. # @炒茄子 2023-07-08
3. """
4. 当前脚本用于对FY4A L1(REGC中国区域)产品进行辐射定标, 基于FY4A的定位配准产品(*GEO*.HDF)进行地理查找表的校正, 最后输出为GeoTIFF文件(WGS84)
5. --批量处理
6. """
8. # 函数
9. # @炒茄子 2023-07-09
10. import glob
11. import os
12. import h5py
13. import numpy as np
14. from numpy import sqrt, sin, cos, power, arctan, pi
15. from osgeo import gdal
16. from scipy.interpolate import griddata
19. # 读取HDF5文件数据集
20. def h5_data_get(hdf_path, dataset_name):
21. """
22. 读取HDF5文件中的数据集
23. :param hdf_path: HDF5文件的路径
24. :param dataset_name: HDF5文件中数据集的名称
25. :return: 数据集
26. """
28. with h5py.File(hdf_path, 'r') as f:
29. dataset = f[dataset_name][:]
31. return dataset
34. # 读取HDF5文件数据集属性
35. def h5_attr_get(hdf_path, dataset_name, attr_name):
36. """
37. 获取HDF5文件数据集的属性
38. :param hdf_path: HDF5文件的路径
39. :param dataset_name: HDF5文件中数据集的名称
40. :param attr_name: HDF5文件中数据集属性的名称
41. :return: 数据集属性
42. """
44. with h5py.File(hdf_path, 'r') as f:
45. attr = f[dataset_name].attrs[attr_name]
47. return attr
50. # 对FY4A数据集进行辐射定标
51. def fy4a_calibration(hdf_path, nom_channel_name, cal_channel_name):
52. """
53. 获取FY4A数据集并对其进行辐射定标
54. :param hdf_path: HDF5文件的路径
55. :param nom_channel_name: 待定标的数据集名称
56. :param cal_channel_name: 用于定标的数据集名称
57. :return: 辐射定标后的数据集
58. """
60. # 读取数据集
61. nom_channel = h5_data_get(hdf_path, nom_channel_name)
62. cal_channel = h5_data_get(hdf_path, cal_channel_name)
64. # 读取数据集属性
65. nom_min, nom_max = h5_attr_get(hdf_path, nom_channel_name, 'valid_range')
66. cal_min, cal_max = h5_attr_get(hdf_path, cal_channel_name, 'valid_range')
67. nom_fill_value = h5_attr_get(hdf_path, nom_channel_name, 'FillValue')[0]
68. cal_fill_value = h5_attr_get(hdf_path, cal_channel_name, 'FillValue')[0]
70. # 数据集掩码和填充值预准备
71. nom_mask = (nom_channel > nom_min) & (nom_channel < nom_max) & (nom_channel != nom_fill_value)
72. cal_mask = (cal_channel > cal_min) | (cal_channel < cal_max)
73. cal_channel[cal_channel == cal_fill_value] = int(cal_min - 10)
75. # 辐射定标
76. target_channel = np.zeros_like(nom_channel, dtype=np.float32)
77. target_channel[nom_mask] = cal_channel[cal_mask][nom_channel[nom_mask]]
79. # 无效值处理(包括不在范围及填充值)
80. target_channel[~nom_mask] = np.nan
81. target_channel[target_channel == int(cal_min - 10)] = np.nan
83. return target_channel
86. def get_lon_lat(dataset):
87. """
88. FY4A数据集经纬度计算(基于FY4A数据格式说明书的公式)
89. :param dataset: FY4A数据集
90. :return: 经纬度
91. """
92. # 行列号 ==> 经纬度, 公式由FY4A数据格式说明书给出
94. # 获取数据集的行列号
95. rows, cols = dataset.shape
96. # 生成行列号矩阵
97. col_mesh, row_mesh = np.meshgrid(np.arange(cols), np.arange(rows))
99. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
100. ea = 6378.137 # 地球长半轴, 单位km
101. eb = 6356.7523 # 地球短半轴, 单位km
102. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
103. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
104. coff = 1373.5 # 列偏移
105. cfac = 10233137 # 列比例因子
106. loff = 1373.5 # 行偏移
107. lfac = 10233137 # 行比例因子
109. # step1. 求x, y
110. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
111. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
113. # step2. 求sd, sn, s1, s2, s4, sxy
114. sd = sqrt(
115. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
116. )
117. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
118. s1 = h - sn * cos(x) * cos(y)
119. s2 = sn * sin(x) * cos(y)
120. s3 = -sn * sin(y)
121. sxy = sqrt(power(s1, 2) + power(s2, 2))
123. # step3. 求lon, lat
124. lon = lon0 + arctan(s2 / s1) * 180 / pi
125. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
127. return lon, lat
130. def read_geo(geo_path, row_name, col_name):
131. """
132. 基于FY4A的GEO文件(定位配准产品中的对地图像定位产品)获取像元级的行列号矩阵, 并通过公式计算获取对应的经纬度数据集
133. :param geo_path: FY4A的GEO文件路径
134. :param row_name: 行号数据集名称
135. :param col_name: 列号数据集名称
136. :return: 经纬度数据集, 形式(lon, lat)
137. """
138. # 读取数据集及其属性
139. original_row_mesh = h5_data_get(geo_path, row_name)
140. original_col_mesh = h5_data_get(geo_path, col_name)
141. row_fill_value = h5_attr_get(geo_path, col_name, 'FillValue')[0]
142. col_fill_value = h5_attr_get(geo_path, row_name, 'FillValue')[0]
143. # 有效值掩码
144. mask = (original_row_mesh != row_fill_value) & (original_col_mesh != col_fill_value)
145. # 有效值
146. row_mesh = original_row_mesh[mask]
147. col_mesh = original_col_mesh[mask]
149. # 基本参数(均只用于FY4A 4000m分辨率, 其它分辨率需要依据官方说明修改参数)
150. ea = 6378.137 # 地球长半轴, 单位km
151. eb = 6356.7523 # 地球短半轴, 单位km
152. h = 42164 # 卫星高度, 单位km, 即地心到卫星质心的距离
153. lon0 = 104.7 # 投影中心经度, 单位度, 也即卫星星下点所在的经度
154. coff = 1373.5 # 列偏移
155. cfac = 10233137 # 列比例因子
156. loff = 1373.5 # 行偏移
157. lfac = 10233137 # 行比例因子
159. # step1. 求x, y
160. x = (pi * (col_mesh - coff)) / (180 * (2 ** (-16)) * cfac)
161. y = (pi * (row_mesh - loff)) / (180 * (2 ** (-16)) * lfac)
163. # step2. 求sd, sn, s1, s2, s4, sxy
164. sd = sqrt(
165. power(h * cos(x) * cos(y), 2) - (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2)) * (h ** 2 - ea ** 2)
166. )
167. sn = (h * cos(x) * cos(y) - sd) / (power(cos(y), 2) + power(ea / eb, 2) * power(sin(y), 2))
168. s1 = h - sn * cos(x) * cos(y)
169. s2 = sn * sin(x) * cos(y)
170. s3 = -sn * sin(y)
171. sxy = sqrt(power(s1, 2) + power(s2, 2))
173. # step3. 求lon, lat
174. lon = lon0 + arctan(s2 / s1) * 180 / pi
175. lat = arctan(power(ea / eb, 2) * s3 / sxy) * 180 / pi
177. # 输出
178. out_lon = np.zeros_like(original_row_mesh, dtype=np.float32)
179. out_lat = np.zeros_like(original_row_mesh, dtype=np.float32)
180. out_lon[mask] = lon
181. out_lat[mask] = lat
182. out_lon[~mask] = np.nan
183. out_lat[~mask] = np.nan
185. return out_lon, out_lat
188. # 读取地理对照表(.raw, FY4A)
189. def read_glt(raw_path, shape):
190. """
191. 该函数用于读取地理对照表
192. :param raw_path: 地理查找表的路径
193. :param shape: 地理查找表的形状
194. :return: 二维数组(lon, lat)
195. """
196. # 读取二进制数据，将其转换为两个2D数组（经度和纬度）
197. with open(raw_path, 'rb') as f:
198. raw_data = np.fromfile(f, dtype=np.float64) # 读取二进制数据, '
199. raw_lat = raw_data[::2].reshape(shape)
200. raw_lon = raw_data[1::2].reshape(shape)
201. raw_lat[(raw_lat < -90.0) | (raw_lat > 90.0)] = np.nan
202. raw_lon[(raw_lon < -180.0) | (raw_lon > 180.0)] = np.nan
204. """
205. 官方文档中的说明：文件从北向南按行（从西到东）填写，每个网格对应 16 字节，前8字节为经度值，后 8 字节为纬度值，double 类型，高字节在前
206. -----------------------------------------------------------------------------------------------------------------
207. | 通过查看输出的lon, lat发现, lon的范围为[-90, 90], lat的范围为[-180, 180];
208. | 通过读取二进制流时，如果以高字节在前(Big-Endian/小端模式)的方式(dtype=‘
209. | 即默认是小端模式)读取数据时，会出现错误的结果;
210. -----------------------------------------------------------------------------------------------------------------
211. 综上所述, 官方文档的说明存在错误: 应该为前8字节为纬度值，后 8 字节为经度值，double 类型,低字节在前;
212. 参考: https://blog.csdn.net/modabao/article/details/83834580
213. """
215. return raw_lon, raw_lat
218. def clip(original_band, original_lon, original_lat, lon_lat_range=None):
219. """
220. 裁剪, 非常规意义的矩形裁剪, 而是基于经纬度范围的裁剪并进行重组
221. :param original_band: 原始数据集
222. :param original_lon: 对应的经度矩阵
223. :param original_lat: 对应的纬度矩阵
224. :param lon_lat_range: 经纬度范围, 默认为中国区域(形式:[lon_min, lon_max, lat_min, lat_max])
225. :return: 裁剪并重组好的数据集, 包括(band, lon, lat)
226. """
227. if lon_lat_range is None:
228. lon_lat_range = [73, 136, 18, 54] # 中国区域(为兼容REGC, 仅南至海南省主岛屿)
229. lon_min, lon_max, lat_min, lat_max = lon_lat_range
231. # 获取经纬度范围内的掩码
232. mask = (original_lon >= lon_min) & (original_lon <= lon_max) & (original_lat >= lat_min) & (original_lat <= lat_max)
233. valid_lon = original_lon[mask]
234. valid_lat = original_lat[mask]
235. valid_band = original_band[mask]
236. valid_num = np.sum(mask) # True为1, False为0
237. valid_num_sqrt = int(sqrt(valid_num))
238. # 重组数据
239. reform_lon = valid_lon[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
240. reform_lat = valid_lat[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
241. reform_band = valid_band[:valid_num_sqrt ** 2].reshape(valid_num_sqrt, valid_num_sqrt)
243. return reform_band, reform_lon, reform_lat
246. def glt_warp(original_dataset, original_lon, original_lat, out_res, method='nearest'):
247. """
248. 基于地理查找表的校正
249. :param original_dataset: 待校正的数据
250. :param original_lon: 原始数据的经度
251. :param original_lat: 原始数据的纬度
252. :param out_res: 输出数据的分辨率
253. :param method: 插值方法, 默认为最近邻插值('linear', 'nearest', 'cubic')
254. :return: 校正后的数据
255. """
257. # 校正数据的最大最小经纬度
258. lon_min, lon_max = np.nanmin(original_lon), np.nanmax(original_lon)
259. lat_min, lat_max = np.nanmin(original_lat), np.nanmax(original_lat)
261. # 生成目标数据的经纬度矩阵
262. grid_lon, grid_lat = np.meshgrid(
263. np.arange(lon_min, lon_max, out_res),
264. np.arange(lat_max, lat_min, -out_res),
265. )
267. # 插值
268. interp_dataset = griddata(
269. (original_lon.ravel(), original_lat.ravel()), # 原始数据的经纬度, ravel()将多维数组转换为一维数组
270. original_dataset.ravel(), # 原始数据
271. (grid_lon, grid_lat), # 目标数据的经纬度
272. method=method, # 插值方法
273. fill_value=np.nan, # 用于填充超出边界的点
274. )
276. return interp_dataset
279. # 输出TIF
280. def write_tiff(out_path, dataset, transform, nodata=np.nan):
281. """
282. 输出TIFF文件
283. :param out_path: 输出文件的路径
284. :param dataset: 待输出的数据
285. :param transform: 坐标转换信息(形式:[左上角经度, 经度分辨率, 旋转角度, 左上角纬度, 旋转角度, 纬度分辨率])
286. :param nodata: 无效值
287. :return: None
288. """
290. # 创建文件
291. driver = gdal.GetDriverByName('GTiff')
292. out_ds = driver.Create(out_path, dataset[0].shape[1], dataset[0].shape[0], len(dataset), gdal.GDT_Float32)
294. # 设置基本信息
295. out_ds.SetGeoTransform(transform)
296. out_ds.SetProjection('WGS84')
298. # 写入数据
299. for i in range(len(dataset)):
300. out_ds.GetRasterBand(i + 1).WriteArray(dataset[i]) # GetRasterBand()传入的索引从1开始, 而非0
301. out_ds.GetRasterBand(i + 1).SetNoDataValue(nodata)
302. out_ds.FlushCache()
305. # 准备
306. in_path = r'===>\REGC'
307. geo_path = r'===>\FY4A-_AGRI--_N_REGC_1047E_L1-_GEO-_MULT_NOM_20200606194500_20200606194917_4000M_V0001.HDF'
308. out_path = r'===>\output'
309. row_name = 'LineNumber' # 像元行号矩阵, 来自*GEO*.hdf(上geo_path)
310. col_name = 'ColumnNumber' # 像元列号矩阵, 同上
311. out_res = 0.036 # 0.036为输出的分辨率(°/度), 约为4km
313. # 检索HDF5文件
314. paths = glob.iglob(in_path + r'\*REGC*FDI*.HDF')
315. # 通过REGC区域的行列号矩阵计算经纬度矩阵
316. lon, lat = read_geo(geo_path, row_name, col_name)
317. for path in paths:
318. # 获取当前循环文件名
319. hdf_name = os.path.basename(path).split('.')[0]
320. # 读取HDF5文件
321. hdf = h5py.File(path, 'r')
322. nom_channels_name = [key for key in hdf.keys() if 'NOMChannel' in key]
323. cal_channels_name = [key for key in hdf.keys() if 'CALChannel' in key]
324. channels_name = zip(nom_channels_name, cal_channels_name)
325. # 箱子
326. bands = []
327. transform = None
328. for nom_channel_name, cal_channel_name in channels_name:
329. # 辐射定标
330. band = fy4a_calibration(path, nom_channel_name, cal_channel_name)
331. # 裁剪
332. clip_band, clip_lon, clip_lat = clip(band, lon, lat)
333. # GLT校正
334. processed_band = glt_warp(clip_band, clip_lon, clip_lat, out_res, method='linear')
335. # 仿射变换参数
336. transform = (
337. np.nanmin(clip_lon), # 左上角经度
338. out_res, # 东西向分辨率(经度分辨率)
339. 0, # 旋转角度
340. np.nanmax(clip_lat), # 左上角纬度
341. 0, # 旋转角度
342. -out_res # 南北向分辨率(纬度分辨率)
343. )
344. # 添加
345. bands.append(processed_band)
346. # 输出
347. write_tiff(os.path.join(out_path, hdf_name + '.tiff'), bands, transform)

3.3 基于行列号计算经纬度数据集(DISK/全圆盘)的FY4A L1产品的单个文件预处理

此处时间精力有限，并没有进行批量处理的改写，函数也是调用的，需要自行修改。 

1. # @炒茄子 2023-07-06
3. """
4. 当前脚本用于对FY4A L1(DISK全圆盘)产品进行辐射定标, 通过公式计算经纬度数据集(由行列数计算, 公式见官方文档)进行地理查找表的校正(未使用官方
5. 基于DISK的地理对照表, 但效果一致),最后输出为GeoTIFF文件(WGS84)
6. """
9. from lib.func import write_tiff, glt_warp, fy4a_calibration, get_lon_lat, clip
11. if __name__ == '__main__':
12. # 准备工作
13. in_path = r'===>\FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20200601010000_20200601011459_4000M_V0001.HDF'
14. out_path = r'===>\Python_DISK_no_raw.tif'
15. nom_channel_name = 'NOMChannel02'
16. cal_channel_name = 'CALChannel02'
17. # 辐射定标
18. band = fy4a_calibration(in_path, nom_channel_name, cal_channel_name)
19. # 依据行列号获取经纬度
20. lon, lat = get_lon_lat(band)
21. # 裁剪
22. reform_band, reform_lon, reform_lat = clip(band, lon, lat)
23. # 校正
24. interp_band = glt_warp(reform_band, reform_lon, reform_lat, 0.036) # 0.036°大约是4km
25. # 输出
26. write_tiff(out_path, [interp_band], reform_lon, reform_lat)

以下其它单个文件单个波段的处理，可能存在错误，因为函数进行了一定的修改。

1. # @炒茄子 2023-07-04
2. """
3. 当前脚本用于处理FY4A数据，包括辐射定标、基于地理查找表的校正、GeoTIFF文件的创建等, 最后输出为GeoTIFF文件(WGS84)
4. """
6. # 准备
7. in_path = ''
8. georaw_path = ''
9. out_path = ''
10. target_channel_name = 'NOMChannel01'
11. cal_channel_name = 'CALChannel01'
12. # 辐射定标
13. band = fy4a_calibration(in_path, target_channel_name, cal_channel_name)
14. # 读取地理查找表
15. glt_shape = band.shape # 地理查找表的shape与band的shape相同,或可从官方文档中获取((2748, 2748))
16. lon, lat = read_glt(georaw_path, glt_shape)
17. # 重投影和重采样区域
18. clip_band, clip_lon, clip_lat = clip(band, lon, lat)
19. # 基于地理查找表的校正(我认为这实际上可能仅仅是重投影和重采样, 因为我只做了这一些)
20. data = glt_warp(clip_band, clip_lon, clip_lat, 0.036) # 0.04为输出的分辨率(°/度)
21. # 将处理后的数据和地理查找表的信息写入GeoTIFF文件
22. write_tiff(out_path, [data], clip_lon, clip_lat)

1. # @炒茄子 2023-07-06
3. """
4. 当前脚本用于对FY4A L1(REGC中国区域)产品进行辐射定标, 基于FY4A的定位配准产品(*GEO*.HDF)进行地理查找表的校正, 最后输出为GeoTIFF文件(WGS84)
5. """
7. # 准备工作
8. in_path = r''
9. geo_path = r'' # 定位配准文件
10. out_path = r'' # 待输出文件路径
11. nom_channel_name = 'NOMChannel02' # 通道名, 来自*FDI*.hdf
12. cal_channel_name = 'CALChannel02' # 通道名, 来自*FDI*.hdf, 用于辐射定标
13. row_name = 'LineNumber' # 像元行号矩阵, 来自*GEO*.hdf
14. col_name = 'ColumnNumber' # 像元列号矩阵, 同上
16. # 辐射定标
17. band = fy4a_calibration(in_path, nom_channel_name, cal_channel_name)
18. # 依据行列号获取经纬度
19. lon, lat = read_geo(geo_path, row_name, col_name)
20. # 裁剪
21. reform_band, reform_lon, reform_lat = clip(band, lon, lat)
22. # 校正
23. interp_band = glt_warp(reform_band, reform_lon, reform_lat, 0.036) # 0.036°大约是4km
24. # 输出
25. write_tiff(out_path, [interp_band], reform_lon, reform_lat)

04 结果展示 

4.1 ArcGIS显示

 （没有明显的偏移）

4.2 ENVI 显示

---

 暂且如此(代码还有很多可以拓展的地方，但时间精力有限)，