MODIS数据批量处理之HDF多栅格文件提取特征信息（未拼接）

本文链接：https://blog.csdn.net/nuliplus/article/details/138604986

对于MOD06_L2、MOD07_L2等Level 2 swath 产品，处理起来比较麻烦。。。。。。

记录一下。。。

目前没尝试：是先把栅格图像拼接后提取所需特征信息还是反着来，哪个效率更高。

python+GDAL批处理方法

参考博文：python+gdal实现MODIS_HDF文件转TIF（可批量） - 知乎 (zhihu.com)

代码：

from osgeo import gdal
import numpy as np
import os
from osgeo import osr

# #  gdal打开hdf数据集
# datasets = gdal.Open(r"F:\MOD06_L2\HDF\MOD06_L2.A2015001.0220.061.2017318204746.hdf")

# #  获取hdf中的子数据集
# SubDatasets = datasets.GetSubDatasets()
# #  获取子数据集的个数
# SubDatasetsNum =  len(datasets.GetSubDatasets())
# #  输出各子数据集的信息
# print("子数据集一共有{0}个: ".format(SubDatasetsNum))
# for i in range(SubDatasetsNum):
#     print(datasets.GetSubDatasets()[i])

# #  获取hdf中的元数据
# Metadata = datasets.GetMetadata()
# #  获取元数据的个数
# MetadataNum = len(Metadata)
# #  输出各子数据集的信息
# print("元数据一共有{0}个: ".format(MetadataNum))
# for key,value in Metadata.items():
#     print('{key}:{value}'.format(key = key, value = value))
# DatasetTS = datasets.GetSubDatasets()[20][0]
# DatasetPS = datasets.GetSubDatasets()[21][0]
# RasterTS = gdal.Open(DatasetTS)
# TS = RasterTS.ReadAsArray()

# print(DatasetTS)
# print(RasterTS)
# print(TS.shape)


#  数组保存为tif
def array2raster(TifName, GeoTransform, array):
    cols = array.shape[1]  # 矩阵列数
    rows = array.shape[0]  # 矩阵行数
    driver = gdal.GetDriverByName('GTiff')
    outRaster = driver.Create(TifName, cols, rows, 1, gdal.GDT_Float32)
    # 括号中两个0表示起始像元的行列号从(0,0)开始
    outRaster.SetGeoTransform(tuple(GeoTransform))
    # 获取数据集第一个波段，是从1开始，不是从0开始
    outband = outRaster.GetRasterBand(1)
    outband.WriteArray(array)
    outRasterSRS = osr.SpatialReference()
    # 代码4326表示WGS84坐标
    outRasterSRS.ImportFromEPSG(4326)
    outRaster.SetProjection(outRasterSRS.ExportToWkt())
    outband.FlushCache()
    
#  hdf批量转tif
def hdf2tif_batch(hdfFolder):
    #  获取文件夹内的文件名
    hdfNameList = os.listdir(hdfFolder)
    for i in range(len(hdfNameList)):
        #  判断当前文件是否为HDF文件
        if(os.path.splitext(hdfNameList[i])[1] == ".hdf"):
            hdfPath = hdfFolder+"/"+hdfNameList[i]
            #  gdal打开hdf数据集
            datasets = gdal.Open(hdfPath)
            
            #  获取hdf中的元数据
            Metadata = datasets.GetMetadata()
            #  获取四个角的维度
            Latitudes = Metadata["GRIN