Last active
July 21, 2016 08:33
Revisions
-
pepijn-devries revised this gist
Jul 21, 2016 . 1 changed file with 1 addition and 1 deletion.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal 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)) } -
pepijn-devries revised this gist
Jul 21, 2016 . 1 changed file with 3 additions and 1 deletion.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal 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. 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 -
pepijn-devries revised this gist
Jul 19, 2016 . 1 changed file with 1 addition and 1 deletion.There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal 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 ratio between the current time index and the total ## number of indices. This serves as a progress indicator. print(i/length(part.results)) } -
pepijn-devries created this gist
Jul 18, 2016 .There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode charactersOriginal 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)) }