如何使用 facet_wrap 将 NSE 和 PBIAS 结果添加到 ggplot？

Question

我是初学者，我创建了一个函数（在下方）来计算模拟数据与观察数据的百分比偏差 (PBIAS) 和 Nash-Sutcliffe 效率 (NSE)。但是我只能使用我的整个数据集来计算这些测试。

model.assess <- function(Sim, Obs) { 
  rmse = sqrt( mean( (Sim - Obs)^2, na.rm = TRUE) ) #Formula to calculate RMSE
  RSR <- rmse / sd(Obs) #object producing RSR test from the RMSE formula
  PBIAS <- 100 *(sum((Sim - Obs)/sum(Obs), na.rm =TRUE)) #object producing PBIAS test
  NSE <- 1 - sum((Obs - Sim)^2)/sum((Obs - mean(Obs))^2, na.rm =TRUE) #object producing NSE test
  stats <- print(paste0("RSR = ", sprintf("%.3f", round(RSR, digits=3)), "    PBIAS = ", sprintf("%.3f",round(PBIAS, digits=3)),"    NSE = ", sprintf("%.3f",round(NSE, digits=3))))  
  return(stats) #returns the results of the tests with 3 decimals and spacing in between

这是我的四个不同站点（SNS、MRC、TLG、SJF）的数据集，每月流量：

StationID Date      Obs_flow    Sim_flow      Month     Year
SNS    1950-10-01   0.010170    0.030687967 October 1950-01-01      
SNS    1950-11-01   0.366260    0.416466741 November 1950-01-01     
SNS    1950-12-01   0.412210    0.496136731 December 1950-01-01     
SNS    1951-01-01   0.119520    0.182072570 January 1951-01-01      
SNS    1951-02-01   0.113480    0.142611192 February 1951-01-01     
SNS    1951-03-01   0.127090    0.176350274 March   1951-01-01  
SNS    1951-04-01   0.175120    0.193221389 April   1951-01-01      
SNS    1951-05-01   0.208940    0.275980903 May     1951-01-01  
SNS    1951-06-01   0.114420    0.144675317 June    1951-01-01      
SNS    1951-07-01   0.032280    0.018057796 July    1951-01-01

要使用我使用的方程和 R 平方绘制 Obs 与 Sim 的散点图：

dataset %>%
  filter(StationID == "SNS") %>%
ggplot(aes(x = Obs_flow, y = Sim_flow)) + 
   geom_point(aes(Obs_flow, Sim_flow), alpha = 0.3)+
     stat_smooth(aes(x = Obs_flow, y = Sim_flow), 
                method = "lm", se = TRUE, colour="#FC4E07", fullrange = TRUE) + 
   stat_poly_eq(formula = "y~x", 
             aes(label = paste0(..eq.label..)),  #adding the equation on the top
             parse = TRUE, label.x.npc = "center", label.y.npc = 0.97, size = 3.45, family= "Times New Roman")+
     stat_poly_eq(formula = "y~x", 
             aes(label = paste0(..rr.label..)), #adding the Rsquared at the bottom
             parse = TRUE, label.x.npc = 0.95, label.y.npc = 0.05, size = 3.45, family= "Times New Roman")+
  
  annotate("text", x = 0, y = 1.3,, label = paste0(model.assess(dataset$Sim_flow, dataset$Obs_flow)),  collapse = "\n", hjust = 0, size=2.4, family= "Times New Roman") +
  
   facet_wrap(~ Month, ncol=4,  labeller = labeller(StationID = c("MRC" = "Merced River", "SJF"= "Upper San Joaquin River", "SNS" = "Stanislaus River", "TLG" = "Tuolumne River")), scales = "fixed")

stat_poly_eq 为每个方面添加了一个方程和 Rsquared，但注释为所有方面添加了相同的数字。有没有办法分别为每个方面添加 NSE 和 PBIAS？我尝试了包 HydroGOF，但我得到了相同的结果。请原谅审美。

Answer 1

示例数据将有助于其他人为您提供帮助。请查看此 link 以备将来查询。

你有一些问题。您的 model.assess() 函数提供一条记录，而您需要每个方面的值。所以，我使用代码

创建了一个虚拟对象

ll <- data.frame(Month=c(),label=c())
nM <- length(Month)
lapply(1:nM, function(i){
  a <- Sim_flow*i*i*0.5
  b <- Obs_flow*i
  m <- model.assess(a,b)
  ll <<- rbind(ll,data.frame(Month=Month[i],label=m))
})
labels <- ll

接下来，您需要使用geom_label而不是中提到的注释。下面的代码

ggplot(data=dataset, aes(x = Obs_flow, y = Sim_flow)) + 
  geom_point(aes(Obs_flow, Sim_flow), alpha = 0.3)+
  stat_smooth(aes(x = Obs_flow, y = Sim_flow), 
              method = "lm", se = TRUE, colour="#FC4E07", fullrange = TRUE) + 
  stat_poly_eq(formula = "y~x", 
               aes(label = paste0(..eq.label..)),  #adding the equation on the top
               parse = TRUE, label.x.npc = "center", label.y.npc = 0.97, size = 3.45, family= "Times New Roman") +
  stat_poly_eq(formula = "y~x", 
               aes(label = paste0(..rr.label..)), #adding the Rsquared at the bottom
               parse = TRUE, label.x.npc = 0.95, label.y.npc = 0.05, size = 3.45, family= "Times New Roman") +
  facet_wrap(~Month, ncol=4,  labeller = labeller(StationID = c("MRC" = "Merced River", "SJF"= "Upper San Joaquin River", "SNS" = "Stanislaus River", "TLG" = "Tuolumne River")), scales = "fixed") +
  geom_label(data = labels, aes(label=label, x = Inf, y = -Inf), 
             hjust=1, vjust=0, size=1.8,
             inherit.aes = FALSE)

给出以下

如何使用 facet_wrap 将 NSE 和 PBIAS 结果添加到 ggplot？

How can I add NSE and PBIAS results to ggplot using facet_wrap?

r

ggplot2

facet-wrap