ptrackr: Tracking Individual Particles Through 2D- and 3D-Space

Description Usage Arguments Details Value Examples

View source: R/loopit_2D3D.R

Wrapper function to increase performance by looping the trackit-functions in small time intervalls.

loopit_2D3D(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_seeded`	matrix of particles with 3 colums (lon, lat, depth)
`romsobject`	list of matrices containing ROMS-model cell-values (lon_u, lat_u, h, i_u, i_v, i_w)
`roms_slices`	number of time-frames to use in the particle-tracking
`start_slice`	determines which roms_slice the particle-tracking starts with
`domain`	either "2D" or "3D"
`trajectories`	TRUE/FALSE statement to define whether to store particle trajectories (model runs much faster without storing trajectories). Default is FALSE.
`speed`	(w_sink) sinking rate m/days
`runtime`	(time) total number fo days to run the model
`looping_time`	default at 0.25 which is equal to the 6h intervall of the ROMS-model
`sedimentation`	TRUE/FALSE statement whether particles should settle on the seafloor depending on current speed and particle density (McCave & Swift 1976). Default is FALSE,
`particle_radius`	radius of the particles; this influences the sedimentation rate with smaller values meaning less sedimentation
`uphill_restricted`	define whether particles are restricted from moving uphill, defined as from how many meters difference particles cannot cross between cells
`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

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)

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)

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))