SWATplusR基本操作②

# はじめに
SWATplusRのバージョンが1年前から大きく変わっていたため改めて書き直す。
用いるRの基本操作は[前回の記事](https://qiita.com/Bulgent/items/8a444a7fc043c5671a37)を見てください。
現バージョンは、前回のバージョンと大きく変わった点はなく細かい不都合が減った印象。

### 環境
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パッケージをインストールすることで使えるようになる。
以下のコードでインストール
```R
 #　remotesをインストールしていなければ使用
install.packages('remotes')

remotes::install_github('chrisschuerz/SWATplusR')
```

## 実行フォルダの構築
SWATplusRはTextInOutに対して実行される。
その際、TextInOutには実行ファイル(.exe)が含まれている必要がある(例: rev60.5.4_64rel.exe)。
実行ファイルが新し過ぎたり、SWATplusRが古すぎたりすると実行したときにエラーが生じる。

実行ファイルは2通りの方法で入手できる。
#### 1.デモファイルからの入手
以下のコードでデモファイルを入手することができる。
そのファイルの実行ファイルを実行したいTxtInOutに使用する。

```R
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に使用。

# 実行
## シミュレーション
シミュレーションの実行は以下のコードで行う。
```R
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種類がある。
また、サブ流域や土地利用を指定して細かくパラメータの変更もできる(詳しくは[公式](https://chrisschuerz.github.io/SWATplusR/articles/SWATplusR.html#model-parameter-alteration))。

::: note info
パラメータを変更しているにも関わらず，モデル実行結果に変動がなければSWATplusRを入れ直すといいかも
:::

#### 1. 指定した値に変更
変更方法はabschg(その値の分パラメータを加減), pctchg(その値の分パラメータを百分率で変更), relchg(その値分パラメータを割合で変更), absval(その値でパラメータを置換)の4種類、pctchgとrelchgはほぼ同じ。
パラメータの後ろの.hruや.solはTxtInOutのcal_params.calで確認可能。
以下のコードで実行可能。
```R
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で決める



```R
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で読み込むとよい。
```R
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を使用(グローバル感度解析です)
ナッシュサトクリフ係数を用いてパラメータの感度を評価。

:::note info
実測値のcsvの形は以下
date, flow_obs
2000/1/1, 13
2000/1/2, 12
.
.

:::
```R
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