海洋初级生产力产品，hdf格式，这个数据比较有特色的一点是只有一层，没有子数据集，就是一个BAND
实例代码
from osgeo import gdal,osr
import numpy as np
import os
输入HDF文件路径
hdf_file = 'E:/data_modis/vgpm2011001.hdf'
输出GeoTIFF文件路径
tif_file = 'E:/output/output111.tif'
打开HDF文件
hdf_dataset = gdal.Open(hdf_file)
if hdf_dataset is None:
print('无法打开HDF文件')
exit(1)
获取HDF数据集的大小和波段数
rows = hdf_dataset.RasterYSize
cols = hdf_dataset.RasterXSize
bands = hdf_dataset.RasterCount
第一个子数据集合,也就是NDVI数据
DatasetNDVI = hdf_dataset.GetRasterBand(1)
RasterNDVI = gdal.Open(DatasetNDVI)
NDVI = DatasetNDVI.ReadAsArray()
图像四个角的地理坐标
GeoCoordinates = np.zeros((4, 2), dtype = "float32")
GeoCoordinates[0] = np.array([float(-180),float(90)])
GeoCoordinates[1] = np.array([float(180),float(90)])
GeoCoordinates[2] = np.array([float(180),float(-90)])
GeoCoordinates[3] = np.array([float(-180),float(-90)])
列数
Columns = float(cols)
行数
Rows = float(rows)
图像四个角的图像坐标
PixelCoordinates = np.array([[0, 0],
[Columns - 1, 0],
[Columns - 1, Rows - 1],
[0, Rows - 1]], dtype = "float32")
计算仿射变换矩阵
from scipy.optimize import leastsq
def func(i):
Transform0, Transform1, Transform2, Transform3, Transform4, Transform5 = i[0], i[1], i[2], i[3], i[4], i[5]
return [Transform0 + PixelCoordinates[0][0] * Transform1 + PixelCoordinates[0][1] * Transform2 - GeoCoordinates[0][0],
Transform3 + PixelCoordinates[0][0] * Transform4 + PixelCoordinates[0][1] * Transform5 - GeoCoordinates[0][1],
Transform0 + PixelCoordinates[1][0] * Transform1 + PixelCoordinates[1][1] * Transform2 - GeoCoordinates[1][0],
Transform3 + PixelCoordinates[1][0] * Transform4 + PixelCoordinates[1][1] * Transform5 - GeoCoordinates[1][1],
Transform0 + PixelCoordinates[2][0] * Transform1 + PixelCoordinates[2][1] * Transform2 - GeoCoordinates[2][0],
Transform3 + PixelCoordinates[2][0] * Transform4 + PixelCoordinates[2][1] * Transform5 - GeoCoordinates[2][1],
Transform0 + PixelCoordinates[3][0] * Transform1 + PixelCoordinates[3][1] * Transform2 - GeoCoordinates[3][0],
Transform3 + PixelCoordinates[3][0] * Transform4 + PixelCoordinates[3][1] * Transform5 - GeoCoordinates[3][1]]
最小二乘法求解
GeoTransform = leastsq(func,np.asarray((1,1,1,1,1,1)))
print(GeoTransform)
保存为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()
TifName = "aa20250113new.tif"
array2raster(TifName, GeoTransform[0], NDVI)
关闭数据集
hdf_dataset = None
tif_dataset = None
