Simulating elevation data in R

I was looking for a quick way to create dummy digital elevation model (DEM) datasets in R, and a kind soul in StackOverload shared the following code. It is an implementation of the Perlin noise algorithm written by Ken Parlin in the early 1980s and is now widely used to create realistic looking CGI landscapes and other natural patterns.

The function:

perlin_noise <- function( 
  n = 5,   m = 7,    # Size of the grid for the vector field
  N = 100, M = 100   # Dimension of the image
) {
  # For each point on this n*m grid, choose a unit 1 vector
  vector_field <- apply(
    array( rnorm( 2 * n * m ), dim = c(2,n,m) ),
    2:3,
    function(u) u / sqrt(sum(u^2))
  )
  f <- function(x,y) {
    # Find the grid cell in which the point (x,y) is
    i <- floor(x)
    j <- floor(y)
    stopifnot( i >= 1 || j >= 1 || i < n || j < m )
    # The 4 vectors, from the vector field, at the vertices of the square
    v1 <- vector_field[,i,j]
    v2 <- vector_field[,i+1,j]
    v3 <- vector_field[,i,j+1]
    v4 <- vector_field[,i+1,j+1]
    # Vectors from the point to the vertices
    u1 <- c(x,y) - c(i,j)
    u2 <- c(x,y) - c(i+1,j)
    u3 <- c(x,y) - c(i,j+1)
    u4 <- c(x,y) - c(i+1,j+1)
    # Scalar products
    a1 <- sum( v1 * u1 )
    a2 <- sum( v2 * u2 )
    a3 <- sum( v3 * u3 )
    a4 <- sum( v4 * u4 )
    # Weighted average of the scalar products
    s <- function(p) 3 * p^2 - 2 * p^3
    p <- s( x - i )
    q <- s( y - j )
    b1 <- (1-p)*a1 + p*a2
    b2 <- (1-p)*a3 + p*a4
    (1-q) * b1 + q * b2
  }
  xs <- seq(from = 1, to = n, length = N+1)[-(N+1)]
  ys <- seq(from = 1, to = m, length = M+1)[-(M+1)]
  outer( xs, ys, Vectorize(f) )
}

The function perlin_noise(xFieldSize, yFieldSize, x_grid, y_grid) produces a matrix of Parlin noise values that can be quickly visualised:

image(1:150, 1:100, perlin_noise(10,7,150,100), col = grey.colors(100), asp=T)

Producing something like this:

————————–

Appendum (14.3.13)

As the aforementioned kind soul (Vincent Zoonekynd) in *StackOverflow* has added, it is possible to apply a fractal dimension to the function above to achieve a better output.

scaler <- .6
fractDim <- 7
m <- perlin_noise(2, 2, 2^fractDim, 2^fractDim)
for( i in 2:fractDim ){
  m <- m + scaler^i * perlin_noise(2^i, 2^i, 2^fractDim, 2^fractDim)
}

Each iteration overlays a new fractal level of noise:

The {rgl} package (here I've used the surface3d function) is an effective way to visualise the data in 3D: