Python空间数据处理3： HDF格式转GDAL支持格式

HDF格式
一种能有条不紊 、完备地保存遥感影像的属性和空间信息数据 , 同时使查询和提取相关数据也很方便容易的数据格式格式。HDF格式是一种较为常见的遥感数据格式（例如MODIS、3B43降水数据），但GDAL暂时未能给HDF格式提供支持，需要PYHDF与GDAL交互，从而将HDF格式转换为GDAL支持的格式。

PyHDF
python处理HDF数据的开源包，笔者安装的pyhdf4(电脑配置WIN64+python2.7)，在我的另一篇文章《Python开源包安装：PYHDF安装》中，有详细介绍PyHDF的安装方法。

下面笔者就以 MODIS 1B 数据为例，根据PyHDF和GDAL，介绍转HDF格式影像为包含红、红外两个波段的GDAL所支持的栅格格式影像。

1、pyhdf 获取待处理影像栅格矩阵
导入需要用到的模块：

import osfrom pyhdf.SD import SDimport numpy as np

将工作空间转到待处理影像所在文件夹：

os.chdir(r'D:\modis_ndvi_classificate\fix_proj\pre_ndvi')

SD打开待处理影像

FILE_NAME = 'MOD09Q1.A2015169.h26v04.006.2015301111712.hdf'hdf = SD(FILE_NAME, SDC.READ)

查看影像内所包含的数据：

print hdf.datasets()

输入“print hdf.datasets()”后，显示：

其中，sur_refl_b02, sur_refl_b01 分别代表MODIS的红外、红两个波段。

获取红、红外两个波段的影像栅格，data1与data2：

data1 = hdf.select('sur_refl_b01')data1 = data1[:]data2 = hdf.select('sur_refl_b02')data2 = data2[:]

将data1与data2叠加：

data = np.concatenate ([data1[np.newaxis,:], data2[np.newaxis,:]], axis=0)

2、获取仿射矩阵
查看影像头文件信息：

attr = hdf.attributes(full=1)attNames = attr.keys()attNames.sort()print attNames

输入”print attNames” 后会显示该影像所包含头文件的目录有[‘ArchiveMetadata.0’, ‘CoreMetadata.0’,’HDFEOSVersion’,’StructMetadata.0’,’identifier_product_doi’,’identifier_product_doi_authority’]，我们打开’ArchiveMetadata.0’，中间包含了影像的起始点坐标和分辨率信息，：

t = attr['ArchiveMetadata.0']print t[0]

输入‘print t[0]’后，我们获得该影像起始点坐标：经度104.43258313171，纬度49.9999999955098；红，红外波段分辨率为231.656358263889，则该影像仿射矩阵为：

geot = (104.43258313171,231.656358263889,0,49.9999999955098,0,-231.656358263889)

3、给影像添加用户所指定的投影：
（以WGS84为例）

from osgeo import osrproj= osr.SpatialReference()proj.ImportFromEPSG(4326)proj.ExportToWkt()

4、重写影像
有了影像的栅格矩阵、仿射矩阵，投影信息后，我们就可以按照《Python空间数据处理1： GDAL读写遥感图像》部分的内容，重新将影像写成GDAL所支持的栅格格式文件了。