In my previous post, I showed how particle movement can be simulated in the marine environment. Now that we have tracked the particle movement, how can we use this information to translate this to something more meaningful? This can be done by determining the particle density on the map.
The most simple way would be to lay a grid over the particles and simply count the number of particles in each grid cells. However, I would like to go for something a little bit complicated that gives some smoother results. I will use kernel density estimates. This means that I will consider the particles as a density function, were the position of the particle gives the highest density for that particle. You can get a density map of all particles by simply adding up the probabilities of the density functions of each individual particle. Fortunately, the KernSmooth package will take care of all this.
Remember that in the original simulation, I released 5,000 particles within a specific time-frame. With the kernel density estimated, I can estimate the number of particles per square kilometer at any time and location. This will become even more interesting if these particles represent something in the real world, let's say oil. If we know how much cubic meter of oil is represented by each particle, the estimates can be transformed in to the thickness of the oil slick. Keep in mind that this is just an initial test simulation that does not include all processes in order to properly simulate oil spills (e.g., oil evaporation, degradation and dispersion is not simulated).
Here's a video created with ffmpeg of the particle density movement over time.
And here's the fully commented script used to create the frames for the video.
| 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)) | |
| } |
