make a point set, put a parameteric movement kernel on it and plot

lscape.R

library(viridis)
library(spatstat)
library(truncdist)

# source: https://github.com/slwu89/PlosCompBio_2013/blob/master/fig2/code/functionsModel.R
points.clustered = function(n, meanParents = 10, clusteredness = .25, ...){
  meanDist = clusteredness / sqrt(meanParents)
  meanChildren = n / meanParents
  ps = rMatClust(meanParents, meanDist, meanChildren, ...)
  while(ps$n != n){
    ps = rMatClust(meanParents, meanDist, meanChildren, ...)
  }
  return(ps)
}

# 2d plane
xy_plane = owin(xrange = c(0,1),yrange = c(0,1))

# poisson scatter, matérn clustering, overdispersed
points = points.clustered(n = 20,meanParents = 5, clusteredness = 0.25, win = xy_plane)
xy_points = cbind(points$x,points$y)

# movement kernel
exp_fit <- function(d,q,up=1){
  f_opt = function(x){
    abs(d - qexp(p = q,rate = x))
  }
  sol = optimise(f = f_opt,lower = 0,upper = up,maximum = FALSE)
  return(sol$minimum)
}

exp_kern <- exp_fit(d = 0.1,q = 0.75)

# movement matrices
dist <- as.matrix(dist(xy_points,diag = TRUE,upper = TRUE))
movement <- apply(X = dist,MARGIN = 1,FUN = function(x){dtrunc(x = x,spec = "exp",a = 1e-12,b = Inf, rate = 1/exp_kern)})
movement <- movement/rowSums(movement)


###############################################################################
# Plot landscape
###############################################################################

movement_quantile = cut(as.vector(movement),breaks=quantile(as.vector(movement),probs=seq(0, 1, 0.2)),include.lowest = TRUE,labels = FALSE)
movement_color <- matrix(data = plasma(n = length(unique(movement_quantile)),alpha=0.5)[movement_quantile],nrow = nrow(movement),ncol = ncol(movement))

par(bg = grey(0.15))
plot.new()
for(i in 1:ncol(movement)){
  for(j in 1:nrow(movement)){
    segments(x0 = xy_points[i,1],y0 = xy_points[i,2],
             x1 = xy_points[j,1],y1 = xy_points[j,2],
             col = movement_color[i,j],lty = 1.15,lwd = 1.15)
  }
}
points(xy_points,pch=21,cex=5,bg=grey(level = 0.95,alpha = 0.85),col="white")
text(xy_points,labels=as.character(1:20),col="black")
par(bg = "white")