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对比

  • 24
    点赞
  • 28
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论
以下是一些示例代码,演示如何将 Sentinel-2 和 Sentinel-1 数据进行融合: 1. 利用Pythonsentinelsat库下载Sentinel-2和Sentinel-1数据: ```python from sentinelsat import SentinelAPI, read_geojson, geojson_to_wkt # 登录sentinelsat账号 api = SentinelAPI('username', 'password', 'https://scihub.copernicus.eu/dhus') # 下载Sentinel-2数据 footprint = geojson_to_wkt(read_geojson('path/to/footprint.geojson')) products = api.query(footprint, platformname='Sentinel-2', cloudcoverpercentage=(0, 30), producttype='S2MSI1C') # 下载Sentinel-1数据 products = api.query(footprint, platformname='Sentinel-1', polarisationmode='VV VH', producttype='GRD', orbitdirection='ASCENDING') ``` 2. 使用Python的gdal库读取和处理Sentinel-1数据: ```python from osgeo import gdal # 读取Sentinel-1数据 s1_vv = gdal.Open('path/to/sentinel1_vv.tif') s1_vh = gdal.Open('path/to/sentinel1_vh.tif') # 将Sentinel-1数据转换为dB单位 s1_vv_db = 10 * np.log10(s1_vv.ReadAsArray()) s1_vh_db = 10 * np.log10(s1_vh.ReadAsArray()) # 对Sentinel-1数据进行滤波和校正 # ... # 将Sentinel-1数据重采样到Sentinel-2的分辨率 # ... # 将Sentinel-1数据和Sentinel-2数据进行融合 # ... ``` 3. 使用Python的scikit-image库将Sentinel-2和Sentinel-1数据进行融合: ```python from skimage import exposure # 将Sentinel-2数据进行拉伸和直方图匹配,使其与Sentinel-1数据的动态范围一致 s2_rgb = exposure.rescale_intensity(s2_rgb, in_range=(0, 0.3), out_range=(0, 1)) s2_rgb_matched = exposure.match_histograms(s2_rgb, s1_vv_db) # 将Sentinel-1数据和Sentinel-2数据进行加权融合 s1_weight = 0.6 s2_weight = 0.4 fused = (s1_weight * s1_vv_db + (1 - s1_weight) * s1_vh_db) * s2_weight + (1 - s2_weight) * s2_rgb_matched ``` 这只是一些示例代码,具体的融合方法和参数需要根据具体的应用场景进行调整和优化。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

炒茄子

不装逼我浑身难受aaa

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值