Skip to content

Instantly share code, notes, and snippets.

@pepijn-devries
Last active July 21, 2016 08:33

Revisions

  1. pepijn-devries revised this gist Jul 21, 2016. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion 2d.particle.tracking.kernel.r
    Original file line number Diff line number Diff line change
    @@ -95,4 +95,4 @@ for (i in 1:length(part.results)) {
    ## print ratio between the current time index and the total
    ## number of indices. This serves as a progress indicator.
    print(i/length(part.results))
    }
    }
  2. pepijn-devries revised this gist Jul 21, 2016. 1 changed file with 3 additions and 1 deletion.
    4 changes: 3 additions & 1 deletion 2d.particle.tracking.kernel.r
    Original file line number Diff line number Diff line change
    @@ -3,7 +3,9 @@ require(KernSmooth)
    require(maptools)

    ## load the simulation results. See the previous blog post
    ## how this file was created.
    ## 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
  3. pepijn-devries revised this gist Jul 19, 2016. 1 changed file with 1 addition and 1 deletion.
    2 changes: 1 addition & 1 deletion 2d.particle.tracking.kernel.r
    Original file line number Diff line number Diff line change
    @@ -90,7 +90,7 @@ for (i in 1:length(part.results)) {
    ## close the graphics device
    dev.off()

    ## print ration between the current time index and the total
    ## print ratio between the current time index and the total
    ## number of indices. This serves as a progress indicator.
    print(i/length(part.results))
    }
  4. pepijn-devries created this gist Jul 18, 2016.
    96 changes: 96 additions & 0 deletions 2d.particle.tracking.kernel.r
    Original file line number Diff line number Diff line change
    @@ -0,0 +1,96 @@
    require(raster)
    require(KernSmooth)
    require(maptools)

    ## load the simulation results. See the previous blog post
    ## how this file was created.
    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 ration between the current time index and the total
    ## number of indices. This serves as a progress indicator.
    print(i/length(part.results))
    }