In this post I will discuss an implementation of sequential Gaussian simulation (SGS) from the field of geostatistics. Geostatistics is simply a statistical consideration of spatially distributed data. Sequential Gaussian simulation is a technique used to “fill in” a grid representing the area of interest using a smattering of observations, and a model of the observed trend. The basic workflow incorporates three steps:

Modeling the measured variation using a semivariogram
Using the semivariogram to perform interpolation by kriging
Running simulations to estimate the spatial distribution of the variable(s) of interest

We will use code from previous posts to perform simple kriging and z-score transformations, which will be available at cjohnson318/geostatistics. The four-inch paint brush version of two dimensional sequential Gaussian simulation is as follows:

Sequential Gaussian Simulation

Define a grid, this may be fine or coarse
Place the z-transformed initial data into the nearest grid cells
Visit the rest of the grid cells at random and perform kriging using all of the data available
Back-transform the data to retrieve the approximate distribution

There are two important things to note here:

SGS expects normally distributed data, so if the data is not normally distributed, we use a z-score transformation before SGS and a back-transformation afterwards
In step 3. above, as we move randomly about the grid, we add each newly kriged estimate to the data used for kriging each time

The second point means that our data set for kriging increases by one data point for each step until we visit all of the cells in our grid. The randomization allows us to run multiple simulations and then take their mean.

Data

For this example, I will be using the 140 data points found in the data set cluster.dat prepared by Clayton V. Deutsch. This data set and others in the GSLIB: Geostatistical Software Library and User’s Guide may be found here.

I have prepared a Python function for reading this data format.

def read_data( fn ):
    f = open( fn, "r" )
    # title of the data set
    title = f.readline()
    # number of variables
    nvar = int( f.readline() )
    # variable names
    columns = [ f.readline().strip() for i in range( nvar ) ]
    # make a list for the data
    data = list()
    # for each line of the data
    while True:
        # read a line
        line = f.readline()
        # if that line is empty
        if line == '':
            # the jump out of the while loop
            break
        # otherwise, append that line to the data list
        else:
            data.append( line )
    # strip off the newlines, and split by whitespace
    data = [ i.strip().split() for i in data ]
    # turn a list of list of strings into an array of floats
    data = np.array( data, dtype=np.float )
    # combine the data with the variable names into a DataFrame
    df = pandas.DataFrame( data, columns=columns,  )
    return df

We’ll import some packages and load our data into our workspace using the following lines.

import time
import numpy as np
import scipy.stats
import scipy.optimize
import scipy.interpolate

import krige
import utilities

import random

clstr = read_data( 'cluster.dat' )

Gridding

The next step is to define a grid. The data in cluster.dat lies in 50 by 50 foot region, so we can form a grid from zero to fifty with thirty steps as follows.

# we want 30 steps between 0 and 50, in both dimensions
nx = 50 ; xsteps = 30
ny = 50 ; ysteps = 30
xrng = np.linspace(0,nx,num=xsteps)
yrng = np.linspace(0,ny,num=ysteps)

Randomization

This part creates a random path through the spatial domain. First I shuffle a vector of indices, attach some data to those indices, and then sort them in order to shuffle the data. The data is the cell address, an ordered pair of integers, and the physical location, which is an ordered pair of floats, describing a coordinate in feet.

# total number of cells
N = xsteps * ysteps

# an array of values 0 to N-1
idx = np.arange( N )

# randomize idx in place
random.shuffle( idx )

# create a list named path
path = list()

# initialize an index for idx
t = 0

# walk through the 2D grid..
for i in range( xsteps ):
    for j in range( ysteps ):
        
        # the first element is the shuffled index,
        # then the cell address, then the physical location
        path.append( [ idx[t], (i,j), (xrng[i],yrng[j]) ] )
        
        # increment the index
        t += 1

# sort the shuffled index, shuffling the cell 
# addresses and the physical locations (in feet)
path.sort()

Simulation

This part transforms the data before it goes into the sgs() function, and bundles the location information with the the transformed variable. The sgs() function then visits all of the cells in the grid at random, calculates a kriging estimate, adds that to the original data, and repeats.

# x and y coordinates in feet
locations = np.array( clstr[['Xlocation','Ylocation']] )

# value of the variable of interest, e.g., porosity or permeability
variable = np.array( clstr['Primary'] )

# the z-score transformation
norm, inv, param, mu, sd = to_norm( variable, 1000 )

# the data to be kriged and appended to
# during the sequential Gaussian simulation
data = np.vstack((locations.T,norm)).T

def sgs( data, bw, path, xsteps, ysteps ):
    '''
    Input:  (data)   <N,3> NumPy array of data
            (bw)     bandwidth of the semivariogram
            (path)   randomized path through region of interest
            (xsteps) number of cells in the x dimension
            (ysteps) number of cells in the y dimension
    Output: (M)      <xsteps,ysteps> NumPy array of data
                     representing the simulated distribution
                     of the variable of interest 
    '''
    # create array for the output
    M = np.zeros((xsteps,ysteps))
    # generate the lags for the semivariogram
    hs = np.arange(0,50,bw)
    # for each cell in the grid..
    for step in path :
        # grab the index, the cell address, and the physical location
        idx, cell, loc = step
        # perform the kriging
        kv = krige( data, spherical, hs, bw, loc, 4 )
        # add the kriging estimate to the output
        M[ cell[0], cell[1] ] = kv
        # add the kriging estimate to a spatial location
        newdata = [ loc[0], loc[1], kv ]
        # add this new point to the data used for kriging
        data = np.vstack(( data, newdata ))
    return M

We can time the execution of our code as follows:

t0 = time.time()
M = sgs( data, 5, path, xsteps, ysteps )
t1 = time.time()
print (t1-t0)/60.0

Visualization

This is the best part. This is what the simulation looks like before we transform it back from the z-score transformation.

matshow( M[:,::-1].T, cmap=my_cmap )
colorbar();
xs, ys = str(xsteps), str(ysteps)
xlabel('X') ; ylabel('Y') ; title('Raw SGS '+xs+'x'+ys+'\n')
savefig( 'sgs_untransformed_'+xs+'_by_'+ys+'.png', fmt='png', dpi=200 )

This is what the simulation looks like after we transform it back from the z-score transform.

z = from_norm( M.ravel(), inv, param, mu, sd )
z = np.reshape( z, ( xsteps, ysteps ) )
matshow( z[:,::-1].T, cmap=my_cmap )
colorbar() ;
xs, ys = str(xsteps), str(ysteps)
xlabel('X') ; ylabel('Y') ; title('SGS '+xs+'x'+ys+'\n')
savefig( 'sgs_transformed_'+xs+'_by_'+ys+'.png', fmt='png', dpi=200 )

Connor Johnson

Two Dimensional Sequential Gaussian Simulation in Python