Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
# make a nice viridis GMRF tesselation for Richard
rm(list = ls())
set.seed(1)
library(INLA)
library(raster)
library(viridis)
library(fields)
# grid sizes for sampling the GRF and for the final image
nrow <- ncol <- 3000
nrow_sim <- ncol_sim <- 100
grid <- list(x = seq(0, 1, length.out = ncol_sim),
y = seq(0, 1, length.out = nrow_sim))
obj <- Exp.image.cov(grid = grid, theta = 0.1, setup = TRUE)
r <- raster(sim.rf(obj))
image <- raster(nrow = nrow, ncol = ncol)
extent(image) <- c(0, 1, 0, 1)
pts <- expand.grid(seq(0, 1, length.out = 10),
seq(0, 1, length.out = 10))
pts <- pts + cbind(rnorm(100, 0, 0.1),
rnorm(100, 0, 0.1))
# make an inla mesh
sp <- as(extent(r), 'SpatialPolygons')
mesh <- inla.mesh.2d(loc = pts,
boundary = inla.sp2segment(sp),
max.edge = 0.1,
offset = 0)
# sample GRF at nodes
z <- extract(r, mesh$loc[, 1:2])
# get projection to raster
image_coords<- xyFromCell(image, 1:ncell(image))
A <- inla.spde.make.A(mesh, loc = image_coords)
# instead of linear interpolation, average the three node values
A2 <- A
A2@x[A2@x > 0] <- 1/3
image[] <- (A2 %*% z)[, 1]
# make a pretty plot
png('triangulation_QESIG_logo.png',
width = 1000,
height = 1000,
pointsize = 30)
par(mar = rep(0, 4))
image(image,
col = viridis(1000),
asp = 1,
axes = FALSE,
xlab = '',
ylab = '')
plot(mesh,
add = TRUE,
edge.color = grey(0.4),
lwd = 1,
draw.segments = FALSE)
dev.off()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment