主體內容來自 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