# 1km土壤表层水分数据预处理
你是否发现这个数据很好下载,但是使用有门槛?
加载到ArcGIS中没有投影,批量设置正确的投影很难?
为了让数据分析更加正确,预处理非常重要,希望可以让自己的初步数据处理帮助大家减少很多麻烦。
# 步骤
#0.下载数据:
其中,官方网站下载HDF5数据投影信息.docx(我因为没有仔细看文件,自己预测投影浪费了1天时间,对不上分辨率、经纬度范围和栅格数量[吐血])
#1.运行 Visual studio code,在终端中运行
A设置环境
B下载必要的一些包,其中gdal 本身很复杂,因此安装使用conda会比较快
#创建兼容的 Python 环境
conda create -n numpy1_env python=3.10 # 使用 Python 3.10 更稳定
conda activate numpy1_env
pip install "numpy<2" # 强制安装 NumPy 1.x 最新版本(如 1.26.4)
pip install h5py pybind11 "gdal>=3.4" # 确保安装支持旧版 NumPy 的依赖
#检查
pip list
conda install gdal
#2.运行重投影代码(利用chatgpt 写的)
import os
import h5py
# 在文件头部添加环境检查
import numpy as np
print(f"当前 NumPy 版本: {np.__version__}") # 确保输出 1.x 版本
assert np.__version__.startswith('1.'), "必须使用 NumPy 1.x 版本!"
from osgeo import gdal
# 设置输入输出路tif_path径
input_root_dir = r"储存h5数据的文件夹" # 输入目录(需要替换)
output_root_dir = r"储存h5数据投影后tif的文件夹" # 输出目录(需要替换)
# 确保输出目录存在
os.makedirs(output_root_dir, exist_ok=True)
# 投影和地理转换信息
geo_transform = [
7087254.893403063, # 左上角X坐标
926.6258333333334, # X向分辨率
0.0, # X向仿射变换角度
5896446.849020582, # 左上角Y坐标
0.0, # Y向仿射变换角度
-926.6254166666668 # Y向分辨率
]
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]]"""
# 遍历2003到2020年的每一年
for year in range(2004, 2021):#[提供了2003-2023年的日数据,按需取用]
# 设置每年输入和输出路径
input_dir = os.path.join(input_root_dir, str(year)) # 输入路径:根据年份指定
output_tif_path = os.path.join(output_root_dir, f"year_{year}_transformed") # 输出路径:根据年份指定
# 确保输出目录存在
os.makedirs(output_tif_path, exist_ok=True)
# 遍历输入目录中的所有HDF5文件
for filename in os.listdir(input_dir):
if filename.endswith('.h5') or filename.endswith('.hdf5'):
# 构建完整输入路径
input_hdf5_path = os.path.join(input_dir, filename)
# 打开HDF5文件
with h5py.File(input_hdf5_path, 'r') as hdf5_file:
# 假设HDF5文件中只有一个数据集
dataset_name = list(hdf5_file.keys())[0]
data = hdf5_file[dataset_name][...]
# 创建输出文件名(保留原文件名)
output_tif = os.path.join(output_tif_path, f"{os.path.splitext(filename)[0]}.tif")
# 创建并写入GeoTIFF
driver = gdal.GetDriverByName('GTiff')
out_dataset = driver.Create(output_tif, data.shape[1], data.shape[0], 1, gdal.GDT_Float32)
out_dataset.SetGeoTransform(geo_transform)
out_dataset.SetProjection(projection)
out_dataset.GetRasterBand(1).WriteArray(data)
out_dataset = None # 关闭数据集
print(f"{year}年的所有HDF5文件已成功转换为TIFF!")
print("所有HDF5文件已成功转换为TIFF!")
#转换完成后加载到ArcGIS 中可以做自己的后期处理了!