はじめに
SWATplusRのバージョンが1年前から大きく変わっていたため改めて書き直す。
用いるRの基本操作は前回の記事を見てください。
現バージョンは、前回のバージョンと大きく変わった点はなく細かい不都合が減った印象。
環境
R version 4.1.1
SWATplusR 0.6
rev60.5.4_64rel.exe
Windows10
SWATplusRとは
RでSWAT+(SWAT2012)が動かせるようになる、Rのライブラリ。
モデルのパラメータやシミュレーション期間の変更が可能であり、SWAT+ToolBoxと違い並列処理も可能。
出力データはRでさらに加工が可能。
SWAT+ToolBoxよりも効率的にキャリブレーションが可能です(当社比)。
SWATplusEditorでTxtInOutを作成した後、SWATplusRでシミュレーションがよろしい。
セットアップ
SWATplusRのインストール
SWATplusRはGithub上のRパッケージをインストールすることで使えるようになる。
以下のコードでインストール
# remotesをインストールしていなければ使用
install.packages('remotes')
remotes::install_github('chrisschuerz/SWATplusR')
実行フォルダの構築
SWATplusRはTextInOutに対して実行される。
その際、TextInOutには実行ファイル(.exe)が含まれている必要がある(例: rev60.5.4_64rel.exe)。
実行ファイルが新し過ぎたり、SWATplusRが古すぎたりすると実行したときにエラーが生じる。
実行ファイルは2通りの方法で入手できる。
1.デモファイルからの入手
以下のコードでデモファイルを入手することができる。
そのファイルの実行ファイルを実行したいTxtInOutに使用する。
library(SWATplusR)
demo_path <- 'C:/user/hoge'
path_plus <- load_demo(dataset = 'project',
path = demo_path,
version = 'plus')
2.SWATplusEditorからの入手
以下にSWATplusEditorで使用している実行ファイルがある。
C:\SWAT\SWATPlus\SWATPlusEditor\resources\app.asar.unpacked\swat_exe
この中のrel.exeを実行したいTxtInOutに使用。
実行
シミュレーション
シミュレーションの実行は以下のコードで行う。
library(SWATplusR)
path_plus = "C:/user/hoge/TextInOut"
q_sim <- run_swatplus(project_path = path_plus,
# 出力する項目を選択(詳しい項目は https://swatplus.gitbook.io/docs/user/io)
output = list(q_sur = define_output(file = "basin_wb",
variable = "surq_gen",
unit = 1),
q_lat = define_output(file = "basin_wb",
variable = "latq",
unit = 1),
Flow_Sim = define_output(file = "channel_sd",
variable = "flo_out",
unit = 1)),
# シミュレーション開始日
start_date = "1999-01-01",
# シミュレーション終了日
end_date = "2005-12-31",
# ウォームアップ期間(初めから何年出力しないか)
years_skip = 1,
# 使用実行ファイルのバージョン(無くてもいい)
revision = 60.5)
パラメータ変更
パラメータの変更は指定した値に変更する方法、与えられた範囲内で値を変更する方法に2種類がある。
また、サブ流域や土地利用を指定して細かくパラメータの変更もできる(詳しくは公式)。
パラメータを変更しているにも関わらず，モデル実行結果に変動がなければSWATplusRを入れ直すといいかも
1. 指定した値に変更
変更方法はabschg(その値の分パラメータを加減), pctchg(その値の分パラメータを百分率で変更), relchg(その値分パラメータを割合で変更), absval(その値でパラメータを置換)の4種類、pctchgとrelchgはほぼ同じ。
パラメータの後ろの.hruや.solはTxtInOutのcal_params.calで確認可能。
以下のコードで実行可能。
library(SWATplusR)
par_chg = c('lat_ttime.hru | change = absval' = 0.5, # lat_ttimeの値を0.5に置換
'lat_len.hru | change = abschg' = 30, # lat_lenの値に30加える
'esco.hru | change = pctchg' = 50) # escoの値を50%増しにする
path_plus = "C:/user/hoge/TextInOut"
q_sim <- run_swatplus(project_path = path_plus,
output = list(q_sur = define_output(file = "basin_wb",
variable = "surq_gen",
unit = 1),
q_lat = define_output(file = "basin_wb",
variable = "latq",
unit = 1),
Flow_Sim = define_output(file = "channel_sd",
variable = "flo_out",
unit = 1)),
start_date = "1999-01-01",
end_date = "2005-12-31",
years_skip = 1,
revision = 60.5,
# ここでパラメータを変更
parameter = par_chg)
2. 範囲を決めて値を変更
範囲を定める場合は試行回数も入力する。
シミュレーション回数が増えるため、並列処理をすることが望ましい。
並列処理数はn_threadで決める
library(SWATplusR)
times = 30 # 変更するパラメータセット数
par_set <- tibble("alpha.aqu|change = absval" = runif(times,0.01, 0.9), # 0.9から0.01の範囲でalphaの値を置換
"awc.sol|change = abschg" = runif(times,-0.1, 0.2),
"cn2.bsn|change = pctchg" = runif(times,-30,30))
path_plus = "C:/user/hoge/TextInOut"
q_sim <- run_swatplus(project_path = path_plus,
output = list(q_sur = define_output(file = "basin_wb",
variable = "surq_gen",
unit = 1),
q_lat = define_output(file = "basin_wb",
variable = "latq",
unit = 1),
Flow_Sim = define_output(file = "channel_sd",
variable = "flo_out",
unit = 1)),
start_date = "1999-01-01",
end_date = "2005-12-31",
years_skip = 1,
revision = 60.5,
# ここでパラメータを変更
parameter = par_set.
# 同時に6回シミュレーションを行う。
n_thread = 6)
モデル解析
モデル評価
モデルの評価にはナッシュサトクリフ係数、決定係数、パーセントバイアスを使用。ライブラリはhydroGOF。
実測値はcsv等をtibbleで読み込むとよい。
library(hydroGOF)
# flowの列名が"Flow_Sim"を選択、flowの列名が"Flow_Obs"を選択しNSEを計算
# Nash-Sutcliffe efficiency, NSE(sim, obs)
NSE(select(.data = flow, "Flow_Sim"), select(.data = flow,"Flow_Obs"))
# Coefficient of determination, br2(sim, obs)
br2(select(.data = flow, "Flow_Sim"), select(.data = flow,"Flow_Obs"))
#Percent Bias, pbias(sim, obs)
pbias(select(.data = flow, "Flow_Sim"), select(.data = flow,"Flow_Obs"))
感度評価
パラメータ感度を確認する際はsensitivityを使用(グローバル感度解析です)
ナッシュサトクリフ係数を用いてパラメータの感度を評価。
実測値のcsvの形は以下
date, flow_obs
2000/1/1, 13
2000/1/2, 12
.
.
library(sensitivity)
library(hydroGOF)
library(ggplot2)
library(SWATplusR)
library(fast)
# モデルのpath
proj_path = "D:/user/SWATplusR/TxtInOut"
# 実測値のpath
path_obs_flow = "D:/user/SWATplusR/ObservedFlow.csv"
# 感度評価したいパラメータ、変更方法、その範囲
par_con <- c("cn2.hru | change = pctchg", -40, 60,
"lat_ttime.hru | change = absval", 0.5, 70,
"lat_len.hru | change = absval", 10, 100)
# 入力パラメータの形を整形
par_names = c()
par_mins = c()
par_maxs = c()
for (i in 1:(length(par_con)/3)){
par_names = append(par_names, par_con[3*(i-1)+1])
par_mins = append(par_mins , as.numeric(par_con[3*(i-1)+2]))
par_maxs = append(par_maxs , as.numeric(par_con[3*(i-1)+3]))
}
par_fast <- fast_parameters(minimum = par_mins,
maximum = par_maxs,
names = par_names)
# シミュレーションを実行
q_fast <- run_swatplus(project_path = proj_path,
output = list(q_out = define_output(file = "channel_sd",
variable = "flo_out",
unit = 1)),
parameter = par_fast,
start_date = "1999-01-01",
end_date = "2014-12-31",
years_skip = 1,
n_thread = 6)
# 実測値の読み取り
q_obs <- read_csv(path_obs_flow,show_col_types = FALSE)
# NSEの計算
nse_fast <- q_fast$simulation$q_out %>%
select(-date) %>%
map_dbl(., ~NSE(.x, q_obs$flow_obs))
#Sensitivityの計算
sens_fast <- sensitivity(nse_fast, length(par_names))
# 結果を整形して出力
result_fast <- tibble(parameter = q_fast$parameter$definition$par_name,
Sensitivity = sens_fast) %>%
mutate(parameter = factor(parameter) %>% fct_reorder(., Sensitivity))
print(result_fast)
# 結果を棒グラフで出力
grap = ggplot(result_fast,
aes(x = Sensitivity, y = parameter)) +
geom_col()
print(grap)
参考文献
https://www.marshallplan.at/images/All-Papers/MP-2018/Sch%C3%BCrz+Christoph+822.pdf
https://chrisschuerz.github.io/SWATplusR/index.html
https://cran.r-project.org/web/packages/hydroGOF/hydroGOF.pdf
https://cran.r-project.org/web/packages/sensitivity/sensitivity.pdf