LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@Bulgent(Aguro Moritaka)

SWATplusR 基本操作

Last updated at Posted at 2021-09-13

はじめに

最近使い始めたのでその備忘録

環境

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

パラメータ変更①

こちらによると、あらかじめ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)

参考文献

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?