歡迎關注博主的微信公眾號：“智慧遙感”。

該公眾號將為您奉上Python地學分析、爬蟲、資料分析、Web開發、機器學習、深度學習等熱門原始碼。

本人的GitHub程式碼資料主頁（持續更新中，多給Star，多Fork）：

https://github.com/xbr2017

CSDN也在同步更新：

https://blog.csdn.net/XBR_2014

“ 本節通過Python + GDAL讀取MODIS HDF格式的資料集。”

作為一名遙感仁，在接觸Python之前，一直使用IDL語言（遙感專業的一門語言）處理遙感資料，從幾何校正、輻射定標、影象拼接、向量裁剪柵格以及幾乎所有的遙感影象運算、資料讀取與常見格式的儲存，無一不是通過IDL來實現的。在讀研期間，甚至一度迷上用IDL製作讀取與展示遙感資料的小軟體介面，在畢業論文裡竟然專門有一小節介紹該軟體，現在想想，是多麼的蒼白無力

IDL是一門小眾語言，冷了不能再冷的語言。牆裂推薦：不管你是否搞科研，還是做開發，對於沒有計算機程式設計基礎（C、C++、Java等）的小夥伴們，包括軟妹子吆，Python是價效比最高的語言，沒有之一，簡單易懂，上手快。來吧，讓我們一起高喊：“人生苦短，我用Python”！

順便多說一句，既然我們使用了Python來處理遙感資料，那就準備拋棄ArcGis、MeteoInfo之類的二次開發環境，所以在今後的遙感資料處理介紹中，爭取不使用arcpy包。因為要使用該模組，自己的本本上需要安裝ArcGis軟體，所以還是儘量不使用arcpy。

GDAL讀取HDF

今天以較為經典的衛星資料MODIS為例，講解如何使用Python+GDAL讀取HDF檔案，以及輸出指定投影的tif檔案，這裡的tif是指Geotiff檔案，其中包含地理資訊。

HDF為分層資料格式的檔案，可以使用GetSubDatasets函式獲取它們的列表，然後使用該資訊開啟你想要的那個資料集。

例如，讓我們開啟MODIS檔案中包含的子資料集。請注意，預設情況下，HDAL驅動程式不包含在GDAL中，因此如果你的GDAL版本不包含HDF支援，此示例將不適用於你。假設你可以使用HDF檔案，第一步是將HDF檔案作為資料集開啟：

from osgeo import gdal

ds = gdal.Open('MYD13Q1.A2014313.h20v11.005.2014330092746.hdf')

現在，你可以獲取此開放資料集中包含的子資料集列表。GetSubDatasets方法返回元組列表，每個子資料集有一個元組。每個元組按順序包含子資料集的名稱和描述。以下程式碼段獲取此列表，然後打印出每個子資料集的名稱和描述：

subdatasets = ds.GetSubDatasets()
print 'Number of subdatasets: {}'.format(len(subdatasets))

for sd in subdatasets:
    print 'Name: {0}\nDescription:{1}\n'.format(*sd)

前幾行輸出看起來像這樣，並顯示第一個子資料集的資訊，即NDVI（歸一化差異植被指數），這裡只展示第一個，剩餘11個就不放在這顯示了：

Number of subdatasets: 12

Name: HDF4_EOS:EOS_GRID:"D:/osgeopy-data/Modis/MYD13Q1.A2014313.h20v11.005.2014330092746.hdf":MODIS_Grid_16DAY_250m_500m_VI:250m 16 days NDVI
Description:[4800x4800] 250m 16 days NDVI MODIS_Grid_16DAY_250m_500m_VI (16-bit integer)

要開啟子資料集，請將其名稱傳遞給gdalOpen。例如，這將獲取與第一個子資料集對應的元組，從元組中獲取第一個項（名稱），然後使用它來開啟子資料集：

ndvi_ds = gdal.Open(subdatasets[0][0])

同樣，使用subdatasets[4] [0]表示開啟第五個子資料集。一旦打開了這樣的子資料集，它就可以像任何其他資料集一樣處理。例如，你可以使用ndvi_ds.GetRasterBand(1)獲取NDVI子資料集中的第一個波段。

如果沒有特殊說明，HDF一般指hdf4，現在很多衛星資料開始使用hdf5格式儲存，對於GDAL都可以去讀取操作。當然，也有專門讀取hdf4與hdf5的包，感興趣的可以pip install安裝一下相關的包就可以直接使用啦。這裡給大家稍微展示一下具體程式碼使用：

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/28 22:47'

import h5py
import numpy as np

# HDF5的讀取：
in_ds = h5py.File('FY4A-_AGRI--_N_DISK_1047E_L1-_FDI-_MULT_NOM_20180101024500_20180101025959_4000M_V0001.HDF','r')   
# 讀取FY4A靜止衛星第一波段資料
array = in_ds['NOMChannel01'][:]  
in_ds.close()

GDAL將資料儲存為tif格式

遙感影象需要包含經緯度、空間解析度、投影座標等地理資訊，這是普通影象所不具備的，一般彩色影象只有RGB三個波段，而遙感影象可以包含很多波段。下面的程式是讀取MODIS歸一化植被指數（NDVI），並將資料儲存為帶有地理資訊的tif檔案，必要的解釋在程式碼中已給出。

# _*_ coding: utf-8 _*_
__author__ = 'xbr'
__date__ = '2018/11/26 17:48'

import gdal, osr


def array2raster(newRasterfn, rasterOrigin, xsize, ysize, array):
    """
     newRasterfn: 輸出tif路徑
     rasterOrigin: 原始柵格資料路徑
     xsize: x方向像元大小
     ysize: y方向像元大小
     array: 計算後的柵格資料
    """
    cols = array.shape[1]  # 矩陣列數
    rows = array.shape[0]  # 矩陣行數
    originX = rasterOrigin[0]  # 起始像元經度
    originY = rasterOrigin[1]  # 起始像元緯度
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(newRasterfn, cols, rows, 1, gdal.GDT_Float32)
    # 括號中兩個0表示起始像元的行列號從(0,0)開始
    outRaster.SetGeoTransform((originX, xsize, 0, originY, 0, ysize))
    # 獲取資料集第一個波段，是從1開始，不是從0開始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    # 程式碼4326表示WGS84座標
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()


ds = gdal.Open('D:/osgeopy-data/Modis/MYD13Q1.A2014313.h20v11.005.2014330092746.hdf')

subdatasets = ds.GetSubDatasets()
print('Number of subdatasets: {}'.format(len(subdatasets)))
for sd in subdatasets:
    print('Name: {0}\nDescription:{1}\n'.format(*sd))

ndvi_ds = gdal.Open(subdatasets[0][0]).ReadAsArray()
dst_filename = "D:/nc/try.tif"
xsize = 0.0025
ysize = 0.0025

array2raster(dst_filename, [90, 75], xsize, ysize, ndvi_ds)

NDVI第一波段的灰色影象展示

NDVI第一波段的彩色影象展示