R中的简单供水需求模型
Simple water demand supply model in R
我正在尝试在 R 中计算出一个简单的土壤水分平衡。这是我需要执行的步骤:
对于给定的年份,从第 200 天开始,我需要确定通过以下公式计算的土壤水分 (SWi)
`SW(i) = SW(i-1) + Rain(i) - ETo(i)`
其中SW(i-1)
是前一天的含水量,Rain(i)
是降雨量,ETo(i)
是第i天的蒸散量
条件是SW(i)不能为负或大于SW(max)即20。
这是一个示例数据:
library(tidyverse)
set.seed(123)
dat <- tibble(
year = rep(1980:2015, each = 100),
day = rep(200:299, times = 36),
rain = sample(0:17, size = 100*36,replace = T),
eto = sample(2:9, size = 100*36,replace = T))
SW.initial <- data.frame(year= 1980:2015, SW.199 = sample(1:10, 36, replace = T))
SW.initial 是每年 doy 199 的含水量
SW.max <- 20
dat$SW.fin <- NA
以1980年为例
dat.1980 <- dat[dat$year == 1980,]
SW.initial.1980 <- SW.initial[SW.initial$year== 1980,"SW.199"]
for(doy in dat.1980$day){
SW <- SW.initial.1980
SW <- SW + dat.1980[dat.1980$day == doy, "rain"] - dat.1980[dat.1980$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.initial.1980 <- SW
}
这个循环会给我从 doy 200 到 299 的每一天的 SW,使用:
`SW(i) = SW(i-1) + Rain[i] + ETo[i]``
对于 doy 200,SW(i-1) 是从 SW.initial
数据帧
中给出的
我可以循环遍历所有年份:
for(years in unique(dat$year)){
test <- dat[dat$year == years,]
SW.in <- SW.initial[SW.initial$year == years,"SW.199"]
for(doy in test$day){
SW <- SW.in
SW <- SW + test[test$day == doy, "rain"] - test[test$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.in <- SW
}}
我真的很想避免这个长循环,并且在想是否有更聪明(更快的方法)。
这是你想要的吗?
编辑 -> 添加按年份分组
dat %>% group_by(year) %>% mutate(sw_oneless = c(NA, day[1:length(day)-1]),
sw_oneless + rain - eto)
# A tibble: 3,600 x 6
# Groups: year [36]
year day rain eto sw_oneless `sw_oneless + rain - eto`
<int> <int> <int> <int> <int> <int>
1 1980 200 5 2 NA NA
2 1980 201 14 6 200 208
3 1980 202 7 4 201 204
4 1980 203 15 5 202 212
5 1980 204 16 5 203 214
6 1980 205 0 8 204 196
7 1980 206 9 9 205 205
8 1980 207 16 6 206 216
9 1980 208 9 9 207 207
10 1980 209 8 4 208 212
# ... with 3,590 more rows
要解决第 200 天的 "problem",为什么不从原始数据中过滤掉第 199-300 天?然后您可以 运行 我的代码和 na.omit() 或再次过滤,第 199 天的行消失了。
或者,如果您不能这样做,请将您的 SW.initial 与您的 dat 数据框
合并
我正在尝试在 R 中计算出一个简单的土壤水分平衡。这是我需要执行的步骤:
对于给定的年份,从第 200 天开始,我需要确定通过以下公式计算的土壤水分 (SWi)
`SW(i) = SW(i-1) + Rain(i) - ETo(i)`
其中SW(i-1)
是前一天的含水量,Rain(i)
是降雨量,ETo(i)
是第i天的蒸散量
条件是SW(i)不能为负或大于SW(max)即20。
这是一个示例数据:
library(tidyverse)
set.seed(123)
dat <- tibble(
year = rep(1980:2015, each = 100),
day = rep(200:299, times = 36),
rain = sample(0:17, size = 100*36,replace = T),
eto = sample(2:9, size = 100*36,replace = T))
SW.initial <- data.frame(year= 1980:2015, SW.199 = sample(1:10, 36, replace = T))
SW.initial 是每年 doy 199 的含水量
SW.max <- 20
dat$SW.fin <- NA
以1980年为例
dat.1980 <- dat[dat$year == 1980,]
SW.initial.1980 <- SW.initial[SW.initial$year== 1980,"SW.199"]
for(doy in dat.1980$day){
SW <- SW.initial.1980
SW <- SW + dat.1980[dat.1980$day == doy, "rain"] - dat.1980[dat.1980$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.initial.1980 <- SW
}
这个循环会给我从 doy 200 到 299 的每一天的 SW,使用:
`SW(i) = SW(i-1) + Rain[i] + ETo[i]``
对于 doy 200,SW(i-1) 是从 SW.initial
数据帧
我可以循环遍历所有年份:
for(years in unique(dat$year)){
test <- dat[dat$year == years,]
SW.in <- SW.initial[SW.initial$year == years,"SW.199"]
for(doy in test$day){
SW <- SW.in
SW <- SW + test[test$day == doy, "rain"] - test[test$day == doy, "eto"]
SW <- ifelse(SW < 0, 0, ifelse(SW >= SW.max, SW.max, SW))
dat[dat$year == years & dat$day == doy,"SW.fin"] <- SW
SW.in <- SW
}}
我真的很想避免这个长循环,并且在想是否有更聪明(更快的方法)。
这是你想要的吗?
编辑 -> 添加按年份分组
dat %>% group_by(year) %>% mutate(sw_oneless = c(NA, day[1:length(day)-1]),
sw_oneless + rain - eto)
# A tibble: 3,600 x 6
# Groups: year [36]
year day rain eto sw_oneless `sw_oneless + rain - eto`
<int> <int> <int> <int> <int> <int>
1 1980 200 5 2 NA NA
2 1980 201 14 6 200 208
3 1980 202 7 4 201 204
4 1980 203 15 5 202 212
5 1980 204 16 5 203 214
6 1980 205 0 8 204 196
7 1980 206 9 9 205 205
8 1980 207 16 6 206 216
9 1980 208 9 9 207 207
10 1980 209 8 4 208 212
# ... with 3,590 more rows
要解决第 200 天的 "problem",为什么不从原始数据中过滤掉第 199-300 天?然后您可以 运行 我的代码和 na.omit() 或再次过滤,第 199 天的行消失了。 或者,如果您不能这样做,请将您的 SW.initial 与您的 dat 数据框
合并