briskaR
This vignette explains how to use the briskaR
package.
The R package briskaR
offers a generic framework for the simulation of spatially-explicit exposure-hazard model for ecotoxicological studies.
The modeling framework has four main components:
A simulation run is divided in 7 steps:
knitr::opts_chunk$set(fig.height = 6, fig.width = 6, fig.align = "center")
# 1. Load the main package library("briskaR") ## Set briskaR working intern projection to LAMBERT 93 (default) # briskaRSetInternProjection(LAMBERT_93) # 2. Load 'ggplot2' for plotting objects. Not include in the package # 3. load 'sf' for extra function on shapefiles # 4. load 'raster' for extra function on raster library("ggplot2") library("sf") library("raster") library("sp") # plotting require this package library("dplyr")
data("sfMaize65") # PLOT ------------------------------------------------------------------------ ggplot() + theme_minimal() + scale_fill_manual(values = c("grey", "orange"), name = "Maize") + geom_sf(data = sfMaize65, aes(fill = as.factor(maize)))
sfMaize65$maize_GM<-sfMaize65$maize*c(0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,1,0,1,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,1,0,0,0,1,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,0,1,0,0,0,0,0,0,0,0,1,0,0,0,0,0,1,0,0,0,0,1,0,1,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,0,0,1,1,0,0,1,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,1,1,0,0,0,0,1,0,0,1,0,0,0,0,0,1,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0)
# Random 20% are GM field sfMaize65$maize_GM <- sfMaize65$maize * rbinom(n = nrow(sfMaize65), size = 1 , prob = 0.2)
# filter table to have only GM fields sfMaize65_GM <- sfMaize65[sfMaize65$maize_GM == 1,] ## PLOT ------------------------------------------------------------------------ plt_GM <- ggplot() + theme_minimal() + scale_fill_manual(values = c("grey", "red"), name = "GM maize") + geom_sf(data = sfMaize65, aes(fill = as.factor(maize_GM))) plt_GM + geom_sf_text(data = sfMaize65_GM, aes(label = label))
Before proceding the dispersal, we can add a buffer around the landscape. A good approach would be to define sources and receptors at this step in order to embedded all area within the square frame. The eposure is going to be computed over this square frame.
squareFrame_sfMaize65 <- st_squared_geometry(list(sfMaize65), buffer = 200) # PLOT ------------------------------------------------------------------------ plt_GM + geom_sf(data = squareFrame_sfMaize65, fill = NA)
stack_dispersal <- brk_dispersal(sfMaize65_GM, size_raster = 2^8, kernel = "geometric", kernel.options = list("a" = -2.63), squared_frame = squareFrame_sfMaize65) # PLOT ------------------------------------------------------------------------ raster::plot(stack_dispersal[[1:6]])
brk_emission <- function(sf, # geometry on which we applied the pollen emission pattern keyTime, key, # name of the column which is going to be created FUN # a function of timeline? ){ if("stackTimeline" %in% colnames(sf)) { stop("Please rename column 'stackTimeline' in sf object.") } # CREATE NEW COLUMN sf[["key_temp"]] <- lapply(1:nrow(sf), FUN) # check length with timeline if(all(sapply(sf[[keyTime]], length) != sapply(sf[["key_temp"]], length))){ stop(paste0("element within list returned by FUN' has not the same length as list element of '", keyTime, "' column")) } # ADD StackTimeLine stackTimeline = sort(unique(do.call("c", sf[[keyTime]]))) sf[[key]] = lapply(1:nrow(sf), function(i){ index_matching = match(sf[[keyTime]][[i]], stackTimeline) res = rep(0,length(stackTimeline)) res[index_matching] = sf[["key_temp"]][[i]] return(res) }) # /!\ REPLACE timeline sf[[keyTime]] = lapply(1:nrow(sf), function(i) stackTimeline) warning(paste("The column variable", keyTime, "may have changed")) # remove "key_temp" sf[["key_temp"]] = NULL return(sf) }
Emission, Exposure, and then individuals developement, migration and so need a time line. This is done using the function brk_timeline
as follow:
sfMaize65_GM_Pollen <- brk_timeline(sf = sfMaize65_GM, key = "timeline", from = as.Date("01-07-2018", format = "%d-%m-%y"), to = as.Date("01-09-2018", format = "%d-%m-%y"), by = "days")
# Profile of emission: data("maize.proportion_pollen") graphics::plot(maize.proportion_pollen, type= "l") funTimePollen <- function(time){ density = runif(1, 7, 11) pollen = rgamma(1, shape = 1.6, scale = 1 / (2 * 10 ^ -7)) nbr_days = length(time) deb = sample(1:(nbr_days - length(maize.proportion_pollen)), 1) end = (deb + length(maize.proportion_pollen) - 1) pollen_emission <- rep(0, nbr_days) pollen_emission[deb:end] <- as.numeric(pollen * density * maize.proportion_pollen) return(pollen_emission) } # ADD Column EMISSION To sfMaize65_GM_Pollen <- brk_emission(sf = sfMaize65_GM_Pollen, keyTime = "timeline", # Length of the reference timeline key = "EMISSION", # Name of the new column # function over each line of sf FUN = function(i){ funTimePollen(sfMaize65_GM_Pollen$timeline[[i]]) })
stackTimeline = seq(from = min(do.call("c", sfMaize65_GM_Pollen[["timeline"]])), to = max(do.call("c", sfMaize65_GM_Pollen[["timeline"]])), by = "days") stack_exposure <- brk_exposure(stack_dispersal, sfMaize65_GM_Pollen, key = "EMISSION", # Name of the new column keyTime = "timeline", # Length of the reference timeline loss = 0.1, beta = 0.2, quiet = TRUE) # PLOT ------------------------------------------------------------------------ raster::plot(stack_exposure[[1:6]])
Plot of exposure out of GM fields
# A MULTIPOLYGON OUT OF GM FIELDS sfMaize65_outGM <- st_sf(geometry = st_difference(x = st_geometry(squareFrame_sfMaize65), y = st_union(st_geometry(sfMaize65_GM)))) # POINTS OUT OF GM FIELDS gridPOINT_squareFrame <- st_make_grid(squareFrame_sfMaize65, n = 30, what = "centers") # grid n x n !! gridPOINT_outGM <- st_sf(geometry = st_intersection(x = st_geometry(gridPOINT_squareFrame), y = st_union(st_geometry(sfMaize65_outGM)))) # PLOT ------------------------------------------------------------------------ ggplot() + theme_minimal() + geom_sf(data = sfMaize65_outGM, fill = "grey30") + geom_sf(data = gridPOINT_outGM) # point out of GM field
# Create raster df_outGM = as.data.frame(raster::extract(x = stack_exposure, y = sf::as_Spatial(gridPOINT_outGM))) #colnames(df_outGM) = paste0("Time_", gsub("-", "\\1", stackTimeline)) sf_outGM = st_as_sf(geometry = st_geometry(gridPOINT_outGM), df_outGM) # PLOT ------------------------------------------------------------------------ ggplot() + theme_minimal() + labs(title = paste("Exposure at", colnames(df_outGM)[40])) + scale_color_continuous(low = "green", high = "red", trans = "log", name = "log scaled") + geom_sf(data = sf_outGM, aes(color = df_outGM[,40]))
First of all is to define the receptor area.
We can take any set of place within the area where exposure has been computed. A good way is to consider both receptor and source at the beginning in order to embed both in squareFrame
# Margins around each fields sfMaize65_receptor = st_multibuffer(sfMaize65_GM, dist = rep(100, nrow(sfMaize65_GM))) # PLOT ------------------------------------------------------------------------ plt_GMreceptor <- ggplot() + theme_minimal() + geom_sf(data = sfMaize65_GM, fill = "red") + geom_sf(data = sfMaize65_receptor, fill = "#669900") plt_GMreceptor
# 1. Number of site for the first generation nbrSite = 100 # 2. Set eggs sfLarvae <- brk_newPoints(sf = sfMaize65_receptor, size = nbrSite) # size = number of sites # PLOT ------------------------------------------------------------------------ plt_GMreceptor + geom_sf(data = sfLarvae)
DateEmergence = sample(seq(as.Date("01-07-2018", format = "%d-%m-%y"), as.Date("01-09-2018", format = "%d-%m-%y"), by = "days"), size = nbrSite, replace = TRUE) sfLarvae = brk_timeline( sf = sfLarvae, key = "Date", from = DateEmergence, to = DateEmergence + 20, by = "days")
First step is to recover the time line of emission.
stackTimelineEMISSION = sort(unique(do.call("c", sfMaize65_GM_Pollen[["timeline"]])))
Then we make matching the emission/deposition profile with host
exposureINDIVIDUAL <- brk_exposureMatch(stackRaster_exposure = stack_exposure, sf = sfLarvae, stackTimeline = stackTimelineEMISSION, keyTime = "Date", key = "EXPOSURE")
Finally we compute a Da
damageLethal = function(x,LC50, slope){ return(1/(1+(x/LC50)^slope)) } LC50DR = 451*10^4 slopeDR = -2.63 damageINDIVIDUAL <- exposureINDIVIDUAL %>% dplyr::mutate(DAMAGE = lapply(EXPOSURE, function(expos){damageLethal(expos, LC50DR, slopeDR)}))
DFdamage = data.frame( DAMAGE = do.call("c", damageINDIVIDUAL[["DAMAGE"]]), Date = do.call("c", damageINDIVIDUAL[["Date"]]) ) %>% dplyr::group_by(Date) %>% dplyr::summarise(mean_DAMAGE = mean(DAMAGE, na.rm = TRUE), q025_DAMAGE = quantile(DAMAGE, probs = 0.025, na.rm = TRUE), q975_DAMAGE = quantile(DAMAGE, probs = 0.975, na.rm = TRUE), min_DAMAGE = min(DAMAGE, na.rm = TRUE), max_DAMAGE = max(DAMAGE, na.rm = TRUE)) minDateDAMAGE = data.frame(Date = do.call("c", damageINDIVIDUAL[["Date"]])) ggplot() + theme_minimal() + labs(x = "Time", y = "Probability Distribution of Damage") + geom_line(data = DFdamage, aes(x = Date, y = mean_DAMAGE), color = "red") + geom_ribbon(data = DFdamage, aes(x = Date, ymin = q025_DAMAGE, ymax = q975_DAMAGE), alpha = 0.5, color = NA, fill = "grey10") + geom_ribbon(data = DFdamage, aes(x = Date, ymin = min_DAMAGE, ymax = max_DAMAGE), alpha = 0.5, color = NA, fill = "grey90")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.