SIG_mesh_logo.R

# 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()