Landscape simulation in S

Various ways of generating 1- and 2-dimensional random landscapes with a specified second-order correlation structure (spatial autocovariance, spatial power spectrum). In broad outline, there are three major classes of methods, which trade off conceptual simplicity and ease of programming for efficiency and power (see Cressie, Statistics for Spatial Data, for more detail on anything related to spatial statistics): I have chosen to go with the happy medium--FFT algorithms.

The S language

I've chosen to write these algorithms in the S language. S is a statistics/graphics language, currently available (as S-PLUS) from Mathsoft for Windows or Unix, for a lot of money. There is also a public domain version called R coming along; at the moment there are good versions for Windows and Unix and a backward version for the Mac (alas, a year behind; I've only tried the SunOS/Unix version). The only thing R doesn't (yet) have that I want is 3-D perspective plotting.

It would be about the same amount of work (i.e. not much) to re-write the algorithms in MATLAB or some other high-level language (you could probably do it in Mathematica but it would almost certainly be excruciatingly slow), and it would be more work but would run much faster to code them in C or C++ (or FORTRAN); FFT routines should be available from netlib or Numerical Recipes (Press et al., Cambridge University Press).

I hope the examples below should be relatively coherent even if you don't know S: a couple of things to note are that "<-" is the assignment symbol in S, and that "outer" is an "outer product" function for generating tables. ("rnorm" and "runif" generate vectors of Gaussian and uniform deviates respectively. "apply" applies a function to a vector or matrix.)

1-D fractal landscapes

Algorithm as prescribed by Keitt or Hastings and Sugihara. Note that I'm specifying the standard deviation (n^{-(H+0.5)}), not the variance (n^{-(2H+1)}). 
  fland1 <- function(H=0.5,n=16) {
     b <- H+0.5
     sdvec <- (1:n)^(-b)
     return(Re(fft(complex(modulus=rnorm(n,mean=0,sd=sdvec),
                           arg=runif(n,-pi,pi)),
               inverse=T)))
  }
Examples with H ranging from 0.1 to 0.9 (results from the following). 
par(mfrow=c(3,3))
  for (i in seq(0.1,0.9,by=0.1))
   plot(fland1(H=i,n=1024),type="l")

2-D fractal landscapes

Sources as above. Landscapes of up to about 256x256 are feasible in S (it takes a few seconds, but certainly not more than a minute); after that it runs out of memory. This suggests that you could probably squeeze out 1024x1024 with no real problem in raw C or C++ if you needed to. 
 fland2 <- function(H=0.5,n=16) {
     b <- H+0.5
     sdvec <- (1:n)^(-b)
     sdmat <- outer(sdvec,sdvec,"*")
     return(Re(fft(matrix(complex(modulus=rnorm(n^2,mean=0,sd=sdmat),
                           arg=runif(n^2,-pi,pi)),
                     nrow=n),
               inverse=T)))
  }
# example
 image(fland2(H=0.7,n=256))
Does anyone know why this looks so "rectangular"? (Actually, I think I do now -- and I'm surprised that Hastings and Sugihara don't think about this, if I'm right. The formula for the power spectrum is S(m,n) = c*m^(-B)*n^(-B) (with B=2*H+1), which says that power falls off linearly in the x and y directions, not radially. I would guess that S(m,n) = c*(m^2+n^2)^(-B/2) would work better.)

This is what it looks like (for example) if I take the modulus instead of the real part of the landscape ...

Non-fractal 2-D landscape

Here's an example of generating a non-fractal landscape with a general (non-fractal) correlation structure. The first bit of code takes any landscape at all, FFTs it, randomizes the phase, and inverse-FFTs it. There are other possibilities for randomizing, e.g. generating a chi-square(2) deviate with the right expected value for the modulus instead of leaving it unchanged. 
 # randomize phase of a 2-D landscape
  rland1 <- function(x) {
    n <- nrow(x)   # linear dimension
    l <- length(x) # total number of nodes
    y <- fft(x)
    z <- matrix(complex(modulus=Mod(y),arg=runif(l,-pi,pi)))
    y <- fft(z,inverse=T)/l
    return(matrix(Mod(y),nrow=n))
  }
Once we have that code, we can generate a "landscape" which is really just the two-dimensional realization of the radial correlation function. This code takes a vector or a function of distance representing the radial ACF, applies it to a 2-D grid, and passes the result to the previous code. 
simland2 <- function(r,n) {
## given a correlation-by-(radial)-distance function r
##  (which can be either a vector or a function)
##  simulate a 2-D landscape with linear size n 
## substitute into a grid of radial distances
 g <- outer(1:n,1:n,function(x,y)sqrt(x^2+y^2))
 if (is.function(r)) {
   gr <- apply(g,c(1,2),r)
  }  else if (is.numeric(r)) {
   gr <- matrix(r[round(g)],nrow=n) }
 rland1(gr)  # pass it to a function which FFTs, randomizes phases,
             # and inverse-FFTs
}
A sample, using correlation function exp(-x/20)+0.1*N(0,1): 
 n <- 128
# pick some correlation function
##  test: exponential + Gaussian noise
   r <- exp(-(1:(2*n))/20)+0.1*rnorm(2*n)
   r[r>1] <- 1
g2 <- simland2(r,n)
With this correlation function:

We get this landscape:

Of course, I should check this by going back and computing the ACF or the variogram from the simulated landscape.

Update, June 1999: playing around with this some more with one-dimensional landscapes, actually checking the covariances, gives the following conclusions: