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

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值