【亲测】Python-GDAL 把MODIS HDF影像写入并导出为GeoTiff

写在前边，昨日鏖战许久未能在Windows平台上搞定基于R语言（主要是rgdal 与gdalUtils）的HDF文件的预处理，主要报错为：

Error in gdal_chooseInstallation(hasDrivers = of) : 
  No installations match.

然后，夜以继日地在网上找解决方案，发现只有问题，没有方案。。。内心崩溃之余，顺手换了平台与语言，以下是比较成功的完整的解决流程。

目标问题

将MODIS HDF 文件转化为R等数据分析语言可以识别并批量处理的GeoTiff格式。

问题分析

1. Windows + R已无法解决且网上没有比较成熟的Solution
2. Mac及GDAL对Python的兼容性具高

具体解决流程

1. Mac 安装 GDAL库

brew install gdal

电脑上没有Homebrew的可以看下文:
Mac 上安装Homebrew 教程

2. Python gdal模块

conda install gdal

电脑上没有Conda的可以看下文:
Mac 上安装Conda教程

3. 基于Python-GDAL将HDF文件转为GeoTiff格式

import gdal,osr

###
def array2raster(newRasterfn, rasterOrigin, xsize, ysize, array):
    """
     newRasterfn: 输出tif路径
     rasterOrigin: 原始栅格数据Extent
     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()

###
from os import listdir
import re

inputf = '/Users/jerseyshen/Documents/Global_NDVI/'
output = '/Users/jerseyshen/Documents/Global_NDVI_new/'
files = listdir(inputf)

pattern = '.hdf$'
timer =  0 
for i in files:
    temp_input = inputf+i
    
    temp_out_name = re.split(pattern,i)[0]+'.tif'
    temp_output = output+temp_out_name
    
    df = ds = gdal.Open(temp_input)
    subdatasets = ds.GetSubDatasets()
    ndvi_ds = gdal.Open(subdatasets[0][0]).ReadAsArray()
    
    xsize = 0.05
    ysize = 0.05    
    #对于全球尺度rasterOrigin 为[0(xmin),-90(ymin)]
    array2raster(temp_output,[0,-90],xsize,ysize,ndvi_ds) 
    timer = timer + 1
    print(timer)

4. GeoTiff文件在R语言中的呈现

image.png