命令行记录-python读取hdf图层，转成tif文件

主体内容来自 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

转载于:https://www.cnblogs.com/vividautumn/p/11548225.html