2d.particle.tracking.r

## load required packages:
require(sp)
require(raster)
require(ncdf4)
require(ncdf.tools)
require(rgeos)
require(maptools)

###############################################################
## STEP 1: load the water currents and wind data
###############################################################

## function that reads netcdf files and puts them in a workable stack object:
read.mover <- function(ncdfile, xvar = "u", yvar = "v") {
  u <- stack(ncdfile, varname = xvar)
  v <- stack(ncdfile, varname = yvar)
  ## Get time indices:
  z_time      <- attr(u, "z")
  if (class(z_time[[1]]) == "character") {
    z_time <- as.POSIXct(z_time[[1]], tz = "UTC")
  } else {
    z_time_info <- tolower(strsplit(tolower(names(z_time)), " since ")[[1]])
    z_time      <- z_time[[1]]*switch(z_time_info[1], seconds = 1, minutes = 60, hours = 60*60, days = 60*60*24)
    z_time      <- as.POSIXct(z_time, origin = as.POSIXct(z_time_info[2], tz = "UTC"), tz = "UTC")
    rm(z_time_info)
  }

  return(list(u = u, v = v, z = z_time))
}

## Read the water current data:
temp       <- read.mover("input/MetO-NWS-PHYS-hi-Agg-all-levels_1468581872238.nc",
                         "vozocrtx", "vomecrty")
water_u    <- temp$u
water_v    <- temp$v
water_time <- temp$z
rm(temp)

## Read the wind data:
temp       <- read.mover("input/CERSAT-GLO-BLENDED_WIND_L4-V5-OBS_FULL_TIME_SERIE_1468508078544.nc",
                         "eastward_wind_rms", "northward_wind_rms")
wind_u     <- temp$u
wind_v     <- temp$v
wind_time  <- temp$z
rm(temp)

###############################################################
## STEP 2: specify particle release details
###############################################################

## time at which the simulated release should start:
release.start         <- as.POSIXct("2016-07-01 08:00", tz = "UTC")
## time at which teh simulated release should end:
release.end           <- as.POSIXct("2016-07-06 08:00", tz = "UTC")
## number of particles to release in total:
release.n             <- 5000
## release lon-coordinate:
release.x             <- 3.7
## release lat-coordinate:
release.y             <- 52
## release depth (not used yet in simulation)
release.z             <- 0

###############################################################
## STEP 3: specify simulation details
###############################################################

## time at which the simulation should start:
simulation.start.time <- as.POSIXct("2016-07-01 00:00", tz = "UTC")
## time step of simulation in seconds:
simulation.time.step  <- 60*60/4
## integer value indicating every n-th time step that should be exported:
simulation.out.step   <- 4
## time at which the simulation should stop
simulation.end.time   <- as.POSIXct("2016-07-10 08:00", tz = "UTC")
## create a time series based on the start and end time and the time step:
simulation.time       <- seq(simulation.start.time, simulation.end.time - simulation.time.step, simulation.time.step)

## projection to be used to calculate particle displacement in meters:
part.proj             <- CRS("+proj=utm +zone=31 +datum=WGS84")

## drift factor indicates with which fraction particles are displaced by the wind speed:
drift.factor    <- 0.03
## the diffusion coefficient, used in random walk simulation
diffusion.coeff <- 1e5

## Countries for which the shorelines are downloaded and plotted
countries <- "NLD"
gd_load <- lapply(as.list(countries), function(x) {
  gadm <- ""
  lev <- 0
  if (x == "NLD") lev <- 1
  try(gadm <- raster::getData("GADM", country = x, level = lev))
  if (class(gadm) == "SpatialPolygonsDataFrame" && lev == 1 && !is.character(gadm)) {
    gadm <- unionSpatialPolygons(gadm, gadm@data$ENGTYPE_1 != "Water body")[2]
  }
  gadm
})

## for each time step determine how many particles are released:
n.particles           <- rep(NA, length(simulation.time))
n.particles[simulation.time <= release.start] <- 0
n.particles[simulation.time >= release.end]   <- release.n
n.particles[is.na(n.particles)]               <- approx(c(0, release.n), n = 2 + sum(is.na(n.particles)))$y[c(-1, -(2 + sum(is.na(n.particles))))]
n.particles           <- round(n.particles)

## create a SpatialPointsDataFrame for the particles:
particles              <- data.frame(pos.x    = rep(release.x, release.n),
                                     pos.y    = rep(release.y, release.n),
                                     pos.z    = rep(release.z, release.n))
coordinates(particles) <- ~ pos.x + pos.y
proj4string(particles) <- CRS("+proj=longlat +datum=WGS84")

## create a list in which results should be stored:
part.results          <- list()
previous.position     <- particles

## function that extrapolates water currents for NA values
extend.data <- function(r, band = 21) {
  ## replace NA values of raster data at coasts with nearest neighbours:
  fill.na <- function(x, i = round(band^2/2)) {
    if(is.na(x)[i]) {
      return(mean(x, na.rm = T))
    } else {
      return(x[i])
    }
  }
  r <- setValues(r, getValues(focal(r, w = matrix(1,band,band), fun = fill.na, 
                                    pad = TRUE, na.rm = FALSE )))
  r
}

