Modelación de zonas buffer
Rafael Navas -
Última actualización: 03 junio 2020
1 Escenarios
Se definen escenarios similares al punto anterior, pero en esta ocasión se varían de forma uniforme en cada zona. En otras palabras, no se diferencia entre tributario y cauce principal.
2 Resultados
La siguiente figura muestra la probabilidad de no excedencia de la concentración de fósforo para el escenario base (situación actual) y los escenarios de zona buffer 40-20 y 10-05 (Principal-tributario en metros). Ambos escenarios de zona buffer producen resultados bastante similares. Este resultado deberá ser interpretado con cuidado.
La siguiente aplicación permite visualizar el efecto de distintas configuraciones de FILTERW en la producción de fósforo.
https://rafaelnavas23.shinyapps.io/app_EscenariosSantaLucia/
rm(list = ls())
library(tibble)
x_con = readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A0_B0_C0.RDS")$PTcon
PTbase <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A0_B0_C0.RDS")$rch_41
PT40 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A40_B40_C40.RDS")$rch_41
PT20 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A20_B20_C20.RDS")$rch_41
PT10 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A10_B10_C10.RDS")$rch_41
PT05 <- readRDS("~/R_projects/SWAT_calibration_ETmodis/resultados/ZonaBuffer_optim22/output_buffer_diarioFW/PTcon_ecdf_A5_B5_C5.RDS")$rch_41
plot(x_con,PTbase$q50, type="l", ylim=c(0.01,1),lty=1,lwd=2,
xlab="concentración (mg/l)", ylab="Probabilidad de no excedencia", log="x")
lines(x_con, PT40$q50,lty=1,col=2,lwd=2)
lines(x_con, PT20$q50,lty=1,col=3,lwd=2)
lines(x_con, PT10$q50,lty=1,col=4,lwd=2)
lines(x_con, PT05$q50,lty=1,col=5,lwd=2)
legend("topleft", c("base", "P40", "P20","P10","P05"),lwd=2,col=1:5)
abline(v=0.025)
Anexo: código para correr los escenarios
rm(list=ls())
library(tibble)
library(SWATplusR)
library(stringi)
library(stringr)
library(emdbook)
# lectura de parámetros optimizados
mypar = readRDS("./input_buffer/best_par.RDS")
# escenarios a modelar
escenarios = read.table("./input_buffer/FILTER_W_escenarios.csv",sep=";", header = T)
# Zonificación
zonaA = c(1,7,8,10,14,15,16,17,18,19,20,22,24,28)
zonaB = c(4,5,25,31,3,27,2,23,29)
zonaC = c(41,12,21,30,9,13,32,38,36,34,6,11,37,35,39,40,26,33)
## Inicio de loop en los escenarios
########################################
for(n in 1:nrow(escenarios)){
print(paste("escenario:" ,n))
par_set = tibble("r1_CN2::CN2.mgt|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_CN2"]*100,
"r1_ESCO::ESCO.hru|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_ESCO"],
"r1_SOL_AWC::SOL_AWC.sol|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_SOL_AWC"]*100,
"r1_ALPHA_BF::ALPHA_BF.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_ALPHA_BF"],
"r1_GWQMN::GWQMN.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_GWQMN"],
"r1_GW_DELAY::GW_DELAY.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_GW_DELAY"],
"r1_REVAPMN::REVAPMN.gw|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_REVAPMN"],
"r1_GW_REVAP::GW_REVAP.gw|change = absval |
sub == c(1, 16, 18, 22,23)" =
mypar[,"r1_GW_REVAP"],
"r1_OV_N::OV_N.hru|change = absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_OV_N"],
"r1_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" =
mypar[,"r1_SLSUBBSN"]*100,
"r2_CN2::CN2.mgt|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_CN2"]*100,
"r2_ESCO::ESCO.hru|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_ESCO"],
"r2_SOL_AWC::SOL_AWC.sol|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_SOL_AWC"]*100,
"r2_ALPHA_BF::ALPHA_BF.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_ALPHA_BF"],
"r2_GWQMN::GWQMN.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_GWQMN"],
"r2_GW_DELAY::GW_DELAY.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_GW_DELAY"],
"r2_REVAPMN::REVAPMN.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_REVAPMN"],
"r2_GW_REVAP::GW_REVAP.gw|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_GW_REVAP"],
"r2_OV_N::OV_N.hru|change = absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_OV_N"],
"r2_SLSUBBSN::SLSUBBSN.hru|change = pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_SLSUBBSN"]*100,
"r1_HRU_SLP::HRU_SLP.hru|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_HRU_SLP"]*100,
"r1_USLE_K::USLE_K.sol|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[,"r1_USLE_K"]*100,
"r1_USLE_P::USLE_P.mgt|change= pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[,"r1_USLE_P"]*100,
"r1_P_UPDIS::P_UPDIS.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_P_UPDIS"],
"r1_FRT_SURFACE::FRT_SURFACE.mgt|change=absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[,"r1_FRT_SURFACE"],
"r1_FRT_KG::FRT_KG.mgt|change=pctchg |
sub == c(1,8,10,14,15,16,17,18,19,20,22)" = mypar[,"r1_FRT_KG"]*100,
"r1_NPERCO::NPERCO.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"=mypar[,"r1_NPERCO"],
"r1_ERORGN::ERORGN.hru|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_ERORGN"],
"r2_HRU_SLP::HRU_SLP.hru|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_HRU_SLP"]*100,
"r2_USLE_K::USLE_K.sol|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)" =
mypar[,"r2_USLE_K"]*100,
"r2_USLE_P::USLE_P.mgt|change= pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_USLE_P"]*100,
"r2_P_UPDIS::P_UPDIS.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_P_UPDIS"],
"r2_FRT_SURFACE::FRT_SURFACE.mgt|change=absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_FRT_SURFACE"],
"r2_FRT_KG::FRT_KG.mgt|change=pctchg |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_FRT_KG"]*100,
"r2_NPERCO::NPERCO.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_NPERCO"],
"r2_ERORGN::ERORGN.hru|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_ERORGN"],
"r1_PPERCO::PPERCO.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_PPERCO"],
"r2_PPERCO::PPERCO.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_PPERCO"],
"r1_CDN::CDN.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_CDN"],
"r2_CDN::CDN.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_CDN"],
"r1_CMN::CMN.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_CMN"],
"r2_CMN::CMN.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_CMN"],
"r1_PSP::PSP.bsn|change= absval |
sub == c(1,8,10,14,15,16,17,18,19,20,22)"= mypar[,"r1_PSP"],
"r2_PSP::PSP.bsn|change= absval |
sub == c(2,3,4,5,6,7,9,11,12,13,21,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41)"=
mypar[,"r2_PSP"],
"FILTERW_A::FILTERW.mgt|change = absval |
sub == c(1,7,8,10,14,15,16,17,18,19,20,22,24,28)" = as.vector(rep(as.numeric(escenarios$FW_A[n]),nrow(mypar))),
"FILTERW_B::FILTERW.mgt|change = absval |
sub == c(4,5,25,31,3,27,2,23,29)" = as.vector(rep(as.numeric(escenarios$FW_B[n]),nrow(mypar))),
"FILTERW_C::FILTERW.mgt|change = absval |
sub == c(41,12,21,30,9,13,32,38,36,34,6,11,37,35,39,40,26,33)" = as.vector(rep(as.numeric(escenarios$FW_C[n]),nrow(mypar)))
)
# ejecucion del modelo
buffer_run = run_swat2012(project_path = "./SWAT_SL_20200401/",
output = list(FLOW = define_output(file = "rch",
variable = "FLOW_OUT",
unit = 1:41),
PT = define_output(file = "rch",
variable = "TOT_Pkg",
unit = 1:41)),
parameter = par_set,
start_date = as.Date("2005-01-01"),
end_date = as.Date("2015-12-31"),
years_skip = 1,
add_parameter = FALSE,
n_thread = 12)
# calculo de la concentración de fósforo:
flow = buffer_run[1:41]
pt = buffer_run[42:82]
# Concentración de fósforo por simulación
fun1 = function(x,y){
out = 1000*unlist(x)/(unlist(y)*60*60*24)
out = ifelse(is.finite(out),out,NA)
}
PTcon = list()
for(i in 1:41){
pt_rch = pt[[i]][,-1]
flow_rch = flow[[i]][,-1]
PTcon[[i]] = mapply(fun1, pt_rch, flow_rch)
}
# ECDF de concentración de fósforo
x_con = lseq(0.0002,2,length.out = 1000)
fun2 = function(x){apply(x,2, function(x){ecdf(x)(x_con)})}
PT_ecdf = lapply(PTcon, fun2)
# extremos y medio de ECDF
PT_ecdf_max = lapply(PT_ecdf, function(x){round(apply(x,1,max),4)})
PT_ecdf_min = lapply(PT_ecdf, function(x){round(apply(x,1,min),4)})
PT_ecdf_q50 = lapply(PT_ecdf, function(x){round(apply(x,1,function(x){quantile(x,0.5)}),4)})
PT_ecdf_q25 = lapply(PT_ecdf, function(x){round(apply(x,1,function(x){quantile(x,0.25)}),4)})
PT_ecdf_q75 = lapply(PT_ecdf, function(x){round(apply(x,1,function(x){quantile(x,0.75)}),4)})
PT_ecdf_out = list()
PT_ecdf_out[["PTcon"]] = x_con
for(i in 1:41){
PT_ecdf_out[[paste0("rch_", str_pad(i, width = 2, pad="0"))]] = data.frame(q50=PT_ecdf_q50[[i]],
q25=PT_ecdf_q25[[i]],
q75=PT_ecdf_q75[[i]],
qmin = PT_ecdf_min[[i]],
qmax = PT_ecdf_max[[i]])
}
saveRDS(PT_ecdf_out,paste0("./output_buffer/PTcon_ecdf",
"_A",paste0(escenarios[n,"FW_A"]),
"_B",paste0(escenarios[n,"FW_B"]),
"_C",paste0(escenarios[n,"FW_C"]),
".RDS"))
}A work by Rafael Navas