主体内容来自 https://www.cnblogs.com/ninicwang/p/11535170.html
1、安装pyhdf包
2、读hdf4文件
#导入包
from pyhdf.SD import *
from osgeo import osr
import numpy as np
(1) #读取文件
file='3B43.20100501.7.HDF'
hdf=SD(file)
#可查看hdf的函数
>>> dir(hdf)
['__class__', '__del__', '__delattr__', '__dict__', '__dir__
', '__doc__', '__eq__', '__format__', '__ge__', '__getattr__
', '__getattribute__', '__gt__', '__hash__', '__init__', '__
init_subclass__', '__le__', '__lt__', '__module__', '__ne__'
, '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__s
etattr__', '__sizeof__', '__str__', '__subclasshook__', '__w
eakref__', '_id', 'attr', 'attributes', 'create', 'datasets'
, 'end', 'info', 'nametoindex', 'reftoindex', 'select', 'set
fillmode']
#获取hdf文件的属性信息
attr=hdf.attributes()
##获取hdf文件的图层以及每个图层所对应的行与列等信息,datasets是一个字典
datasets=hdf.datasets()
attr.keys()
#输出为 dict_keys(['FileHeader', 'FileInfo', 'GridHeader'])
#其中只有GridHeader中的信息较为可用
>>> print(attr['GridHeader'])
BinMethod=ARITHMETIC_MEAN;
Registration=CENTER;
LatitudeResolution=0.25;
LongitudeResolution=0.25;
NorthBoundingCoordinate=50;
SouthBoundingCoordinate=-50;
EastBoundingCoordinate=180;
WestBoundingCoordinate=-180;
Origin=SOUTHWEST;
#读某一图层(我这里的图层是降雨图层),图层名可通过datasets看到。并获取图层数据
rainfall_array=hdf.select("precipitation").get()
#图层的旋转与转置(这里我要做旋转转置是由于我的图层需求,可通过查看自己的图层是否需要做这一步),依据Origin=SOUTHWEST以及data.shape行列号(行1440,列400)获取原始栅格数据的方向
#对矩阵进行左右翻转
data=np.fliplr(rainfall_array)
#对矩阵进行转置
data=np.transpose(data)
#从图层中读取到的数据是5月份逐小时数据,我们这里转换成逐月数据
data=data*24*31
#将数据转换成两个字节的整型数据,注意是uint而非unit
d=data.astype(np.uint16)
(2)转tif文件
#创建一个空的投影对象,实例化
ref=osr.SpatialReference()
#定义投影(导入投影参数)
#设置坐标系统为WGS84
ref.ImportFromEPSG(4326)
#查看参考系信息
s=ref.ExportToWkt()
#写tif文件
from osgeo import gdal
driver=gdal.GetDriverByName("GTiff")
#通过data.shape判读数组im_bands,im_height,im_width波段数,行列号
#通过data.dtype.name获取栅格数据类型,gdal.GDT_Uint16而非uint16同一个东西不同的说法
dataset=driver.Create("prec_new.tif",1440,400,1,gdal.GDT_UInt16)
#仿射变换参数是自己编的
im_geotrans=(-180,0.25,0.0,50,0.0,-0.25)
dataset.SetGeoTransform(im_geotrans)
dataset.SetProjection(s)
dataset.GetRasterBand(1).WriteArray(d)
del dataset