[Python知识库]Python读取hdf5文件

# 将代码和数据文件（火灾和土地使用文件）都放在同一个文件夹
# 读取网格数据的面积
import numpy as np
import h5py
from scipy import interpolate
import netCDF4 as nc
import pandas as pd
import pylab as pl
import matplotlib as mpl


# 步骤1 ------------------------ 提高两个矩阵的分辨率
# 读火灾数据hdf5文件
string='GFED4.1s_2015.hdf5'
f = h5py.File(string, 'r')

basis_regions = f['/ancill/basis_regions'][:]  #各个区域的标识矩阵
grid_area     = f['/ancill/grid_cell_area'][:] #每个网格面积的矩阵

# 读一下火灾原始经纬度720*1440
lat=f['lat'][:]
lon=f['lon'][:]

# 重新定义更高分辨率的经纬度
lon_new=np.linspace(-180+0.25/3/2,180-0.25/3/2,1440*3)
lat_new=np.linspace(90-0.25/3/2,-90+0.25/3/2,720*3)

# 使用插值函数把分辨率提高
grid_area_newfunc = interpolate.interp2d(lon[0,:], lat[:,0], grid_area , kind='linear')
grid_area_new = grid_area_newfunc(lon_new ,lat_new) #应该是正确的完成了
grid_area_new

# 现在对那个大洲的标识矩阵也进行相同的操作
basis_regions_newfunc = interpolate.interp2d(lon[0,:], lat[:,0], basis_regions , kind='linear')
basis_regions_new = basis_regions_newfunc(lon_new ,lat_new) 
basis_regions_new = np.fix(basis_regions_new) # 需要注意现在是上下颠倒的，需要调回来
basis_regions_new =np.flipud(basis_regions_new)

# 步骤2------------------------------------计算某一年的土地使用面积
# 对于某一年为例
# 某个区域的耕地面积计算=sum(grid_fraction*grid_area_new/9.0*mask)
# 需要注意这里的grid_fraction是上下颠倒的

# 读.nc文件
file = 'KK10.nc'
dataset = nc.Dataset(file)

time=7901   #时间点的维度
mask_values=14  #1-14，区域的数目

#新建一个数据框来储存最后的计算结果     
data_frame=np.ones((time,mask_values))
data_frame=pd.DataFrame(data_frame)

# 循环计算
for t in range(time): # range(2,time)     0,1
    print(t+1)
    land_use_frac = dataset.variables['land_use'][t,:,:]
    land_use_frac =np.flipud(land_use_frac)
    
    for mask_val in [1]:    # [1,5]
        mask=(basis_regions_new==mask_val)
        print(np.sum(land_use_frac*grid_area_new/9.0*mask))
        data_frame.iloc[t,mask_val-1]=np.sum(land_use_frac*grid_area_new/9.0*mask)

# 将数据输出为csv文件
data_frame.to_csv('land_use.csv')