#' Loopit 2D/3D
#'
#' Wrapper function to increase performance by looping the trackit-functions in small time intervalls.
#'
#' Function to run the functions loopit_trackit_2D/loopit_trackit_3D to follow particles through different consecutive ROMS-sclices. Looping can also increase performance when using very large number of particles by looping through shorter time steps.
#' Loops are set to run in half day intervals. If no runtime is defined, the function will loop depending on the depth of the deepest cell and the sinking speed to allow each particle to possibly sink to the seafloor (2*max(h)/speed)
#'
#' @param pts_seeded matrix of particles with 3 colums (lon, lat, depth)
#' @param romsobject list of matrices containing ROMS-model cell-values (lon_u, lat_u, h, i_u, i_v, i_w)
#' @param roms_slices number of time-frames to use in the particle-tracking
#' @param start_slice determines which roms_slice the particle-tracking starts with
#' @param domain either "2D" or "3D"
#' @param trajectories TRUE/FALSE statement to define whether to store particle trajectories (model runs much faster without storing trajectories). Default is FALSE.
#' @param speed (w_sink) sinking rate m/days
#' @param runtime (time) total number fo days to run the model
#' @param looping_time default at 0.25 which is equal to the 6h intervall of the ROMS-model
#' @param sedimentation TRUE/FALSE statement whether particles should settle on the seafloor depending on current speed and particle density (McCave & Swift 1976). Default is FALSE,
#' @param particle_radius radius of the particles; this influences the sedimentation rate with smaller values meaning less sedimentation
#' @param uphill_restricted define whether particles are restricted from moving uphill, defined as from how many meters difference particles cannot cross between cells
#' @param sed_at_max_speed particles will settle at all times only depending on the highest current speed given in any of the ROMS slices. Currently this only work when 4 sclices are available. Default is FALSE
#'
#' @return list(pts=pts, pend=pend, stopindex=obj$stopindex, ptrack=obj$ptrack, lon_list=lon_list, idx_list=idx_list, idx_list_2D=idx_list_2D, id_list=id_list)
#' @export
#' @examples
#' data(surface_chl)
#' data(toyROMS)
#'
#' ########## 3D-tracking:
#' pts_seeded <- create_points_pattern(surface_chl, multi=100)
#' run <- loopit_2D3D(pts_seeded = pts_seeded, romsobject = toyROMS, roms_slices = 4, speed = 100, runtime = 50, domain = "3D", trajectories = TRUE)
#'
#' ## testing the output
#' library(rasterVis)
#' library(rgdal)
#' library(rgl)
#' ra <- raster(nrow=50,ncol=50,ext=extent(surface_chl))
#' r_roms <- rasterize(x = cbind(as.vector(toyROMS$lon_u), as.vector(toyROMS$lat_u)), y= ra, field = as.vector(-toyROMS$h))
#' pr <- projectRaster(r_roms, crs = "+proj=laea +lon_0=137 +lat_0=-66") #get the right projection (through the centre)
#'
#' plot3D(pr, adjust = FALSE, zfac = 50) # plot bathymetry with 50x exaggerated depth
#' pointsxy <- project(as.matrix(run$pend[,1:2]), projection(pr)) #projection on Tracking-points
#' points3d(pointsxy[,1],pointsxy[,2],run$pend[,3]*50)#,xlim=xlim,ylim=ylim)
#'
#'
#' ########## 2D-tracking:
#' pts_seeded <- create_points_pattern(surface_chl, multi=100)
#' run <- loopit_2D3D(pts_seeded = pts_seeded, roms_slices = 4, romsobject = toyROMS, speed = 100, runtime = 50, sedimentation = TRUE)
#'
#' plot(pts_seeded)
#' points(run$pend, col="red", cex=0.6)
#' points(run$pts , col="blue", cex=0.6)
#'
#'
#' ########## 2D-tracking with storing trajectories:
#' pts_seeded <- create_points_pattern(surface_chl, multi=100)
#' run <- loopit_2D3D(pts_seeded = pts_seeded, roms_slices = 4, particle_radius = 0.00001, romsobject = toyROMS, speed = 100, runtime = 50, sedimentation = TRUE, trajectories = TRUE)
#'
#' plot(pts_seeded)
#' points(run$pend, col="red", cex=0.6)
#' points(run$pts , col="blue", cex=0.6)
#'
#' ## looking at the horizontal flux: this should be abother function to handle the output
#' ra <- raster(nrow=50,ncol=50,ext=extent(surface_chl))
#' mat_list <- list()
#' for(islices in 1:length(run$idx_list_2D)){
#' mat_list[[islices]] <- matrix(unlist(run$idx_list_2D[[islices]]),ncol=12)
#' }
#' testmatrix <- do.call(rbind, mat_list)
#' testid <- unlist(run$id_list)
#' flux_list <- split(testmatrix,testid)
#' for(k in 1:length(flux_list)){
#' ## cells visited by a particle ("presence-only")
#' flux_list[[k]] <- unique(flux_list[[k]])
#' ## drop first and last value (input and setting cell)
#' flux_list[[k]] <- flux_list[[k]][-c(1,length(flux_list[[k]]))]
#' }
#' flux <- as.vector(unlist(flux_list))
#'
#' xlim <- c(xmin(ra),xmax(ra))
#' ylim <- c(ymin(ra),ymax(ra))
#' df <- data.frame(cbind(toyROMS$lon_u[flux],toyROMS$lat_u[flux]))
#' df$cell <- cellFromXY(ra, df)
#' ra[] <- tabulate(df$cell, ncell(ra))
#' plot(ra)
#'
#' ## looking at the current-slices
#' roms_list <- list()
#' par(mfrow=c(2,2))
#' for(i in 1:4){
#' roms_list[[i]] <- rasterize(x = cbind(as.vector(toyROMS$lon_u), as.vector(toyROMS$lat_u)), y= ra, field = as.vector(sqrt((toyROMS$i_u[,,,i]^2)+(toyROMS$i_v[,,,i])^2)))
#' plot(roms_list[[i]])
#' }
#' par(mfrow=c(1,1))
loopit_2D3D <- function(pts_seeded, romsobject, roms_slices = 1, start_slice = 1, domain = "2D", trajectories = FALSE,
speed, runtime = 10, looping_time = 0.25, sedimentation=FALSE, particle_radius=0.00016,
time_steps_in_s=1800, uphill_restricted=NULL, sed_at_max_speed=FALSE, mean_move=FALSE){
pts <- pts_seeded
loop_length <- looping_time*24*2
# h <<- romsobject$h
all_i_u <- romsobject$i_u
all_i_v <- romsobject$i_v
all_i_w <- romsobject$i_w
## setup kdtree
sknn <- with(romsobject, setup_knn(lon_u, lat_u, hh)) # (lon_roms=lon_u, lat_roms=lat_u, depth_roms=hh)
romsobject$kdtree <- sknn$kdtree
romsobject$kdxy <- sknn$kdxy
## create lists to store all particles that settled at the end of each tracking-loop
lon_list <- list()
lat_list <- list()
if(domain == "3D")
depth_list <- list()
## create lists to store the positions of each particle in each time-step
if(trajectories == TRUE){
idx_list <- list()
idx_list_2D <- list()
id_list <- list()
id_vec <- seq_len(nrow(pts_seeded))
}
# if(missing(runtime)){
# runtime <- ceiling(max(h)/speed) ## no runtime defined
# } else
runtime <- runtime ## runtime defined
curr_vector <- rep(1:roms_slices,runtime)
## allow for different starting ROMS-slices (re-arrange the vector)
sliced_vector <- curr_vector[c(start_slice:length(curr_vector),1:(start_slice-1))]
# sedimentationparams <- NULL
if(domain == "2D") sedimentationparams <- buildparams(-speed/(60*60*24), time_step_in_s = time_steps_in_s, r=particle_radius) ## sinking speed transformation into m/sec
## prepare list with parameters of the ROMS
romsparams <- list()
romsparams$h <- romsobject$h
## assign current-speed/direction to the cells
romsparams$i_u <- all_i_u
romsparams$i_v <- all_i_v
romsparams$i_w <- all_i_w
## boundaries of the ROMS-area
romsparams$roms_ext <- c(min(romsobject$lon_u), max(romsobject$lon_u), min(romsobject$lat_u), max(romsobject$lat_u))
## for a differnet approach to calculate sedimentation:
if(sed_at_max_speed==TRUE & roms_slices==4){
nrow <- dim(romsobject$i_u)[1]
ncol <- dim(romsobject$i_u)[2]
nlayer <- dim(romsobject$i_u)[3]
i_u_max <- array(NA, dim = c(nrow,ncol,nlayer))
i_v_max <- array(NA, dim = c(nrow,ncol,nlayer))
for(irow in 1:nrow){
for(icol in 1:ncol){
for(ilayer in 1:nlayer){
i_u_max[irow,icol,ilayer] <- max(abs(romsobject$i_u[irow,icol,ilayer,]))
i_v_max[irow,icol,ilayer] <- max(abs(romsobject$i_v[irow,icol,ilayer,]))
}
}
}
romsparams$i_u_max <- i_u_max
romsparams$i_v_max <- i_v_max
}
runtime <- roms_slices*runtime ## counting full days
## loop over different time-slices
for(irun in 1:runtime){
if(irun==1) message(paste0("starting # of particles: ",dim(pts)[1]))
print(paste0(irun,".loop"))
## re-assign current-speed/direction to the cells for mulitple slices
if(roms_slices>1){
romsparams$i_u <- all_i_u[,,,sliced_vector[irun]]
romsparams$i_v <- all_i_v[,,,sliced_vector[irun]]
romsparams$i_w <- all_i_w[,,,sliced_vector[irun]]
}
## save an id for each particle to follow its path
if(trajectories) id_list[[irun]] <- id_vec
## run the particle-tracking for all floating particles
if(domain == "3D"){ obj <- trackit_3D(pts=pts, romsobject=romsobject, w_sink=speed, time=looping_time, romsparams=romsparams,
loop_trackit=TRUE, time_steps_in_s=time_steps_in_s)
}else{ obj <- trackit_2D(pts=pts, romsobject=romsobject, w_sink=speed, time=looping_time, romsparams=romsparams,
loop_trackit=TRUE, time_steps_in_s=time_steps_in_s,
sedimentationparams=sedimentationparams, sedimentation=sedimentation,
particle_radius=particle_radius, uphill_restricted=uphill_restricted, sed_at_max_speed=sed_at_max_speed, mean_move=mean_move)
}
## store the particles that stopped (settled)
lon_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 1, obj$stopindex)])
lat_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 2, obj$stopindex)])
## only for 3D
if(domain == "3D") depth_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 3, obj$stopindex)])
if(trajectories){
## store the cell-indices of each pts from each time-slice
idx_list[[irun]] <- obj$indices
idx_list_2D[[irun]] <- obj$indices_2D
## reduce the id_vec to new number of pts
id_vec <- id_vec[obj$stopindex==0]
## create vector to check if list has some NULL in following if-statement
NULL_test <- as.character(idx_list[[irun]])
if(any(NULL_test == "NULL") == TRUE){
fill_up_seq <- which(NULL_test == "NULL")
for(ifill in fill_up_seq){
idx_list[[irun]][[ifill]] <- matrix(NA, nrow = nrow(idx_list[[irun]][[1]]))
idx_list_2D[[irun]][[ifill]] <- matrix(NA, nrow = nrow(idx_list_2D[[irun]][[1]]))
}
}
}
##re-assign coordinates of floating particles to re-run in "trackit"-function
if (length(unlist(lon_list))!=nrow(pts_seeded) ## check if all particles are settled
& !is.null(nrow(obj$ptrack[obj$stopindex==0,,dim(obj$ptrack)[3]])) ## if there's only one particle left it bugs around...
& length(lon_list)!=runtime){ ## if the run stops before all particles are settled, pts should not be overwritten
if(domain == "3D"){
pts <- matrix(obj$ptrack[obj$stopindex==0,,dim(obj$ptrack)[3]],ncol=3)
} else pts <- matrix(obj$ptrack[obj$stopindex==0,,dim(obj$ptrack)[3]],ncol=3)
} else break
message(paste0(dim(pts)[1]," particles floating"))
}
## store stopping-locations of particles
if(domain == "3D"){
pend <- cbind(unlist(lon_list), unlist(lat_list), unlist(depth_list))
}else pend <- cbind(unlist(lon_list), unlist(lat_list))
if(nrow(pts)==1) pend <- pend[-nrow(pend),]
message(paste0((dim(pts_seeded)[1]-dim(pend)[1])," particle(s) still floating"))
## store the output
if(trajectories==TRUE){
list(pts=pts, pend=pend, pnow=obj$pnow, stopindex=obj$stopindex, ptrack=obj$ptrack, lon_list=lon_list, idx_list=idx_list, idx_list_2D=idx_list_2D, id_list=id_list)
}else {
list(pts=pts, pend=pend, pnow=obj$pnow, stopindex=obj$stopindex, ptrack=obj$ptrack, lon_list=lon_list)
}
}
#############################
## for the settlingmodel??
## this function allows to track a greater number of particles when an error in trackit occurs due to limited RAM
## this is done by looping a number of short trackings, each using the output of the previous short tracking as their input
## "time" needs to be defined as a fraction of the total "runtime" of the tracking
# loopit <- function(vertical_movement=0.001,time=86400,time_step=3600,in_days=F,runtime=4320000){
# if(in_days==TRUE){
# vertical_movement <- vertical_movement*60*60*24
# time <- time*60*60*24
# time_step <- time_step*60*60*24
# runtime <- runtime*60*60*24
# }
# ## create lists to store all particles that settled at the end of each tracking-loop
# long_list <- list()
# lat_list <- list()
# vertical_list <- list()
# #if (missing(runtime)){
# # runtime <- ceiling(max(h)/vertical_movement) ## no runtime defined
# #} else runtime <- runtime ## runtime defined
#
# for(irun in 1:(runtime/time)){
# message(paste0(irun,".loop"))
# ## run the particle-tracking for all floating particles
# obj <- trackit(vertical_movement,time,time_step)
# ## store the particles that stopped (settled)
# long_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 1, obj$stopindex)])
# lat_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 2, obj$stopindex)])
# vertical_list[irun] <- list(obj$ptrack[cbind(seq(nrow(obj$ptrack)), 3, obj$stopindex)])
# ##re-assign coordinates of floating particles to re-run in "trackit"-function
# if (length(unlist(long_list))!=nrow(check) ## check if all particles are settled
# & !is.null(nrow(obj$ptrack[obj$stopindex==0,,dim(obj$ptrack)[3]])) ## if there's only one particle left it bugs around...
# & length(long_list)!=runtime){ ## if the run stops before all particles are settled, pts should not be overwritten
# ## Mike doesn't like this!! But it works... How to do it better?
# pts <<- matrix(obj$ptrack[obj$stopindex==0,,dim(obj$ptrack)[3]],ncol=3)
# } else break
# }
# ## store stopping-locations of particles
# pend <- cbind(unlist(long_list),unlist(lat_list),unlist(vertical_list))
# ## store remaining unsettled particles
# list(pts=pts,pend=pend,stopindex=obj$stopindex,ptrack=obj$ptrack,long_list=long_list)
# }
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.