R/functions.R
In drunita/mr-model: Risk prediction model Mr
Defines functions
weather_data_prototipo
weather_data_NC
LW_satate
LWD_daily_estimation
LWD_estimation
weather_data
ploting_the_Risk
ploting_the_disease
barplot_model_results
plots_model_results
running_simulation
Estimating_control_action_days
daily_weather
harvest_remotion
risk_factor
prob_infect
Documented in
barplot_model_results
daily_weather
Estimating_control_action_days
harvest_remotion
LWD_daily_estimation
LWD_estimation
LW_satate
ploting_the_disease
ploting_the_Risk
plots_model_results
prob_infect
risk_factor
running_simulation
weather_data
weather_data_NC
weather_data_prototipo
#' Infection progress function
#'
#' This function is the main function and uses all the other function to estimates
#' the infection progress for each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param Temp number given by daily_weather function
#' @param Hum number given by daily_weather function
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param modeling_days increasing number
#' @param R factor given by risk_factor function
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return the new stimates for mature,sporulated,healthy,latent,modeling_days,harvest,latent_removed,sporulated_removed,p ,R,pod_data,control_daysfactor
#' @export
modelo_pod_infection <- function (pod_data,plot_data , Temp , WH , x, mature, sporulated,modeling_days,R, y){
control_days <- Estimating_control_action_days(x,y,modeling_days)
p<-prob_infect( x["cultivar" ], Temp, WH,plot_data,modeling_days,control_days)
pod_data["maturity"] <- growth_factor(Temp , pod_data)
pod_data["Healthy"] <- new_healthy(p , pod_data)
pod_data["sporulate_level"] <- sporulation_grade(Temp , pod_data)
#############
new_pods<-new_pods_function(x,modeling_days)
pod_data <- new_data_functtion(p,pod_data , new_pods )
mature_sporulated <-mature_sporulated_healty_latent_pod(pod_data,mature,sporulated,modeling_days)
mature <- as.numeric(mature_sporulated[1])
sporulated<- as.numeric(mature_sporulated[2])
modeling_days<- as.numeric(mature_sporulated[3])
pod_data <- as.data.frame(mature_sporulated[4])
harvest_data<-harvest_remotion(pod_data , mature , sporulated, modeling_days , x, control_days)
mature <- as.numeric(harvest_data[1])
sporulated<- as.numeric(harvest_data[2])
harvest<-as.numeric(harvest_data[3])
healthy<- as.numeric(harvest_data[4])
latent <- as.numeric(harvest_data[5])
latent_removed<- as.numeric(harvest_data[6])
sporulated_removed<- as.numeric(harvest_data[7])
pod_data <- as.data.frame(harvest_data[8])
R<-risk_factor(p,R,pod_data,healthy,sporulated,latent)
return( list(mature,sporulated,healthy,latent,modeling_days,harvest,latent_removed,sporulated_removed,p ,R,pod_data,control_days))
}
#' Growing pod_data file function
#'
#' This function adds a row at the end of the pod_data data frame with the new information
#' of the simulated day
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @param new_pods number specified by the user in the crop characteristic dataframe
#' @return the pod_data data frame with the new row
#' @export
new_data_functtion <- function (p,pod_data,new_pods){
pod_data <- pod_data %>% add_row(Healthy = new_pods, maturity = 0,Disease = new_infected (p, pod_data),sporulate_level=0,.before = 1)
return(pod_data)
}
#' new pods function
#'
#' This function calculates the number of pods entering in the crop base on the initiation of the harvest day
#' asuming trees are more active in the first two month. Then, tree are less active until the thir month and do not produce new pods
#'
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param modeling_days increasing number
#' @return the number of new pods that enter to the crop that day
#' @export
new_pods_function<-function (x,modeling_days){
if(modeling_days<60){return( round(as.numeric(x["new_pods" ])))}
if(modeling_days>=60&modeling_days<90){return( round(as.numeric(x["new_pods" ])/2))}
return(round(as.numeric(x["new_pods" ])/3))
}
#' Grow factor function
#'
#' This function calculates the maturity status of the pods (grow factor)
#' for each simulation day using the Zuidema et al 2003 model
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param Temp number given by daily_weather function
#' @return a vector with the new maturity status for the pods in the crop which can be added to the pod_data data frame
#' @export
growth_factor <- function (Temp,pod_data){
temperature_growth_factor <- c(19.99 , 20.98 , 21.97 , 22.96 , 23.99 , 24.48 , 24.95 , 25.97 , 26.96 , 27.97)
index_growth_factor <- c(0.0049 , 0.0053 , .0056 , .006 , .0064 , .0066 , .0067 , .0069 , .0071 , .0073)
mod <- lm(index_growth_factor ~ temperature_growth_factor)
gi <- predict(mod,newdata = data.frame(temperature_growth_factor = Temp))
growing <- function(x){x+gi}
a <- map(pod_data["maturity"],growing)
return(a)
}
#' Probability of Infection function
#'
#' This function estimates the probability of a new infection to develop
#' for each simulation day.
#'
#' @param cultivar character specified by the user in the crop characteristic dataframe
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param modeling_days increasing number
#' @param Temp number given by daily_weather function
#' @param Hum number given by daily_weather function
#' @param control_days list with a vector and a number given by Estimating_control_action_days function
#' @return the probability of a new infection the simulated day
#' @export
prob_infect<- function(cultivar, Temp, WH,plot_data,modeling_days,control_days){
incidence<-tail(plot_data$Sporulated, n=1)/(tail(plot_data$Healthy, n=1)+tail(plot_data$Sporulated, n=1)+tail(plot_data$Latent, n=1))
cd<-as.numeric(control_days[2])
if (cd <15){return(0)}### asuminedo que fumigacion protege mazorcas 15 dias
if (WH >8){
if (incidence>0.1){return(0.05)} #incidence of 60% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time) when no cultural sanitation is made
else{return(0.02)}#incidence of 40% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
}
if (incidence>0.1){return(0.01)}#incidence of 20% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
else{return(0)}#incidence of 5% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
}
#' Risk factor function
#'
#' This function estimates the risk that a new epidemic develops for each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param healthy number given by mature_sporulated_healty_latent_pod function
#' @param latent number given by mature_sporulated_healty_latent_pod function
#' @param R factor given by risk_factor function for the day before
#' @param p probability of a new infection given by prob_infect function
#' @return the risk of epidemy development
#' @export
risk_factor<-function(p,R,pod_data,healthy,sporulated,latent){
#young_pod<-sum(pod_data$Healthy[pod_data$maturity<0.33])/(healthy+sporulated+latent)
#daily_risk<-p*young_pod
#if(p==0.001){R<-R+daily_risk}
if(p==0.05){R<-R+10}
else if (p==0.02) {R<-R+8}
else if (p==0.01) {R<-R+5}
else {R<-R-5}
if(R<0){R<-0}
return(R)
}
#' new infected function
#'
#' This function estimates the number of pod infected in each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @return a number showing how many pods got infected in each day
#' @export
new_infected <- function (p, pod_data) {
new_inf <- sum(pod_data$Healthy[pod_data["maturity"]<0.6]*p)
return(new_inf)
}
#' new healthy function
#'
#' This function estimates the number of pod that remained healthy in each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @return a number showing how many pods remained infected in each day
#' @export
new_healthy <- function (p , pod_data) {
a<-pod_data$Healthy[pod_data["maturity"]<0.6]*(1-p)
b <- pod_data$Healthy[pod_data["maturity"]>=0.6]
c<-c(a,b)
return(c)
}
#' Sporulation factor function
#'
#' This function calculates the sporulation status of the pods (grow factor)
#' for each simulation day
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param Temp number given by daily_weather function
#' @return a vector with the new sporulation status for the pods in the crop which can be added to the pod_data data frame
#' @export
#'
sporulation_grade <- function (Temp,pod_data){
si <- 1/65
sporulating <- function(x){x+si}
a <- map(pod_data["sporulate_level"],sporulating)
return(a)
}
#' new mature and sporulated function
#'
#' This function estimates the number of pod that grew to maturity and sporulation
#' in each simulation day. Then uses this information to actualize the pod_data
#' data file (remove last rows in which the pods completed growth and sporulation)
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param modeling_days increasing number
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @return a list with numbers representing the new mature and sporulating pods, the increased modeling day and the updated pod_data
#' @export
#'
mature_sporulated_healty_latent_pod <- function (pod_data , mature , sporulated , modeling_days) {
if (length(subset(pod_data,pod_data["maturity"]>=1)["Healthy"][,1])>0){
mature<-mature+sum(subset(pod_data,pod_data["maturity"]>=1)["Healthy"])
pod_data[pod_data$maturity>=1,]["Healthy"]<-0
}
if (length(subset(pod_data,pod_data["sporulate_level"]>=1)["Disease"][,1]) > 0){
sporulated<-sporulated+sum(subset(pod_data,pod_data["sporulate_level"]>=1)["Disease"])
pod_data[pod_data$sporulate_level>=1,]["Disease"]<-0}
rows_to_drop<-pod_data$Healthy==0 & pod_data$Disease==0
pod_data<-pod_data[!rows_to_drop,]
modeling_days<-modeling_days+1
return( list(mature,sporulated, modeling_days , pod_data))
}
#' pods harvested and removed function
#'
#' This function estimates the number of healthy pods that were harvested
#' or diseased pods removed from the crop in each simulation day.Then uses this
#' information to actualize the pod_data and the pods with each disease status
#' remaining in the status afetr the remotions
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param modeling_days increasing number
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param control_days list with a vector and a number given by Estimating_control_action_days function
#' @return a list with numbers representing the new mature,sporulated,harvest,healthy,latent,disease_remove, sporulated_remove and the updated pod_data
#' @export
#'
harvest_remotion <- function(pod_data , mature , sporulated, modeling_days , x , control_days){
if (modeling_days %% as.numeric(x["croping_frequency"]) == 0) {
harvest<-mature
mature<-0}else{harvest<-0}
if (control_days[1]==0 & length(pod_data["Disease"])!=0) {
disease_remove<-sum(pod_data["Disease"])*as.numeric(x["removing_latent_eficacy"])
sporulated_remove<-sporulated*as.numeric(x["removing_sporulated_eficacy"])
sporulated<-sporulated*(1-as.numeric(x["removing_sporulated_eficacy"]))
pod_data["Disease"] <- pod_data["Disease"]*(1-as.numeric(x["removing_latent_eficacy"]))
}else{
disease_remove<-0
sporulated_remove<-0
}
if (length(pod_data["Healthy"])!=0){
healthy<- sum(pod_data["Healthy"])
}else(healthy<-0)
if (length(pod_data["Disease"])!=0){
latent<- sum(pod_data["Disease"])
}else(latent<-0)
return(list(mature,sporulated,harvest,healthy,latent,disease_remove, sporulated_remove,pod_data))
}
#' daily weather function
#'
#' This function estimates the temperature and presipotation for each day.
#'
#' @param md data frame with the wheather
#' @param modeling_days increasing number
#' @return a list with numbers representing the day temperature, presipitation, RH and hours of wetennes
#' @export
#'
daily_weather<-function(md,modeling_days){
Temp<-as.numeric(md$Temperature[modeling_days])
Hum<-as.numeric(md$Humidity[modeling_days])
Rain<-as.numeric(md$Rain[modeling_days])
WH_CART<-as.numeric(md$CART_WH[modeling_days])
WH_TH<-as.numeric(md$TH_WH[modeling_days])
return(list(Temp,Hum,WH_CART,WH_TH,Rain))
}
#' Control action function
#'
#' This function estimates the day when control actions (remotios and funmigations)
#' were done. base on the information provided by the user
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param modeling_days increasing number
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return a list especifying the days in which remotions were done and the last fumigation day
#' @export
Estimating_control_action_days<-function(x,y,modeling_days){
if (length(na.omit(y$Fumigation_dates))==0 & length(na.omit(y$Remotion_dates))==0){
return(list(100,100))
}
if (length(na.omit(y$Fumigation_dates))!=0 & length(na.omit(y$Remotion_dates))!=0){
Days_fumigation<-as.numeric(as.Date(as.character(y$Fumigation_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
Days_remotion<-as.numeric(as.Date(as.character(y$Remotion_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_fumigation[1]>modeling_days & Days_remotion[1]>modeling_days){
return(list(100,100))
}
if (Days_fumigation[1]<=modeling_days){
last_fumigation<- max(na.omit(Days_fumigation[Days_fumigation<=modeling_days]))
Days_from_last_fumigation<-modeling_days-last_fumigation
}else{Days_from_last_fumigation<-100}
if (Days_remotion[1]<=modeling_days){
last_remotion<- max(na.omit(Days_remotion[Days_remotion<=modeling_days]))
Days_from_last_remotion<-modeling_days-last_remotion
}else{Days_from_last_remotion<-100}
return(list(as.numeric(Days_from_last_remotion),as.numeric(Days_from_last_fumigation)))
}
if (length(na.omit(y$Fumigation_dates))!=0){
Days_fumigation<-as.numeric(as.Date(as.character(y$Fumigation_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_fumigation[1]>modeling_days){
return(list(100,100))
}
if (Days_fumigation[1]<=modeling_days){
last_fumigation<- max(na.omit(Days_fumigation[Days_fumigation<=modeling_days]))
Days_from_last_fumigation<-modeling_days-last_fumigation
return(list(100,as.numeric(Days_from_last_fumigation)))
}
}
if (length(na.omit(y$Remotion_dates))!=0){
Days_remotion<-as.numeric(as.Date(as.character(y$Remotion_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_remotion[1]>modeling_days){
return(list(100,100))
}
if (Days_remotion[1]<=modeling_days){
last_remotion<- max(na.omit(Days_remotion[Days_remotion<=modeling_days]))
Days_from_last_remotion<-modeling_days-last_remotion
return(list(as.numeric(Days_from_last_remotion),100))
}}
}
#' Running the simulation function
#'
#' This function runs the simulation for a specified number of days
#' using the modelo_pod_infection
#' @param weather_data data frame with the wheather
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param modeling_days increasing number
#' @param R factor given by risk_factor function
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return the plot_data with the results of the simulation
#' @export
#'
running_simulation<-function(weather_data,modeling_days,pod_data,plot_data,x,mature, sporulated, R,y,model_WH ){
for (i in 1:length(weather_data$Temperature)) {
Temp<-as.numeric(daily_weather(weather_data,modeling_days)[1])
Hum<-as.numeric(daily_weather(weather_data,modeling_days)[2])
if (model_WH=="CART"){WH<-as.numeric(daily_weather(weather_data,modeling_days)[3])}
else if (model_WH=="TH"){ WH<-as.numeric(daily_weather(weather_data,modeling_days)[4])}
##############
##########
else if (model_WH=="WH prototipo"){ WH<-
as.numeric(daily_weather(weather_data,modeling_days)[4])}
##########
############
Rain<-as.numeric(daily_weather(weather_data,modeling_days)[5])
a<-modelo_pod_infection(pod_data,plot_data,Temp,WH,x,mature, sporulated , modeling_days, R,y)
mature<-round(as.numeric(a[1]))
sporulated<-round(as.numeric(a[2]))
healthy<-round(as.numeric(a[3]))
latent<-round(as.numeric(a[4]))
modeling_days<-as.numeric(a[5])
harvest<-round(as.numeric(a[6]))
latent_removed<- round(as.numeric(a[7]))
sporulated_removed<-round(as.numeric(a[8]))
p<-round(as.numeric(a[9]),3)
R<-round(as.numeric(a[10]),3)
pod_data<-as.data.frame(a[11])
plot_data <- plot_data %>% add_row(Days = modeling_days, Healthy=healthy,Mature = mature,
Latent = latent,Sporulated=sporulated,Harvest=harvest ,
Latent_removed= latent_removed, Sporulated_removed=sporulated_removed,
Risk_factor= R)
}
return(plot_data)
}
#' line plot function
#'
#' This function generates a line plot showing represent the simulation results
#'
#' @param x a numeric vector especifing the x factor
#' @param y a numeric vector especifing the y factor
#' @param name_x a character specifing the x axis name
#' @param name_y a character specifing the y axis name
#' @param group a character vector specifying the grouping factor
#' @param group_name a character specifying the group factor name
#' @param levels_order a character vector specifying the grouping factor order for the legend
#' @param y_max a number especifying the maximum of the y axis
#' @param y_min a number especifying the manimum of the y axis
#' @param title_legend a character specifying the legend title
#' @return a line plot
#' @export
#' @examples
#' y<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
#' x<-rep(plot_data$Days,3)
#' n<-length(plot_data$Days)
#' group<-c(rep("Healthy",n),rep("Latent",n),rep("Sporulated",n))
#' plots_model_results(x,"Days",y,"Pods proportion",group,"Disease status",c("Healthy","Latent","Sporulated"),0,1, "Disease status")
plots_model_results<-function(x,name_x,y,name_y,group,group_name,levels_order,y_min,y_max,title_legend){
plot_data<-data.frame(name_x=x,name_y=y,group_name=group)
p<-ggplot(plot_data,aes(as.numeric(x),as.numeric(y),fill = factor(group, levels=levels_order)))
p<- p + theme_classic()
p<- p + geom_line(aes(color=group),size=1.1)+ scale_color_manual(values=c("#69b3a2", "#56B4E9", "#E69F00"))
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + xlab(paste("\n",name_x,sep="")) + ylab(paste(name_y,"\n",sep=""))+ylim(y_min, y_max)
p<- p +theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(0.5, 1.1, 0.5, 1.3, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
p <- p + labs(color=title_legend)
p<- p+ theme(legend.text = element_text(size=12, angle=0),
legend.title = element_text(face="bold",size=12, angle=0))
return(p)
}
#' bar plot function
#'
#' This function generates a bar plot showing represent the simulation results
#'
#' @param x a numeric vector especifing the x factor
#' @param y a numeric vector especifing the y factor
#' @param name_x a character specifing the x axis name
#' @param name_y a character specifing the y axis name
#' @param group a character vector specifying the grouping factor
#' @param group_name a character specifying the group factor name
#' @param levels_order a character vector specifying the grouping factor order for the legend
#' @param y_max a number especifying the maximum of the y axis
#' @param y_min a number especifying the manimum of the y axis
#' @param title_legend a character specifying the legend title
#' @return a bar plot
#' @export
#' @examples
#' y<-c(sum(plot_data$Harvest),sum(plot_data$Latent_removed),sum(plot_data$Sporulated_removed))
#' x<-c("Healthy","Latent","Sporulated")
#' group<-x
#' p4<-barplot_model_results(x,"Disease status",y,"Total pods harvest",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
barplot_model_results<-function(x,name_x,y,name_y,group,group_name,levels_order,y_min,y_max,title_legend){
plot_data<-data.frame(name_x=x,name_y=y,group_name=group)
p<-ggplot(plot_data,aes(factor(x),as.numeric(y),fill = factor(group, levels=levels_order)))
p<- p + theme_classic()
p<- p + geom_bar(stat="identity")+ scale_fill_manual(values=c("#69b3a2", "#56B4E9", "#E69F00"))
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + xlab(paste("\n",name_x,sep="")) + ylab(paste(name_y,"\n",sep=""))+ylim(y_min, y_max)
p<- p +theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(0.5, 1.1, 0.5, 1.3, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
#p <- p + labs(fill=title_legend)
p<- p+ theme(legend.position="none")
return(p)
}
#' Plot the simulation function
#'
#' This function plots the results of the simulation
#' @param plot_data final .csv given by running_simulation function
#' @return A plot with for panels showing diseases progres, the amount of pods harves and remove total and during the simulation
#' @export
#'
ploting_the_disease<-function(plot_data){
y<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
x<-rep(plot_data$Days,3)
n<-length(plot_data$Days)
group<-c(rep("Healthy",n),rep("Latent",n),rep("Sporulated",n))
p1<-plots_model_results(x,"Days",y,"Pods",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease status")
y1<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
tot<-plot_data$Healthy+plot_data$Latent+plot_data$Sporulated
y<-y1/rep(tot,3)
p2<-plots_model_results(x,"Days",y,"Pods proportion",group,"Disease status",c("Healthy","Latent","Sporulated"),0,1, "Disease status")
y<-c(plot_data$Harvest,plot_data$Latent_removed,plot_data$Sporulated_removed)
p3<-plots_model_results(x,"Days",y,"Pods harvest",group,"Pods harvest",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
y<-c(sum(plot_data$Harvest),sum(plot_data$Latent_removed),sum(plot_data$Sporulated_removed))
x<-c("Healthy","Latent","Sporulated")
group<-x
p4<-barplot_model_results(x,"Disease status",y,"Total pods harvest",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
require(gridExtra)
plot<-grid.arrange(p1,p2,p3,p4, ncol=2,nrow=2)
return(plot)
}
#' Plot risk function
#'
#' This function plots the evolution of the risk factoer and incidence of sporulated
#' pods during the simulation and gives control resomentations
#'
#' @param plot_data final .csv given by running_simulation function
#' @return A plot with for panels showing risk progres and control recomendations
#' @export
#'
ploting_the_Risk<-function(plot_data){
df<-data.frame("Inoculum_pressure"=plot_data$Sporulated/(plot_data$Healthy+plot_data$Latent+plot_data$Sporulated),
"New_infection"=plot_data$Risk_factor, "Days"= plot_data$Days)
scaleFactor <-1
if(max(df$New_infection)!=0 & max(df$Inoculum_pressure)!=0){
scaleFactor <- max(df$Inoculum_pressure) / max(df$New_infection)}
p <- ggplot(df, aes(x=Days)) +theme_classic()
p<- p + geom_line(aes(y=Inoculum_pressure), col="#E69F00",size=1.1) +
geom_line(aes(y=New_infection * scaleFactor), col="#56B4E9",size=1.1) +
scale_y_continuous(name="Inoculum Pressure Risk", sec.axis=sec_axis(~./scaleFactor, name="New Infection Risk"))
p<- p+theme(axis.title.y.left=element_text(color="#E69F00"),axis.title.y.right=element_text(color="#56B4E9"))
if(tail(df,1)$Inoculum_pressure>0.095){
p <- p +geom_line(aes(y = 0.1), color = "#E69F00", linetype = "dotted",size=1.5)
}
p <- p +geom_line(aes(y = 30*scaleFactor), color = "#56B4E9", linetype = "dotted",size=1.5)
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(1.6, 1.5,0.8, 1.5, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
p<- p+ theme(legend.text = element_text(size=12, angle=0),
legend.title = element_text(face="bold",size=12, angle=0))
if(tail(df,1)$New_infection>30 & tail(df,1)$Inoculum_pressure<0.1){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$New_infection*scaleFactor,label="Recomiendaciones:\n1. Fumigar", hjust = "inward", vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure>0.1 & tail(df,1)$New_infection<30){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$Inoculum_pressure,label="Recomiendaciones:\n1. Hacer remocion", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure>0.1 & tail(df,1)$New_infection>30){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$Inoculum_pressure,label="Recomiendaciones:\n1. Hacer remocion\n2. Fumigar", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure<0.1 & tail(df,1)$New_infection<30){
p<- p + annotate(geom="text",x=0, y=30*scaleFactor,label="Recomiendaciones:\n1. No accion de control necesaria", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
p
return(p)
}
#' Getting the wheather from the IDEAM
#'
#' This function fetches data from the IDAEM API
#'
#' @param location the zipcode of the city for wich the data is required"=,"Remotion_dates"=)
#' @return The data frame with the weather
#' @export
#'
weather_data <- function(location){
####Inspect/Network/Headers
###location = "68001000", para bucaramanga
get_weather<-GET(paste("http://mipronostico.ideam.gov.co/IdeamWebApp2/Ideam/getData/reverse/reverse.php?cod=",location, sep = ""))
get_weather_text<-content(get_weather, "text")
get_weather_json<-fromJSON(get_weather_text,flatten = TRUE)
df <- as.data.frame(get_weather_json)
get_weather_df<-data.frame(Fecha=df$pronostico.fecha, Hora=df$pronostico.hora, Temperatura=df$pronostico.temperature,Humedad=df$pronostico.humidity,
DewPoint=df$pronostico.dewPoint,Nubosidad=df$pronostico.totalCloudCover,Presipitacion=df$pronostico.precipitation1h,wind_speed=df$pronostico.windSpeedMS)
return(get_weather_df)
}
#' Daily wet hours for the simulation dates
#'
#' This function uses the results from the LWD_daily_estimation function to generate a data frame indicating the number of hours
#' that the cocoa pod has a water film on the surface (WH) and the longest period of dry hours for the day (maxDLH)
#'
#' @param df a Data frame with the wheather parameter for each hour of each day. This df is provided by the fuction whather data
#' @return a data frame with the date, the CART_WH, TH_WH, maxDH, Temperature, Humidity and Rain
#' @export
#'
LWD_estimation <- function(df1) {
days<-names(table(df1$Fecha))
days <- format(days, format="%Y-%m-%d")
LWD_CART <- rep(0, length(days))
LWD_TH <- rep(0, length(days))
max_D <- rep(0, length(days))
Temp <- rep(0, length(days))
Hum <- rep(0, length(days))
Rain <- rep(0, length(days))
for (i in 1:length(days)) {
df <- subset(df1, df1$Fecha == days[i])
val <- LWD_daily_estimation(df)
LWD_CART[i] <- as.numeric(val[1])
LWD_TH[i]<- as.numeric(val[2])
Temp[i]<- signif(mean(as.numeric(df$Temperatura)),3)
Hum[i]<-signif(mean(as.numeric(df$Humedad)),3)
Rain[i]<-sum(as.numeric(df$Presipitacion))
}
df2 <- data.frame(Fecha = days ,
CART_WH = LWD_CART,
TH_WH= LWD_TH,
Temperature=Temp,
Humidity=Hum,
Rain=Rain)
return(df2)
}
#' Daily wet hours
#'
#' This function calculates the number of hours that the cocoa pod has a water film on the surface on a day and the longest period of dry hours for the day using the results from the
#' LW_satate function
#'
#' @param df a Data frame with the wheather parameter for each hour of each day. This df is provided by the fuction whather data
#' @return a list of three elements with the wet hours and and the longest period of dry hours for the day acording to the CART model and the wet hours acording to TH model
#' @export
#'
LWD_daily_estimation <- function(df) {
DPT <- as.numeric(df$DewPoint)
RH <- as.numeric(df$Humedad)
Temp <- as.numeric(df$Temperatura)
wind <- as.numeric(df$wind_speed)
wet_hr_CART <- rep(0, length(DPT))
wet_hr_TH <- rep(0, length(DPT))
lws<- LW_satate(DPT[1], RH[1], Temp[1], wind[1],RH[1])
wet_hr_CART[1] <- as.numeric(lws[1])
wet_hr_TH[1]<-as.numeric(lws[2])
for (i in 2:length(DPT)) {
lws<- LW_satate(DPT[i], RH[i], Temp[i], wind[i],RH[i-1])
wet_hr_CART[i] <- as.numeric(lws[1])
wet_hr_TH[i]<-as.numeric(lws[2])
}
#max(df$Temperature) & min(df$Temperature)
return(list(sum(wet_hr_CART),sum(wet_hr_TH)))
}
#' Wetness state of the cocoa pod
#'
#' This function uses the CART and TH models (Kim et al.,2002; Bregaglio; 2011) to
#' estimate whether the cocoa pod has a water film on the surface depending on the whather parameters
#'
#' @param DPT dewpoint
#' @param RH relative humidity
#' @param Temp Temperature
#' @param Wind Wind speed_km_h
#' @return a boolean list of two elements showing the results of the CART and TH models with 0:dry, 1:wet
#' @export
#'
LW_satate <- function(DPT, RH, Temp, wind, RH_ant) {
DPD <- Temp - DPT
In1 <-1.6064 * sqrt(Temp) + 0.0036 * Temp ^ 2 + 0.1531 * RH - 0.4599 * wind *DPD - 0.0035 * Temp * RH
In2 <-0.7921 * sqrt(Temp) + 0.0046 * RH - 2.3889 * wind - 0.039 * Temp * wind +1.061 * wind * DPD
wet_CART<-0
wet_TH<-0
if (RH <= 70){wet_TH <- 0} else if(RH >= 87){wet_TH<-1} else {
if (RH-RH_ant > 2.5){wet_TH<-1}else{wet_TH<-0}
}
if (DPD >= 3.7) {
return(list(wet_CART,wet_TH))
}
if (wind < 2.5) {
return(list(as.numeric(In1 > 14.4674),wet_TH)) # when false return 0 (dry variable)
}
if (RH <= 87) {
return(list(as.numeric(In2 > 37),wet_TH)) # when false return 0 (dry variable)
}
return(list(wet_CART,wet_TH))
}
#'
#' Getting the wheather from the NC
#'
#' This function fetches data from the IDAEM API
#'
#' @param q the query for the weather api in qls; q = 'SELECT mean("temp") as temp,mean("hum") as hum, difference(mean("rain")) as rain, mean("windv") as wind FROM /^DL00011$/ WHERE time >= now() - 24d GROUP BY time(1h) fill(none)'
#' @return The data frame with the weather
#' @export
#'
weather_data_NC <- function(q, user, password){
encoded <- URLencode(q)
get_weather<-GET(paste("https://login.daimob.co/api/datasources/proxy/4/query?db=agro&q=", encoded, sep = ""),
authenticate(user, password))
get_weather_text<-content(get_weather, "text")
get_weather_json<-fromJSON(get_weather_text,flatten = TRUE)
colum.names<-unlist(get_weather_json$results$series[[1]]$columns)
df1<-as.data.frame(get_weather_json$results$series[[1]]$values)
colnames(df1)<-colum.names
df<-data.frame(Fecha=format(as.Date(df1$time,format="%Y-%m-%d")), Hora=format(parse_datetime(df1$time), format="%H:%M:%S", tz="America/Bogota"), Temperatura=as.numeric(df1$temp),Humedad=as.numeric(df1$hum),
DewPoint=as.numeric(df1$temp)-(100-as.numeric(df1$hum))/5,Presipitacion=as.numeric(df1$rain),wind_speed=as.numeric(df1$wind))
return(na.omit(df))
}
#' Getting the wheather from the prototype
#'
#' This function gets the weather from prototype
#'
#' @param df a data frame of the prototype measurements
#' @return The data frame with the weather
#' @export
#'
weather_data_prototipo <- function(df){
df1<-data.frame(Fecha=format(as.Date(df$Fecha,format="%m/%d/%Y")), Hora=format(parse_datetime(paste(df$Fecha, df$Hora)," %m/%d/%Y %I:%M:%S %p"),format="%I"),Temperatura=as.numeric(df$Temperatura),Humedad=as.numeric(df$Humedad.relativa),
DewPoint=as.numeric(df$Temperatura)-(100-as.numeric(df$Humedad.relativa))/5,Presipitacion=0,wind_speed=0)
df2<-subset( df1, df1$Hora==names(table(df1$Hora))[1]&df1$Fecha==names(table(df1$Fecha))[1])
df3<-data.frame(Fecha=df2$Fecha[1], Hora=df2$Hora[1],Temperatura=mean(df2$Temperatura),Humedad=mean(df2$Humedad),
DewPoint=mean(df2$DewPoint),Presipitacion=sum(df2$Presipitacion),wind_speed=mean(df2$wind_speed))
for(j in 1:length(names(table(df1$Hora)))){
for(i in 2:length(names(table(df1$Hora)))){
df2<-subset( df1, df1$Hora==names(table(df1$Hora))[i]&df1$Fecha==names(table(df1$Fecha))[j])
df3 <- df3 %>% add_row(Fecha=df2$Fecha[1], Hora=df2$Hora[1],Temperatura=mean(df2$Temperatura),Humedad=mean(df2$Humedad),
DewPoint=mean(df2$DewPoint),Presipitacion=sum(df2$Presipitacion),wind_speed=mean(df2$wind_speed))
}}
return(na.omit(df3))
}
drunita/mr-model documentation built on Dec. 20, 2021, 1:20 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions weather_data_prototipo weather_data_NC LW_satate LWD_daily_estimation LWD_estimation weather_data ploting_the_Risk ploting_the_disease barplot_model_results plots_model_results running_simulation Estimating_control_action_days daily_weather harvest_remotion risk_factor prob_infect
Documented in barplot_model_results daily_weather Estimating_control_action_days harvest_remotion LWD_daily_estimation LWD_estimation LW_satate ploting_the_disease ploting_the_Risk plots_model_results prob_infect risk_factor running_simulation weather_data weather_data_NC weather_data_prototipo
#' Infection progress function
#'
#' This function is the main function and uses all the other function to estimates
#' the infection progress for each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param Temp number given by daily_weather function
#' @param Hum number given by daily_weather function
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param modeling_days increasing number
#' @param R factor given by risk_factor function
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return the new stimates for mature,sporulated,healthy,latent,modeling_days,harvest,latent_removed,sporulated_removed,p ,R,pod_data,control_daysfactor
#' @export
modelo_pod_infection <- function (pod_data,plot_data , Temp , WH , x, mature, sporulated,modeling_days,R, y){
control_days <- Estimating_control_action_days(x,y,modeling_days)
p<-prob_infect( x["cultivar" ], Temp, WH,plot_data,modeling_days,control_days)
pod_data["maturity"] <- growth_factor(Temp , pod_data)
pod_data["Healthy"] <- new_healthy(p , pod_data)
pod_data["sporulate_level"] <- sporulation_grade(Temp , pod_data)
#############
new_pods<-new_pods_function(x,modeling_days)
pod_data <- new_data_functtion(p,pod_data , new_pods )
mature_sporulated <-mature_sporulated_healty_latent_pod(pod_data,mature,sporulated,modeling_days)
mature <- as.numeric(mature_sporulated[1])
sporulated<- as.numeric(mature_sporulated[2])
modeling_days<- as.numeric(mature_sporulated[3])
pod_data <- as.data.frame(mature_sporulated[4])
harvest_data<-harvest_remotion(pod_data , mature , sporulated, modeling_days , x, control_days)
mature <- as.numeric(harvest_data[1])
sporulated<- as.numeric(harvest_data[2])
harvest<-as.numeric(harvest_data[3])
healthy<- as.numeric(harvest_data[4])
latent <- as.numeric(harvest_data[5])
latent_removed<- as.numeric(harvest_data[6])
sporulated_removed<- as.numeric(harvest_data[7])
pod_data <- as.data.frame(harvest_data[8])
R<-risk_factor(p,R,pod_data,healthy,sporulated,latent)
return( list(mature,sporulated,healthy,latent,modeling_days,harvest,latent_removed,sporulated_removed,p ,R,pod_data,control_days))
}
#' Growing pod_data file function
#'
#' This function adds a row at the end of the pod_data data frame with the new information
#' of the simulated day
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @param new_pods number specified by the user in the crop characteristic dataframe
#' @return the pod_data data frame with the new row
#' @export
new_data_functtion <- function (p,pod_data,new_pods){
pod_data <- pod_data %>% add_row(Healthy = new_pods, maturity = 0,Disease = new_infected (p, pod_data),sporulate_level=0,.before = 1)
return(pod_data)
}
#' new pods function
#'
#' This function calculates the number of pods entering in the crop base on the initiation of the harvest day
#' asuming trees are more active in the first two month. Then, tree are less active until the thir month and do not produce new pods
#'
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param modeling_days increasing number
#' @return the number of new pods that enter to the crop that day
#' @export
new_pods_function<-function (x,modeling_days){
if(modeling_days<60){return( round(as.numeric(x["new_pods" ])))}
if(modeling_days>=60&modeling_days<90){return( round(as.numeric(x["new_pods" ])/2))}
return(round(as.numeric(x["new_pods" ])/3))
}
#' Grow factor function
#'
#' This function calculates the maturity status of the pods (grow factor)
#' for each simulation day using the Zuidema et al 2003 model
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param Temp number given by daily_weather function
#' @return a vector with the new maturity status for the pods in the crop which can be added to the pod_data data frame
#' @export
growth_factor <- function (Temp,pod_data){
temperature_growth_factor <- c(19.99 , 20.98 , 21.97 , 22.96 , 23.99 , 24.48 , 24.95 , 25.97 , 26.96 , 27.97)
index_growth_factor <- c(0.0049 , 0.0053 , .0056 , .006 , .0064 , .0066 , .0067 , .0069 , .0071 , .0073)
mod <- lm(index_growth_factor ~ temperature_growth_factor)
gi <- predict(mod,newdata = data.frame(temperature_growth_factor = Temp))
growing <- function(x){x+gi}
a <- map(pod_data["maturity"],growing)
return(a)
}
#' Probability of Infection function
#'
#' This function estimates the probability of a new infection to develop
#' for each simulation day.
#'
#' @param cultivar character specified by the user in the crop characteristic dataframe
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param modeling_days increasing number
#' @param Temp number given by daily_weather function
#' @param Hum number given by daily_weather function
#' @param control_days list with a vector and a number given by Estimating_control_action_days function
#' @return the probability of a new infection the simulated day
#' @export
prob_infect<- function(cultivar, Temp, WH,plot_data,modeling_days,control_days){
incidence<-tail(plot_data$Sporulated, n=1)/(tail(plot_data$Healthy, n=1)+tail(plot_data$Sporulated, n=1)+tail(plot_data$Latent, n=1))
cd<-as.numeric(control_days[2])
if (cd <15){return(0)}### asuminedo que fumigacion protege mazorcas 15 dias
if (WH >8){
if (incidence>0.1){return(0.05)} #incidence of 60% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time) when no cultural sanitation is made
else{return(0.02)}#incidence of 40% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
}
if (incidence>0.1){return(0.01)}#incidence of 20% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
else{return(0)}#incidence of 5% (defiend as number of esporulated pods per total pods in the crop at the end of simulating time)
}
#' Risk factor function
#'
#' This function estimates the risk that a new epidemic develops for each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param healthy number given by mature_sporulated_healty_latent_pod function
#' @param latent number given by mature_sporulated_healty_latent_pod function
#' @param R factor given by risk_factor function for the day before
#' @param p probability of a new infection given by prob_infect function
#' @return the risk of epidemy development
#' @export
risk_factor<-function(p,R,pod_data,healthy,sporulated,latent){
#young_pod<-sum(pod_data$Healthy[pod_data$maturity<0.33])/(healthy+sporulated+latent)
#daily_risk<-p*young_pod
#if(p==0.001){R<-R+daily_risk}
if(p==0.05){R<-R+10}
else if (p==0.02) {R<-R+8}
else if (p==0.01) {R<-R+5}
else {R<-R-5}
if(R<0){R<-0}
return(R)
}
#' new infected function
#'
#' This function estimates the number of pod infected in each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @return a number showing how many pods got infected in each day
#' @export
new_infected <- function (p, pod_data) {
new_inf <- sum(pod_data$Healthy[pod_data["maturity"]<0.6]*p)
return(new_inf)
}
#' new healthy function
#'
#' This function estimates the number of pod that remained healthy in each simulation day.
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param p probability of a new infection given by prob_infect function
#' @return a number showing how many pods remained infected in each day
#' @export
new_healthy <- function (p , pod_data) {
a<-pod_data$Healthy[pod_data["maturity"]<0.6]*(1-p)
b <- pod_data$Healthy[pod_data["maturity"]>=0.6]
c<-c(a,b)
return(c)
}
#' Sporulation factor function
#'
#' This function calculates the sporulation status of the pods (grow factor)
#' for each simulation day
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param Temp number given by daily_weather function
#' @return a vector with the new sporulation status for the pods in the crop which can be added to the pod_data data frame
#' @export
#'
sporulation_grade <- function (Temp,pod_data){
si <- 1/65
sporulating <- function(x){x+si}
a <- map(pod_data["sporulate_level"],sporulating)
return(a)
}
#' new mature and sporulated function
#'
#' This function estimates the number of pod that grew to maturity and sporulation
#' in each simulation day. Then uses this information to actualize the pod_data
#' data file (remove last rows in which the pods completed growth and sporulation)
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param modeling_days increasing number
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @return a list with numbers representing the new mature and sporulating pods, the increased modeling day and the updated pod_data
#' @export
#'
mature_sporulated_healty_latent_pod <- function (pod_data , mature , sporulated , modeling_days) {
if (length(subset(pod_data,pod_data["maturity"]>=1)["Healthy"][,1])>0){
mature<-mature+sum(subset(pod_data,pod_data["maturity"]>=1)["Healthy"])
pod_data[pod_data$maturity>=1,]["Healthy"]<-0
}
if (length(subset(pod_data,pod_data["sporulate_level"]>=1)["Disease"][,1]) > 0){
sporulated<-sporulated+sum(subset(pod_data,pod_data["sporulate_level"]>=1)["Disease"])
pod_data[pod_data$sporulate_level>=1,]["Disease"]<-0}
rows_to_drop<-pod_data$Healthy==0 & pod_data$Disease==0
pod_data<-pod_data[!rows_to_drop,]
modeling_days<-modeling_days+1
return( list(mature,sporulated, modeling_days , pod_data))
}
#' pods harvested and removed function
#'
#' This function estimates the number of healthy pods that were harvested
#' or diseased pods removed from the crop in each simulation day.Then uses this
#' information to actualize the pod_data and the pods with each disease status
#' remaining in the status afetr the remotions
#'
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param modeling_days increasing number
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param control_days list with a vector and a number given by Estimating_control_action_days function
#' @return a list with numbers representing the new mature,sporulated,harvest,healthy,latent,disease_remove, sporulated_remove and the updated pod_data
#' @export
#'
harvest_remotion <- function(pod_data , mature , sporulated, modeling_days , x , control_days){
if (modeling_days %% as.numeric(x["croping_frequency"]) == 0) {
harvest<-mature
mature<-0}else{harvest<-0}
if (control_days[1]==0 & length(pod_data["Disease"])!=0) {
disease_remove<-sum(pod_data["Disease"])*as.numeric(x["removing_latent_eficacy"])
sporulated_remove<-sporulated*as.numeric(x["removing_sporulated_eficacy"])
sporulated<-sporulated*(1-as.numeric(x["removing_sporulated_eficacy"]))
pod_data["Disease"] <- pod_data["Disease"]*(1-as.numeric(x["removing_latent_eficacy"]))
}else{
disease_remove<-0
sporulated_remove<-0
}
if (length(pod_data["Healthy"])!=0){
healthy<- sum(pod_data["Healthy"])
}else(healthy<-0)
if (length(pod_data["Disease"])!=0){
latent<- sum(pod_data["Disease"])
}else(latent<-0)
return(list(mature,sporulated,harvest,healthy,latent,disease_remove, sporulated_remove,pod_data))
}
#' daily weather function
#'
#' This function estimates the temperature and presipotation for each day.
#'
#' @param md data frame with the wheather
#' @param modeling_days increasing number
#' @return a list with numbers representing the day temperature, presipitation, RH and hours of wetennes
#' @export
#'
daily_weather<-function(md,modeling_days){
Temp<-as.numeric(md$Temperature[modeling_days])
Hum<-as.numeric(md$Humidity[modeling_days])
Rain<-as.numeric(md$Rain[modeling_days])
WH_CART<-as.numeric(md$CART_WH[modeling_days])
WH_TH<-as.numeric(md$TH_WH[modeling_days])
return(list(Temp,Hum,WH_CART,WH_TH,Rain))
}
#' Control action function
#'
#' This function estimates the day when control actions (remotios and funmigations)
#' were done. base on the information provided by the user
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param modeling_days increasing number
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return a list especifying the days in which remotions were done and the last fumigation day
#' @export
Estimating_control_action_days<-function(x,y,modeling_days){
if (length(na.omit(y$Fumigation_dates))==0 & length(na.omit(y$Remotion_dates))==0){
return(list(100,100))
}
if (length(na.omit(y$Fumigation_dates))!=0 & length(na.omit(y$Remotion_dates))!=0){
Days_fumigation<-as.numeric(as.Date(as.character(y$Fumigation_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
Days_remotion<-as.numeric(as.Date(as.character(y$Remotion_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_fumigation[1]>modeling_days & Days_remotion[1]>modeling_days){
return(list(100,100))
}
if (Days_fumigation[1]<=modeling_days){
last_fumigation<- max(na.omit(Days_fumigation[Days_fumigation<=modeling_days]))
Days_from_last_fumigation<-modeling_days-last_fumigation
}else{Days_from_last_fumigation<-100}
if (Days_remotion[1]<=modeling_days){
last_remotion<- max(na.omit(Days_remotion[Days_remotion<=modeling_days]))
Days_from_last_remotion<-modeling_days-last_remotion
}else{Days_from_last_remotion<-100}
return(list(as.numeric(Days_from_last_remotion),as.numeric(Days_from_last_fumigation)))
}
if (length(na.omit(y$Fumigation_dates))!=0){
Days_fumigation<-as.numeric(as.Date(as.character(y$Fumigation_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_fumigation[1]>modeling_days){
return(list(100,100))
}
if (Days_fumigation[1]<=modeling_days){
last_fumigation<- max(na.omit(Days_fumigation[Days_fumigation<=modeling_days]))
Days_from_last_fumigation<-modeling_days-last_fumigation
return(list(100,as.numeric(Days_from_last_fumigation)))
}
}
if (length(na.omit(y$Remotion_dates))!=0){
Days_remotion<-as.numeric(as.Date(as.character(y$Remotion_dates), format="%Y-%m-%d")-as.Date(as.character(x$Crop_starting_date), format="%Y-%m-%d"))
if (Days_remotion[1]>modeling_days){
return(list(100,100))
}
if (Days_remotion[1]<=modeling_days){
last_remotion<- max(na.omit(Days_remotion[Days_remotion<=modeling_days]))
Days_from_last_remotion<-modeling_days-last_remotion
return(list(as.numeric(Days_from_last_remotion),100))
}}
}
#' Running the simulation function
#'
#' This function runs the simulation for a specified number of days
#' using the modelo_pod_infection
#' @param weather_data data frame with the wheather
#' @param pod_data data frame("Healthy"=,"maturity"=,"Disease"=,"sporulate_level"= )
#' @param plot_data data.frame("Days"=,"Healthy"=,"Mature"=,"Latent"=,"Sporulated"=, "Harvest"=,"Latent_removed"=, "Sporulated_removed"=, "Risk_factor"=)
#' @param x data frame with crop characteristics data.frame("Crop_starting_date"=""%Y-%d-%m"","new_pods"=, "cultivar"=, "croping_frequency" =, "removing_sporulated_eficacy"=, "removing_latent_eficacy"=)
#' @param mature number given by mature_sporulated_healty_latent_pod function
#' @param sporulated number given by mature_sporulated_healty_latent_pod function
#' @param modeling_days increasing number
#' @param R factor given by risk_factor function
#' @param y data frame with control actions data.frame("Fumigation_dates"=,"Remotion_dates"=)
#' @return the plot_data with the results of the simulation
#' @export
#'
running_simulation<-function(weather_data,modeling_days,pod_data,plot_data,x,mature, sporulated, R,y,model_WH ){
for (i in 1:length(weather_data$Temperature)) {
Temp<-as.numeric(daily_weather(weather_data,modeling_days)[1])
Hum<-as.numeric(daily_weather(weather_data,modeling_days)[2])
if (model_WH=="CART"){WH<-as.numeric(daily_weather(weather_data,modeling_days)[3])}
else if (model_WH=="TH"){ WH<-as.numeric(daily_weather(weather_data,modeling_days)[4])}
##############
##########
else if (model_WH=="WH prototipo"){ WH<-
as.numeric(daily_weather(weather_data,modeling_days)[4])}
##########
############
Rain<-as.numeric(daily_weather(weather_data,modeling_days)[5])
a<-modelo_pod_infection(pod_data,plot_data,Temp,WH,x,mature, sporulated , modeling_days, R,y)
mature<-round(as.numeric(a[1]))
sporulated<-round(as.numeric(a[2]))
healthy<-round(as.numeric(a[3]))
latent<-round(as.numeric(a[4]))
modeling_days<-as.numeric(a[5])
harvest<-round(as.numeric(a[6]))
latent_removed<- round(as.numeric(a[7]))
sporulated_removed<-round(as.numeric(a[8]))
p<-round(as.numeric(a[9]),3)
R<-round(as.numeric(a[10]),3)
pod_data<-as.data.frame(a[11])
plot_data <- plot_data %>% add_row(Days = modeling_days, Healthy=healthy,Mature = mature,
Latent = latent,Sporulated=sporulated,Harvest=harvest ,
Latent_removed= latent_removed, Sporulated_removed=sporulated_removed,
Risk_factor= R)
}
return(plot_data)
}
#' line plot function
#'
#' This function generates a line plot showing represent the simulation results
#'
#' @param x a numeric vector especifing the x factor
#' @param y a numeric vector especifing the y factor
#' @param name_x a character specifing the x axis name
#' @param name_y a character specifing the y axis name
#' @param group a character vector specifying the grouping factor
#' @param group_name a character specifying the group factor name
#' @param levels_order a character vector specifying the grouping factor order for the legend
#' @param y_max a number especifying the maximum of the y axis
#' @param y_min a number especifying the manimum of the y axis
#' @param title_legend a character specifying the legend title
#' @return a line plot
#' @export
#' @examples
#' y<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
#' x<-rep(plot_data$Days,3)
#' n<-length(plot_data$Days)
#' group<-c(rep("Healthy",n),rep("Latent",n),rep("Sporulated",n))
#' plots_model_results(x,"Days",y,"Pods proportion",group,"Disease status",c("Healthy","Latent","Sporulated"),0,1, "Disease status")
plots_model_results<-function(x,name_x,y,name_y,group,group_name,levels_order,y_min,y_max,title_legend){
plot_data<-data.frame(name_x=x,name_y=y,group_name=group)
p<-ggplot(plot_data,aes(as.numeric(x),as.numeric(y),fill = factor(group, levels=levels_order)))
p<- p + theme_classic()
p<- p + geom_line(aes(color=group),size=1.1)+ scale_color_manual(values=c("#69b3a2", "#56B4E9", "#E69F00"))
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + xlab(paste("\n",name_x,sep="")) + ylab(paste(name_y,"\n",sep=""))+ylim(y_min, y_max)
p<- p +theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(0.5, 1.1, 0.5, 1.3, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
p <- p + labs(color=title_legend)
p<- p+ theme(legend.text = element_text(size=12, angle=0),
legend.title = element_text(face="bold",size=12, angle=0))
return(p)
}
#' bar plot function
#'
#' This function generates a bar plot showing represent the simulation results
#'
#' @param x a numeric vector especifing the x factor
#' @param y a numeric vector especifing the y factor
#' @param name_x a character specifing the x axis name
#' @param name_y a character specifing the y axis name
#' @param group a character vector specifying the grouping factor
#' @param group_name a character specifying the group factor name
#' @param levels_order a character vector specifying the grouping factor order for the legend
#' @param y_max a number especifying the maximum of the y axis
#' @param y_min a number especifying the manimum of the y axis
#' @param title_legend a character specifying the legend title
#' @return a bar plot
#' @export
#' @examples
#' y<-c(sum(plot_data$Harvest),sum(plot_data$Latent_removed),sum(plot_data$Sporulated_removed))
#' x<-c("Healthy","Latent","Sporulated")
#' group<-x
#' p4<-barplot_model_results(x,"Disease status",y,"Total pods harvest",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
barplot_model_results<-function(x,name_x,y,name_y,group,group_name,levels_order,y_min,y_max,title_legend){
plot_data<-data.frame(name_x=x,name_y=y,group_name=group)
p<-ggplot(plot_data,aes(factor(x),as.numeric(y),fill = factor(group, levels=levels_order)))
p<- p + theme_classic()
p<- p + geom_bar(stat="identity")+ scale_fill_manual(values=c("#69b3a2", "#56B4E9", "#E69F00"))
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + xlab(paste("\n",name_x,sep="")) + ylab(paste(name_y,"\n",sep=""))+ylim(y_min, y_max)
p<- p +theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(0.5, 1.1, 0.5, 1.3, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
#p <- p + labs(fill=title_legend)
p<- p+ theme(legend.position="none")
return(p)
}
#' Plot the simulation function
#'
#' This function plots the results of the simulation
#' @param plot_data final .csv given by running_simulation function
#' @return A plot with for panels showing diseases progres, the amount of pods harves and remove total and during the simulation
#' @export
#'
ploting_the_disease<-function(plot_data){
y<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
x<-rep(plot_data$Days,3)
n<-length(plot_data$Days)
group<-c(rep("Healthy",n),rep("Latent",n),rep("Sporulated",n))
p1<-plots_model_results(x,"Days",y,"Pods",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease status")
y1<-c(plot_data$Healthy,plot_data$Latent,plot_data$Sporulated)
tot<-plot_data$Healthy+plot_data$Latent+plot_data$Sporulated
y<-y1/rep(tot,3)
p2<-plots_model_results(x,"Days",y,"Pods proportion",group,"Disease status",c("Healthy","Latent","Sporulated"),0,1, "Disease status")
y<-c(plot_data$Harvest,plot_data$Latent_removed,plot_data$Sporulated_removed)
p3<-plots_model_results(x,"Days",y,"Pods harvest",group,"Pods harvest",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
y<-c(sum(plot_data$Harvest),sum(plot_data$Latent_removed),sum(plot_data$Sporulated_removed))
x<-c("Healthy","Latent","Sporulated")
group<-x
p4<-barplot_model_results(x,"Disease status",y,"Total pods harvest",group,"Disease status",c("Healthy","Latent","Sporulated"),0,max(y)+0.01*max(y), "Disease statust")
require(gridExtra)
plot<-grid.arrange(p1,p2,p3,p4, ncol=2,nrow=2)
return(plot)
}
#' Plot risk function
#'
#' This function plots the evolution of the risk factoer and incidence of sporulated
#' pods during the simulation and gives control resomentations
#'
#' @param plot_data final .csv given by running_simulation function
#' @return A plot with for panels showing risk progres and control recomendations
#' @export
#'
ploting_the_Risk<-function(plot_data){
df<-data.frame("Inoculum_pressure"=plot_data$Sporulated/(plot_data$Healthy+plot_data$Latent+plot_data$Sporulated),
"New_infection"=plot_data$Risk_factor, "Days"= plot_data$Days)
scaleFactor <-1
if(max(df$New_infection)!=0 & max(df$Inoculum_pressure)!=0){
scaleFactor <- max(df$Inoculum_pressure) / max(df$New_infection)}
p <- ggplot(df, aes(x=Days)) +theme_classic()
p<- p + geom_line(aes(y=Inoculum_pressure), col="#E69F00",size=1.1) +
geom_line(aes(y=New_infection * scaleFactor), col="#56B4E9",size=1.1) +
scale_y_continuous(name="Inoculum Pressure Risk", sec.axis=sec_axis(~./scaleFactor, name="New Infection Risk"))
p<- p+theme(axis.title.y.left=element_text(color="#E69F00"),axis.title.y.right=element_text(color="#56B4E9"))
if(tail(df,1)$Inoculum_pressure>0.095){
p <- p +geom_line(aes(y = 0.1), color = "#E69F00", linetype = "dotted",size=1.5)
}
p <- p +geom_line(aes(y = 30*scaleFactor), color = "#56B4E9", linetype = "dotted",size=1.5)
p<- p + theme(plot.title=element_text(size=18,face="bold",hjust = 0.5))
p<- p + theme(axis.title=element_text(size=14,face="bold"))
p<- p + theme(axis.text.x = element_text(face="bold",size=12, angle=0),
axis.text.y = element_text(face="bold",size=12, angle=0))
p<- p + theme(panel.background = element_rect(fill = "white"),plot.margin = margin(1.6, 1.5,0.8, 1.5, "cm"),
plot.background = element_rect(fill = "white",colour = "black",size = 0.3))
p<- p+ theme(legend.text = element_text(size=12, angle=0),
legend.title = element_text(face="bold",size=12, angle=0))
if(tail(df,1)$New_infection>30 & tail(df,1)$Inoculum_pressure<0.1){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$New_infection*scaleFactor,label="Recomiendaciones:\n1. Fumigar", hjust = "inward", vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure>0.1 & tail(df,1)$New_infection<30){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$Inoculum_pressure,label="Recomiendaciones:\n1. Hacer remocion", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure>0.1 & tail(df,1)$New_infection>30){
p<- p + annotate(geom="text",x=0, y=tail(df,1)$Inoculum_pressure,label="Recomiendaciones:\n1. Hacer remocion\n2. Fumigar", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
if(tail(df,1)$Inoculum_pressure<0.1 & tail(df,1)$New_infection<30){
p<- p + annotate(geom="text",x=0, y=30*scaleFactor,label="Recomiendaciones:\n1. No accion de control necesaria", hjust = "inward" , vjust = "inward" , color="Black",size=8)
}
p
return(p)
}
#' Getting the wheather from the IDEAM
#'
#' This function fetches data from the IDAEM API
#'
#' @param location the zipcode of the city for wich the data is required"=,"Remotion_dates"=)
#' @return The data frame with the weather
#' @export
#'
weather_data <- function(location){
####Inspect/Network/Headers
###location = "68001000", para bucaramanga
get_weather<-GET(paste("http://mipronostico.ideam.gov.co/IdeamWebApp2/Ideam/getData/reverse/reverse.php?cod=",location, sep = ""))
get_weather_text<-content(get_weather, "text")
get_weather_json<-fromJSON(get_weather_text,flatten = TRUE)
df <- as.data.frame(get_weather_json)
get_weather_df<-data.frame(Fecha=df$pronostico.fecha, Hora=df$pronostico.hora, Temperatura=df$pronostico.temperature,Humedad=df$pronostico.humidity,
DewPoint=df$pronostico.dewPoint,Nubosidad=df$pronostico.totalCloudCover,Presipitacion=df$pronostico.precipitation1h,wind_speed=df$pronostico.windSpeedMS)
return(get_weather_df)
}
#' Daily wet hours for the simulation dates
#'
#' This function uses the results from the LWD_daily_estimation function to generate a data frame indicating the number of hours
#' that the cocoa pod has a water film on the surface (WH) and the longest period of dry hours for the day (maxDLH)
#'
#' @param df a Data frame with the wheather parameter for each hour of each day. This df is provided by the fuction whather data
#' @return a data frame with the date, the CART_WH, TH_WH, maxDH, Temperature, Humidity and Rain
#' @export
#'
LWD_estimation <- function(df1) {
days<-names(table(df1$Fecha))
days <- format(days, format="%Y-%m-%d")
LWD_CART <- rep(0, length(days))
LWD_TH <- rep(0, length(days))
max_D <- rep(0, length(days))
Temp <- rep(0, length(days))
Hum <- rep(0, length(days))
Rain <- rep(0, length(days))
for (i in 1:length(days)) {
df <- subset(df1, df1$Fecha == days[i])
val <- LWD_daily_estimation(df)
LWD_CART[i] <- as.numeric(val[1])
LWD_TH[i]<- as.numeric(val[2])
Temp[i]<- signif(mean(as.numeric(df$Temperatura)),3)
Hum[i]<-signif(mean(as.numeric(df$Humedad)),3)
Rain[i]<-sum(as.numeric(df$Presipitacion))
}
df2 <- data.frame(Fecha = days ,
CART_WH = LWD_CART,
TH_WH= LWD_TH,
Temperature=Temp,
Humidity=Hum,
Rain=Rain)
return(df2)
}
#' Daily wet hours
#'
#' This function calculates the number of hours that the cocoa pod has a water film on the surface on a day and the longest period of dry hours for the day using the results from the
#' LW_satate function
#'
#' @param df a Data frame with the wheather parameter for each hour of each day. This df is provided by the fuction whather data
#' @return a list of three elements with the wet hours and and the longest period of dry hours for the day acording to the CART model and the wet hours acording to TH model
#' @export
#'
LWD_daily_estimation <- function(df) {
DPT <- as.numeric(df$DewPoint)
RH <- as.numeric(df$Humedad)
Temp <- as.numeric(df$Temperatura)
wind <- as.numeric(df$wind_speed)
wet_hr_CART <- rep(0, length(DPT))
wet_hr_TH <- rep(0, length(DPT))
lws<- LW_satate(DPT[1], RH[1], Temp[1], wind[1],RH[1])
wet_hr_CART[1] <- as.numeric(lws[1])
wet_hr_TH[1]<-as.numeric(lws[2])
for (i in 2:length(DPT)) {
lws<- LW_satate(DPT[i], RH[i], Temp[i], wind[i],RH[i-1])
wet_hr_CART[i] <- as.numeric(lws[1])
wet_hr_TH[i]<-as.numeric(lws[2])
}
#max(df$Temperature) & min(df$Temperature)
return(list(sum(wet_hr_CART),sum(wet_hr_TH)))
}
#' Wetness state of the cocoa pod
#'
#' This function uses the CART and TH models (Kim et al.,2002; Bregaglio; 2011) to
#' estimate whether the cocoa pod has a water film on the surface depending on the whather parameters
#'
#' @param DPT dewpoint
#' @param RH relative humidity
#' @param Temp Temperature
#' @param Wind Wind speed_km_h
#' @return a boolean list of two elements showing the results of the CART and TH models with 0:dry, 1:wet
#' @export
#'
LW_satate <- function(DPT, RH, Temp, wind, RH_ant) {
DPD <- Temp - DPT
In1 <-1.6064 * sqrt(Temp) + 0.0036 * Temp ^ 2 + 0.1531 * RH - 0.4599 * wind *DPD - 0.0035 * Temp * RH
In2 <-0.7921 * sqrt(Temp) + 0.0046 * RH - 2.3889 * wind - 0.039 * Temp * wind +1.061 * wind * DPD
wet_CART<-0
wet_TH<-0
if (RH <= 70){wet_TH <- 0} else if(RH >= 87){wet_TH<-1} else {
if (RH-RH_ant > 2.5){wet_TH<-1}else{wet_TH<-0}
}
if (DPD >= 3.7) {
return(list(wet_CART,wet_TH))
}
if (wind < 2.5) {
return(list(as.numeric(In1 > 14.4674),wet_TH)) # when false return 0 (dry variable)
}
if (RH <= 87) {
return(list(as.numeric(In2 > 37),wet_TH)) # when false return 0 (dry variable)
}
return(list(wet_CART,wet_TH))
}
#'
#' Getting the wheather from the NC
#'
#' This function fetches data from the IDAEM API
#'
#' @param q the query for the weather api in qls; q = 'SELECT mean("temp") as temp,mean("hum") as hum, difference(mean("rain")) as rain, mean("windv") as wind FROM /^DL00011$/ WHERE time >= now() - 24d GROUP BY time(1h) fill(none)'
#' @return The data frame with the weather
#' @export
#'
weather_data_NC <- function(q, user, password){
encoded <- URLencode(q)
get_weather<-GET(paste("https://login.daimob.co/api/datasources/proxy/4/query?db=agro&q=", encoded, sep = ""),
authenticate(user, password))
get_weather_text<-content(get_weather, "text")
get_weather_json<-fromJSON(get_weather_text,flatten = TRUE)
colum.names<-unlist(get_weather_json$results$series[[1]]$columns)
df1<-as.data.frame(get_weather_json$results$series[[1]]$values)
colnames(df1)<-colum.names
df<-data.frame(Fecha=format(as.Date(df1$time,format="%Y-%m-%d")), Hora=format(parse_datetime(df1$time), format="%H:%M:%S", tz="America/Bogota"), Temperatura=as.numeric(df1$temp),Humedad=as.numeric(df1$hum),
DewPoint=as.numeric(df1$temp)-(100-as.numeric(df1$hum))/5,Presipitacion=as.numeric(df1$rain),wind_speed=as.numeric(df1$wind))
return(na.omit(df))
}
#' Getting the wheather from the prototype
#'
#' This function gets the weather from prototype
#'
#' @param df a data frame of the prototype measurements
#' @return The data frame with the weather
#' @export
#'
weather_data_prototipo <- function(df){
df1<-data.frame(Fecha=format(as.Date(df$Fecha,format="%m/%d/%Y")), Hora=format(parse_datetime(paste(df$Fecha, df$Hora)," %m/%d/%Y %I:%M:%S %p"),format="%I"),Temperatura=as.numeric(df$Temperatura),Humedad=as.numeric(df$Humedad.relativa),
DewPoint=as.numeric(df$Temperatura)-(100-as.numeric(df$Humedad.relativa))/5,Presipitacion=0,wind_speed=0)
df2<-subset( df1, df1$Hora==names(table(df1$Hora))[1]&df1$Fecha==names(table(df1$Fecha))[1])
df3<-data.frame(Fecha=df2$Fecha[1], Hora=df2$Hora[1],Temperatura=mean(df2$Temperatura),Humedad=mean(df2$Humedad),
DewPoint=mean(df2$DewPoint),Presipitacion=sum(df2$Presipitacion),wind_speed=mean(df2$wind_speed))
for(j in 1:length(names(table(df1$Hora)))){
for(i in 2:length(names(table(df1$Hora)))){
df2<-subset( df1, df1$Hora==names(table(df1$Hora))[i]&df1$Fecha==names(table(df1$Fecha))[j])
df3 <- df3 %>% add_row(Fecha=df2$Fecha[1], Hora=df2$Hora[1],Temperatura=mean(df2$Temperatura),Humedad=mean(df2$Humedad),
DewPoint=mean(df2$DewPoint),Presipitacion=sum(df2$Presipitacion),wind_speed=mean(df2$wind_speed))
}}
return(na.omit(df3))
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.