博客已同步微信公众号:GIS茄子
;若博客出现纰漏或有更多问题交流欢迎关注GIS茄子
,或者邮箱联系(推荐-见主页).
Python:(Sentinel-1)如何解析SNAP输出的HDF5文件并输出为GeoTIFF?
最近在了解sentinel-1的预处理过程,但是由于影响太大了,常规的GeoTIFF无法输出预处理结果,BigTIFF输出时似乎也遇到了一些问题(好在后面解决了,所以正好做一下HDF5文件输出的TIFF文件与BigTIFF文件的对比),对于输出的HDF5文件则完全没有问题。但是问题在于HDF5文件的结构尚不了解,因此对于其中的地理信息如何提取很关键(当然你可以使用ArcGIS或者ENVI打开其中的VV和VH波段,但是都无法自动读取到其中的地理信息或者坐标系信息)。
由于我处理的Sentinel-1时IW的VV和VH,因此输出的HDF5文件存在两个波段:
下方是关于这个地理信息的参数(ps:找了我好久,里面的属性信息真的太多了,而且官方文档似乎对于这个HDF5文件的结构并没有说明,真的象拔蚌了🌿):
未必准确,但是目前从得到的结果与BigTIFF对比是几乎完全一致的地理位置(如果有更详细的文档或者准确信息,请微信公众号或者邮箱联系我,这对我帮助很大)。
first
表示第一行,last
表示最后一行,near
表示扫描线的起点,far
表示扫描线的终点。
其实这里搞不懂为什么要有四个点位的信息?一般的角点信息只需要左上和右下两个点位就足够了,算了我不是这个方向的多说无益。
那么,其实说到这里其实已经搞定了,WGS84坐标系有了,仿射参数也已经有了,VV和VH波段数据也有了。
# @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)
基本上一致