使用Fortran读写风云卫星HDF5文件

Note: 本节介绍的安装方法可能不再适用，请参考最新的NETCDF/HDF5安装说明 — 2023-11-15

最近需要大批量读写一些FY卫星观测产品，决定使用fortran来读写，相对于NCL和Python会快不少，但是需要库依赖。
网上的指导很少，也不全面，自己参照着摸索了半天，也算可以满足使用要求了，决定记录一下。

系统：CentOS release 6.7 (Final)
权限：普通用户
gfortran 版本： gcc version 4.4.7 20120313 (Red Hat 4.4.7-18) (GCC)
HDF5版本: hdf5-1.8.22.tar.gz

HDF5依赖库

和依赖NetCDF库来读取NC文件一样，Fortran在读取HDF文件时需要HDF库的支持。
HDFx库为诸多语言比如C,C++,Pyhton,Java,Fortran,IDL等提供HDF格式文件操作扩展。
这里选择的是HDF5库，官网下载 https://portal.hdfgroup.org/display/support/Downloads.
我选择的版本是hdf5-1.8.22.tar.gz。

安装HDF5

1 从官网下载hdf5，将下载的tar.gz文件/home/hjh/hdf5/
2 解压并进入解压文件夹

bash

[hjh@node01] /home/hjh/hdf5$ tar -xvf hdf5-1.8.22.tar.gz 
[hjh@node01] /home/hjh/hdf5$ cd hdf5-1.8.22

3 编译安装

bash

[hjh@node01] /home/hjh/hdf5/hdf5-1.8.22$ ./configure --prefix=/home/hjh/hdf5 --enable-fortran

这里我们没有root权限，指定安装在自己的文件夹下/home/hjh/hdf5。
安装支持Fortran的依赖，--enable-fortran会编译生成h5fc编译命令并且会在lib文件夹下生成xxFortran.a等一些列Fortran扩展库。
默认情况下不会生成适用于Fortran的扩展，所以必须要要指定一下。

继续执行：

plaintext

$ make
$ make check   
$ make install

检查/home/hjh/hdf5内容
bin/ : h5cc(C语言编译器) h5dump (类似于ncdump) h5fc (fortran编译器) 等
lib/ ： libhdf5_fortran.so libhdf5hl_fortran.la libhdf5hl_fortran.so.10 等
include/ ：hdf5.mod 等

4 添加到环境变量
~/.bashrc 文件添加以下内容

bash

#HDF5
export HDF5=/home/hjh/hdf5
export PATH=$HDF5/bin:$PATH

保存退出，命令行运行 source ~/.bashrc

尝试一下h5dump,h5fc,h5cc,h5copy等命令。

编译、运行

bash

gfortran read.fengyun3b.hdf5.f90  -I/home/hjh/hdf5/include -L/home/hjh/hdf5/lib -lhdf5_fortran

默认生成a.out可执行文件,使用-o executable.exe来指定生成可执行文件的文件名。

或者h5fc编译Fortran脚本，h5cc编译c脚本。

bash

h5fc read.fengyun3b.hdf5.f90

后一种方式会同时生成read.fengyun3b.hdf5.o库文件和a.out可执行文件。
运行很简单：./a.out。

示例：打开风云3B MWRI数据并读取亮温

获取风云卫星的HDF结构

shell

h5dump -H  /data04/0/MWRI/L1_Recalibrated/FY3B/L1/ASCEND/2011/20110102/FY3B_MWRIA_GBAL_L1_20110102_0140_010KM_MS.HDF

HDF文件结构如下，关注一些DEM地形和EARTH_OBSERVE_BT_10_to_89GHz亮温两个变量，同位于calibration组：

xml

