Python：（Sentinel-1）如何解析SNAP输出的HDF5文件并输出为GeoTIFF？

博客已同步微信公众号：GIS茄子；若博客出现纰漏或有更多问题交流欢迎关注GIS茄子，或者邮箱联系(推荐-见主页).

Python：（Sentinel-1）如何解析SNAP输出的HDF5文件并输出为GeoTIFF？

01 前言

最近在了解sentinel-1的预处理过程，但是由于影响太大了，常规的GeoTIFF无法输出预处理结果，BigTIFF输出时似乎也遇到了一些问题（好在后面解决了，所以正好做一下HDF5文件输出的TIFF文件与BigTIFF文件的对比），对于输出的HDF5文件则完全没有问题。但是问题在于HDF5文件的结构尚不了解，因此对于其中的地理信息如何提取很关键（当然你可以使用ArcGIS或者ENVI打开其中的VV和VH波段，但是都无法自动读取到其中的地理信息或者坐标系信息）。

02 解析HDF5文件

由于我处理的Sentinel-1时IW的VV和VH，因此输出的HDF5文件存在两个波段：

VV和VH相关波段信息

下方是关于这个地理信息的参数（ps：找了我好久，里面的属性信息真的太多了，而且官方文档似乎对于这个HDF5文件的结构并没有说明，真的象拔蚌了🌿）：

元数据

那么我们来解释一下其中关键的8个参数：

first_near_lat = 30.710711909958782; // double

first_near_long = 106.20485428671394; // double

first_far_lat = 30.710711909958782; // double

first_far_long = 109.12878070499457; // double

last_near_lat = 28.79451557740343; // double

last_near_long = 106.20485428671394; // double

last_far_lat = 28.79451557740343; // double

last_far_long = 109.12878070499457; // double

未必准确，但是目前从得到的结果与BigTIFF对比是几乎完全一致的地理位置（如果有更详细的文档或者准确信息，请微信公众号或者邮箱联系我，这对我帮助很大）。

first 表示第一行，last表示最后一行，near表示扫描线的起点，far表示扫描线的终点。

其实这里搞不懂为什么要有四个点位的信息？一般的角点信息只需要左上和右下两个点位就足够了，算了我不是这个方向的多说无益。

那么，其实说到这里其实已经搞定了，WGS84坐标系有了，仿射参数也已经有了，VV和VH波段数据也有了。

03 代码

# @Author   : ChaoQiezi
# @Time     : 2023/12/18  8:40
# @Email    : chaoqiezi.one@qq.com"""
This script is used to 读取HDF5、BigTIFF文件
"""import os.path
import h5py
from osgeo import gdal, osr# 准备
h5_path = r'H:\Datasets\Objects\TobacooLeafRecognition\Data\HDF5\S1A_IW_GRDH_1SDV_20220602T103546_20220602T103611_043483_05311F_8F62_NR_Orb_Cal_Spk_TC_dB.h5'
tiff_path = r'H:\Datasets\Objects\TobacooLeafRecognition\Data\BigTIFF\S1A_IW_GRDH_1SDV_20220602T103546_20220602T103611_043483_05311F_8F62_NR_Orb_Cal_Spk_TC_dB.tif'
out_dir = r'H:\Datasets\Objects\TobacooLeafRecognition\Data'
out_path = os.path.join(out_dir, 'vv_vh.tiff')
vh_name = 'bands/Sigma0_VH_db'
vv_name = 'bands/Sigma0_VV_db'
metadata_name = 'metadata/Abstracted_Metadata'
lon_min_name = 'first_near_long'
lon_max_name = 'last_far_long'
lat_min_name = 'last_far_lat'
lat_max_name = 'first_near_lat'
lon_res_name = 'lon_pixel_res'
lat_res_name = 'lat_pixel_res'# 探索HDF5文件
with h5py.File(h5_path) as h5:vh, vv = h5[vh_name][:], h5[vv_name][:]metadata = h5[metadata_name]lon_min = metadata.attrs[lon_min_name]lon_max = metadata.attrs[lon_max_name]lat_min = metadata.attrs[lat_min_name]lat_max = metadata.attrs[lat_max_name]lon_res = metadata.attrs[lon_res_name]lat_res = metadata.attrs[lat_res_name]
# 提取栅格信息
rows, cols = vv.shape
transform = [lon_min, lon_res, 0, lat_max, 0, -lon_res]
# 定义地理信息(WGS84)
srs = osr.SpatialReference()
srs.ImportFromEPSG(4326)  # WGS84
# 输出
driver = gdal.GetDriverByName('GTiff')
ds = driver.Create(out_path, cols, rows, 2, gdal.GDT_Float32)
ds.SetProjection(srs.ExportToWkt())  # 设置坐标系
ds.SetGeoTransform(transform)  # 设置仿射参数
[ds.GetRasterBand(_ix+1).WriteArray(_band) for _ix, _band in enumerate([vv, vh])]  # 写入数据
ds.FlushCache()
ds = None
# 探索BigTIFF文件
ds = gdal.Open(tiff_path)
bands = ds.ReadAsArray()
proj = ds.GetProjection()
tiff_transform = ds.GetGeoTransform()
print('HDF5的proj: {}'.format(srs.ExportToWkt()))
print('BigTIFF的proj: {}'.format(proj))
print('HDF5的仿射变换参数: {}'.format(transform))
print('BigTIFF的proj: {}'.format(tiff_transform))

输出：

HDF5的proj: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]BigTIFF的proj: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]HDF5的仿射变换参数: [106.20485428671394, 8.983152841195215e-05, 0, 30.710711909958782, 0, -8.983152841195215e-05]BigTIFF的proj: (106.20485428671394, 8.983152841195215e-05, 0.0, 30.710711909958782, 0.0, -8.983152841195215e-05)

基本上一致

HDF5输出与BigTIFF对比