NCL读取hdf5文件——以AMSR2海冰密集度数据为例
由于大部分计算都是在服务器进行。因此考虑增进NCL和Shell脚本水平，之前的博文Python读取HDF5介绍了HDF5文件及其读写，但是由于服务器本身系统版本较老，且属于内网原因，决定直接使用自带的NCL脚本对其进行撰写记录。
读取
NCL6版本以后便支持HDF5文件的读写，官方文档也给出了读写HDF5文件的实例：
;----------------------------------------------------------------------
; hdf5_1.ncl
;
; Concepts illustrated:
; - Reading group data off an HDF5 file using two methods
; - Creating a color map using RGB triplets
; - Drawing raster contours for faster results
;----------------------------------------------------------------------
; These files are loaded by default in NCL V6.2.0 and newer
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
; load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;----------------------------------------------------------------------begin
;---It is necessary to use advanced file structure to read group datasetfileoption("nc", "FileStructure", "Advanced")fname = "MSG3-SEVI-MSG15-0100-NA-20130521001244.164000000Z-1074164.h5"f = addfile(fname,"r");------------------------------------------------------------
; There are two ways to read variables off a files
; containing nested groups.
;
; Method 1:
; - Open the group with the "=>" syntax
; - Read the data using the usual "->" syntax
;
; Method 2:
; - Open the variable directly using the full path
; to the file
;
; Note: there's an inconsistency with how the group names
; are constructed for the two methods. We have opened a
; trouble ticket for this.
;------------------------------------------------------------USE_FULL_PATH = Trueif(.not.USE_FULL_PATH) then
;---First method; open group, then read data from it.group_name ="/U_MARF/MSG/Level1_5/DATA/Channel_07"iname = "IMAGE_DATA"pname = "Palette"gp = f=>$group_name$fdata = tofloat(gp->$iname$)palette = tofloat(gp->$pname$)else
;---Second method; use full path to variable.group_name = "/U-MARF/MSG/Level1.5/DATA/Channel 07"iname = group_name + "/IMAGE_DATA"pname = group_name + "/Palette"fdata = tofloat(f->$iname$)palette = tofloat(f->$pname$)end ifwks = gsn_open_wks("png","hdf5") ; send graphics to PNG fileres = True ; plot mods desiredres@gsnMaximize = True ; maximize plot in frameres@cnFillOn = True ; color fill res@cnFillMode = "RasterFill" ; Raster mode is much faster; and uses less memory.res@cnFillPalette = palette/255.res@cnLinesOn = False ; turn off contour linesres@cnLevelSelectionMode = "ManualLevels" ; set manual contour levelsres@cnMinLevelValF = min(fdata) ; set min contour levelres@cnMaxLevelValF = max(fdata) ; set max contour levelres@cnLevelSpacingF = 10.0 ; set contour spacingres@lbOrientation = "vertical" ; vertical labelbarres@lbBoxLinesOn = Falseres@gsnTickMarksOn = Falseif(USE_FULL_PATH)res@tiMainString = "Variable read using full group path"elseres@tiMainString = "Variable read by opening group first"end ifplot = gsn_csm_contour(wks,fdata, res) ; create plot
end
主要有两种读法，但不论哪一种，都需要对HDF5文件本身的结构有所了解，在上一篇博文中。我主要使用了matlab的h5info函数来直接查看，而NCL也提供了对应的函数：ncl_filedump，
其用法也非常简单：
ncl_filedump K1VHR_15NOV2013_1200_L02_OLR.h5 | less
或者直接：
setfileoption("h5", "FileStructure", "Advanced")
f = addfile("ncl_wrt_uvt.h5", "r")
print(f)
也可查看
值得注意的是：
matlab h5info返回的信息与NCL返回的信息并不一致，这导致很多时候我们并不能像matlab或者python、NCL官方示例一样，从不同的GROUP中读取数据，NCL对于HDF5数据的处理很多时候和nc一样，直接给出变量名。
举个实例，同样是AMSR2文件，matlab给出的结构为：
Group '/HDFEOS' Group '/HDFEOS/ADDITIONAL' Group '/HDFEOS/ADDITIONAL/FILE_ATTRIBUTES' Group '/HDFEOS/GRIDS' Group '/HDFEOS/GRIDS/NpPolarGrid12km' Dataset 'XDim' Size: 608MaxSize: 608Datatype: H5T_IEEE_F32LE (single)ChunkSize: []Filters: noneFillValue: 0.000000Attributes:'CLASS': 'DIMENSION_SCALE''NAME': 'This is a netCDF dimension but not a netCDF variable.''axis': 'X''long_name': 'x-coordinate in Cartesian system''REFERENCE_LIST': H5T_COMPOUNDDataset 'YDim' Size: 896MaxSize: 896Datatype: H5T_IEEE_F32LE (single)ChunkSize: []Filters: noneFillValue: 0.000000Attributes:'CLASS': 'DIMENSION_SCALE''NAME': 'This is a netCDF dimension but not a netCDF variable.''axis': 'Y''long_name': 'y-coordinate in Cartesian system''REFERENCE_LIST': H5T_COMPOUNDDataset 'lat' Size: 608x896MaxSize: 608x896Datatype: H5T_IEEE_F32LE (single)ChunkSize: []Filters: noneFillValue: 0.000000Attributes:'DIMENSION_LIST': H5T_VLEN'long_name': 'Latitude''units': 'degrees_north''standard_name': 'latitude'Dataset 'lon' Size: 608x896MaxSize: 608x896Datatype: H5T_IEEE_F32LE (single)ChunkSize: []Filters: noneFillValue: 0.000000Attributes:'DIMENSION_LIST': H5T_VLEN'long_name': 'Longitude''standard_name': 'longitude''units': 'degrees_east'Group '/HDFEOS/GRIDS/NpPolarGrid12km/Data Fields' Dataset 'SI_12km_NH_18H_ASC' Size: 608x896MaxSize: 608x896Datatype: H5T_STD_I32LE (int32)ChunkSize: []Filters: noneFillValue: 0Attributes:'DIMENSION_LIST': H5T_VLEN'coordinates': 'lon lat''long_name': '18.7 GHz horizontal daily average ascending Tbs''units': 'degree_kelvin''standard_name': 'brightness_temperature''packing_convention': 'netCDF''packing_convention_description': 'unpacked = scale_factor x packed + add_offset''scale_factor': 0.100000 'add_offset': 0.000000 '_FillValue': 0
算是标准的HDF5结构了，而NCL print出来则是：
dimensions:YDim_NpPolarGrid12km = 896XDim_NpPolarGrid12km = 608YDim_SpPolarGrid12km = 664XDim_SpPolarGrid12km = 632variables:double GridLat_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )projection : Polar Stereographiccorners : ( -39.29786148540076, -39.29786148540076, -41.51518496442763, -41.51518496442763 )long_name : latitudeunits : degrees_northdouble GridLon_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )projection : Polar Stereographiccorners : ( -42.23673723653989, 42.23673723653989, 135, -135 )long_name : longitudeunits : degrees_eastinteger SI_12km_SH_SNOWDEPTH_5DAY_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : 110 -- missing/not calculated, 120 -- Land, 130 -- Open water, 140 -- multiyear sea ice, 150 -- variability in snow depth, 160 -- snow meltunits : cmlong_name : 5-day snow depthcoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : 5-day snow depthprojection : Polar Stereographicinteger SI_12km_SH_ICEDIFF_DAY_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 1 to 100 or -1 to -100 -- Percent difference between algorithms, 120 -- Land, 200 to 300 -- Missing NT2 value (200+BBA), -200 to -300 -- Missing BBA value (-200 - NT2)units : percentlong_name : Sea ice concentration daily average using the difference between BBA and NT2 algorithmscoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily average using the difference between BBA and NT2 algorithmsprojection : Polar Stereographicinteger SI_12km_SH_ICEDIFF_DSC_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 1 to 100 or -1 to -100 -- Percent difference between algorithms, 120 -- Land, 200 to 300 -- Missing NT2 value (200+BBA), -200 to -300 -- Missing BBA value (-200 - NT2)units : percentlong_name : Sea ice concentration daily descending average using the difference between BBA and NT2 algorithmscoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily descending average using the difference between BBA and NT2 algorithmsprojection : Polar Stereographicinteger SI_12km_SH_ICEDIFF_ASC_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 1 to 100 or -1 to -100 -- Percent difference between algorithms, 120 -- Land, 200 to 300 -- Missing NT2 value (200+BBA), -200 to -300 -- Missing BBA value (-200 - NT2)units : percentlong_name : Sea ice concentration daily ascending average using the difference between BBA and NT2 algorithmscoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily ascending average using the difference between BBA and NT2 algorithmsprojection : Polar Stereographicinteger SI_12km_SH_ICECON_DAY_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 110 -- missing/not calculated, 120 -- Landunits : percentlong_name : Sea ice concentration daily averagecoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily averageprojection : Polar Stereographicinteger SI_12km_SH_ICECON_DSC_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 110 -- missing/not calculated, 120 -- Landunits : percentlong_name : Sea ice concentration daily descending averagecoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily descending averageprojection : Polar Stereographicinteger SI_12km_SH_ICECON_ASC_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmcomment : data value meaning: 0 -- Open Water, 110 -- missing/not calculated, 120 -- Landunits : percentlong_name : Sea ice concentration daily ascending averagecoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : Sea ice concentration daily ascending averageprojection : Polar Stereographicinteger SI_12km_SH_89H_DAY_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12km_FillValue : 0add_offset : 0scale_factor : 0.1packing_convention_description : unpacked = scale_factor x packed + add_offsetpacking_convention : netCDFstandard_name : brightness_temperatureunits : degree_kelvinlong_name : 89.0 GHz horizontal daily average Tbscoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : 89.0 GHz horizontal daily average Tbsprojection : Polar Stereographicinteger SI_12km_SH_89H_DSC_SpPolarGrid12km ( YDim_SpPolarGrid12km, XDim_SpPolarGrid12km )coordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12km_FillValue : 0add_offset : 0scale_factor : 0.1packing_convention_description : unpacked = scale_factor x packed + add_offsetpacking_convention : netCDFstandard_name : brightness_temperatureunits : degree_kelvinlong_name : 89.0 GHz horizontal daily average descending Tbscoordinates : GridLat_SpPolarGrid12km, GridLon_SpPolarGrid12kmlong_name : 89.0 GHz horizontal daily average descending Tbsprojection : Polar Stereographic
没有了Grouo信息，而是直接将变量提取出来。
因此，在使用NCL读取HDF5时，应当注意其NCL对其的处理方式。
实例AMSR2海冰数据
接下来介绍实例，用ncl读取AMSR2 12.5Km的海冰密集度数据：
; - Reading am AMSR2 h5 file containing "sea ice"
; - Creating the appropriate grid and time coordinates
; - Chang format to WPS_int fileload "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/esmf/ESMF_regridding.ncl"beginsetfileoption("h5", "FileStructure", "Advanced")ymdd1="20191201"ymdd2="2019-12-01"
; ymdd1=getenv("basedate1")
; ymdd2=getenv("basedate2")data_filename="/public/home/zhangzilu/AMSR/2019/AMSR_U2_L3_SeaIce12km_B04_"+ymdd1+".he5"print(data_filename)f = addfile(data_filename, "r")lat=f->GridLat_NpPolarGrid12kmlon=f->GridLon_NpPolarGrid12kmice_con=tofloat(f->SI_12km_NH_ICECON_DAY_NpPolarGrid12km)printVarSummary(ice_con)
; print(lat)
end便可，批量读取只要与shell脚本结合便可高效批量完成。
读取HDF5数据唯一的难点便在于对于数据结构的理解。
