# 将代码和数据文件(火灾和土地使用文件)都放在同一个文件夹
# 读取网格数据的面积
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')
Python读取hdf5文件
最新推荐文章于 2023-10-28 20:40:24 发布