C++的HDF5库将h5图像转为tif格式：szip压缩的图像也可转换

  本文介绍基于C++ 语言的hdf5库与gdal库，将.h5格式的多波段HDF5图像批量转换为.tif格式的方法；其中，本方法支持对szip压缩的HDF5图像（例如高分一号卫星遥感影像）加以转换。

  将HDF5图像批量转换为.tif格式，在部分场景下操作并不难——在我们之前的文章Python中ArcPy实现栅格图像文件由HDF格式批量转换为TIFF格式（https://blog.csdn.net/zhebushibiaoshifu/article/details/124200816）中，就介绍过基于Python中的arcpy模块实现这一需求的方法。但是，正如我们在文章Windows打开HDF5图像：HDFView软件的下载、安装（https://blog.csdn.net/zhebushibiaoshifu/article/details/142642691）中提到的那样，由于szip这个压缩模块不再受到hdf5库的支持，导致用szip程序压缩的HDF5图像，比如高分系列遥感影像的.h5文件，就没办法在Windows中通过Python的h5py、gdal等库直接打开了。

  那么在这里，我们就介绍一下基于C++ 语言的hdf5库，打开.h5格式图像（包括那些用到szip压缩程序的HDF5图像）的方法。不过需要注意，我这里是在Linux的Ubuntu系统中操作的，至少可以保证这个代码在Linux下可以正常运行；但能否在Windows中的C++ 环境下也正常运行，我暂时还没试过——按道理应该也是可行的，大家如果有需要的话可以尝试一下。

  本文所用代码如下。

#include <iostream>
#include <sstream>
#include <vector>
#include <filesystem>
#include <gdal.h>
#include <gdal_priv.h>
#include "hdf5.h"
#include "ogr_spatialref.h" int main(int argc, char *argv[]) {const std::string h5_path = "/home/ctj/data/H5/";const std::string tif_path = "/home/ctj/data/TIFF_48SUB/";// const std::string h5_path = argv[1];// const std::string tif_path = argv[2];const char *dataset_0 = "/Cloud_Mask/cloudmask";const char *dataset_1 = "/GeometricCorrection/DataSet_16_1";const char *dataset_2 = "/GeometricCorrection/DataSet_16_2";const char *dataset_3 = "/GeometricCorrection/DataSet_16_3";const char *dataset_4 = "/GeometricCorrection/DataSet_16_4";const char *projection_para = "ProjectionPara";const char *projection_str = "ProjectionStr";hid_t file_id;hid_t dataset_id;hid_t attr_id;hid_t attr_dtype;herr_t status;hid_t mem_type_id = H5T_NATIVE_UINT16;int size = 6863;int band_num = 5;// namespace fs = filesystem;status = H5open();GDALAllRegister();for (const auto& entry : std::filesystem::directory_iterator(h5_path)) {if (entry.path().extension() == ".h5") {std::string filename = entry.path().filename().string();std::cout << filename << std::endl;std::string baseName = filename.substr(0, filename.find_last_of('.'));const std::string output_filename = tif_path + baseName + ".tif";file_id = H5Fopen((h5_path + filename).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);attr_id = H5Aopen(file_id, projection_para, H5P_DEFAULT);attr_dtype = H5Aget_type(attr_id);size_t string_length = H5Tget_size(attr_dtype);char *attr_data = new char[1000];status = H5Aread(attr_id, attr_dtype, attr_data);std::istringstream iss(attr_data);std::vector<double> transform(6);int i = 0;std::string str;while (getline(iss, str, ',')) {transform[i] = stod(str);++i;}attr_id = H5Aopen(file_id, projection_str, H5P_DEFAULT);attr_dtype = H5Aget_type(attr_id);char *attr_data_str = new char[1000];status = H5Aread(attr_id, attr_dtype, attr_data_str);dataset_id = H5Dopen1(file_id, dataset_0);std::vector<u_int16_t> data_0(size * size);status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_0.data());dataset_id = H5Dopen1(file_id, dataset_1);std::vector<u_int16_t> data_1(size * size);status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_1.data());dataset_id = H5Dopen1(file_id, dataset_2);std::vector<u_int16_t> data_2(size * size);status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_2.data());dataset_id = H5Dopen1(file_id, dataset_3);std::vector<u_int16_t> data_3(size * size);status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_3.data());dataset_id = H5Dopen1(file_id, dataset_4);std::vector<u_int16_t> data_4(size * size);status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_4.data());status = H5Fclose(file_id);GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");GDALDataset *poDstDS = poDriver->Create(output_filename.c_str(), size, size, band_num, GDT_UInt16, nullptr);u_int16_t *band_data_0 = &data_0[0];poDstDS->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, size, size, band_data_0, size, size, GDT_UInt16, 0, 0);u_int16_t *band_data_1 = &data_1[0];poDstDS->GetRasterBand(2)->RasterIO(GF_Write, 0, 0, size, size, band_data_1, size, size, GDT_UInt16, 0, 0);u_int16_t *band_data_2 = &data_2[0];poDstDS->GetRasterBand(3)->RasterIO(GF_Write, 0, 0, size, size, band_data_2, size, size, GDT_UInt16, 0, 0);u_int16_t *band_data_3 = &data_3[0];poDstDS->GetRasterBand(4)->RasterIO(GF_Write, 0, 0, size, size, band_data_3, size, size, GDT_UInt16, 0, 0);u_int16_t *band_data_4 = &data_4[0];poDstDS->GetRasterBand(5)->RasterIO(GF_Write, 0, 0, size, size, band_data_4, size, size, GDT_UInt16, 0, 0);for (int i = 1; i <= band_num; ++i) {GDALRasterBand *poBand = poDstDS->GetRasterBand(i);if (poBand != nullptr) {poBand->SetNoDataValue(0);}}poDstDS->SetGeoTransform(transform.data());poDstDS->SetProjection(attr_data_str);GDALClose(poDstDS);}}status = H5close();return 0;
}

  上述是本文完整代码。接下来，就分段介绍一下每段代码的具体含义。

  首先，需要包含必要的头文件。在这里，包括标准输入输出、字符串流、向量、文件系统等功能，以及hdf5库与gdal库。同时，定义了两个常量字符串h5_path与tif_path，分别指向转换前的HDF5图像和转换后的TIFF图像的目录。

