C++使用HDF5庫實(shí)現(xiàn)將h5圖像轉(zhuǎn)為tif格式

更新時(shí)間：2024年12月09日 11:44:57 作者：瘋狂學(xué)習(xí)GIS

這篇文章主要為大家詳細(xì)介紹了C++ 語言如何通過hdf5庫與gdal庫實(shí)現(xiàn)將.h5格式的多波段HDF5圖像批量轉(zhuǎn)換為.tif格式,有需要的可以參考一下

本文介紹基于C++ 語言的hdf5庫與gdal庫，將.h5格式的多波段HDF5圖像批量轉(zhuǎn)換為.tif格式的方法；其中，本方法支持對szip壓縮的HDF5圖像（例如高分一號(hào)衛(wèi)星遙感影像）加以轉(zhuǎn)換。

將HDF5圖像批量轉(zhuǎn)換為.tif格式，在部分場景下操作并不難——在我們之前的文章Python中ArcPy實(shí)現(xiàn)柵格圖像文件由HDF格式批量轉(zhuǎn)換為TIFF格式中，就介紹過基于Python中的arcpy模塊實(shí)現(xiàn)這一需求的方法。但是，正如我們在文末補(bǔ)充內(nèi)容中提到的那樣，由于szip這個(gè)壓縮模塊不再受到hdf5庫的支持，導(dǎo)致用szip程序壓縮的HDF5圖像，比如高分系列遙感影像的.h5文件，就沒辦法在Windows中通過Python的h5py、gdal等庫直接打開了。

那么在這里，我們就介紹一下基于C++ 語言的hdf5庫，打開.h5格式圖像（包括那些用到szip壓縮程序的HDF5圖像）的方法。不過需要注意，我這里是在Linux的Ubuntu系統(tǒng)中操作的，至少可以保證這個(gè)代碼在Linux下可以正常運(yùn)行；但能否在Windows中的C++ 環(huán)境下也正常運(yùn)行，我暫時(shí)還沒試過——按道理應(yīng)該也是可行的，大家如果有需要的話可以嘗試一下。

本文所用代碼如下。

#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;
}

上述是本文完整代碼。接下來，就分段介紹一下每段代碼的具體含義。

首先，需要包含必要的頭文件。在這里，包括標(biāo)準(zhǔn)輸入輸出、字符串流、向量、文件系統(tǒng)等功能，以及hdf5庫與gdal庫。同時(shí)，定義了兩個(gè)常量字符串h5_path與tif_path，分別指向轉(zhuǎn)換前的HDF5圖像和轉(zhuǎn)換后的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/";

隨后，設(shè)定要讀取的HDF5圖像的數(shù)據(jù)集（波段）的路徑，以及空間參考信息的屬性名稱；這些參數(shù)大家就按照自己HDF5圖像的實(shí)際情況來修改即可。

接下來，初始化hdf5庫的狀態(tài)變量，這些變量是hdf5庫操作需要的。同時(shí)，用size表示圖像的寬度和高度，因?yàn)槲疫@里HDF5圖像是正方形，所以只需指定1個(gè)值。此外，band_num表示待轉(zhuǎn)換遙感影像的波段數(shù)。

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驅(qū)動(dòng)程序。

status = H5open();
GDALAllRegister();

隨后，使用std::filesystem::directory_iterator遍歷指定目錄中的所有文件，并只處理擴(kuò)展名為.h5的文件；對于這些文件，構(gòu)建輸出文件名——基于原始文件名，去掉擴(kuò)展名并添加.tif擴(kuò)展名。

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的屬性，讀取其內(nèi)容到attr_data中；隨后，解析attr_data為一個(gè)包含6個(gè)元素的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);

到這里，我們就可以對每個(gè)數(shù)據(jù)集調(diào)用H5Dopen1將其打開，并使用H5Dread將數(shù)據(jù)讀入向量中

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());

???????// ... 重復(fù)上述步驟讀取其他數(shù)據(jù)集 ...

隨后，記得關(guān)閉HDF5圖像以釋放資源。

status = H5Fclose(file_id);

