SWATplusR 基本操作

#はじめに
 最近使い始めたのでその備忘録

##環境
R version 4.1.1
SWATplusR 0.3.5.1
rev60.5.2_64rel.exe
##R基本操作
####カレントディレクトリを出力

    getwd()

####ディレクトリを移動
pathの書き方注意(例 C:/Users/tokachi/sample.txt)

      setwd("C:/Users/tokachi/sample.txt")

####パッケージインストール
パッケージはまず初めにインストールしないとlibrary()出来ない

      install.packages("SWATplusR")

####.rファイルを読み込む

    source("hoge.r")

####引数確認

    help(run_swatplus)

####tibbleをcsvに保存

    write_csv(hoge,"C:/Users/tokachi/hoge.csv")

####tibbleの表示数変更

    #tibbleの表示上限を1000に変更
    options("tibble.print_max" = 1000)

####tibbleの表示範囲選択

    #2行目、5列目のデータ取得
    q_sim[1,5]

    #3行目のデータを取得
    q_sim[3,]

    #3~6列目のデータ取得
    q_sim[,3:6]

####データ型
 - integer 
 
 整数型(1, -1, 210...)

 - double
 
実数型(3.12, 2.7181 ...)

 - character

文字型("Cattle", "Soil" ...)

 - logical

論理型(TRUE, FALSE, T, F, NA ...)

##SWATモデルを動かす
project_pathはTxtInOutの絶対パス。
TxtInOutにはswatplus_rev60.5.2.exeを入れておく。

###libraryの読み込み

    library(SWATplusR)

>以下のメッセージが出力されたが動いてはいるっぽい?
> >警告メッセージ: 
>>Missing column names filled in: 'X6' [6] 



    path_plus = "C:/Users/tokachi/TxtInOut"    

    wb_sim <- run_swatplus(project_path = path_plus,
                           output = list(precip = define_output(file = "basin_wb",
                                                                variable = "precip",
                                                                unit = 1),
                                         q_sur  = define_output(file = "basin_wb",
                                                                variable = "surq_gen",
                                                                unit = 1),
                                         q_lat  = define_output(file = "basin_wb",
                                                                variable = "latq",
                                                                unit = 1),
                                         eta    = define_output(file = "basin_wb",
                                                                variable = "et",
                                                                unit = 1),
                                        Flow_Sim = define_output(file = "channel_sd",
                                                                 variable = "flo_out",
                                                                 unit = 5)),
                           revision = 60.5)

####output
fileはTxtInOutのファイル名(例 basin_wb, channel等), variableはそのfile内の変数, unitのそのfile内にあるunit値を参照

##シミュレーション結果の整理

シミュレーション結果と実測値を同じグラフ上に書く。

###libraryの読み込み


    library(tidyverse)
    library(dplyr)
    library(ggplot2)

###実測値の読み込み＆データ編集
tibbleで読み込んだ実測値(.csv)の実測値部分(Flow_Obs)をシミュレーション結果の表に結合する。

    #tibbleで読み込む
    obs <- read_csv("C:/Users/tokachi/observe.csv")
    #推定値と実測値の表を組み合わせる
    flow = mutate(wb_sim, select(.data = obs, Flow_Obs))

###グラフの表示

    graph <- ggplot(data = flow, aes(date)) +
             geom_line(aes(y = Flow_Obs, colour = "Flow_Obs")) +
             geom_line(aes(y = Flow_Sim, colour = "Flow_sim")) +
             ylab("discharge")

##評価
###libraryの読み込み

    library(hydroGOF)

###NSE, R2, PBIASを用いた評価

    #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"))

##パラメータ変更①

[こちら](https://groups.google.com/g/swatplus/c/PMcC_zKsQMg/m/HXp7HhDWEwAJ)によると、あらかじめfile.cioの22行目を以下の文に置換しないとパラメータ変更が推定値に反映されないらしい。

    chg               cal_parms.cal     calibration.cal   codes.sft         wb_parms.sft      water_balance.sft  null              null              null              null     

pctchgは％分変更, relchgは割合分変更, abschgはその値の分変更, absvalはその値に置換する。


    #HRUレベルで定義されているcn2の初期値を-5%し, GWレベルで定義されているalphaの値を0.5に置き換える。
    par_single <- c("cn2.hru|change = pctchg" = -5,
                    "alpha.aqu|change = absval" = 0.5)
    #変更したパラメータでモデルを動かす。
    q_sim <- run_swatplus(project_path = path_plus,
                          output       = list(Flow_sim = define_output(file = "channel_sd",
                                                       variable = "flo_out",
                                                       unit = 1)),
                          parameter    = par_single)
    

##パラメータ変更②

以下のコードは, cn2とalphanの値を決められた範囲(cn2なら-15~10の間の値のabchg)のランダムな値に8回変化させることを意味する。

    par_set <- tibble("cn2.hru|change = abschg" = runif(8,-15,10),
                  "alpha.aqu|change = absval" = runif(8, 0, 1))

このpar_setを用いて同様にrun_swatplusできる。n_threadは並列で処理する量を意味する。

    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 = 5)),
                      parameter = par_set,
                      start_date  = "2005-01-01",
                      end_date = "2014-12-31",
                      years_skip = 2,
                      n_thread = 6)

どんなパラメータをとったかは, 以下のコードで確認できる。

    q_sim$parameter

シミュレーション結果は以下のコードで確認できる。

    q_sim$simulation

##ソースコード
パラメータをランダムに変えて, その評価をcsvで出力してくれるやつ
参考までに


    #ライブラリのインポート
    library(SWATplusR)
    library(tidyverse)
    library(dplyr)
    library(ggplot2)
    library(hydroGOF)

    ModelRun <- function(){
        #SWATモデルを2005-2014の期間でまわし、surq_gen, latq, flo_outを出力する。
        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 = 5)),
                    parameter = par_set,
                    #並列処理数
                    n_thread = 6,
                    start_date  = "2005-01-01",
                    end_date = "2014-12-31",
                    years_skip = 2,
                    revision = 60.5)
        return(q_sim)
    }

    Evaluation <- function(x){
        #flo_outのsimulationとobservationを評価する。
        names <- c()
        nse <- c()
        r2 <- c()
        pbias <- c()
        obs <- read_csv(path_obs)
        num <- length(x$simulation$Flow_Sim)
        for (i in 2:num){
            sim <- x$simulation$Flow_Sim[i]
            obs <- select(.data = obs, Flow_Obs)
            name <- paste("run_", i-1, sep = "")
            n <- NSE(sim, obs)
            r <- br2(sim, obs)
            p <- pbias(sim, obs)
            nse = append(nse,n,after = i)
            r2 = append(r2,r,after = i)
            pbias = append(pbias,p,after = i)
            names = append(names,name,after = i)
        }
        score <- data.frame(
        Name = names,
        NSE = nse,
        PBIAS = pbias,
        R2 = r2)

        score <- mutate(.data = score,q_sim$parameter$value)
        return(score)
    }

    WriteCsv <- function(name){
        #flo_outと評価をそれぞれcsvに書き出す
        write_csv(q_sim$simulation$Flow_Sim,paste("C:/User/",name,"_flow.csv",sep = ""))
        write_csv(score,paste("C:/User/",name,"_score.csv",sep = ""))
    }

    #TxtInOutにはrev.exeを入れておく。
    path_plus <- "C:/User/SWATplusR/TxtInOut"

    #observationはsimilationと期間を合わせておく
    path_obs <- "C:/User/SWATplusR/obs_flow2007_2014.csv"

    #csv名決定
    value = "caliblation"

    #シミュレーション回数
    fois = 54

    #パラメータ設定
    par_set <- tibble("awc.sol|change = abschg" = runif(fois,-0.3,-0.1),
                    "cn2.bsn|change = pctchg" = runif(fois,-20,-10),
                    "cn3_swf.hru|change = pctchg" = runif(fois,-10.5,-10),
                    "z.sol|change = pctchg" = runif(fois,45,60),
                    "bd.sol|change = abschg" = runif(fois,-0.5,-0.3),
                    "alpha.aqu|change = abschg" = runif(fois,0.2,0.21),
                    "k.sol|change = pctchg" = runif(fois,20,30),
                    "snomelt_tmp.hru|change = abschg" = runif(fois,-4,-1),
                    "snofall_tmp.hru|change = abschg" = runif(fois,0,6),
                    "snomelt_max.hru|change = abschg" = runif(fois,-5,0),
                    "snomelt_lag.hru|change = abschg" = runif(fois,-0.5,-0.1))

    start.time<- proc.time()

    q_sim <- ModelRun()
    score <- Evaluation(q_sim)
    print(score)
    WriteCsv(value)

    print(proc.time()-start.time)

##参考文献

https://chrisschuerz.github.io/SWATplusR/articles/SWATplusR.html#model-parameter-alteration

https://cran.r-project.org/web/packages/hydroGOF/hydroGOF.pdf

https://www.marshallplan.at/images/All-Papers/MP-2018/Sch%C3%BCrz+Christoph+822.pdf