###############################################################
## STEP 4: loop all simulation time steps
###############################################################

## start looping the simulation times
for (i in 1:length(simulation.time)) {
  print(i)
  ## Get the water currents that are closest to the current simulation time
  tdif   <- abs(as.numeric(water_time) - as.numeric(simulation.time[i]))
  tsel_c <- which(tdif == min(tdif))[[1]]
  twat_u <- extend.data(raster(water_u, tsel_c))
  twat_v <- extend.data(raster(water_v, tsel_c))
  tdif   <- abs(as.numeric(wind_time) - as.numeric(simulation.time[i]))
  tsel_w <- which(tdif == min(tdif))[[1]]
  twin_u <- extend.data(raster(wind_u, tsel_w))
  twin_v <- extend.data(raster(wind_v, tsel_w))

  ## Get the last known particle locations and transform coordinates into meters:
  part.coords <- coordinates(spTransform(previous.position, part.proj))

  ## displace particles with with water_currents
  part.coords[0:n.particles[i],] <- part.coords[0:n.particles[i],] +
    simulation.time.step*cbind(raster::extract(twat_u, previous.position, "bilinear"),
                               raster::extract(twat_v, previous.position, "bilinear"))[0:n.particles[i],]
  ## displace particles  with wind
  part.coords[0:n.particles[i],] <- part.coords[0:n.particles[i],] +
    simulation.time.step*drift.factor*cbind(raster::extract(twin_u, previous.position, "bilinear"),
                                            raster::extract(twin_v, previous.position, "bilinear"))[0:n.particles[i],]
  ## displace with random walk (diffusion):
  part.coords[0:n.particles[i],] <- part.coords[0:n.particles[i],] +
    cbind(runif(n.particles[i], -1, 1),
          runif(n.particles[i], -1, 1))*sqrt(diffusion.coeff*simulation.time.step*6/1e4)

  ## turn particle information into a SpatialPointsDataFrame object:
  part.coords <- as.data.frame(part.coords)
  coordinates(part.coords) <- ~ pos.x + pos.y
  proj4string(part.coords) <- part.proj
  ## backtransform coordinates to original projection:
  part.coords <- spTransform(part.coords, CRS(proj4string(previous.position)))

  ## store the new locations in the results list:
  previous.position     <- part.coords

  ## for each n-th simulation step store the results (where n = simulation.out.step):
  if ((i - 1)%%4 == 0) {
    ## add an indicator of the release state of the particles:
    part.coords <- SpatialPointsDataFrame(part.coords, data.frame(released = rep(F, release.n)))
    part.coords@data$released[0:n.particles[i]] <- T
    part.results[[(i - 1)/4 + 1]] <- part.coords

    ## name the list element after the time index it represents:
    names(part.results)[(i - 1)/4 + 1] <- paste(format(simulation.time[i], "%Y-%m-%d %H:%M:%S"),
                                                attr(as.POSIXlt(simulation.time[i]),"tzone"))

    png(sprintf("output\\out%03d.png", (i - 1)/4), 640, 480, res = 75, type = "cairo")
    plot.ylim <- c(51.9, 53)
    plot(part.results[[(i - 1)/4 + 1]], pch =".", xlim = c(3, 5), ylim = plot.ylim,
         main = paste(format(simulation.time[i], "%Y-%m-%d %H:%M:%S"),
                      attr(as.POSIXlt(simulation.time[i]),"tzone")))
    
    ## draw water current arrows
    wat_arrow <- as.data.frame(coordinates(raster(water_u, tsel_c)))
    wat_arrow$u <- as.vector(t(as.matrix(raster(water_u, tsel_c))))
    wat_arrow$v <- as.vector(t(as.matrix(raster(water_v, tsel_c))))
    scale = 0.05
    aspect.ratio <- 1/cos(mean(plot.ylim)*pi/180)
    arrows(wat_arrow$x, wat_arrow$y,
           wat_arrow$x + scale*aspect.ratio*wat_arrow$u,
           wat_arrow$y + scale*wat_arrow$v, length = 0.01, col = "lightpink")
    
    ## draw wind current arrows
    win_arrow <- as.data.frame(coordinates(raster(wind_u, tsel_w)))
    win_arrow$u <- as.vector(t(as.matrix(raster(wind_u, tsel_w))))
    win_arrow$v <- as.vector(t(as.matrix(raster(wind_v, tsel_w))))
    scale <- 0.05
    aspect.ratio <- 1/cos(mean(plot.ylim)*pi/180)
    arrows(win_arrow$x, win_arrow$y,
           win_arrow$x + scale*aspect.ratio*win_arrow$u,
           win_arrow$y + scale*win_arrow$v, length = 0.01, col = "lightblue")
    
    ## draw coastlines:
    lapply(gd_load, plot, add = T, col = "lightgrey")
    text(2.5, 52.9, "r-coders-anonymous.blogspot.com", adj = c(0, 0.5), cex = 1.5)
    axis(1)
    axis(2)
    dev.off()
  }
}

## save the resulting particle information:
save(part.results, file = "output/out.RData", compress = T)