当前位置：首页 > news >正文

（Arcgis）Python编程批量将HDF5文件转换为TIFF格式并应用地理转换和投影信息

news 来源：原创 2024/4/27 22:25:10

国家青藏高原科学数据中心下载中国1千米分辨率逐日全天候地表土壤水分数据集（2003-2022）
问题：数据在arcgis打开特别大，无法和矢量数据重合，没有设置地理坐标系
数据在网站上提供了投影信息，提示可以进行python编程

环境：pycharm2020python3.8（编译器用的是arcgis pro自带的python3.8）

一、打开HDF5数据投影信息（文档）如下

1.地理转换六参数：
(1)	（左上角X坐标）7087254.893403063,
(2)（X向分辨率）926.6258333333334,
(3)（X向仿射变换角度）0.0,
(4)（左上角Y坐标）5896446.849020582,
(5)（Y向仿射变换角度）0.0,
(6)（Y向分辨率）-926.6254166666668

2.该投影信息对应于GDAL开源包中GDALDataset接口的GetProjection()方法返回参数如下：
PROJCS["AMSRE_Coordinate_System",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Cylindrical_Equal_Area"],PARAMETER["standard_parallel_1",30],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]

二、Python编程将HDF5文件转换为TIFF格式并应用地理转换和投影信息

修改这两个就行
1.input_folder 输入文件夹
2.output_folder 输出文件夹

from osgeo import gdal, osr
import os# 定义输入文件夹和输出文件夹路径
input_folder = r'C:\Users\53446\Desktop\1\1'
output_folder = r'C:\Users\53446\Desktop\1\2'# 遍历输入文件夹中的所有HDF5文件
for file_name in os.listdir(input_folder):if file_name.endswith('.h5'):# 构建输入文件路径input_file = os.path.join(input_folder, file_name)# 打开现有的HDF5文件dataset = gdal.Open(input_file)if dataset is not None:# 设置地理转换六参数transform = (7087254.893403063, 926.6258333333334, 0.0, 5896446.849020582, 0.0, -926.6254166666668)dataset.SetGeoTransform(transform)# 设置投影信息projection = 'PROJCS["AMSRE_Coordinate_System",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"]],AUTHORITY["EPSG","4326"]],PROJECTION["Cylindrical_Equal_Area"],PARAMETER["standard_parallel_1",30],PARAMETER["central_meridian",0],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]'srs = osr.SpatialReference()srs.ImportFromWkt(projection)dataset.SetProjection(srs.ExportToWkt())# 构建输出文件路径output_file = os.path.join(output_folder, file_name.replace('.h5', '.tif'))# 使用Translate函数将HDF5文件转换为TIFF格式并应用地理转换和投影信息gdal.Translate(output_file, dataset, format='GTiff')# 关闭文件dataset = Noneelse:print(f"无法打开文件: {input_file}")