HDF5 "/data04/0/MWRI/L1_Recalibrated/FY3B/L1/ASCEND/2011/20110102/FY3B_MWRIA_GBAL_L1_20110102_0140_010KM_MS.HDF" {
   GROUP "Calibration" {
      DATASET "DEM" {
         DATATYPE  H5T_STD_I16LE
         DATASPACE  SIMPLE { ( 1826, 240 ) / ( 1826, 240 ) }
		 ....
	  }
	....
	   DATASET "EARTH_OBSERVE_BT_10_to_89GHz" {
         DATATYPE  H5T_STD_I16LE
         DATASPACE  SIMPLE { ( 10, 1826, 240 ) / ( 10, 1826, 240 ) }
		 ....
	  }
   }

由于不同轨道的行数是不同的，所以需要用可变数组来读取我们需要的亮温等，并且根据打开的文件来确定数组的大小。

代码：
read.fengyun3b.hdf5.f90

fortran

  USE HDF5 ! This module contains all necessary modules 
  
  CHARACTER filename*100
   
  ! CHARACTER(LEN=8), PARAMETER :: filename = "dsetf.h5" ! File name
  CHARACTER(LEN=3) dsetname      ! Dataset name 

  INTEGER(HID_T) :: file_id       ! File identifier  文件句柄
  INTEGER(HID_T) :: dset_id       ! Dataset identifier 变量句柄
  INTEGER(HID_T) :: grp_id        ! Dataset identifier  group句柄
  INTEGER(HID_T) :: dspace_id     ! Dataset identifier  只可意会，和维数和大小相关
  INTEGER     	 :: error       ! Error flag - success：0

  INTEGER(HID_T) :: dem_id  
  INTEGER, DIMENSION(1826,240) :: DEM ! Data buffers 对于事先指导变量维数的可以直接声明

  INTEGER, DIMENSION(:,:,:),allocatable :: tb_out ! Data buffers  本例中对于大小不定的数据使用可变数组

  INTEGER(HSIZE_T), DIMENSION(3) :: tb_dims,maxdims
  
  INTEGER(HSIZE_T), DIMENSION(2) :: dem_dims
  
   ! folder="/data04/0/MWRI/L1_Recalibrated/FY3B/L1/ASCEND/2011/20110102/"
  filename="FY3B_MWRIA_GBAL_L1_20110102_0140_010KM_MS.HDF"

  CALL h5open_f(error) ! Initialize FORTRAN interface.
 
  CALL h5fopen_f (filename, H5F_ACC_RDWR_F, file_id, error) ! Open an existing file.

  CALL h5gopen_f (file_id, "Calibration", grp_id, error)  ! Open an existing group.

  CALL h5dopen_f(grp_id, "EARTH_OBSERVE_BT_10_to_89GHz", dset_id, error)  !open dataset 就是变量
   
  ! Figure out the size of the EARTH_OBSERVE_BT_10_to_89GHz.
  CALL h5dget_space_f(dset_id,dspace_id,error)  ! Get the dataspace ID

  ! Get dims from dataspace
  CALL h5sget_simple_extent_dims_f(dspace_id, tb_dims, maxdims, error)   

  ! Allocate memory for the array.
  ALLOCATE(tb_out(tb_dims(1),tb_dims(2),tb_dims(3)))  
  
  ! Read the dataset.	 
  CALL h5dread_f(dset_id, H5T_NATIVE_INTEGER, tb_out, tb_dims, error)

  !! do something
  PRINT*,tb_out(1,1)*0.01+327.68  !! record(short型) * scalefactor+ offset
 
  CALL h5dclose_f(dset_id, error)   ! Close the dataset.
  DEALLOCATE(tb_out)  
  
!! 读取 DEM数据，这个是固定已知大小的读取
  CALL h5dopen_f(grp_id, "DEM", dem_id, error)  !open dem 
 
  dem_dims=(/1826,240/)
  CALL h5dread_f(dem_id, H5T_NATIVE_INTEGER, DEM, dem_dims, error)  ! Read the dem
  PRINT*,DEM(1,1)

  CALL h5dclose_f(dem_id, error)   ! Close the dataset.

!!!!!
  CALL h5gclose_f(grp_id, error) 	! Close the dataset.
  CALL h5fclose_f(file_id, error)  	! Close the file.
  CALL h5close_f(error)				! Close FORTRAN interface.

end

写出HDF5文件

主程序调用

f90

integer,parameter					:: nlat25=720+1, nlon25=1440+1
real,dimension(nlon25,nlat25,2) :: grid25_tb89s
character  dsetname*20, filein*80
integer retcode
....
dsetname = "TB_89VH_GRID25"
filein="TBs/grid25_tb89s."//yyyymmdd//".hdf5"
call write_hdf5(filein,nlon25, nlat25, 2, grid25_tb89s, dsetname,retcode) 

....

子程序如下：

f90

SUBROUTINE write_hdf5(filename,dim1, dim2,ch, dset_data, dsetname,retcode)

      USE HDF5   
      IMPLICIT NONE
	  Integer :: dim1, dim2,ch
      CHARACTER(*):: filename
      CHARACTER(*):: dsetname != "TB_89VH_GRID25" ! Dataset name   
      ! Identifiers
      INTEGER(HID_T) :: file_id       ! File identifier
      ! INTEGER(HID_T) :: group_id      ! Group identifier
      INTEGER(HID_T) :: dset_id      ! Dataset 1 identifier
      INTEGER(HID_T) :: dspace_id    ! Dataspace 1 identifier
	  
      ! FP array
      INTEGER :: rank ,retcode                ! Dataset rank
      INTEGER(HSIZE_T), DIMENSION(3) :: data_dims
      REAL, DIMENSION(dim1,dim2,ch) :: dset_data   ! Data buffers
 
      INTEGER :: error 

  ! =====================================================================
 
      ! Initialize the dset_data array 
      data_dims = (/dim1,dim2,ch/)
      rank = 3

    
      ! Initialize Fortran interface
      CALL h5open_f(error)   
      ! Create a new file
      CALL h5fcreate_f(trim(adjustl(filename)), H5F_ACC_TRUNC_F, file_id, error)
 
      CALL h5screate_simple_f(rank, data_dims, dspace_id, error)

      CALL h5dcreate_f(file_id, trim(adjustl(dsetname)), H5T_NATIVE_REAL, dspace_id,dset_id, error)

      CALL h5dwrite_f(dset_id, H5T_NATIVE_REAL, dset_data, data_dims, error)

      CALL h5dclose_f(dset_id, error)

      CALL h5sclose_f(dspace_id, error)
      
     ! Close the file
     CALL h5fclose_f(file_id, error)
     ! Close FORTRAN interface
     CALL h5close_f(error)
	 
END SUBROUTINE write_hdf5