#include <iostream>
#include <sstream>
#include <vector>
#include <filesystem>
#include <gdal.h>
#include <gdal_priv.h>
#include "hdf5.h"
#include "ogr_spatialref.h" int main(int argc, char *argv[]) {const std::string h5_path = "/home/ctj/data/H5/";const std::string tif_path = "/home/ctj/data/TIFF_48SUB/";

  随后，设定要读取的HDF5图像的数据集（波段）的路径，以及空间参考信息的属性名称；这些参数大家就按照自己HDF5图像的实际情况来修改即可。

  接下来，初始化hdf5库的状态变量，这些变量是hdf5库操作需要的。同时，用size表示图像的宽度和高度，因为我这里HDF5图像是正方形，所以只需指定1个值。此外，band_num表示待转换遥感影像的波段数。

const char *dataset_0 = "/Cloud_Mask/cloudmask";
const char *dataset_1 = "/GeometricCorrection/DataSet_16_1";
// ... 省略部分代码 ...
const char *projection_para = "ProjectionPara";
const char *projection_str = "ProjectionStr";hid_t file_id;
hid_t dataset_id;
hid_t attr_id;
hid_t attr_dtype;
herr_t status;
hid_t mem_type_id = H5T_NATIVE_UINT16;
int size = 6863;
int band_num = 5;

  紧接着，初始化hdf5库，注册所有可用的GDAL驱动程序。

status = H5open();
GDALAllRegister();

  随后，使用std::filesystem::directory_iterator遍历指定目录中的所有文件，并只处理扩展名为.h5的文件；对于这些文件，构建输出文件名——基于原始文件名，去掉扩展名并添加.tif扩展名。

for (const auto& entry : std::filesystem::directory_iterator(h5_path)) {if (entry.path().extension() == ".h5") {std::string filename = entry.path().filename().string();std::cout << filename << std::endl;std::string baseName = filename.substr(0, filename.find_last_of('.'));const std::string output_filename = tif_path + baseName + ".tif";

  随后，使用H5Fopen打开HDF5图像，在这里选择以只读模式访问。

file_id = H5Fopen((h5_path + filename).c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

  随后，需要读取原本HDF5图像的空间参考信息。在这里，首先打开名为projection_para的属性，读取其内容到attr_data中；随后，解析attr_data为一个包含6个元素的double向量transform——这些元素用于地理变换。

attr_id = H5Aopen(file_id, projection_para, H5P_DEFAULT);
attr_dtype = H5Aget_type(attr_id);
size_t string_length = H5Tget_size(attr_dtype);
char *attr_data = new char[1000];
status = H5Aread(attr_id, attr_dtype, attr_data);std::istringstream iss(attr_data);
std::vector<double> transform(6);
int i = 0;
std::string str;
while (getline(iss, str, ',')) {transform[i] = stod(str);++i;
}

  类似地，读取名为projection_str的属性，该属性包含投影信息的WKT字符串。

attr_id = H5Aopen(file_id, projection_str, H5P_DEFAULT);
attr_dtype = H5Aget_type(attr_id);
char *attr_data_str = new char[1000];
status = H5Aread(attr_id, attr_dtype, attr_data_str);

  到这里，我们就可以对每个数据集调用H5Dopen1将其打开，并使用H5Dread将数据读入向量中

dataset_id = H5Dopen1(file_id, dataset_0);
std::vector<u_int16_t> data_0(size * size);
status = H5Dread(dataset_id, mem_type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT, data_0.data());// ... 重复上述步骤读取其他数据集 ...

  随后，记得关闭HDF5图像以释放资源。

status = H5Fclose(file_id);

  接下来，就该gdal库登场了。使用gdal库创建一个新的TIFF文件，并使用RasterIO方法将每个波段的数据写入到TIFF文件中。

GDALDriver *poDriver = GetGDALDriverManager()->GetDriverByName("GTiff");GDALDataset *poDstDS = poDriver->Create(output_filename.c_str(), size, size, band_num, GDT_UInt16, nullptr);
u_int16_t *band_data_0 = &data_0[0];
poDstDS->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, size, size, band_data_0, size, size, GDT_UInt16, 0, 0);
// ... 写入其他波段 ...

  同时，设置每个波段的NoData值为0，同时按照前述从HDF5图像中读取到的信息，设置TIFF图像的地理变换参数和投影信息。

for (int i = 1; i <= band_num; ++i) {GDALRasterBand *poBand = poDstDS->GetRasterBand(i);if (poBand != nullptr) {poBand->SetNoDataValue(0);}
}poDstDS->SetGeoTransform(transform.data());
poDstDS->SetProjection(attr_data_str);GDALClose(poDstDS);

  最后，不要忘记关闭hdf5库以释放资源。

status = H5close();

  至此，大功告成。

欢迎关注：疯狂学习GIS