python处理1km土壤表层水分数据【适合python小白】

# 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 中可以做自己的后期处理了！