ouyangwenyu / hydromodel Goto Github PK
View Code? Open in Web Editor NEW新安江水文模型
License: GNU General Public License v3.0
新安江水文模型
License: GNU General Public License v3.0
您好,我使用仓库中的样例数据运行了 hydromodel.models.xaj 方法,返回值q_sim是一个数组。他代表的是什么意思呢,是每个网格的预测流量吗
运行calibrate_xaj.py时,无论如何修改时间长度总是报错IndexError,调试后发现是xaj.py第880行索引j超出范围,但Traceback中错误信息显示为axis 0
将xaj.py中j范围改为range(len(l)-1)后此bug不出现,变为RMSE计算中出现空值
Please add documentation
作者您好,可以大概给介绍一下代码当中用于训练的数据的构成么?感激不尽!
您好,很感谢您能分享这么优秀且有意义的项目。但是实际应用中,非源头水域水库的入库计算的需求感觉还是比较大的。为了适用不是源头的流域,需要考虑上游的汇流情况。想着是不是可以将马斯京根算法集成到这个项目,以适应复杂水域的情况,以下是自己的一个简单想法及简单实现,不知道对不对。
Describe the feature (e.g., new functions/tutorials) you would like to propose.
Tell us what can be achieved with this new feature and what's the expected outcome.
##马斯京根算法的计算流程。
class MuskingumReach:
def __init__(self, k, x):
self.k = k
self.x = x
def route_hydrograph(self, inflows, dt):
outflows = list()
outflows.append((inflows[0]))
c0 = ((dt / self.k) - (2 * self.x)) / ((2 * (1 - self.x)) + (dt / self.k))
c1 = ((dt / self.k) + (2 * self.x)) / ((2 * (1 - self.x)) + (dt / self.k))
c2 = ((2 * (1 - self.x)) - (dt / self.k)) / ((2 * (1 - self.x)) + (dt / self.k))
for i in range(len(inflows) - 1):
q_out = (c0 * inflows[i + 1]) + (c1 * inflows[i]) + (c2 * outflows[i])
q_out = max(min(inflows), q_out)
outflows.append(q_out)
return outflows
这个是仓库代码的XAJ模型部分,将马斯京根算法集成到xaj里面,汇总得到最终的汇流返回。
def xaj(
p_and_e,
params: Union[np.array, list],
return_state=False,
warmup_length=365,
**kwargs,
) -> Union[tuple, np.array]:
"""
run XAJ model
Parameters
----------
p_and_e
prcp and pet; sequence-first (time is the first dim) 3-d np array: [time, basin, feature=2]
params
parameters of XAJ model for basin(s);
2-dim variable -- [basin, parameter]:
the parameters are B IM UM LM DM C SM EX KI KG A THETA CI CG (notice the sequence)
return_state
if True, return state values, mainly for warmup periods
warmup_length
hydro models need a warm-up period to get good initial state values
kwargs
route_method
now we provide two ways: "CSL" (recession constant + lag time) and "MZ" (method from mizuRoute)
source_type
default is "sources" and it will call "sources" function; the other is "sources5mm",
and we will divide the runoff to some <5mm pieces according to the books in this case
source_book
When source_type is "sources5mm" there are two implementions for dividing sources,
as the methods in "ShuiWenYuBao" and "GongChengShuiWenXue"" are different.
Hence, both are provided, and the default is the former.
Returns
-------
Union[np.array, tuple]
streamflow or (streamflow, states)
"""
# default values for some function parameters
model_name = kwargs.get("name", "xaj")
source_type = kwargs.get("source_type", "sources")
source_book = kwargs.get("source_book", "HF")
pr_file = kwargs.get("param_range_file", None)
model_param_dict = read_model_param_dict(pr_file)
# params
param_ranges = model_param_dict[model_name]["param_range"]
if model_name == "xaj":
route_method = "CSL"
elif model_name == "xaj_mz":
route_method = "MZ"
else:
raise NotImplementedError(
"We don't provide this route method now! Please use 'CSL' or 'MZ'!"
)
if np.isnan(params).any():
raise ValueError(
"Parameters contain NaN values. Please check your opt algorithm"
)
xaj_params = [
(value[1] - value[0]) * params[:, i] + value[0]
for i, (key, value) in enumerate(param_ranges.items())
]
k = xaj_params[0]
b = xaj_params[1]
im = xaj_params[2]
um = xaj_params[3]
lm = xaj_params[4]
dm = xaj_params[5]
c = xaj_params[6]
sm = xaj_params[7]
ex = xaj_params[8]
ki_ = xaj_params[9]
kg_ = xaj_params[10]
# ki+kg should be smaller than 1; if not, we scale them
ki = np.where(ki_ + kg_ < 1.0, ki_, (1.0 - PRECISION) / (ki_ + kg_) * ki_)
kg = np.where(ki_ + kg_ < 1.0, kg_, (1.0 - PRECISION) / (ki_ + kg_) * kg_)
if route_method == "CSL":
cs = xaj_params[11]
l = xaj_params[12]
elif route_method == "MZ":
# we will use routing method from mizuRoute -- http://www.geosci-model-dev.net/9/2223/2016/
a = xaj_params[11]
theta = xaj_params[12]
# make it as a parameter
kernel_size = int(xaj_params[15])
else:
raise NotImplementedError(
"We don't provide this route method now! Please use 'CSL' or 'MZ'!"
)
ci = xaj_params[13]
cg = xaj_params[14]
# initialize state values
if warmup_length > 0:
p_and_e_warmup = p_and_e[0:warmup_length, :, :]
_q, _e, *w0, s0, fr0, qi0, qg0 = xaj(
p_and_e_warmup,
params,
return_state=True,
warmup_length=0,
**kwargs,
)
else:
w0 = (0.5 * um, 0.5 * lm, 0.5 * dm)
s0 = 0.5 * sm
fr0 = np.full(ex.shape, 0.1)
qi0 = np.full(ci.shape, 0.1)
qg0 = np.full(cg.shape, 0.1)
# state_variables
inputs = p_and_e[warmup_length:, :, :]
runoff_ims_ = np.full(inputs.shape[:2], 0.0)
rss_ = np.full(inputs.shape[:2], 0.0)
ris_ = np.full(inputs.shape[:2], 0.0)
rgs_ = np.full(inputs.shape[:2], 0.0)
es_ = np.full(inputs.shape[:2], 0.0)
for i in range(inputs.shape[0]):
if i == 0:
'''
这个函数实现了新安江模型中的单步径流产生过程。在水文学中,径流是指雨水通过地表流向河流、湖泊等水体的过程,是水文循环中重要的组成部分之一。
'''
(r, rim, e, pe), w = generation(
inputs[i, :, :], k, b, im, um, lm, dm, c, *w0
) ###Runoff(径流) Runoff from impermeable areas(不透水面径流) Evapotranspiration(蒸发蒸腾量) Net Precipitation(净降水)
if source_type == "sources":
(rs, ri, rg), (s, fr) = sources(
pe, r, sm, ex, ki, kg, s0, fr0, book=source_book
)
elif source_type == "sources5mm":
(rs, ri, rg), (s, fr) = sources5mm(
pe, r, sm, ex, ki, kg, s0, fr0, book=source_book
) ## 总的地表径流 总的壤中流 总的地下径流 模拟得到的最终的自由水蓄水容量深度 模拟得到的最终的初始面积。
else:
raise NotImplementedError("No such divide-sources method")
else:
(r, rim, e, pe), w = generation(
inputs[i, :, :], k, b, im, um, lm, dm, c, *w
)
if source_type == "sources":
(rs, ri, rg), (s, fr) = sources(
pe, r, sm, ex, ki, kg, s, fr, book=source_book
)
elif source_type == "sources5mm":
(rs, ri, rg), (s, fr) = sources5mm(
pe, r, sm, ex, ki, kg, s, fr, book=source_book
)
else:
raise NotImplementedError("No such divide-sources method")
# impevious part is pe * im
runoff_ims_[i, :] = rim
# so for non-imprvious part, the result should be corrected
rss_[i, :] = rs * (1 - im) ### im 是透水率(impermeability)
ris_[i, :] = ri * (1 - im)
rgs_[i, :] = rg * (1 - im)
es_[i, :] = e
# seq, batch, feature
runoff_im = np.expand_dims(runoff_ims_, axis=2)
rss = np.expand_dims(rss_, axis=2)
es = np.expand_dims(es_, axis=2)
qs = np.full(inputs.shape[:2], 0.0)
if route_method == "CSL": ###qs 是输出的流量(汇流后的径流)。qt 是通过简单叠加 qs_、qi 和 qg 计算得到的流量。lag 是汇流时段的滞后时间。cs 是汇流过程中的汇流常数。
qt = np.full(inputs.shape[:2], 0.0)
for i in range(inputs.shape[0]):
if i == 0:
qi = linear_reservoir(ris_[i], ci, qi0) #地表径流
qg = linear_reservoir(rgs_[i], cg, qg0) #地下径流
else:
qi = linear_reservoir(ris_[i], ci, qi)
qg = linear_reservoir(rgs_[i], cg, qg)
qs_ = rss_[i]
qt[i, :] = qs_ + qi + qg
for j in range(len(l)):
lag = int(l[j])
for i in range(lag):
qs[i, j] = qt[i, j]
for i in range(lag, inputs.shape[0]):
qs[i, j] = cs[j] * qs[i - 1, j] + (1 - cs[j]) * qt[i - lag, j]
elif route_method == "MZ":
rout_a = a.repeat(rss.shape[0]).reshape(rss.shape)
rout_b = theta.repeat(rss.shape[0]).reshape(rss.shape)
conv_uh = uh_gamma(rout_a, rout_b, kernel_size)
qs_ = uh_conv(runoff_im + rss, conv_uh)
for i in range(inputs.shape[0]):
if i == 0:
qi = linear_reservoir(ris_[i], ci, qi0)
qg = linear_reservoir(rgs_[i], cg, qg0)
else:
qi = linear_reservoir(ris_[i], ci, qi)
qg = linear_reservoir(rgs_[i], cg, qg)
qs[i, :] = qs_[i, :, 0] + qi + qg
else:
raise NotImplementedError(
"We don't provide this route method now! Please use 'CS' or 'MZ'!"
)
# seq, batch, feature
q_sim = np.expand_dims(qs, axis=2)
if return_state:
return q_sim, es, *w, s, fr, qi, qg
******马斯京根算法集成在这里进行上游的汇流演算******
q_sim = q_sim + q_mk1 + q_mk2 + q_mk3
*************************************************************
return q_sim, es
不知道这种方式是否可行?不知道作者大大什么时候可以实现相关的半分布式模型?可以加您个联系方式吗?想请教一下您?
Ouyang你好,我目前也从事新安江水文模型的相关研究,有幸学习到您的代码。在研读wiki中实测径流分析模块时,发现缺少退水曲线和洪水分割这两部分的代码。而我自身编程水平有限,还在学习作者的思路中,因此想请问作者什么时候更新这部分的code,或者我们可以交流探讨一下。
你好,我看示例数据中时段是天为单位的,我想问下如果我想跑3h时段长的数据是否需要修改些内容?我制作输入文件试了下是可以跑通,但是不太确定计算是否可以直接用于小时时段的,参数需要处理下吗?
在运行test.py 的时候程序报错:没有hydroeval 和 platypus 模块 。 请问这两个是什么库?
在终端中pip install 对应的库过后,任然找不到对应的子模块
(hydromodel) C:\Users\Harry\OneDrive\hydromodel file\PRBxaj>python calibrate_xaj.py --data_type owndata --data_dir "C:/Users/Harry/OneDrive /hydromodel file/PRBxaj" --exp bada001 --cv_fold 1 --warmup 720 --period "2010-01-01 00:00" "2024-12-31 23:00" --calibrate_period "2012-01-01 00:00" "2020-12-31 23:00" --test_period "2021-01-01 00:00" "2023-12-31 23:00" --basin_id bada --model "{"name": "xaj", "source_type": "sources5mm", "source_book": "HF"}" --param_range_file "C:/Users/14992/OneDrive/hydromodel file/hydromodel-master/hydromodel/models/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}"
[18:53:19] DEBUG No module named 'forge' signatures.py:40
DEBUG No module named 'forge' signatures.py:40
'charmap' codec can't decode byte 0x8d in position 115: character maps to
Traceback (most recent call last):
File "C:\Users\Harry\OneDrive\hydromodel file\PRBxaj\calibrate_xaj.py", line 23, in
from hydromodel.datasets.data_preprocess import (
File "C:\Users\Harry\OneDrive\hydromodel file\PRBxaj\hydromodel\datasets\data_preprocess.py", line 23, in
from hydromodel import CACHE_DIR, SETTING
ImportError: cannot import name 'SETTING' from 'hydromodel' (C:\Users\Harry\OneDrive\hydromodel file\PRBxaj\hydromodel_init_.py)
我发现模型里面没有用到二维数据,比如代表地形的dem数据,代表流域边界的shp数据等。模型中没有用到这些数据吗
Traceback (most recent call last):
File "G:\yanling\bysj\xinanj\hydro-model-xaj-master\hydro-model-xaj-master\hydromodel\data\data_postprocess.py", line 263, in
read_and_save_et_ouputs(one_model_one_hyperparam_setting_dir, fold=0)
File "G:\yanling\bysj\xinanj\hydro-model-xaj-master\hydro-model-xaj-master\hydromodel\data\data_postprocess.py", line 226, in read_and_save_et_ouputs
_, e_train = xaj(
File "G:\yanling\bysj\xinanj\hydro-model-xaj-master\hydro-model-xaj-master\hydromodel\models\xaj.py", line 785, in xaj
_q, _e, *w0, s0, fr0, qi0, qg0 = xaj(
File "G:\yanling\bysj\xinanj\hydro-model-xaj-master\hydro-model-xaj-master\hydromodel\models\xaj.py", line 808, in xaj
(r, rim, e, pe), w = generation(
File "G:\yanling\bysj\xinanj\hydro-model-xaj-master\hydro-model-xaj-master\hydromodel\models\xaj.py", line 222, in generation
eu, el, ed = calculate_evap(lm, c, wu0, wl0, prcp, pet)
File "D:\miniconda\envs\xaj\lib\site-packages\numba\np\arraymath.py", line 3481, in where_impl
raise ValueError("all inputs should have the same shape")
ValueError: all inputs should have the same shape
Process finished with exit code 1
Please Check your directory:
ROOT_DIR of the repo: D:\hydro-model-xaj-master\hydro-model-xaj-master
DATASET_DIR of the repo: D:\data
Compare evaluation results of different calibrated models in an exp directory
All settings in an exp directory: 0it [00:00, ?it/s]
Traceback (most recent call last):
File "datapostprocess4calibrate.py", line 156, in
statistics(the_args)
File "datapostprocess4calibrate.py", line 112, in statistics
mean_df_test = pd.concat(mean_lst_test, axis=1).T
File "D:\ProgramData\Anaconda3\envs\xaj\lib\site-packages\pandas\util_decorators.py", line 311, in wrapper
return func(*args, **kwargs)
File "D:\ProgramData\Anaconda3\envs\xaj\lib\site-packages\pandas\core\reshape\concat.py", line 304, in concat
sort=sort,
File "D:\ProgramData\Anaconda3\envs\xaj\lib\site-packages\pandas\core\reshape\concat.py", line 351, in init
raise ValueError("No objects to concatenate")
ValueError: No objects to concatenate
Appreciate for your work.Here I have some questions. Are the parameters "a and theta" just same as "n and k"? If not ,how to convert?
你好,我看到您的example中给出的训练集测试集中他们是将4个流域的数据放到一起,并且他们有着相同的时间长度,我现在想要使用自己的数据集进行训练,我的数据是一个流域的多个场次洪水,他们不是连续的时间序列,那么我该如何进行处理
暂时解决,感谢!
我看到readme文件中提到模型的数据是一个三维数组[时间,流域,变量]。那这里的变量是指什么,我看到代码中变量应该是降雨和蒸散发数据,对吗。
换句话说,预报模型中,需要降雨数据/蒸散发数据和率定出来的16个变量这三类数据对吗。
你好,对differentiable version of XAJ的实现十分感兴趣,但在hydromodel/models/路径下的文件中没有找到相关代码,可以提供下相关访问方式吗
用pycharm运行calibrate_xaj.py时提示出现该问题,但是所有包都安装完毕,也没有修改程序。请教是什么原因
不同的新安江模型好像有不同的参数,我有调研到一种实现的参数如下:
WM
WUM
WLM
KC
B
SMD
EX
SKIGD
KID
IMP
C
CID
CGD
N
NK
您知道这些参数与该项目参数的区别与联系吗?是否可以与本项目的参数相互等价替换。
以上的模型实现中,汇流阶段好像采用的是纳什瞬时单位线,而不是采用CS,L实现的(你能给出对应实现吗?)。
我的简单实现如下:
if route_method == "CSL":
qt = np.full(inputs.shape[:2], 0.0)
for i in range(inputs.shape[0]):
if i == 0:
qi = linear_reservoir(ris_[i], ci, qi0)
qg = linear_reservoir(rgs_[i], cg, qg0)
else:
qi = linear_reservoir(ris_[i], ci, qi)
qg = linear_reservoir(rgs_[i], cg, qg)
qs_ = rss_[i]
qt[i, :] = qs_ + qi + qg
for j in range(inputs.shape[1]):
qs[:, j] = nash_iu(qt[:, j], N, NK, inputs.shape[0])
# for j in range(len(l)):
# lag = int(l[j])
# for i in range(lag):
# qs[i, j] = qt[i, j]
# for i in range(lag, inputs.shape[0]):
# qs[i, j] = cs[j] * qs[i - 1, j] + (1 - cs[j]) * qt[i - lag, j]
def nash_iu(r, k, n, time_steps):
"""
计算纳什瞬时单位线输出
r: 输入径流 (输入时间序列)
k: 储水库时间常数
n: 储水库数量
time_steps: 总时间步数
"""
from scipy.special import gamma, gammaincc
t = np.arange(time_steps)
unit_hydrograph = (t ** (n - 1) * np.exp(-t / k)) / (k ** n * gamma(n))
outflow = np.convolve(r, unit_hydrograph)[:time_steps]
return outflow
还有个问题想请教一下是,我已经初步将马斯京根算法集成到该项目,但是率定效果不是很理想,率定的NSE才80出头,想请教一下您有哪些率定的技巧吗?我目前主要是调整sce_ua的超参数。好像有种做法是通过场次洪水率定,这种情况可以直接获取场次洪水数据进行输入率定还是输入所有数据但是仅计算洪水区间的loss?
default={
"name": "SCE_UA",
"random_seed": 2048,
# these params are just for test
"rep": 1500,
"ngs": 22,
"kstop": 5,
"peps": 0.1,
"pcento": 0.1,
},
期待您的回复。
我想问一下关于数据camels_us用于新安江率定的参数初始值是多少?从哪个文件看呢
水平有限,还没跑起来。wiki里面对模型原理的讲解非常详细,有很多探究和思考。比较惊叹贡献者的知识体系和数学功底。
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.