hydrodataset数据集制作

原创已于 2025-08-19 17:50:52 修改 · 199 阅读

CC 4.0 BY-SA版权

文章标签：

于 2025-08-18 14:41:55 首次发布

conda activate hydromodel

python /home/dell4/workspace/ldn/hydromodel-master/data/Getdata.py

/home/dell4/workspace/ldn/hydromodel-master/data/Getdata.py

import os
import pandas as pd
import numpy as np
import xarray as xr
from datetime import datetime

# 配置参数
output_dir = "/home/dell4/workspace/ldn/hydromodel-master/data"  # 输出目录
basin_ids = ['xxx']  # 流域ID列表
start_date = "2012-06-10"  # 起始日期
end_date = "2022-12-31"    # 结束日期


# 创建输出目录
os.makedirs(output_dir, exist_ok=True)

# 定义气候类型
climate_zones = {
    'temperate': {'prcp_range': (1, 15), 'pet_range': (0.5, 5), 'season_shift': 0},
    'mediterranean': {'prcp_range': (0.1, 25), 'pet_range': (1, 8), 'season_shift': 0.5},
    'tropical': {'prcp_range': (3, 40), 'pet_range': (2, 10), 'season_shift': 0.25}
}

# 定义土壤类型
soil_types = {
    'sandy': {'infiltration': 0.7, 'baseflow': 0.3},
    'clay': {'infiltration': 0.3, 'baseflow': 0.1},
    'loamy': {'infiltration': 0.5, 'baseflow': 0.2}
}

# 生成流域属性表
attributes_data = []
for basin_id in basin_ids:
    # 随机分配气候区和土壤类型
    climate = np.random.choice(list(climate_zones.keys()))
    soil = np.random.choice(list(soil_types.keys()))
    
    # 随机生成流域属性
    area = np.random.randint(50, 2000)  # 流域面积(km²)
    elevation = np.random.randint(50, 3000)  # 海拔高度(m)
    
    attributes_data.append({
        'id': basin_id,
        'name': f'Basin {basin_id.upper()}',
        'area': area,  # 简化列名，程序需要'area'字段
        'climate_zone': climate,
        'soil_type': soil,
        'elevation': elevation,  # 简化列名
        'slope': np.random.uniform(1, 30),  # 简化列名
        'drainage_density': np.random.uniform(0.5, 2.5)  # 简化列名，去掉斜杠
    })

# 保存流域属性并转换为NetCDF格式
attributes_df = pd.DataFrame(attributes_data)
# 保存为CSV（可选）
attributes_csv_path = os.path.join(output_dir, "basin_attributes.csv")
attributes_df.to_csv(attributes_csv_path, index=False)
# 转换为NetCDF格式
attrs_ds = xr.Dataset.from_dataframe(attributes_df.set_index('id'))
attrs_ds.to_netcdf(os.path.join(output_dir, "attributes.nc"))

# 生成每个流域的时间序列数据 - 改为日尺度，仅保留每天零点
date_range = pd.date_range(start=start_date, end=end_date, freq='D')  # 每日数据
# 强制时间格式为YYYY-MM-DD 00:00:00
date_range = pd.to_datetime(date_range.strftime('%Y-%m-%d 00:00:00'))
t = np.arange(len(date_range))  # 时间索引按天计算
all_basin_data = []  # 存储所有流域的时序数据

for _, basin in attributes_df.iterrows():
    basin_id = basin['id']
    climate = basin['climate_zone']
    soil = basin['soil_type']
    area = basin['area']
    elevation = basin['elevation']
    
    # 获取气候参数
    prcp_range = climate_zones[climate]['prcp_range']
    pet_range = climate_zones[climate]['pet_range']
    season_shift = climate_zones[climate]['season_shift']
    
    # 获取土壤参数
    infiltration = soil_types[soil]['infiltration']
    baseflow_coeff = soil_types[soil]['baseflow']
    
    # 1. 生成降水(prcp) - 日尺度
    seasonal = np.sin(2 * np.pi * (t/365.25 + season_shift))  # 按年周期(365.25天)
    if climate == 'tropical':
        seasonal += 0.5 * np.sin(4 * np.pi * t/365.25)  # 热带地区有两个雨季
    elevation_factor = 1 + (elevation/3000) * 0.5
    base_prcp = prcp_range[0] + (prcp_range[1] - prcp_range[0]) * (0.5 + 0.5 * seasonal)
    rain_events = np.random.poisson(0.15, len(t))  # 日降水事件频率
    rain_intensity = np.random.exponential(5, len(t))
    prcp = base_prcp + rain_events * rain_intensity
    prcp = np.maximum(prcp, 0)  # 单位：mm/day
    
    # 2. 生成潜在蒸散发(pet) - 日尺度
    pet_seasonal = np.sin(2 * np.pi * t/365.25)
    base_pet = pet_range[0] + (pet_range[1] - pet_range[0]) * (0.5 + 0.5 * pet_seasonal)
    elevation_pet_factor = 1 - (elevation/3000) * 0.3
    pet = base_pet * elevation_pet_factor
    pet = np.maximum(pet, 0.1)  # 单位：mm/day
    
    # 3. 生成实际蒸散发(et) - 日尺度
    soil_moisture = np.zeros(len(t))
    soil_capacity = 200  # mm
    for i in range(1, len(t)):
        recharge = min(prcp[i] * infiltration, soil_capacity - soil_moisture[i-1])
        et_loss = min(pet[i], soil_moisture[i-1])
        soil_moisture[i] = soil_moisture[i-1] + recharge - et_loss
    et = pet * (soil_moisture / soil_capacity)  # 单位：mm/day
    
    # 4. 生成径流(flow) - 日尺度，使用面流量单位避免转换问题
    slope = basin['slope']
    runoff_coeff = 0.7 * (slope/30)
    surface_runoff = prcp * (1 - infiltration) * runoff_coeff  # 单位：mm/day
    baseflow = baseflow_coeff * (0.5 + 0.5 * np.sin(2 * np.pi * t/365.25 - np.pi/4))  # 单位：mm/day
    runoff_depth = surface_runoff + baseflow  # 总径流深度，单位：mm/day
    
    # 使用面流量单位(mm/day)，与降水单位一致，避免维度转换问题
    flow = runoff_depth  # 单位：mm/day
    
    # 5. 生成节点流量(可选) - 日尺度延迟
    if np.random.rand() > 0.3:
        delay_days = int(area ** 0.3)  # 汇流延迟（天）
        node1_flow = np.zeros(len(flow))
        for i in range(delay_days, len(flow)):
            node1_flow[i] = 0.6 * flow[i-delay_days] + 0.4 * baseflow[i]  # 单位：mm/day
    else:
        node1_flow = None
    
    # 创建数据框
    data = {
        'time': date_range,
        'pet': pet,
        'prcp': prcp,
        'flow': flow,
        'et': et,
        'basin': basin_id  # 添加流域ID
    }
    
    if node1_flow is not None:
        data['node1_flow'] = node1_flow
    
    basin_df = pd.DataFrame(data)
    all_basin_data.append(basin_df)

# 合并所有流域数据并保存为NetCDF格式
combined_df = pd.concat(all_basin_data, ignore_index=True)
# 设置索引以便转换为xarray Dataset
combined_df = combined_df.set_index(['basin', 'time'])
combined_ds = xr.Dataset.from_dataframe(combined_df)

# 为变量添加单位属性（关键：确保单位与日尺度匹配）
combined_ds['prcp'].attrs['units'] = 'mm/day'
combined_ds['pet'].attrs['units'] = 'mm/day'
combined_ds['flow'].attrs['units'] = 'mm/day'
combined_ds['et'].attrs['units'] = 'mm/day'
if 'node1_flow' in combined_ds:
    combined_ds['node1_flow'].attrs['units'] = 'mm/day'

# 保存为timeseries.nc
combined_ds.to_netcdf(os.path.join(output_dir, "timeseries.nc"))

# 数据验证函数
def validate_data(df, basin_id):
    """验证生成的数据是否合理"""
    # 检查负值
    if (df.select_dtypes(include=np.number) < 0).any().any():
        print(f"警告! {basin_id} 包含负值")
    
    # 检查缺失值
    if df.isnull().any().any():
        print(f"警告! {basin_id} 包含缺失值")
    
    # 物理一致性检查
    if 'et' in df.columns and 'pet' in df.columns:
        if (df['et'] > df['pet']).any():
            print(f"警告! {basin_id} 实际ET大于潜在ET")

# 对最后一个流域进行验证
if all_basin_data:
    validate_data(all_basin_data[-1], basin_id)

print(f"水文数据集已生成至目录: {output_dir}")
print(f"包含 {len(basin_ids)} 个流域的数据")
print("生成的文件:")
print("- attributes.nc: 流域属性数据")
print("- timeseries.nc: 所有流域的时序数据（每日零点）")

python /home/dell4/workspace/ldn/hydromodel-master/scripts/calibrate_xaj.py --data_type owndata --data_dir "/home/dell4/workspace/ldn/hydromodel-master/data" --exp expbiliuhe002 --cv_fold 1 --warmup 720 --period "2012-06-10 00:00" "2022-08-31 23:00" --calibrate_period "2012-06-10 00:00" "2017-08-31 23:00" --test_period "2017-09-01 00:00" "2022-08-31 23:00" --basin_id xxx --model "{\"name\": \"xaj\", \"source_type\": \"sources5mm\", \"source_book\": \"HF\"}" --param_range_file "/home/dell4/workspace/ldn/hydromodel-master/XAJ/changdian_61561/param_range.yaml" --algorithm "{\"name\": \"SCE_UA\", \"random_seed\": 1234, \"rep\": 10, \"ngs\": 10, \"kstop\": 5, \"peps\": 0.1, \"pcento\": 0.1}" --loss "{\"type\": \"time_series\", \"obj_func\": \"RMSE\", \"events\": null}"

python /home/dell4/workspace/ldn/hydromodel-master/scripts/evaluate_xaj.py --exp expbiliuhe002

python /home/dell4/workspace/ldn/hydromodel-master/scripts/visualize.py --exp expbiliuhe002

代码的主要功能是生成多个流域的水文数据集，包括流域属性信息和时间序列数据。让我为你总结其核心功能和实现逻辑：

功能概述
- 生成多个流域的属性信息（面积、气候区、土壤类型等）
- 为每个流域生成连续多年的水文时间序列数据
- 包含降水、潜在蒸散发、实际蒸散发、流量等关键水文变量
- 具有数据验证功能，确保生成数据的合理性
技术实现
- 使用 pandas 处理和存储数据
- 利用 numpy 生成符合物理规律的随机数据
- 考虑了气候类型、土壤类型、地形等因素对水文过程的影响
- 实现了简化的水文模型来模拟水流过程
数据生成逻辑
- 流域属性：随机分配气候区和土壤类型，生成面积、海拔等属性
- 降水：考虑季节性变化、海拔影响和随机降水事件
- 蒸散发：基于气候特征和海拔生成，受土壤水分限制
- 径流：结合地表径流和地下径流，考虑坡度等地形因素
- 节点流量：模拟上游节点对流域出口流量的贡献
数据验证
- 检查负值和缺失值
- 验证实际蒸散发不大于潜在蒸散发的物理一致性
输出结果
- 流域属性表（CSV 格式）
- 每个流域的时间序列数据（CSV 格式）
- 数据保存在指定目录中