From 9865afe3167d1c539a26f881e36f1d168fd17e36 Mon Sep 17 00:00:00 2001 From: ouyangwenyu Date: Sat, 28 Oct 2023 09:49:14 +0800 Subject: [PATCH 1/5] add bim package for env --- environment-dev.yml | 5 +++++ environment.yml | 3 ++- setup.py | 19 ++++++++++++------- test/test_data.py | 2 +- 4 files changed, 20 insertions(+), 9 deletions(-) diff --git a/environment-dev.yml b/environment-dev.yml index 81764f5..288a934 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -13,9 +13,14 @@ dependencies: - spotpy=1.5.14 - seaborn - tqdm + - yaml + - bmipy + # for dev - sphinx - pytest - black + - flake8 + # pip - pip - pip: - hydrodataset diff --git a/environment.yml b/environment.yml index 02544ac..da831b9 100644 --- a/environment.yml +++ b/environment.yml @@ -13,7 +13,8 @@ dependencies: - spotpy=1.5.14 - seaborn - tqdm - - pytest + - yaml + - bmipy - pip - pip: - hydrodataset diff --git a/setup.py b/setup.py index d019c08..e143dd7 100644 --- a/setup.py +++ b/setup.py @@ -1,11 +1,16 @@ -#!/usr/bin/env python -# -*- coding:utf-8 -*- - +""" +Author: Wenyu Ouyang +Date: 2023-10-28 09:16:46 +LastEditTime: 2023-10-28 09:27:22 +LastEditors: Wenyu Ouyang +Description: setup.py for hydromodel package +FilePath: \hydro-model-xaj\setup.py +Copyright (c) 2023-2024 Wenyu Ouyang. All rights reserved. +""" +import pathlib from setuptools import setup, find_packages -with open("README.md") as readme_file: - readme = readme_file.read() - +readme = pathlib.Path("README.md").read_text() setup( name="hydromodel", # 输入项目名称 version="0.0.1", # 输入版本号 @@ -21,5 +26,5 @@ include_package_data=True, platforms="any", install_requires=[""], # 输入项目所用的包 - python_requires='>= 3.6 ', # Python版本要求 + python_requires=">= 3.7 ", # Python版本要求 ) diff --git a/test/test_data.py b/test/test_data.py index 6844fc7..1bbcb2a 100644 --- a/test/test_data.py +++ b/test/test_data.py @@ -1,7 +1,7 @@ """ Author: Wenyu Ouyang Date: 2022-10-25 21:16:22 -LastEditTime: 2022-11-28 14:50:30 +LastEditTime: 2023-10-28 09:18:48 LastEditors: Wenyu Ouyang Description: Test for data preprocess FilePath: \hydro-model-xaj\test\test_data.py From 7a7ad304778dae3472655170eafdd89690c89548 Mon Sep 17 00:00:00 2001 From: ouyangwenyu Date: Sat, 28 Oct 2023 15:56:11 +0800 Subject: [PATCH 2/5] delete xaj dir --- {xaj => hydromodel/app}/dockerfile | 0 .../calibrate}/calibrate_ga_xaj_bmi.py | 2 +- .../calibrate}/calibrate_sceua_xaj_bmi.py | 2 +- {xaj => hydromodel/models}/configuration.py | 2 +- hydromodel/models/xaj.py | 173 +++ {xaj => hydromodel/models}/xaj_bmi.py | 154 ++- {xaj => hydromodel/utils}/constant_unit.py | 0 test/test_xaj_bmi.py | 4 +- xaj/__init__.py | 0 xaj/params.py | 122 -- xaj/xajmodel.py | 1060 ----------------- 11 files changed, 268 insertions(+), 1251 deletions(-) rename {xaj => hydromodel/app}/dockerfile (100%) rename {xaj => hydromodel/calibrate}/calibrate_ga_xaj_bmi.py (99%) rename {xaj => hydromodel/calibrate}/calibrate_sceua_xaj_bmi.py (99%) rename {xaj => hydromodel/models}/configuration.py (92%) rename {xaj => hydromodel/models}/xaj_bmi.py (67%) rename {xaj => hydromodel/utils}/constant_unit.py (100%) delete mode 100644 xaj/__init__.py delete mode 100644 xaj/params.py delete mode 100644 xaj/xajmodel.py diff --git a/xaj/dockerfile b/hydromodel/app/dockerfile similarity index 100% rename from xaj/dockerfile rename to hydromodel/app/dockerfile diff --git a/xaj/calibrate_ga_xaj_bmi.py b/hydromodel/calibrate/calibrate_ga_xaj_bmi.py similarity index 99% rename from xaj/calibrate_ga_xaj_bmi.py rename to hydromodel/calibrate/calibrate_ga_xaj_bmi.py index 555b671..71dd54a 100644 --- a/xaj/calibrate_ga_xaj_bmi.py +++ b/hydromodel/calibrate/calibrate_ga_xaj_bmi.py @@ -18,7 +18,7 @@ from hydromodel.utils import stat from hydromodel.utils.stat import statRmse from hydromodel.visual.hydro_plot import plot_sim_and_obs, plot_train_iteration -from xaj.xaj_bmi import xajBmi +from hydromodel.models.xaj_bmi import xajBmi def evaluate(individual, x_input, y_true, warmup_length, model): diff --git a/xaj/calibrate_sceua_xaj_bmi.py b/hydromodel/calibrate/calibrate_sceua_xaj_bmi.py similarity index 99% rename from xaj/calibrate_sceua_xaj_bmi.py rename to hydromodel/calibrate/calibrate_sceua_xaj_bmi.py index ac7cf9b..404cd4b 100644 --- a/xaj/calibrate_sceua_xaj_bmi.py +++ b/hydromodel/calibrate/calibrate_sceua_xaj_bmi.py @@ -4,7 +4,7 @@ from spotpy.parameter import Uniform, ParameterSet from spotpy.objectivefunctions import rmse from hydromodel.models.model_config import MODEL_PARAM_DICT -from xaj.xaj_bmi import xajBmi +from hydromodel.models.xaj_bmi import xajBmi class SpotSetup(object): diff --git a/xaj/configuration.py b/hydromodel/models/configuration.py similarity index 92% rename from xaj/configuration.py rename to hydromodel/models/configuration.py index 708f736..219abd8 100644 --- a/xaj/configuration.py +++ b/hydromodel/models/configuration.py @@ -1,6 +1,6 @@ import yaml import numpy as np -from xaj.xajmodel import xaj_state as xaj_state +from hydromodel.models.xaj import xaj_state as xaj_state def read_config(config_file): diff --git a/hydromodel/models/xaj.py b/hydromodel/models/xaj.py index 614d3c0..8aaf9ba 100644 --- a/hydromodel/models/xaj.py +++ b/hydromodel/models/xaj.py @@ -8,6 +8,7 @@ from scipy.special import gamma from hydromodel.models.model_config import MODEL_PARAM_DICT +from xaj.params import MODEL_PARAM_DICT PRECISION = 1e-5 @@ -885,3 +886,175 @@ def xaj( if return_state: return q_sim, es, *w, s, fr, qi, qg return q_sim, es + + +def xaj_state( + p_and_e, + params: Union[np.array, list], + return_state=True, + warmup_length=365, + **kwargs, +) -> Union[tuple, np.array]: + # default values for some function parameters + model_name = kwargs["name"] if "name" in kwargs else "xaj" + source_type = kwargs["source_type"] if "source_type" in kwargs else "sources" + source_book = kwargs["source_book"] if "source_book" in kwargs else "HF" + # 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 'CS' or 'MZ'!" + ) + 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 'CS' 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_state( + # 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:, :, :] + inputs = p_and_e[:, :, :] + 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 + ) + 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) + 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": + 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)): + # print(range(len(l))) + lag = int(l[j]) + for i in range(0,lag): + qs[i, j] = qt[i, j] + if i == lag: + break + for i in range(lag, inputs.shape[0]): + # print(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 'CSL' 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 diff --git a/xaj/xaj_bmi.py b/hydromodel/models/xaj_bmi.py similarity index 67% rename from xaj/xaj_bmi.py rename to hydromodel/models/xaj_bmi.py index f8fc628..94f1d02 100644 --- a/xaj/xaj_bmi.py +++ b/hydromodel/models/xaj_bmi.py @@ -2,8 +2,8 @@ from bmipy import Bmi import numpy as np -from xaj import configuration -from xaj.xajmodel import xaj_route, xaj_runoff +from hydromodel.models import configuration +from hydromodel.models.xajmodel import xaj_route, xaj_runoff import datetime import logging @@ -14,46 +14,62 @@ class xajBmi(Bmi): """Empty model wrapped in a BMI interface.""" + name = "hydro-model-xaj" input_var_names = ("precipitation", "ETp") output_var_names = ("ET", "discharge") - var_units = {"precipitation": "mm/day", "ETp": "mm/day", "discharge": "mm/day", "ET": "mm/day"} + var_units = { + "precipitation": "mm/day", + "ETp": "mm/day", + "discharge": "mm/day", + "ET": "mm/day", + } def __init__(self): - """Create a BmiHeat model that is ready for initialization.""" + """Create a Bmi model that is ready for initialization.""" self.time_step = 0 def initialize(self, config_file, params, p_and_e): try: logger.info("xaj: initialize_model") config = configuration.read_config(config_file) - model_name = config['model_name'] - source_type = config['source_type'] - source_book = config['source_book'] + model_name = config["model_name"] + source_type = config["source_type"] + source_book = config["source_book"] # forcing_data = pd.read_csv(config['forcing_data']) model_info = { - 'model_name': model_name, - 'source_type': source_type, - 'source_book': source_book, - + "model_name": model_name, + "source_type": source_type, + "source_book": source_book, } # p_and_e_df, p_and_e = configuration.extract_forcing(forcing_data) - p_and_e_warmup = p_and_e[0:config['warmup_length'], :, :] - self.q_sim_state, self.es_state, self.w0, self.w1, self.w2, self.s0, self.fr0, self.qi0, self.qg0 = configuration.warmup( - p_and_e_warmup, params, config['warmup_length'], model_info) + p_and_e_warmup = p_and_e[0 : config["warmup_length"], :, :] + ( + self.q_sim_state, + self.es_state, + self.w0, + self.w1, + self.w2, + self.s0, + self.fr0, + self.qi0, + self.qg0, + ) = configuration.warmup( + p_and_e_warmup, params, config["warmup_length"], model_info + ) self.params = params - self.warmup_length = config['warmup_length'] - self._start_time_str = config['start_time_str'] - self._end_time_str = config['end_time_str'] - self._time_units = config['time_units'] + self.warmup_length = config["warmup_length"] + self._start_time_str = config["start_time_str"] + self._end_time_str = config["end_time_str"] + self._time_units = config["time_units"] # self.p_and_e_df = p_and_e_df self.p_and_e = p_and_e self.config = config self.model_info = model_info - train_period = config['train_period'] - test_period = config['test_period'] + train_period = config["train_period"] + test_period = config["test_period"] self.train_start_time_str = train_period[0] self.train_end_time_str = train_period[1] self.test_start_time_str = test_period[0] @@ -61,6 +77,7 @@ def initialize(self, config_file, params, p_and_e): except: import traceback + traceback.print_exc() raise @@ -68,28 +85,40 @@ def update(self): """Update model for a single time step.""" self.time_step += 1 - p_and_e_sim = self.p_and_e[self.warmup_length:self.time_step + self.warmup_length] + p_and_e_sim = self.p_and_e[ + self.warmup_length : self.time_step + self.warmup_length + ] # p_and_e_sim = self.p_and_e[0:self.time_step] - self.runoff_im, self.rss_, self.ris_, self.rgs_, self.es_runoff, self.rss = xaj_runoff(p_and_e_sim, - w0=self.w0, s0=self.s0, - fr0=self.fr0, - params_runoff=self.params, - return_state=False, - model_info=self.model_info, - ) + ( + self.runoff_im, + self.rss_, + self.ris_, + self.rgs_, + self.es_runoff, + self.rss, + ) = xaj_runoff( + p_and_e_sim, + w0=self.w0, + s0=self.s0, + fr0=self.fr0, + params_runoff=self.params, + return_state=False, + model_info=self.model_info, + ) if self.time_step + self.warmup_length >= self.p_and_e.shape[0]: - q_sim, es = xaj_route(p_and_e_sim, - params_route=self.params, - model_info=self.model_info, - runoff_im=self.runoff_im, - rss_=self.rss_, - ris_=self.ris_, - rgs_=self.rgs_, - rss=self.rss, - qi0=self.qi0, - qg0=self.qg0, - es=self.es_runoff, - ) + q_sim, es = xaj_route( + p_and_e_sim, + params_route=self.params, + model_info=self.model_info, + runoff_im=self.runoff_im, + rss_=self.rss_, + ris_=self.ris_, + rgs_=self.rgs_, + rss=self.rss, + qi0=self.qi0, + qg0=self.qg0, + es=self.es_runoff, + ) self.p_sim = p_and_e_sim[:, :, 0] self.e_sim = p_and_e_sim[:, :, 1] # q_sim = convert_unit( @@ -128,7 +157,7 @@ def get_var_grid(self, name: str) -> int: raise NotImplementedError() def get_var_type(self, name: str) -> str: - return 'float64' + return "float64" def get_var_units(self, name: str) -> str: return self.var_units[name] @@ -143,26 +172,26 @@ def get_var_location(self, name: str) -> str: raise NotImplementedError() def get_start_time(self, period): - if period == 'train': + if period == "train": return self.start_Time(self.train_start_time_str) - if period == 'test': + if period == "test": return self.start_Time(self.test_start_time_str) # return self.start_Time(self._start_time_str) def get_current_time(self): # return self.start_Time(self._start_time_str) + datetime.timedelta(self.time_step+self.warmup_length) - if self._time_units == 'hours': + if self._time_units == "hours": time_step = datetime.timedelta(hours=self.time_step) - elif self._time_units == 'days': + elif self._time_units == "days": time_step = datetime.timedelta(days=self.time_step) return self.start_Time(self._start_time_str) + time_step def get_end_time(self, period): - if period == 'train': + if period == "train": return self.end_Time(self.train_end_time_str) - if period == 'test': + if period == "test": return self.end_Time(self.test_end_time_str) - # return self.end_Time(self._end_time_str) + # return self.end_Time(self._end_time_str) def get_time_units(self) -> str: return self._time_units @@ -175,24 +204,20 @@ def get_value(self, name: str) -> None: return self.get_value_ptr(name).flatten() def get_value_ptr(self, name: str) -> np.ndarray: - if name == 'discharge': + if name == "discharge": return self.q_sim - elif name == 'ET': + elif name == "ET": return self.es - def get_value_at_indices( - self, name: str, inds: np.ndarray - ) -> np.ndarray: - + def get_value_at_indices(self, name: str, inds: np.ndarray) -> np.ndarray: return self.get_value_ptr(name).take(inds) def set_value(self, name: str, src: np.ndarray): - val = self.get_value_ptr(name) val[:] = src.reshape(val.shape) def set_value_at_indices( - self, name: str, inds: np.ndarray, src: np.ndarray + self, name: str, inds: np.ndarray, src: np.ndarray ) -> None: val = self.get_value_ptr(name) val.flat[inds] = src @@ -246,12 +271,11 @@ def get_grid_face_nodes(self, grid: int, face_nodes: np.ndarray) -> np.ndarray: raise NotImplementedError() def get_grid_nodes_per_face( - self, grid: int, nodes_per_face: np.ndarray + self, grid: int, nodes_per_face: np.ndarray ) -> np.ndarray: raise NotImplementedError() def start_Time(self, _start_time_str): - try: if " " in _start_time_str: date, time = _start_time_str.split(" ") @@ -263,18 +287,18 @@ def start_Time(self, _start_time_str): if time: hour, minute, second = time.split(":") - # self._startTime = self._startTime.replace(hour=int(hour), - # minute=int(minute), + # self._startTime = self._startTime.replace(hour=int(hour), + # minute=int(minute), # second=int(second)) - self._startTime = datetime.datetime(int(year), int(month), int(day), int(hour), int(minute), - int(second)) + self._startTime = datetime.datetime( + int(year), int(month), int(day), int(hour), int(minute), int(second) + ) except ValueError: raise ValueError("Invalid start date format!") return self._startTime def end_Time(self, _end_time_str): - try: if " " in _end_time_str: date, time = _end_time_str.split(" ") @@ -286,7 +310,9 @@ def end_Time(self, _end_time_str): if time: hour, minute, second = time.split(":") - self._endTime = datetime.datetime(int(year), int(month), int(day), int(hour), int(minute), int(second)) + self._endTime = datetime.datetime( + int(year), int(month), int(day), int(hour), int(minute), int(second) + ) except ValueError: raise ValueError("Invalid start date format!") return self._endTime diff --git a/xaj/constant_unit.py b/hydromodel/utils/constant_unit.py similarity index 100% rename from xaj/constant_unit.py rename to hydromodel/utils/constant_unit.py diff --git a/test/test_xaj_bmi.py b/test/test_xaj_bmi.py index 0ee679b..c61808c 100644 --- a/test/test_xaj_bmi.py +++ b/test/test_xaj_bmi.py @@ -1,8 +1,8 @@ import logging import definitions -from xaj.configuration import read_config -from xaj.xaj_bmi import xajBmi +from hydromodel.models.configuration import read_config +from hydromodel.models.xaj_bmi import xajBmi logging.basicConfig(level=logging.INFO) diff --git a/xaj/__init__.py b/xaj/__init__.py deleted file mode 100644 index e69de29..0000000 diff --git a/xaj/params.py b/xaj/params.py deleted file mode 100644 index 49f3393..0000000 --- a/xaj/params.py +++ /dev/null @@ -1,122 +0,0 @@ -""" -Author: Wenyu Ouyang -Date: 2022-10-25 21:16:22 -LastEditTime: 2022-12-25 16:06:05 -LastEditors: Wenyu Ouyang -Description: some basic config for hydro-model-xaj models -FilePath: \hydro-model-xaj\hydromodel\models\model_config.py -Copyright (c) 2021-2022 Wenyu Ouyang. All rights reserved. -""" -from collections import OrderedDict - -# NOTE: Don't change the parameter settings - -MODEL_PARAM_DICT = { - "xaj": { - "param_name": [ - # Allen, R.G., L. Pereira, D. Raes, and M. Smith, 1998. - # Crop Evapotranspiration, Food and Agriculture Organization of the United Nations, - # Rome, Italy. FAO publication 56. ISBN 92-5-104219-5. 290p. - "K", # ratio of potential evapotranspiration to reference crop evaporation generally from Allen, 1998 - "B", # The exponent of the tension water capacity curve - "IM", # The ratio of the impervious to the total area of the basin - "UM", # Tension water capacity in the upper layer - "LM", # Tension water capacity in the lower layer - "DM", # Tension water capacity in the deepest layer - "C", # The coefficient of deep evapotranspiration - "SM", # The areal mean of the free water capacity of surface soil layer - "EX", # The exponent of the free water capacity curve - "KI", # Outflow coefficients of interflow - "KG", # Outflow coefficients of groundwater - "CS", # The recession constant of channel system - "L", # Lag time - "CI", # The recession constant of the lower interflow - "CG", # The recession constant of groundwater storage - ], - "param_range": OrderedDict( - { - "K": [0.1, 0.1], - "B": [0.1, 0.1], - "IM": [0.1, 0.1], - "UM": [10, 10.0], - "LM": [70.0, 70.0], - "DM": [80.0, 80.0], - "C": [0.1, 0.1], - "SM": [50, 50], - "EX": [1.2, 1.2], - "KI": [0.3, 0.3], - "KG": [0.3, 0.3], - "CS": [0.5, 0.5], - "L": [5.0, 5.0], # unit is day - "CI": [0.5, 0.5], - "CG": [0.98, 0.98], - } - ), - }, - "xaj_mz": { - "param_name": [ - # Allen, R.G., L. Pereira, D. Raes, and M. Smith, 1998. - # Crop Evapotranspiration, Food and Agriculture Organization of the United Nations, - # Rome, Italy. FAO publication 56. ISBN 92-5-104219-5. 290p. - "K", # ratio of potential evapotranspiration to reference crop evaporation generally from Allen, 1998 - "B", # The exponent of the tension water capacity curve - "IM", # The ratio of the impervious to the total area of the basin - "UM", # Tension water capacity in the upper layer - "LM", # Tension water capacity in the lower layer - "DM", # Tension water capacity in the deepest layer - "C", # The coefficient of deep evapotranspiration - "SM", # The areal mean of the free water capacity of surface soil layer - "EX", # The exponent of the free water capacity curve - "KI", # Outflow coefficients of interflow - "KG", # Outflow coefficients of groundwater - "A", # parameter of mizuRoute - "THETA", # parameter of mizuRoute - "CI", # The recession constant of the lower interflow - "CG", # The recession constant of groundwater storage - "KERNEL", # kernel size of mizuRoute unit hydrograph when using convolution method - ], - "param_range": OrderedDict( - { - "K": [0.1, 1.0], - "B": [0.1, 0.4], - "IM": [0.01, 0.1], - "UM": [0.0, 20.0], - "LM": [60.0, 90.0], - "DM": [60.0, 120.0], - "C": [0.0, 0.2], - "SM": [1.0, 100.0], - "EX": [1.0, 1.5], - "KI": [0.0, 0.7], - "KG": [0.0, 0.7], - "A": [0.0, 2.9], - "THETA": [0.0, 6.5], - "CI": [0.0, 0.9], - "CG": [0.98, 0.998], - "KERNEL": [1, 15], - } - ), - }, - "gr4j": { - "param_name": ["x1", "x2", "x3", "x4"], - "param_range": OrderedDict( - { - "x1": [100.0, 1200.0], - "x2": [-5.0, 3.0], - "x3": [20.0, 300.0], - "x4": [1.1, 2.9], - } - ), - }, - "hymod": { - "param_name": ["cmax", "bexp", "alpha", "ks", "kq"], - "param_range": OrderedDict( - { - "cmax": [1.0, 500.0], - "bexp": [0.1, 2.0], - "alpha": [0.1, 0.99], - "ks": [0.001, 0.10], - "kq": [0.1, 0.99], - } - ), - }, -} diff --git a/xaj/xajmodel.py b/xaj/xajmodel.py deleted file mode 100644 index 2027955..0000000 --- a/xaj/xajmodel.py +++ /dev/null @@ -1,1060 +0,0 @@ -""" -Core code for XinAnJiang model -""" -import logging -from typing import Union -import numpy as np -from numba import jit -from scipy.special import gamma - -from xaj.params import MODEL_PARAM_DICT - -PRECISION = 1e-5 - - -@jit -def calculate_evap(lm, c, wu0, wl0, prcp, pet) -> tuple[np.array, np.array, np.array]: - """ - Three-layers evaporation model from "Watershed Hydrologic Simulation" written by Prof. RenJun Zhao. - The book is Chinese, and its name is 《流域水文模拟》; - The three-layers evaporation model is described in Page 76; - The method is same with that in Page 22-23 in "Hydrologic Forecasting (5-th version)" written by Prof. Weimin Bao. - This book's Chinese name is 《水文预报》 - - Parameters - ---------- - lm - average soil moisture storage capacity of lower layer - c - coefficient of deep layer - wu0 - initial soil moisture of upper layer; update in each time step - wl0 - initial soil moisture of lower layer; update in each time step - prcp - basin mean precipitation - pet - potential evapotranspiration - - Returns - ------- - tuple[np.array,np.array,np.array] - eu/el/ed are evaporation from upper/lower/deeper layer, respectively - """ - eu = np.where(wu0 + prcp >= pet, pet, wu0 + prcp) - ed = np.where((wl0 < c * lm) & (wl0 < c * (pet - eu)), c * (pet - eu) - wl0, 0.0) - el = np.where( - wu0 + prcp >= pet, - 0.0, - np.where( - wl0 >= c * lm, - (pet - eu) * wl0 / lm, - np.where(wl0 >= c * (pet - eu), c * (pet - eu), wl0), - ), - ) - return eu, el, ed - - -@jit -def calculate_prcp_runoff(b, im, wm, w0, pe) -> tuple[np.array, np.array]: - """ - Calculates the amount of runoff generated from rainfall after entering the underlying surface. - - Same in "Watershed Hydrologic Simulation" and "Hydrologic Forecasting (5-th version)" - - Parameters - ---------- - b - B; exponent coefficient - im - IMP; imperiousness coefficient - wm - average soil moisture storage capacity - w0 - initial soil moisture - pe - net precipitation - - Returns - ------- - tuple[np.array,np.array] - r -- runoff; r_im -- runoff of impervious part - """ - wmm = wm * (1.0 + b) - a = wmm * (1.0 - (1.0 - w0 / wm) ** (1.0 / (1.0 + b))) - if np.isnan(a).any(): - raise ArithmeticError("Please check if w0>wm or b is a negative value!") - r_cal = np.where( - pe > 0.0, - np.where( - pe + a < wmm, - # 1e-5 is a precision which we set to guarantee float's calculation is correct - pe - (wm - w0) + wm * (1.0 - np.minimum(a + pe, wmm) / wmm) ** (1.0 + b), - pe - (wm - w0), - ), - np.full(pe.shape, 0.0), - ) - r = np.maximum(r_cal, 0.0) - # separate impervious part with the other - r_im_cal = pe * im - r_im = np.maximum(r_im_cal, 0.0) - return r, r_im - - -def calculate_w_storage( - um, lm, dm, wu0, wl0, wd0, eu, el, ed, pe, r -) -> tuple[np.array, np.array, np.array]: - """ - Update the soil moisture values of the three layers. - - According to the equation 2.60 in the book《水文预报》 - - Parameters - ---------- - um - average soil moisture storage capacity of the upper layer - lm - average soil moisture storage capacity of the lower layer - dm - average soil moisture storage capacity of the deep layer - wu0 - initial values of soil moisture in upper layer - wl0 - initial values of soil moisture in lower layer - wd0 - initial values of soil moisture in deep layer - eu - evaporation of the upper layer; it isn't used in this function - el - evaporation of the lower layer - ed - evaporation of the deep layer - pe - net precipitation; it is able to be negative value in this function - r - runoff - - Returns - ------- - tuple[np.array,np.array,np.array] - wu,wl,wd -- soil moisture in upper, lower and deep layer - """ - # pe>0: the upper soil moisture was added firstly, then lower layer, and the final is deep layer - # pe<=0: no additional water, just remove evapotranspiration, - # but note the case: e >= p > 0 - # (1) if wu0 + p > e, then e = eu (2) else, wu must be zero - wu = np.where( - pe > 0.0, - np.where(wu0 + pe - r < um, wu0 + pe - r, um), - np.where(wu0 + pe > 0.0, wu0 + pe, 0.0), - ) - # calculate wd before wl because it is easier to cal using where statement - wd = np.where( - pe > 0.0, - np.where(wu0 + wl0 + pe - r > um + lm, wu0 + wl0 + wd0 + pe - r - um - lm, wd0), - wd0 - ed, - ) - # water balance (equation 2.2 in Page 13, also shown in Page 23) - # if wu0 + p > e, then e = eu; else p must be used in upper layer, - # so no matter what the case is, el didn't include p, neither ed - wl = np.where(pe > 0.0, wu0 + wl0 + wd0 + pe - r - wu - wd, wl0 - el) - # the water storage should be in reasonable range - wu_ = np.clip(wu, a_min=0.0, a_max=um) - wl_ = np.clip(wl, a_min=0.0, a_max=lm) - wd_ = np.clip(wd, a_min=0.0, a_max=dm) - return wu_, wl_, wd_ - - -def generation(p_and_e, k, b, im, um, lm, dm, c, wu0=None, wl0=None, wd0=None) -> tuple: - """ - Single-step runoff generation in XAJ. - - Parameters - ---------- - p_and_e - precipitation and potential evapotranspiration - k - ratio of potential evapotranspiration to reference crop evaporation - b - exponent parameter - um - average soil moisture storage capacity of the upper layer - lm - average soil moisture storage capacity of the lower layer - dm - average soil moisture storage capacity of the deep layer - im - impermeability coefficient - c - coefficient of deep layer - wu0 - initial values of soil moisture in upper layer - wl0 - initial values of soil moisture in lower layer - wd0 - initial values of soil moisture in deep layer - - Returns - ------- - tuple[tuple, tuple] - (r, rim, e, pe), (wu, wl, wd); all variables are np.array - """ - # make sure physical variables' value ranges are correct - prcp = np.maximum(p_and_e[:, 0], 0.0) - # get potential evapotranspiration - pet = np.maximum(p_and_e[:, 1] * k, 0.0) - # wm - wm = um + lm + dm - if wu0 is None: - # just an initial value - wu0 = 0.6 * um - if wl0 is None: - wl0 = 0.6 * lm - if wd0 is None: - wd0 = 0.6 * dm - w0_ = wu0 + wl0 + wd0 - - # w0 need locate in correct range so that following calculation could be right - # To make sure float data's calculation is correct, we'd better minus a precision (1e-5) - w0 = np.minimum(w0_, wm - 1e-5) - - # Calculate the amount of evaporation from storage - eu, el, ed = calculate_evap(lm, c, wu0, wl0, prcp, pet) - e = eu + el + ed - - # Calculate the runoff generated by net precipitation - prcp_difference = prcp - e - pe = np.maximum(prcp_difference, 0.0) - r, rim = calculate_prcp_runoff(b, im, wm, w0, pe) - # Update wu, wl, wd - wu, wl, wd = calculate_w_storage( - um, lm, dm, wu0, wl0, wd0, eu, el, ed, prcp_difference, r - ) - - return (r, rim, e, pe), (wu, wl, wd) - - -def sources(pe, r, sm, ex, ki, kg, s0=None, fr0=None, book="HF") -> tuple: - """ - Divide the runoff to different sources - - We use the initial version from the paper of the inventor of the XAJ model -- Prof. Renjun Zhao: - "Analysis of parameters of the XinAnJiang model". Its Chinese name is <<新安江模型参数的分析>>, - which could be found by searching in "Baidu Xueshu". - The module's code can also be found in "Watershed Hydrologic Simulation" (WHS) Page 174. - It is nearly same with that in "Hydrologic Forecasting" (HF) Page 148-149 - We use the period average runoff as input and the unit period is day so we don't need to difference it as books show - - We also provide code for formula from《水文预报》"Hydrologic Forecasting" (HF) the fifth version. Page 40-41 and 150-151; - the procedures in 《工程水文学》"Engineering Hydrology" (EH) the third version are different we also provide. - they are in the "sources5mm" function. - - Parameters - ------------ - pe - net precipitation - r - runoff from xaj_generation - sm - areal mean free water capacity of the surface layer - ex - exponent of the free water capacity curve - ki - outflow coefficients of the free water storage to interflow relationships - kg - outflow coefficients of the free water storage to groundwater relationships - s0 - free water capacity of last period - fr0 - runoff area of last period - - Return - ------------ - tuple[tuple, tuple] - rs -- surface runoff; ri-- interflow runoff; rg -- groundwater runoff; - s1 -- final free water capacity; - all variables are numpy array - - """ - # maximum free water storage capacity in a basin - ms = sm * (1.0 + ex) - if fr0 is None: - fr0 = 0.1 - if s0 is None: - s0 = 0.5 * sm - # For free water storage, because s is related to fr and s0 and fr0 are both values of last period, - # we have to trans the initial value of s from last period to this one. - # both WHS（流域水文模拟）'s sample code and HF（水文预报） use s = fr0 * s0 / fr. - # I think they both think free water reservoir as a cubic tank. Its height is s and area of bottom rectangle is fr - # we will have a cubic tank with varying bottom and height, - # and fixed boundary (in HF sm is fixed) or none-fixed boundary (in EH smmf is not fixed) - # but notice r's list like" [1,0] which 1 is the 1st period's runoff and 0 is the 2nd period's runoff - # after 1st period, the s1 could not be zero, but in the 2nd period, fr=0, then we cannot set s=0, because some water still in the tank - # fr's formula could be found in Eq. 9 in "Analysis of parameters of the XinAnJiang model", - # Here our r doesn't include rim, so there is no need to remove rim from r; this is also the method in 《水文预报》（HF） - - # NOTE: when r is 0, fr should be 0, however, s1 may not be zero and it still hold some water, - # then fr can not be 0, otherwise when fr is used as denominator it lead to error, - # so we have to deal with this case later, for example, when r=0, we cannot use pe * fr to replace r - # because fr get the value of last period, and it is not 0 - fr = np.copy(fr0) - if any(fr == 0.0): - raise ArithmeticError( - "Please check fr's value, fr==0.0 will cause error in the next step!" - ) - # r=0, then r/pe must be 0 - fr_mask = r > 0.0 - fr[fr_mask] = r[fr_mask] / pe[fr_mask] - if np.isnan(fr).any(): - raise ArithmeticError("Please check pe's data! there may be 0.0") - - # if fr=0, then we cannot get ss, but ss should not be 0, because s1 of last period may not be 0 and it still hold some water - ss = np.copy(s0) - # same initialization for s - s = np.copy(s0) - - ss[fr_mask] = fr0[fr_mask] * s0[fr_mask] / fr[fr_mask] - - if book == "HF": - ss = np.minimum(ss, sm) - au = ms * (1.0 - (1.0 - ss / sm) ** (1.0 / (1.0 + ex))) - if np.isnan(au).any(): - raise ValueError( - "Error: NaN values detected. Try set clip function or check your data!!!" - ) - # the first condition should be r > 0.0, when r=0, rs must be 0, fig 6-12 in EH or fig 5-4 in HF - # so we have to use "pe" carefully!!! when r>0.0, we use pe, otherwise we don't use it!!! - rs = np.full(r.shape, 0.0) - rs[fr_mask] = np.where( - pe[fr_mask] + au[fr_mask] < ms[fr_mask], - # equation 2-85 in HF - fr[fr_mask] - * ( - pe[fr_mask] - - sm[fr_mask] - + ss[fr_mask] - + sm[fr_mask] - * ( - ( - 1 - - np.minimum(pe[fr_mask] + au[fr_mask], ms[fr_mask]) - / ms[fr_mask] - ) - ** (1 + ex[fr_mask]) - ) - ), - # equation 2-86 in HF - fr[fr_mask] * (pe[fr_mask] + ss[fr_mask] - sm[fr_mask]), - ) - rs = np.minimum(rs, r) - - # ri's mask is not same as rs's, because last period's s may not be 0 - # and in this time, ri and rg could be larger than 0 - # we need firstly calculate the updated s, s's mask is same as fr_mask, - # when r==0, then s will be equal to last period's - # equation 2-87 in HF, some free water leave or save, so we update free water storage - s[fr_mask] = ss[fr_mask] + (r[fr_mask] - rs[fr_mask]) / fr[fr_mask] - s = np.minimum(s, sm) - if np.isnan(s).any(): - raise ArithmeticError("Please check fr's data! there may be 0.0") - - elif book == "EH": - # smmf should be with correpond with s - smmf = ms * (1 - (1 - fr) ** (1 / ex)) - smf = smmf / (1 + ex) - ss = np.minimum(ss, smf) - au = smmf * (1 - (1 - ss / smf) ** (1 / (1 + ex))) - if np.isnan(au).any(): - raise ValueError( - "Error: NaN values detected. Try set clip function or check your data!!!" - ) - - # rs's mask is keep with fr_mask - rs = np.full(r.shape, 0.0) - rs[fr_mask] = np.where( - pe[fr_mask] + au[fr_mask] < smmf[fr_mask], - ( - pe[fr_mask] - - smf[fr_mask] - + ss[fr_mask] - + smf[fr_mask] - * ( - 1 - - np.minimum(pe[fr_mask] + au[fr_mask], smmf[fr_mask]) - / smmf[fr_mask] - ) - ** (ex[fr_mask] + 1) - ) - * fr[fr_mask], - (pe[fr_mask] + ss[fr_mask] - smf[fr_mask]) * fr[fr_mask], - ) - rs = np.minimum(rs, r) - s[fr_mask] = ss[fr_mask] + (r[fr_mask] - rs[fr_mask]) / fr[fr_mask] - s = np.minimum(s, smf) - else: - raise ValueError("Please set book as 'HF' or 'EH'!") - # the following part is same for both HF and EH. Even the formula is different, but their meaning is same - # equation 2-88 in HF, next interflow and ground water will be released from the updated free water storage - # We use the period average runoff as input and the general unit period is day. - # Hence, we directly use ki and kg rather than ki_{Δt} in books. - ri = ki * s * fr - rg = kg * s * fr - # ri = np.where(s < PRECISION, np.full(r.shape, 0.0), ki * s * fr) - # rg = np.where(s < PRECISION, np.full(r.shape, 0.0), kg * s * fr) - # equation 2-89 in HF; although it looks different with that in WHS, they are actually same - # Finally, calculate the final free water storage - s1 = s * (1 - ki - kg) - # s1 = np.where(s1 < PRECISION, np.full(s1.shape, 0.0), s1) - return (rs, ri, rg), (s1, fr) - - -def sources5mm( - pe, - runoff, - sm, - ex, - ki, - kg, - s0=None, - fr0=None, - time_interval_hours=24, - book="HF", -): - """ - Divide the runoff to different sources according to books -- 《水文预报》HF 5th edition and 《工程水文学》EH 3rd edition - - Parameters - ---------- - pe - net precipitation - runoff - runoff from xaj_generation - sm - areal mean free water capacity of the surface layer - ex - exponent of the free water capacity curve - ki - outflow coefficients of the free water storage to interflow relationships - kg - outflow coefficients of the free water storage to groundwater relationships - s0 - initial free water capacity - fr0 - initial area of generation - time_interval_hours - 由于Ki、Kg、Ci、Cg都是以24小时为时段长定义的,需根据时段长转换 - book - the methods in 《水文预报》HF 5th edition and 《工程水文学》EH 3rd edition are different, - hence, both are provided, and the default is the former -- "ShuiWenYuBao"; - the other one is "GongChengShuiWenXue" - - Returns - ------- - tuple[tuple, tuple] - rs_s -- surface runoff; rss_s-- interflow runoff; rg_s -- groundwater runoff; - (fr_ds[-1], s_ds[-1]): state variables' final value; - all variables are numpy array - """ - # 由于Ki、Kg都是以24小时为时段长定义的，需根据时段长转换 - hours_per_day = 24 - # 非整除情况，时段+1 - residue_temp = hours_per_day % time_interval_hours - if residue_temp != 0: - residue_temp = 1 - period_num_1d = int(hours_per_day / time_interval_hours) + residue_temp - # 当kss+kg>1时，根式为偶数运算时，kss_period会成为复数，这里会报错；另外注意分母可能为0，kss不可取0 - # 对kss+kg的取值进行限制 - kss_period = (1 - (1 - (ki + kg)) ** (1 / period_num_1d)) / (1 + kg / ki) - kg_period = kss_period * kg / ki - - # 流域最大点自由水蓄水容量深 - smm = sm * (1 + ex) - if s0 is None: - s0 = 0.50 * sm - if fr0 is None: - fr0 = 0.1 - # don't use np.where here, because it will cause some warning - fr = np.copy(fr0) - fr_mask = runoff > 0.0 - fr[fr_mask] = runoff[fr_mask] / pe[fr_mask] - if np.all(runoff < 5): - n = 1 - else: - # when modeling multiple basins, the number of divides is not the same, so we use the maximum number - r_max = np.max(runoff) - residue_temp = r_max % 5 - if residue_temp != 0: - residue_temp = 1 - n = int(r_max / 5) + residue_temp - rn = runoff / n - pen = pe / n - kss_d = (1 - (1 - (kss_period + kg_period)) ** (1 / n)) / ( - 1 + kg_period / kss_period - ) - kg_d = kss_d * kg_period / kss_period - # kss_d = kss_period - # kg_d = kg_period - rs = np.full(runoff.shape, 0.0) - rss = np.full(runoff.shape, 0.0) - rg = np.full(runoff.shape, 0.0) - - s_ds = [] - fr_ds = [] - s_ds.append(s0) - fr_ds.append(fr0) - - for j in range(n): - fr0_d = fr_ds[j] - s0_d = s_ds[j] - # equation 5-32 in HF, but strange, cause each period, rn/pen is same, fr_d should be same - # fr_d = 1 - (1 - fr) ** (1 / n) - fr_d = fr - - ss_d = np.copy(s0_d) - s_d = np.copy(s0_d) - s1_d = np.copy(s0_d) - - ss_d[fr_mask] = fr0_d[fr_mask] * s0_d[fr_mask] / fr_d[fr_mask] - - if book == "HF": - ss_d = np.minimum(ss_d, sm) - # ms = smm - au = smm * (1.0 - (1.0 - ss_d / sm) ** (1.0 / (1.0 + ex))) - if np.isnan(au).any(): - raise ValueError( - "Error: NaN values detected. Try set clip function or check your data!!!" - ) - rs_j = np.full(rn.shape, 0.0) - rs_j[fr_mask] = np.where( - pen[fr_mask] + au[fr_mask] < smm[fr_mask], - # equation 5-26 in HF - fr_d[fr_mask] - * ( - pen[fr_mask] - - sm[fr_mask] - + ss_d[fr_mask] - + sm[fr_mask] - * ( - ( - 1 - - np.minimum(pen[fr_mask] + au[fr_mask], smm[fr_mask]) - / smm[fr_mask] - ) - ** (1 + ex[fr_mask]) - ) - ), - # equation 5-27 in HF - fr_d[fr_mask] * (pen[fr_mask] + ss_d[fr_mask] - sm[fr_mask]), - ) - rs_j = np.minimum(rs_j, rn) - s_d[fr_mask] = ss_d[fr_mask] + (rn[fr_mask] - rs_j[fr_mask]) / fr_d[fr_mask] - s_d = np.minimum(s_d, sm) - - elif book == "EH": - smmf = smm * (1 - (1 - fr_d) ** (1 / ex)) - smf = smmf / (1 + ex) - ss_d = np.minimum(ss_d, smf) - au = smmf * (1 - (1 - ss_d / smf) ** (1 / (1 + ex))) - if np.isnan(au).any(): - raise ValueError( - "Error: NaN values detected. Try set clip function or check your data!!!" - ) - rs_j = np.full(rn.shape, 0.0) - rs_j[fr_mask] = np.where( - pen[fr_mask] + au[fr_mask] < smmf[fr_mask], - ( - pen[fr_mask] - - smf[fr_mask] - + ss_d[fr_mask] - + smf[fr_mask] - * ( - 1 - - np.minimum(pen[fr_mask] + au[fr_mask], smmf[fr_mask]) - / smmf[fr_mask] - ) - ** (ex[fr_mask] + 1) - ) - * fr_d[fr_mask], - (pen[fr_mask] + ss_d[fr_mask] - smf[fr_mask]) * fr_d[fr_mask], - ) - rs_j = np.minimum(rs_j, rn) - s_d[fr_mask] = ss_d[fr_mask] + (rn[fr_mask] - rs_j[fr_mask]) / fr_d[fr_mask] - s_d = np.minimum(s_d, smf) - - else: - raise NotImplementedError( - "We don't have this implementation! Please chose 'HF' or 'EH'!!" - ) - rss_j = s_d * kss_d * fr_d - rg_j = s_d * kg_d * fr_d - s1_d = s_d * (1 - kss_d - kg_d) - - rs = rs + rs_j - rss = rss + rss_j - rg = rg + rg_j - # 赋值s_d和fr_d到数组中，以给下一段做初值 - s_ds.append(s1_d) - fr_ds.append(fr_d) - - return (rs, rss, rg), (s_ds[-1], fr_ds[-1]) - - -@jit -def linear_reservoir(x, weight, last_y=None) -> np.array: - """ - Linear reservoir's release function - - Parameters - ---------- - x - the input to the linear reservoir - weight - the coefficient of linear reservoir - last_y - the output of last period - - Returns - ------- - np.array - one-step forward result - """ - weight1 = 1 - weight - if last_y is None: - last_y = np.full(weight.shape, 0.001) - y = weight * last_y + weight1 * x - return y - - -def uh_conv(x, uh_from_gamma): - """ - Function for 1d-convolution calculation - - Parameters - ---------- - x - x is a sequence-first variable; the dim of x is [seq, batch, feature=1]; - feature must be 1 - uh_from_gamma - unit hydrograph from uh_gamma; the dim: [len_uh, batch, feature=1]; - feature must be 1 - - Returns - ------- - np.array - convolution - """ - outputs = np.full(x.shape, 0.0) - time_length, batch_size, feature_size = x.shape - if feature_size > 1: - logging.error("We only support one-dim convolution now!!!") - for i in range(batch_size): - uh = uh_from_gamma[:, i, 0] - inputs = x[:, i, 0] - outputs[:, i, 0] = np.convolve(inputs, uh)[:time_length] - return outputs - - -def uh_gamma(a, theta, len_uh=15): - """ - A simple two-parameter Gamma distribution as a unit-hydrograph to route instantaneous runoff from a hydrologic model - The method comes from mizuRoute -- http://www.geosci-model-dev.net/9/2223/2016/ - - Parameters - ---------- - a - shape parameter - theta - timescale parameter - len_uh - the time length of the unit hydrograph - Returns - ------- - torch.Tensor - the unit hydrograph, dim: [seq, batch, feature] - """ - # dims of a: time_seq (same all time steps), batch, feature=1 - m = a.shape - if len_uh > m[0]: - raise RuntimeError( - "length of unit hydrograph should be smaller than the whole length of input" - ) - # aa > 0, here we set minimum 0.1 (min of a is 0, set when calling this func); First dimension of a is repeat - aa = np.maximum(0.0, a[0:len_uh, :, :]) + 0.1 - # theta > 0, here set minimum 0.5 - theta = np.maximum(0.0, theta[0:len_uh, :, :]) + 0.5 - # len_f, batch, feature - t = np.expand_dims( - np.swapaxes(np.tile(np.arange(0.5, len_uh * 1.0), (m[1], 1)), 0, 1), axis=-1 - ) - denominator = gamma(aa) * (theta**aa) - # [len_f, m[1], m[2]] - w = 1 / denominator * (t ** (aa - 1)) * (np.exp(-t / theta)) - w = w / w.sum(0) # scale to 1 for each UH - return w - - -def xaj_runoff( - p_and_e, - params_runoff: Union[np.array, list], - *w0, s0, fr0, - **kwargs, -) -> Union[tuple, np.array]: - model_name = kwargs["name"] if "name" in kwargs else "xaj" - source_type = kwargs["source_type"] if "source_type" in kwargs else "sources" - source_book = kwargs["source_book"] if "source_book" in kwargs else "HF" - # 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 'CS' or 'MZ'!" - ) - xaj_params = [ - (value[1] - value[0]) * params_runoff[:, 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_) - - - # w0 = (0.5 * um, 0.5 * lm, 0.5 * dm) - inputs = p_and_e[:, :, :] - 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 - ) - 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) - 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) - return runoff_im, rss_,ris_, rgs_, es,rss - -def xaj_route( - p_and_e, - params_route: Union[np.array, list], - ris_, - rgs_, - rss_, - rss, - runoff_im, - qi0, - qg0, - es, - **kwargs, -)-> Union[tuple, np.array]: - - model_name = kwargs["name"] if "name" in kwargs else "xaj" - source_type = kwargs["source_type"] if "source_type" in kwargs else "sources" - source_book = kwargs["source_book"] if "source_book" in kwargs else "HF" - # 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 'CS' or 'MZ'!" - ) - xaj_params = [ - (value[1] - value[0]) * params_route[:, 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 'CS' or 'MZ'!" - ) - ci = xaj_params[13] - cg = xaj_params[14] - - inputs = p_and_e[:, :, :] - qs = np.full(inputs.shape[:2], 0.0) - 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(len(l)): - lag = int(l[j]) - for i in range(0,lag): - qs[i, j] = qt[i, j] - if i == lag: - break - 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 'CSL' or 'MZ'!" - ) - - # seq, batch, feature - q_sim = np.expand_dims(qs, axis=2) - return q_sim, es - -def xaj_state( - p_and_e, - params: Union[np.array, list], - return_state=True, - warmup_length=365, - **kwargs, -) -> Union[tuple, np.array]: - # default values for some function parameters - model_name = kwargs["name"] if "name" in kwargs else "xaj" - source_type = kwargs["source_type"] if "source_type" in kwargs else "sources" - source_book = kwargs["source_book"] if "source_book" in kwargs else "HF" - # 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 'CS' or 'MZ'!" - ) - 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 'CS' 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_state( - # 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:, :, :] - inputs = p_and_e[:, :, :] - 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 - ) - 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) - 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": - 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)): - # print(range(len(l))) - lag = int(l[j]) - for i in range(0,lag): - qs[i, j] = qt[i, j] - if i == lag: - break - for i in range(lag, inputs.shape[0]): - # print(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 'CSL' 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 From c222649b076e92dcde2492bf9a0b15aaf1620f49 Mon Sep 17 00:00:00 2001 From: ouyangwenyu Date: Sat, 28 Oct 2023 16:21:53 +0800 Subject: [PATCH 3/5] refactor xajbim-related funcs, but not fully tested yet --- hydromodel/models/xaj.py | 213 ++++++++++++++++++++++++++++++++-- hydromodel/models/xaj_bmi.py | 3 +- test/test_rr_event_iden.py | 96 +++++++++++++++ test/test_session_division.py | 44 ------- test/test_xaj_bmi.py | 203 +++++++++++++++++--------------- 5 files changed, 413 insertions(+), 146 deletions(-) create mode 100644 test/test_rr_event_iden.py delete mode 100644 test/test_session_division.py diff --git a/hydromodel/models/xaj.py b/hydromodel/models/xaj.py index 8aaf9ba..b780d9b 100644 --- a/hydromodel/models/xaj.py +++ b/hydromodel/models/xaj.py @@ -8,7 +8,6 @@ from scipy.special import gamma from hydromodel.models.model_config import MODEL_PARAM_DICT -from xaj.params import MODEL_PARAM_DICT PRECISION = 1e-5 @@ -888,6 +887,203 @@ def xaj( return q_sim, es +def xaj_runoff( + p_and_e, + params_runoff: Union[np.array, list], + *w0, + s0, + fr0, + **kwargs, +) -> Union[tuple, np.array]: + model_name = kwargs.get("name", "xaj") + source_type = kwargs.get("source_type", "sources") + source_book = kwargs.get("source_book", "HF") + # 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 'CS' or 'MZ'!" + ) + xaj_params = [ + (value[1] - value[0]) * params_runoff[:, 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_) + + # w0 = (0.5 * um, 0.5 * lm, 0.5 * dm) + inputs = p_and_e[:, :, :] + 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 + ) + 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) + 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) + return runoff_im, rss_, ris_, rgs_, es, rss + + +def xaj_route( + p_and_e, + params_route: Union[np.array, list], + ris_, + rgs_, + rss_, + rss, + runoff_im, + qi0, + qg0, + es, + **kwargs, +) -> Union[tuple, np.array]: + model_name = kwargs.get("name", "xaj") + source_type = kwargs.get("source_type", "sources") + source_book = kwargs.get("source_book", "HF") + # 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 'CS' or 'MZ'!" + ) + xaj_params = [ + (value[1] - value[0]) * params_route[:, 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 'CS' or 'MZ'!" + ) + ci = xaj_params[13] + cg = xaj_params[14] + + inputs = p_and_e[:, :, :] + qs = np.full(inputs.shape[:2], 0.0) + 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(len(l)): + lag = int(l[j]) + for i in range(lag): + qs[i, j] = qt[i, j] + if i == lag: + break + 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 'CSL' or 'MZ'!" + ) + + # seq, batch, feature + q_sim = np.expand_dims(qs, axis=2) + return q_sim, es + + def xaj_state( p_and_e, params: Union[np.array, list], @@ -944,7 +1140,6 @@ def xaj_state( 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, :, :] @@ -962,8 +1157,8 @@ def xaj_state( qi0 = np.full(ci.shape, 0.1) qg0 = np.full(cg.shape, 0.1) - #state_variables - # inputs = p_and_e[warmup_length:, :, :] + # state_variables + # inputs = p_and_e[warmup_length:, :, :] inputs = p_and_e[:, :, :] runoff_ims_ = np.full(inputs.shape[:2], 0.0) rss_ = np.full(inputs.shape[:2], 0.0) @@ -1028,7 +1223,7 @@ def xaj_state( for j in range(len(l)): # print(range(len(l))) lag = int(l[j]) - for i in range(0,lag): + for i in range(lag): qs[i, j] = qt[i, j] if i == lag: break @@ -1049,10 +1244,10 @@ def xaj_state( 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 'CSL' or 'MZ'!" - ) + else: + raise NotImplementedError( + "We don't provide this route method now! Please use 'CSL' or 'MZ'!" + ) # seq, batch, feature q_sim = np.expand_dims(qs, axis=2) diff --git a/hydromodel/models/xaj_bmi.py b/hydromodel/models/xaj_bmi.py index 94f1d02..435ab67 100644 --- a/hydromodel/models/xaj_bmi.py +++ b/hydromodel/models/xaj_bmi.py @@ -3,10 +3,11 @@ import numpy as np from hydromodel.models import configuration -from hydromodel.models.xajmodel import xaj_route, xaj_runoff import datetime import logging +from hydromodel.models.xaj import xaj_route, xaj_runoff + logger = logging.getLogger(__name__) PRECISION = 1e-5 diff --git a/test/test_rr_event_iden.py b/test/test_rr_event_iden.py new file mode 100644 index 0000000..0034e38 --- /dev/null +++ b/test/test_rr_event_iden.py @@ -0,0 +1,96 @@ +""" +Author: Wenyu Ouyang +Date: 2023-10-28 09:23:22 +LastEditTime: 2023-10-28 16:19:22 +LastEditors: Wenyu Ouyang +Description: Test for rainfall-runoff event identification +FilePath: \hydro-model-xaj\test\test_rr_event_iden.py +Copyright (c) 2023-2024 Wenyu Ouyang. All rights reserved. +""" +import os +import pandas as pd +import definitions +from hydromodel.utils.dmca_esr import ( + step1_step2_tr_and_fluctuations_timeseries, + step3_core_identification, + step4_end_rain_events, + step5_beginning_rain_events, + step6_checks_on_rain_events, + step7_end_flow_events, + step8_beginning_flow_events, + step9_checks_on_flow_events, + step10_checks_on_overlapping_events, +) + + +def test_rainfall_runoff_event_identify(): + rain = pd.read_csv( + os.path.join(definitions.ROOT_DIR, "hydromodel/example/division_rain.csv") + )["rain"] + flow = pd.read_csv( + os.path.join(definitions.ROOT_DIR, "hydromodel/example/division_flow.csv") + )["flow"] + time = rain.index.to_numpy() + rain = rain.to_numpy() / 24 + flow = flow.to_numpy() / 24 + rain_min = 0.02 + max_window = 100 + ( + Tr, + fluct_rain_Tr, + fluct_flow_Tr, + fluct_bivariate_Tr, + ) = step1_step2_tr_and_fluctuations_timeseries(rain, flow, rain_min, max_window) + beginning_core, end_core = step3_core_identification(fluct_bivariate_Tr) + end_rain = step4_end_rain_events( + beginning_core, end_core, rain, fluct_rain_Tr, rain_min + ) + beginning_rain = step5_beginning_rain_events( + beginning_core, end_rain, rain, fluct_rain_Tr, rain_min + ) + ( + beginning_rain_checked, + end_rain_checked, + beginning_core, + end_core, + ) = step6_checks_on_rain_events( + beginning_rain, end_rain, rain, rain_min, beginning_core, end_core + ) + end_flow = step7_end_flow_events( + end_rain_checked, + beginning_core, + end_core, + rain, + fluct_rain_Tr, + fluct_flow_Tr, + Tr, + ) + beginning_flow = step8_beginning_flow_events( + beginning_rain_checked, + end_rain_checked, + rain, + beginning_core, + fluct_rain_Tr, + fluct_flow_Tr, + ) + beginning_flow_checked, end_flow_checked = step9_checks_on_flow_events( + beginning_rain_checked, + end_rain_checked, + beginning_flow, + end_flow, + fluct_flow_Tr, + ) + ( + BEGINNING_RAIN, + END_RAIN, + BEGINNING_FLOW, + END_FLOW, + ) = step10_checks_on_overlapping_events( + beginning_rain_checked, + end_rain_checked, + beginning_flow_checked, + end_flow_checked, + time, + ) + print(BEGINNING_RAIN, END_RAIN, BEGINNING_FLOW, END_FLOW) + print(len(BEGINNING_RAIN), len(END_RAIN), len(BEGINNING_FLOW), len(END_FLOW)) diff --git a/test/test_session_division.py b/test/test_session_division.py deleted file mode 100644 index cc00831..0000000 --- a/test/test_session_division.py +++ /dev/null @@ -1,44 +0,0 @@ -import os - -import pandas as pd - -import definitions -from hydromodel.utils.dmca_esr import step1_step2_tr_and_fluctuations_timeseries, step3_core_identification, \ - step4_end_rain_events, \ - step5_beginning_rain_events, step6_checks_on_rain_events, step7_end_flow_events, step8_beginning_flow_events, \ - step9_checks_on_flow_events, step10_checks_on_overlapping_events - - -def test_session_division_new(): - rain = pd.read_csv(os.path.join(definitions.ROOT_DIR, 'hydromodel/example/division_rain.csv'))['rain'] - flow = pd.read_csv(os.path.join(definitions.ROOT_DIR, 'hydromodel/example/division_flow.csv'))['flow'] - time = rain.index.to_numpy() - rain = rain.to_numpy() / 24 - flow = flow.to_numpy() / 24 - rain_min = 0.02 - max_window = 100 - Tr, fluct_rain_Tr, fluct_flow_Tr, fluct_bivariate_Tr = step1_step2_tr_and_fluctuations_timeseries(rain, flow, - rain_min, - max_window) - beginning_core, end_core = step3_core_identification(fluct_bivariate_Tr) - end_rain = step4_end_rain_events(beginning_core, end_core, rain, fluct_rain_Tr, rain_min) - beginning_rain = step5_beginning_rain_events(beginning_core, end_rain, rain, fluct_rain_Tr, rain_min) - beginning_rain_checked, end_rain_checked, beginning_core, end_core = step6_checks_on_rain_events(beginning_rain, - end_rain, rain, - rain_min, - beginning_core, - end_core) - end_flow = step7_end_flow_events(end_rain_checked, beginning_core, end_core, rain, fluct_rain_Tr, fluct_flow_Tr, Tr) - beginning_flow = step8_beginning_flow_events(beginning_rain_checked, end_rain_checked, rain, beginning_core, - fluct_rain_Tr, - fluct_flow_Tr) - beginning_flow_checked, end_flow_checked = step9_checks_on_flow_events(beginning_rain_checked, end_rain_checked, - beginning_flow, - end_flow, fluct_flow_Tr) - BEGINNING_RAIN, END_RAIN, BEGINNING_FLOW, END_FLOW = step10_checks_on_overlapping_events(beginning_rain_checked, - end_rain_checked, - beginning_flow_checked, - end_flow_checked, - time) - print(BEGINNING_RAIN, END_RAIN, BEGINNING_FLOW, END_FLOW) - print(len(BEGINNING_RAIN), len(END_RAIN), len(BEGINNING_FLOW), len(END_FLOW)) diff --git a/test/test_xaj_bmi.py b/test/test_xaj_bmi.py index c61808c..0b65331 100644 --- a/test/test_xaj_bmi.py +++ b/test/test_xaj_bmi.py @@ -3,9 +3,6 @@ import definitions from hydromodel.models.configuration import read_config from hydromodel.models.xaj_bmi import xajBmi - -logging.basicConfig(level=logging.INFO) - import pandas as pd import os from pathlib import Path @@ -19,8 +16,8 @@ cross_valid_data, split_train_test, ) -from xaj.calibrate_sceua_xaj_bmi import calibrate_by_sceua -from xaj.calibrate_ga_xaj_bmi import ( +from hydromodel.calibrate.calibrate_sceua_xaj_bmi import calibrate_by_sceua +from hydromodel.calibrate.calibrate_ga_xaj_bmi import ( calibrate_by_ga, show_ga_result, ) @@ -34,9 +31,11 @@ ) from hydromodel.utils import hydro_constant +logging.basicConfig(level=logging.INFO) + def test_bmi(): - ''' + """ model = xajBmi() print(model.get_component_name()) model.initialize("runxaj.yaml") @@ -64,23 +63,23 @@ def test_bmi(): }) results.to_csv('/home/wangjingyi/code/hydro-model-xaj/scripts/xaj.csv') model.finalize() - ''' + """ # 模型率定 config = read_config(os.path.relpath("runxaj.yaml")) - forcing_data = Path(str(definitions.ROOT_DIR) + str(config['forcing_data'])) - train_period = config['train_period'] - test_period = config['test_period'] + forcing_data = Path(str(definitions.ROOT_DIR) + str(config["forcing_data"])) + train_period = config["train_period"] + test_period = config["test_period"] # period = config['period'] - json_file = Path(str(definitions.ROOT_DIR) + str(config['json_file'])) - npy_file = Path(str(definitions.ROOT_DIR) + str(config['npy_file'])) - cv_fold = config['cv_fold'] - warmup_length = config['warmup_length'] + json_file = Path(str(definitions.ROOT_DIR) + str(config["json_file"])) + npy_file = Path(str(definitions.ROOT_DIR) + str(config["npy_file"])) + cv_fold = config["cv_fold"] + warmup_length = config["warmup_length"] # model_info - model_name = config['model_name'] - source_type = config['source_type'] - source_book = config['source_book'] + model_name = config["model_name"] + source_type = config["source_type"] + source_book = config["source_book"] # algorithm - algorithm_name = config['algorithm_name'] + algorithm_name = config["algorithm_name"] if not (cv_fold > 1): # no cross validation @@ -93,7 +92,7 @@ def test_bmi(): split_train_test(json_file, npy_file, train_period, test_period) kfold = [ - int(f_name[len("data_info_fold"): -len("_test.json")]) + int(f_name[len("data_info_fold") : -len("_test.json")]) for f_name in os.listdir(os.path.dirname(forcing_data)) if fnmatch.fnmatch(f_name, "*_fold*_test.json") ] @@ -104,19 +103,21 @@ def test_bmi(): os.path.dirname(forcing_data), f"data_info_fold{str(fold)}_train.json" ) train_data_file = os.path.join( - os.path.dirname(forcing_data), f"basins_lump_p_pe_q_fold{str(fold)}_train.npy" + os.path.dirname(forcing_data), + f"basins_lump_p_pe_q_fold{str(fold)}_train.npy", ) test_data_info_file = os.path.join( os.path.dirname(forcing_data), f"data_info_fold{str(fold)}_test.json" ) test_data_file = os.path.join( - os.path.dirname(forcing_data), f"basins_lump_p_pe_q_fold{str(fold)}_test.npy" + os.path.dirname(forcing_data), + f"basins_lump_p_pe_q_fold{str(fold)}_test.npy", ) if ( - os.path.exists(train_data_info_file) is False - or os.path.exists(train_data_file) is False - or os.path.exists(test_data_info_file) is False - or os.path.exists(test_data_file) is False + os.path.exists(train_data_info_file) is False + or os.path.exists(train_data_file) is False + or os.path.exists(test_data_info_file) is False + or os.path.exists(test_data_file) is False ): raise FileNotFoundError( "The data files are not found, please run datapreprocess4calibrate.py first." @@ -127,14 +128,17 @@ def test_bmi(): data_info_test = hydro_utils.unserialize_json_ordered(test_data_info_file) current_time = datetime.now().strftime("%b%d_%H-%M-%S") # one directory for one model + one hyperparam setting and one basin - save_dir = os.path.join(os.path.dirname(forcing_data), current_time + "_" + socket.gethostname() + "_fold" + str(fold)) + save_dir = os.path.join( + os.path.dirname(forcing_data), + current_time + "_" + socket.gethostname() + "_fold" + str(fold), + ) if algorithm_name == "SCE_UA": - random_seed = config['random_seed'] - rep = config['rep'] - ngs = config['ngs'] - kstop = config['kstop'] - peps = config['peps'] - pcento = config['pcento'] + random_seed = config["random_seed"] + rep = config["rep"] + ngs = config["ngs"] + kstop = config["kstop"] + peps = config["peps"] + pcento = config["pcento"] for i in range(len(data_info_train["basin"])): basin_id = data_info_train["basin"][i] basin_area = data_info_train["area"][i] @@ -148,23 +152,23 @@ def test_bmi(): os.makedirs(spotpy_db_dir) db_name = os.path.join(spotpy_db_dir, "SCEUA_" + model_name) sampler = calibrate_by_sceua( - data_train[:, i: i + 1, 0:2], - data_train[:, i: i + 1, -1:], + data_train[:, i : i + 1, 0:2], + data_train[:, i : i + 1, -1:], db_name, warmup_length=warmup_length, model={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book + "name": model_name, + "source_type": source_type, + "source_book": source_book, }, algorithm={ - 'name': algorithm_name, - 'random_seed': random_seed, - 'rep': rep, - 'ngs': ngs, - 'kstop': kstop, - 'peps': peps, - 'pcento': pcento + "name": algorithm_name, + "random_seed": random_seed, + "rep": rep, + "ngs": ngs, + "kstop": kstop, + "peps": peps, + "pcento": pcento, }, ) @@ -182,7 +186,9 @@ def test_bmi(): ) model = xajBmi() - model.initialize(os.path.relpath("runxaj.yaml"), params, data_test[:, i: i + 1, 0:2]) + model.initialize( + os.path.relpath("runxaj.yaml"), params, data_test[:, i : i + 1, 0:2] + ) j = 0 while j <= len(data_info_test["time"]): model.update() @@ -197,7 +203,7 @@ def test_bmi(): ) qobs = hydro_constant.convert_unit( - data_test[warmup_length:, i: i + 1, -1:], + data_test[warmup_length:, i : i + 1, -1:], unit_now="mm/day", unit_final=hydro_constant.unit["streamflow"], basin_area=basin_area, @@ -219,12 +225,12 @@ def test_bmi(): basin_id, test_date, qsim, qobs, save_dir=spotpy_db_dir ) elif algorithm_name == "GA": - random_seed = config['random_seed'] - run_counts = config['run_counts'] - pop_num = config['pop_num'] - cross_prob = config['cross_prob'] - mut_prob = config['mut_prob'] - save_freq = config['save_freq'] + random_seed = config["random_seed"] + run_counts = config["run_counts"] + pop_num = config["pop_num"] + cross_prob = config["cross_prob"] + mut_prob = config["mut_prob"] + save_freq = config["save_freq"] for i in range(len(data_info_train["basin"])): basin_id = data_info_train["basin"][i] basin_area = data_info_train["area"][i] @@ -234,36 +240,36 @@ def test_bmi(): if not os.path.exists(deap_db_dir): os.makedirs(deap_db_dir) calibrate_by_ga( - data_train[:, i: i + 1, 0:2], - data_train[:, i: i + 1, -1:], + data_train[:, i : i + 1, 0:2], + data_train[:, i : i + 1, -1:], deap_db_dir, warmup_length=warmup_length, model={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book + "name": model_name, + "source_type": source_type, + "source_book": source_book, }, ga_param={ - 'name': algorithm_name, - 'random_seed': random_seed, - 'run_counts': run_counts, - 'pop_num': pop_num, - 'cross_prob': cross_prob, - 'mut_prob': mut_prob, - 'save_freq': save_freq + "name": algorithm_name, + "random_seed": random_seed, + "run_counts": run_counts, + "pop_num": pop_num, + "cross_prob": cross_prob, + "mut_prob": mut_prob, + "save_freq": save_freq, }, ) show_ga_result( deap_db_dir, warmup_length=warmup_length, basin_id=basin_id, - the_data=data_train[:, i: i + 1, :], + the_data=data_train[:, i : i + 1, :], the_period=data_info_train["time"], basin_area=basin_area, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book + "name": model_name, + "source_type": source_type, + "source_book": source_book, }, train_mode=True, ) @@ -271,13 +277,13 @@ def test_bmi(): deap_db_dir, warmup_length=warmup_length, basin_id=basin_id, - the_data=data_test[:, i: i + 1, :], + the_data=data_test[:, i : i + 1, :], the_period=data_info_test["time"], basin_area=basin_area, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book + "name": model_name, + "source_type": source_type, + "source_book": source_book, }, train_mode=False, ) @@ -285,24 +291,37 @@ def test_bmi(): raise NotImplementedError( "We don't provide this calibrate method! Choose from 'SCE_UA' or 'GA'!" ) - summarize_parameters(save_dir, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book - }) - renormalize_params(save_dir, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book - }) - summarize_metrics(save_dir, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book - }) - save_streamflow(save_dir, model_info={ - 'name': model_name, - 'source_type': source_type, - 'source_book': source_book, - }, fold=fold) + summarize_parameters( + save_dir, + model_info={ + "name": model_name, + "source_type": source_type, + "source_book": source_book, + }, + ) + renormalize_params( + save_dir, + model_info={ + "name": model_name, + "source_type": source_type, + "source_book": source_book, + }, + ) + summarize_metrics( + save_dir, + model_info={ + "name": model_name, + "source_type": source_type, + "source_book": source_book, + }, + ) + save_streamflow( + save_dir, + model_info={ + "name": model_name, + "source_type": source_type, + "source_book": source_book, + }, + fold=fold, + ) print(f"Finish calibrating the {fold}-th fold") From 51f650d6f703a1c6ae8e2ff95e5882162e3caf91 Mon Sep 17 00:00:00 2001 From: ouyangwenyu Date: Sat, 28 Oct 2023 16:43:03 +0800 Subject: [PATCH 4/5] delete some unused lines in xaj --- hydromodel/models/xaj.py | 35 ----------------------------------- hydromodel/models/xaj_bmi.py | 10 +++++----- test/test_xaj_bmi.py | 5 ++--- 3 files changed, 7 insertions(+), 43 deletions(-) diff --git a/hydromodel/models/xaj.py b/hydromodel/models/xaj.py index b780d9b..665fb4b 100644 --- a/hydromodel/models/xaj.py +++ b/hydromodel/models/xaj.py @@ -900,14 +900,6 @@ def xaj_runoff( source_book = kwargs.get("source_book", "HF") # 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 'CS' or 'MZ'!" - ) xaj_params = [ (value[1] - value[0]) * params_runoff[:, i] + value[0] for i, (key, value) in enumerate(param_ranges.items()) @@ -991,8 +983,6 @@ def xaj_route( **kwargs, ) -> Union[tuple, np.array]: model_name = kwargs.get("name", "xaj") - source_type = kwargs.get("source_type", "sources") - source_book = kwargs.get("source_book", "HF") # params param_ranges = MODEL_PARAM_DICT[model_name]["param_range"] if model_name == "xaj": @@ -1007,20 +997,6 @@ def xaj_route( (value[1] - value[0]) * params_route[:, 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] @@ -1140,17 +1116,6 @@ def xaj_state( 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_state( - # 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) diff --git a/hydromodel/models/xaj_bmi.py b/hydromodel/models/xaj_bmi.py index 435ab67..e94f620 100644 --- a/hydromodel/models/xaj_bmi.py +++ b/hydromodel/models/xaj_bmi.py @@ -76,7 +76,7 @@ def initialize(self, config_file, params, p_and_e): self.test_start_time_str = test_period[0] self.test_end_time_str = test_period[1] - except: + except Exception: import traceback traceback.print_exc() @@ -294,8 +294,8 @@ def start_Time(self, _start_time_str): self._startTime = datetime.datetime( int(year), int(month), int(day), int(hour), int(minute), int(second) ) - except ValueError: - raise ValueError("Invalid start date format!") + except ValueError as e: + raise ValueError("Invalid start date format!") from e return self._startTime @@ -314,6 +314,6 @@ def end_Time(self, _end_time_str): self._endTime = datetime.datetime( int(year), int(month), int(day), int(hour), int(minute), int(second) ) - except ValueError: - raise ValueError("Invalid start date format!") + except ValueError as e: + raise ValueError("Invalid start date format!") from e return self._endTime diff --git a/test/test_xaj_bmi.py b/test/test_xaj_bmi.py index 0b65331..875a9d7 100644 --- a/test/test_xaj_bmi.py +++ b/test/test_xaj_bmi.py @@ -81,14 +81,13 @@ def test_bmi(): # algorithm algorithm_name = config["algorithm_name"] - if not (cv_fold > 1): + if cv_fold <= 1: # no cross validation periods = np.sort( [train_period[0], train_period[1], test_period[0], test_period[1]] ) - if cv_fold > 1: - cross_valid_data(json_file, npy_file, periods, warmup_length, cv_fold) else: + cross_valid_data(json_file, npy_file, periods, warmup_length, cv_fold) split_train_test(json_file, npy_file, train_period, test_period) kfold = [ From 3caac4a9e2aa14f3db20c9c24264ebd84fcf59e0 Mon Sep 17 00:00:00 2001 From: ouyangwenyu Date: Sat, 28 Oct 2023 17:41:05 +0800 Subject: [PATCH 5/5] add files from cookiecutter template --- .editorconfig | 21 +++++++ .github/ISSUE_TEMPLATE/bug_report.md | 25 ++++++++ .github/ISSUE_TEMPLATE/config.yml | 10 ++++ .github/ISSUE_TEMPLATE/feature_request.md | 18 ++++++ .github/workflows/docs-build.yml | 51 ++++++++++++++++ .github/workflows/docs.yml | 28 +++++++++ .github/workflows/macos.yml | 38 ++++++++++++ .github/workflows/publish-to-pypi.yml | 38 ------------ .github/workflows/pypi.yml | 30 ++++++++++ .github/workflows/ubuntu.yml | 48 +++++++++++++++ .github/workflows/windows.yml | 33 +++++++++++ AUTHORS.rst | 13 +++++ MANIFEST.in | 8 ++- environment-dev.yml | 4 ++ requirements.txt | 6 +- requirements_dev.txt | 14 +++++ setup.cfg | 26 +++++++++ setup.py | 71 ++++++++++++++++++----- 18 files changed, 424 insertions(+), 58 deletions(-) create mode 100644 .editorconfig create mode 100644 .github/ISSUE_TEMPLATE/bug_report.md create mode 100644 .github/ISSUE_TEMPLATE/config.yml create mode 100644 .github/ISSUE_TEMPLATE/feature_request.md create mode 100644 .github/workflows/docs-build.yml create mode 100644 .github/workflows/docs.yml create mode 100644 .github/workflows/macos.yml delete mode 100644 .github/workflows/publish-to-pypi.yml create mode 100644 .github/workflows/pypi.yml create mode 100644 .github/workflows/ubuntu.yml create mode 100644 .github/workflows/windows.yml create mode 100644 AUTHORS.rst create mode 100644 requirements_dev.txt create mode 100644 setup.cfg diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 0000000..d4a2c44 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,21 @@ +# http://editorconfig.org + +root = true + +[*] +indent_style = space +indent_size = 4 +trim_trailing_whitespace = true +insert_final_newline = true +charset = utf-8 +end_of_line = lf + +[*.bat] +indent_style = tab +end_of_line = crlf + +[LICENSE] +insert_final_newline = false + +[Makefile] +indent_style = tab diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md new file mode 100644 index 0000000..c89fc62 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -0,0 +1,25 @@ +--- +name: Bug Report +about: Create a bug report to help us improve +labels: bug +--- + + + +### Environment Information + +- hydromodel version: +- Python version: +- Operating System: + +### Description + +Describe what you were trying to get done. +Tell us what happened, what went wrong, and what you expected to happen. + +### What I Did + +``` +Paste the command(s) you ran and the output. +If there was a crash, please include the traceback here. +``` diff --git a/.github/ISSUE_TEMPLATE/config.yml b/.github/ISSUE_TEMPLATE/config.yml new file mode 100644 index 0000000..de46af5 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/config.yml @@ -0,0 +1,10 @@ +contact_links: + - name: Ask questions + url: https://github.com/OuyangWenyu/hydro-model-xaj/discussions/categories/q-a + about: Please ask and answer questions here. + - name: Ideas + url: https://github.com/OuyangWenyu/hydro-model-xaj/discussions/categories/ideas + about: Please share your ideas here. + - name: Ask questions from the GIS community + url: https://gis.stackexchange.com + about: To get answers from questions in the GIS community, please ask and answer questions here. diff --git a/.github/ISSUE_TEMPLATE/feature_request.md b/.github/ISSUE_TEMPLATE/feature_request.md new file mode 100644 index 0000000..cf7acb9 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/feature_request.md @@ -0,0 +1,18 @@ +--- +name: Feature Request +about: Submit a feature request to help us improve +labels: Feature Request +--- + + + +### Description + +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. + +### Source code + +``` +Paste your source code here if have sample code to share. +``` diff --git a/.github/workflows/docs-build.yml b/.github/workflows/docs-build.yml new file mode 100644 index 0000000..876561c --- /dev/null +++ b/.github/workflows/docs-build.yml @@ -0,0 +1,51 @@ +name: docs-build +on: + pull_request: + branches: + - master + +jobs: + deploy: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + with: + fetch-depth: 0 + - uses: actions/setup-python@v4 + with: + python-version: "3.10" + - name: Install GDAL + run: | + python -m pip install --upgrade pip + pip install --find-links=https://girder.github.io/large_image_wheels --no-cache GDAL pyproj + - name: Test GDAL installation + run: | + python -c "from osgeo import gdal" + gdalinfo --version + - name: Install dependencies + run: | + pip install --no-cache-dir Cython + pip install -r requirements.txt -r requirements_dev.txt -r requirements_docs.txt + pip install . + - name: Discover typos with codespell + run: codespell --skip="*.csv,*.geojson,*.json,*.js,*.html,*cff,*.pdf,./.git" --ignore-words-list="aci,acount,acounts,fallow,hart,hist,nd,ned,ois,wqs" + - name: PKG-TEST + run: | + python -m unittest discover tests/ + - name: Build docs + run: | + mkdocs build + # - name: Deploy to Netlify + # uses: nwtgck/actions-netlify@v2.0 + # with: + # publish-dir: "./site" + # production-branch: master + # github-token: ${{ secrets.GITHUB_TOKEN }} + # deploy-message: "Deploy from GitHub Actions" + # enable-pull-request-comment: true + # enable-commit-comment: false + # overwrites-pull-request-comment: true + # env: + # NETLIFY_AUTH_TOKEN: ${{ secrets.NETLIFY_AUTH_TOKEN }} + # NETLIFY_SITE_ID: ${{ secrets.NETLIFY_SITE_ID }} + # timeout-minutes: 10 diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml new file mode 100644 index 0000000..6ea88e5 --- /dev/null +++ b/.github/workflows/docs.yml @@ -0,0 +1,28 @@ +name: docs +on: + push: + branches: + - master +jobs: + deploy: + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: actions/setup-python@v4 + with: + python-version: 3.9 + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install --user --no-cache-dir Cython + pip install --user -r requirements.txt + pip install . + - name: Discover typos with codespell + run: | + pip install codespell + codespell --skip="*.csv,*.geojson,*.json,*.js,*.html,*cff,./.git" --ignore-words-list="aci,acount,acounts,fallow,hart,hist,nd,ned,ois,wqs,watermask" + - name: PKG-TEST + run: | + python -m unittest discover tests/ + - run: pip install -r requirements_docs.txt + - run: mkdocs gh-deploy --force diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml new file mode 100644 index 0000000..4b3ed55 --- /dev/null +++ b/.github/workflows/macos.yml @@ -0,0 +1,38 @@ +on: + push: + branches: + - master + pull_request: + branches: + - master + +name: macOS build +jobs: + test-macOS: + runs-on: ${{ matrix.os }} + name: ${{ matrix.os }} (${{ matrix.python-version}}) + strategy: + fail-fast: false + matrix: + os: ["macOS-latest"] + python-version: ["3.10"] + steps: + - name: Checkout code + uses: actions/checkout@v3 + + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version}} + - name: Install GDAL + run: | + brew install gdal + - name: Test GDAL installation + run: | + gdalinfo --version + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install --no-cache-dir Cython + pip install -r requirements.txt + pip install . diff --git a/.github/workflows/publish-to-pypi.yml b/.github/workflows/publish-to-pypi.yml deleted file mode 100644 index 2bc7304..0000000 --- a/.github/workflows/publish-to-pypi.yml +++ /dev/null @@ -1,38 +0,0 @@ -name: Publish Python distributions to PyPI - -on: - push: - workflow_dispatch: - -jobs: - build-n-publish: - name: Build and publish Python distributions to PyPI - runs-on: ubuntu-latest - - steps: - - uses: actions/checkout@master - - name: Set up Python 3.9 - uses: actions/setup-python@v1 - with: - python-version: 3.9 - - - name: Install pypa/build - run: >- - python -m - pip install - build - --user - - name: Build a binary wheel and a source tarball - run: >- - python -m - build - --sdist - --wheel - --outdir dist/ - . - - - name: Publish distribution to PyPI - uses: pypa/gh-action-pypi-publish@release/v1 - with: - user: __token__ - password: ${{ secrets.HYDRO_XAJ_TOKEN }} diff --git a/.github/workflows/pypi.yml b/.github/workflows/pypi.yml new file mode 100644 index 0000000..2bbb0c0 --- /dev/null +++ b/.github/workflows/pypi.yml @@ -0,0 +1,30 @@ +# This workflows will upload a Python Package using Twine when a release is created +# For more information see: https://help.github.com/en/actions/language-and-framework-guides/using-python-with-github-actions#publishing-to-package-registries + +name: pypi + +on: + release: + types: [created] + +jobs: + deploy: + runs-on: ubuntu-latest + + steps: + - uses: actions/checkout@v3 + - name: Set up Python + uses: actions/setup-python@v4 + with: + python-version: "3.x" + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install setuptools wheel twine + - name: Build and publish + env: + TWINE_USERNAME: ${{ secrets.PYPI_USERS }} + TWINE_PASSWORD: ${{ secrets.PYPI_PASSWORD }} + run: | + python setup.py sdist bdist_wheel + twine upload dist/* diff --git a/.github/workflows/ubuntu.yml b/.github/workflows/ubuntu.yml new file mode 100644 index 0000000..6cb7c0c --- /dev/null +++ b/.github/workflows/ubuntu.yml @@ -0,0 +1,48 @@ +on: + push: + branches: + - master + - dev + pull_request: + branches: + - master + - dev + +name: Linux build +jobs: + py-check: + runs-on: ${{ matrix.config.os }} + name: ${{ matrix.config.os }} (${{ matrix.config.py }}) + strategy: + fail-fast: false + matrix: + config: + - { os: ubuntu-latest, py: "3.8" } + - { os: ubuntu-latest, py: "3.9" } + - { os: ubuntu-latest, py: "3.10" } + - { os: ubuntu-latest, py: "3.11" } + steps: + - name: Checkout Code + uses: actions/checkout@v3 + - name: Setup Python + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.config.py }} + - name: Install GDAL + run: | + python -m pip install --upgrade pip + pip install --no-cache-dir Cython + pip install --find-links=https://girder.github.io/large_image_wheels --no-cache GDAL + - name: Test GDAL installation + run: | + python -c "from osgeo import gdal" + gdalinfo --version + - name: Install dependencies + run: | + pip install --user -r requirements.txt + pip install --user -r requirements_dev.txt + pip install . + - name: PKG-TEST + run: | + python -m unittest discover tests/ + diff --git a/.github/workflows/windows.yml b/.github/workflows/windows.yml new file mode 100644 index 0000000..9e9deaa --- /dev/null +++ b/.github/workflows/windows.yml @@ -0,0 +1,33 @@ +on: + push: + branches: + - master + - dev + pull_request: + branches: + - master + - dev + +name: Windows build +jobs: + test-windows: + runs-on: windows-latest + steps: + - uses: actions/checkout@v3 + - name: Install miniconda + uses: conda-incubator/setup-miniconda@v2 + with: + auto-activate-base: true + python-version: "3.10" + - name: Install GDAL + run: conda install -c conda-forge gdal --yes + - name: Test GDAL installation + run: | + python -c "from osgeo import gdal" + gdalinfo --version + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install --no-cache-dir Cython + pip install -r requirements.txt + pip install . diff --git a/AUTHORS.rst b/AUTHORS.rst new file mode 100644 index 0000000..a199192 --- /dev/null +++ b/AUTHORS.rst @@ -0,0 +1,13 @@ +======= +Credits +======= + +Development Lead +---------------- + +* Wenyu Ouyang + +Contributors +------------ + +None yet. Why not be the first? diff --git a/MANIFEST.in b/MANIFEST.in index 57fdc18..89411aa 100644 --- a/MANIFEST.in +++ b/MANIFEST.in @@ -1 +1,7 @@ -recursive-include XXX *.csv *.txt # 打包需包含csv、txt为后缀的文件；XXX为包名 \ No newline at end of file +include LICENSE +include README.md +include requirements.txt + +recursive-exclude * __pycache__ +recursive-exclude * *.py[co] + diff --git a/environment-dev.yml b/environment-dev.yml index 288a934..d6be0cb 100644 --- a/environment-dev.yml +++ b/environment-dev.yml @@ -23,4 +23,8 @@ dependencies: # pip - pip - pip: + - setuptools + - wheel + - twine + - bump2version - hydrodataset diff --git a/requirements.txt b/requirements.txt index 74513e9..245c186 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,16 +1,14 @@ -python=3.10 ipykernel numpy numba pandas scikit-learn deap -spotpy=1.5.14 +spotpy==1.5.14 seaborn tqdm pytest hydrodataset bmipy pyyaml -requests -spotpy \ No newline at end of file +requests \ No newline at end of file diff --git a/requirements_dev.txt b/requirements_dev.txt new file mode 100644 index 0000000..c4e98f7 --- /dev/null +++ b/requirements_dev.txt @@ -0,0 +1,14 @@ +black +black[jupyter] +pip +bump2version +wheel +watchdog +flake8 +tox +coverage +Sphinx +twine + +pytest +pytest-runner \ No newline at end of file diff --git a/setup.cfg b/setup.cfg new file mode 100644 index 0000000..2e09fd3 --- /dev/null +++ b/setup.cfg @@ -0,0 +1,26 @@ +[bumpversion] +current_version = 0.0.1 +commit = True +tag = True + +[bumpversion:file:setup.py] +search = version='{current_version}' +replace = version='{new_version}' + +[bumpversion:file:hydromodel/__init__.py] +search = __version__ = '{current_version}' +replace = __version__ = '{new_version}' + +[bdist_wheel] +universal = 1 + +[flake8] +exclude = docs + +[aliases] +# Define setup.py command aliases here +test = pytest + +[tool:pytest] +collect_ignore = ['setup.py'] + diff --git a/setup.py b/setup.py index e143dd7..bf0dedf 100644 --- a/setup.py +++ b/setup.py @@ -1,30 +1,71 @@ +#!/usr/bin/env python """ Author: Wenyu Ouyang Date: 2023-10-28 09:16:46 -LastEditTime: 2023-10-28 09:27:22 +LastEditTime: 2023-10-28 17:36:12 LastEditors: Wenyu Ouyang -Description: setup.py for hydromodel package +Description: The setup script FilePath: \hydro-model-xaj\setup.py Copyright (c) 2023-2024 Wenyu Ouyang. All rights reserved. """ +import io import pathlib +from os import path as op from setuptools import setup, find_packages -readme = pathlib.Path("README.md").read_text() +readme = pathlib.Path("README.md").read_text(encoding='utf-8') +here = op.abspath(op.dirname(__file__)) + +# get the dependencies and installs +with io.open(op.join(here, "requirements.txt"), encoding="utf-8") as f: + all_reqs = f.read().split("\n") + +install_requires = [x.strip() for x in all_reqs if "git+" not in x] +dependency_links = [x.strip().replace("git+", "") for x in all_reqs if "git+" not in x] + +requirements = [] + +setup_requirements = [ + "pytest-runner", +] + +test_requirements = [ + "pytest>=3", +] + setup( - name="hydromodel", # 输入项目名称 - version="0.0.1", # 输入版本号 - keywords=[""], # 输入关键词 - description="", # 输入概述 + author="Wenyu Ouyang", + author_email="wenyuouyang@outlook.com", + python_requires=">=3.9", + classifiers=[ + "Intended Audience :: Developers", + "License :: OSI Approved :: GNU General Public License v3 (GPLv3)", + "Natural Language :: English", + "Programming Language :: Python :: 3", + "Programming Language :: Python :: 3.8", + "Programming Language :: Python :: 3.9", + "Programming Language :: Python :: 3.10", + "Programming Language :: Python :: 3.11", + ], + description="hydrological models starting from XinAnJiang", + entry_points={ + "console_scripts": [ + "hydromodel=hydromodel.cli:main", + ], + }, + install_requires=install_requires, + dependency_links=dependency_links, + license="GNU General Public License v3", long_description=readme, long_description_content_type="text/markdown", - url="https://github.com/iHeadWater/hydro-model-xaj", # 输入项目Github仓库的链接 - author="iHeadWater", # 输入作者名字 - author_email="", # 输入作者邮箱 - license="MIT_license", # 此为声明文件，一般填写 MIT_license - packages=find_packages(), include_package_data=True, - platforms="any", - install_requires=[""], # 输入项目所用的包 - python_requires=">= 3.7 ", # Python版本要求 + keywords="hydromodel", + name="hydromodel", + packages=find_packages(include=["hydromodel", "hydromodel.*"]), + setup_requires=setup_requirements, + test_suite="tests", + tests_require=test_requirements, + url="https://github.com/iHeadWater/hydro-model-xaj", + version="0.0.1", + zip_safe=False, )