接下來，就該gdal庫登場了。使用gdal庫創(chuàng)建一個(gè)新的TIFF文件，并使用RasterIO方法將每個(gè)波段的數(shù)據(jù)寫入到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);
// ... 寫入其他波段 ...

同時(shí)，設(shè)置每個(gè)波段的NoData值為0，同時(shí)按照前述從HDF5圖像中讀取到的信息，設(shè)置TIFF圖像的地理變換參數(shù)和投影信息。

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);

最后，不要忘記關(guān)閉hdf5庫以釋放資源。

status = H5close();

至此，大功告成。

知識(shí)補(bǔ)充

Windows打開HDF5圖像：HDFView軟件的下載、安裝

下面為大家介紹在Windows電腦中，下載、安裝用以查看HDF5圖像數(shù)據(jù)的軟件HDFView的方法。

HDF5（Hierarchical Data Format 5）是一種用于存儲(chǔ)和管理大量科學(xué)數(shù)據(jù)的文件格式，其由HDF Group開發(fā)和維護(hù)，廣泛應(yīng)用于科學(xué)計(jì)算、工程、金融和醫(yī)學(xué)等領(lǐng)域。談及HDF5圖像數(shù)據(jù)在Windows中的打開方式，主要包括基于HDF Group開發(fā)的HDFView軟件來打開，以及用C++、Python來打開等2種方式。

在之前，我很少選擇用HDFView軟件來打開HDF5，因?yàn)樵缧r(shí)候這個(gè)軟件的安裝比較麻煩，還需要修改一下環(huán)境變量什么的，不如在Python中配置對應(yīng)的庫（比如h5py、gdal等）然后用代碼讀取來的容易。但是，后來發(fā)現(xiàn)由于szip這個(gè)壓縮模塊不再受到hdf5等庫的支持（我看網(wǎng)上說好像是因?yàn)檫@個(gè)庫不再是非盈利的了還是怎么），導(dǎo)致那些用到szip壓縮的HDF5圖像（比如高分系列遙感影像數(shù)據(jù)的.h5文件），就沒辦法在Windows中通過Python的h5py、gdal等方便地打開了（Linux下C++ 的hdf5庫我試了，還是可以正常打開的，但是Windows中C++ 的hdf5庫是否能打開我還沒試過）。所以，在Windows中，如果只是需要打開、查看一下數(shù)據(jù)的話（不需要代碼執(zhí)行一些分析或批處理），通過HDFView軟件來打開HDF5還是很方便的。

這里就介紹一下HDFView軟件的下載、安裝方法。

首先，我們打開HDFView軟件的Github下載網(wǎng)站（https://github.com/HDFGroup/hdfview/releases）。當(dāng)然，也可以先進(jìn)入官方下載網(wǎng)站（https://portal.hdfgroup.org/downloads/index.html），找到其中的HDFView軟件下載位置，如下圖所示。

隨后，在彈出的窗口中，點(diǎn)擊Download下的here字樣，如下圖所示。

隨后，進(jìn)入的就是前面提到的Github網(wǎng)站。選擇適合自己的軟件版本，如下圖所示。

下載完畢后，將壓縮包放在一個(gè)自己指定的路徑中，并解壓壓縮包，雙擊打開其中的.exe文件，如下圖所示。

隨后，將彈出安裝窗口，如下圖所示。

其中，需要注意的就是配置一下軟件的安裝路徑，如下圖所示。

完成安裝后，可以在開始菜單中看到HDFView軟件的圖標(biāo)，如下圖所示。

雙擊圖標(biāo)，即可打開軟件，如下圖所示。

新版本的HDFView軟件也不需要再額外配置環(huán)境變量了，按照以上方法完成安裝后，直接打開就可以使用。

到此這篇關(guān)于C++使用HDF5庫實(shí)現(xiàn)將h5圖像轉(zhuǎn)為tif格式的文章就介紹到這了,更多相關(guān)C++ HDF5實(shí)現(xiàn)h5轉(zhuǎn)tif內(nèi)容請搜索腳本之家以前的文章或繼續(xù)瀏覽下面的相關(guān)文章希望大家以后多多支持腳本之家！

您可能感興趣的文章: