2d.particle.tracking.kernel.r

require(raster)
require(KernSmooth)
require(maptools)

## load the simulation results. See the previous blog post
## how this file was created. This data file contains
## an object named 'part.results', which is a list
## of SpatialPointsDataFrame objects.
load("output/out.RData")

## Create a SpatialPoints object representing the release location
release.location <- SpatialPoints(t(c(3.7, 52)), proj4string = CRS("+proj=longlat +datum=WGS84"))
## everything will be plotted in UTM zone 31, so we need to transform all spatial data:
release.location <- spTransform(release.location, CRS("+proj=utm +zone=31 +datum=WGS84"))

## 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
})
## Transform to UTM:
gd_load <- lapply(gd_load, spTransform, CRSobj = CRS("+proj=utm +zone=31 +datum=WGS84"))

## x range for the grid in UTM coordinates (m)
x.rast <- c(475000, 635000)
## y range for the grid in UTM coordinates (m)
y.rast <- c(5730000, 5890000)
## dimensions for the grid
grid.size <- c(480, 480)

## grid cell surface area in sqaure kilometers
grid.area <- ((diff(x.rast)/grid.size[1])*(diff(y.rast)/grid.size[2]))/1e6

## note that 'part.results' is a list of SpatialPointsDataFrame objects.
## each element represents a time index and is named after the
## corresponding simulation date and time. Here we loop each time index:
for (i in 1:length(part.results)) {
  ## select data at time index i:
  sel.points <- part.results[[i]]
  ## only select particles that have been released:
  sel.points <- sel.points[sel.points$released,]
  ## only continue if particles have been released:
  if (length(sel.points) > 0) {
    ## project to UTM such that the coordinates are in meters:
    sel.points <- spTransform(sel.points, CRS("+proj=utm +zone=31 +datum=WGS84"))
    ## calculate the kernel density using a 2x2 km bandwidth:
    est <- bkde2D(coordinates(sel.points),
                  bandwidth = c(2000, 2000),
                  gridsize  = grid.size,
                  range.x = list(x.rast, y.rast))
    ## multiply density estimated with number of released particles and divide by
    ## grid cell surface area multiplied with the total density, to get the density
    ## in number of particles per square kilometer:
    est.raster = raster(list(x = est$x1, y = est$x2, z = est$fhat*length(sel.points)/(sum(est$fhat)*grid.area)))
  } else  {
    est.raster = raster(list(x = seq(x.rast[1], x.rast[2], length.out = grid.size[1]),
                             y = seq(y.rast[1], y.rast[2], length.out = grid.size[2]),
                             z = matrix(0, grid.size[1], grid.size[2])))
  }

  projection(est.raster) <- CRS("+proj=utm +zone=31 +datum=WGS84")

  ## open a graphics device to plot the density estimates
  png(sprintf("output\\kernout%03d.png", i), 640, 480, res = 100, type = "cairo")

  ## use print as otherwise lattice based graphics won't be drawn
  ## onto the opened device
  print(
    ## I use "spplot" here for nicer plots. using "plot" here will be faster.
    spplot(est.raster,
           ylab.right = "# / km^2", 
           par.settings = list(layout.widths = list(axis.key.padding = 0, 
                                                    ylab.right = 2)),
           col.regions  = colorRampPalette(c("white", rainbow(15, start = 1/6), "black"), bias = 1)(100),
           at           = seq(0, 20, length.out = 100),
           main         = names(part.results)[i],
           maxpixels    = prod(grid.size),
           panel        = function(...) {
             panel.gridplot(...)
             lapply(gd_load, sp.polygons, fill = "lightgrey")
             sp.points(release.location, col = "red")
             panel.text(x.rast[1] + 10000, y.rast[2] - 10000, "r-coders-anonymous.blogspot.com", adj = c(0, 0.5))
           })
  )
  ## close the graphics device
  dev.off()

  ## print ratio between the current time index and the total
  ## number of indices. This serves as a progress indicator.
  print(i/length(part.results))
}