SgFm - Sediment Transport @ Geological Scale

Post Processing Tips Exploring Lecode Output

Lecode output structure

For each <outputTime> step, Lecode generates a series of output files. These files are stored on a local folder which name is defined in the <struct-output> structure within the XmL Input file.
The <OutputDirectory> contains two subfolders:

runfiles which contains information for checkpointing and restart.

outputs which contains information about stratigraphy, flow walkers and TIN.

The outputs directory mostly contains binary HdF5 files written at each time step and for each processor. For each time step, these Hdf5 files are bundle within a XmF file containing the collection of attributes (both for nodes and cells). On top of the output structure, a Xdmf file is used to wrap time-series informations.

To visualise Lecode results a visualisation platform able to read these types of files is required. We recommend using Paraview or Visit .

To start looking at your outputs and the model evolution through time, it is only necessary to load the _series.xdmf files. Depending on the set of parameters defined in your simulation, the following time series files will exist:

Stratal_series.xdmf file containing a combination of the Stratigraphic information (mean grain size, sediment proportion, cell thicknesses, layer time...) over all processors and for all time steps;

TIN_series.xdmf file containing the Triangular Irregular Network grid information for all time steps and the associated flow accumulation dataset;

FW_series.xdmf file containing the Flow Walkers information (flow velocity, stream width and depth ...) over all processors and for all time steps;

Ocean_series.xdmf file representing Sea Level fluctuations (if any) for all time steps.

Using Paraview scripts

To look at simulation results such as time series, XmF and XdmF files we recommend to use Paraview. ParaView is an open-source, multi-platform data analysis and visualisation application built on top of VTK.

In addition to the built-in filters and functionalities embedded in this visualisation software, we have developed several post-processing scripts based on NumPy and VTK libraries.
These scripts help to generate additional informations from Lecode simulation (e.g., consecutive time stepping elevation differences to estimate erosion - deposition surfaces through time, or cumulative thickness computation for specific types of materials).

Feel free to contact us if you want to get some of these scripts, we will be happy to share ...

Cumulative change of elevation programmable filter example

<SeismicLine> Plugin

Lecode Seismic Plugin makes it possible to extract a given cross section within the modelled stratigraphic record along the main axes of the domain. Two files are generated by Lecode in the runfiles output folder.
The first records for each node along the desired seismic line: its position on the stratigraphic layer, its mean grain size and its porosity value.
The second gives the envelop (bottom and top locations) of the stratigraphic record.

Below we provide a seismic script which shows a really simple example of workflow that can be used to extract the seismic attributes from the model for further processing in a seismic modelling framework such Fatiando or RSF Madagascar.

Click here to see a Seismic script example

Our future goals are to develop some algorithms based on geophysical forward modelling, uncertainty quantification and inversion to minimise global mismatch between stratigraphic predictions and observations.

<RMSModel> Plugin

The RMSModel Plugin is designed to facilitate the export of Lecode outputs into a geo-cellular model. This later can then be used to perform reservoir and fluid flow modelling.
Lecode created geo-cellular models have been successfully exported and tested with:

Roxar RMS geological modelling package and,

TemisFlow a multi-dimensional basin modelling package from BeicipFranlab.

"""
Transforming Lecode irregularly sampled cross section dataset to a regularly spaced grid.
python lecode_gridding.py -s seismic.csv -e envelop.csv -l 54375 -r 108750 -b -2740 -t -470 -d 5
"""
import sys, math, numpy, getopt
from scipy import interpolate
from numpy import *
from scipy import *
import matplotlib 
from fatiando import gridder, utils
from fatiando.vis import mpl

#def main(argv):

# Define font size and type
font = {'family' : 'normal',
        'size'   : 6}
matplotlib.rc('font', **font)

total = len(sys.argv)

if total == 1:
    print 'Some required fields are missing'
    print 'Use: python lecode_gridding.py -h for parameter declaration'
    sys.exit(2)
elif total > 2 and total < 15:
    print 'Some required fields are missing'
    print 'Use: python lecode_gridding.py -h for parameter declaration'
    sys.exit(2)

try:
    opts, args = getopt.getopt(sys.argv[1:],"hs:e:l:r:b:t:d:")
except getopt.GetoptError:
    print 'Usage: python lecode_gridding.py -s <seismic> -e <envelop> -l <hmin> -r <hmax> -b <vmin> -t <vmax> -d <dx>'
    print '-s : name of the seismic file exported from Lecode simulation'
    print '-e : name of the envelop file exported from Lecode simulation'
    print '-l : regular grid minimal horizontal value (left)'
    print '-r : regular grid maximal horizontal value (right)'
    print '-b : regular grid minimal vertical value (bottom)'
    print '-t : regular grid maximal vertical value (top)'
    print '-d : regular grid spacing'
    print 'Note: you need to ensure that the vertical and horizontal are both multiple of the specify spacing'
    sys.exit(2)
 
for opt, arg in opts:
    if opt == '-h':
    	print 'Usage: python lecode_gridding.py -s <seismic> -e <envelop> -l <hmin> -r <hmax> -b <vmin> -t <vmax> -d <dx>'
    	print '-s : name of the seismic file exported from LECODE simulation'
    	print '-e : name of the envelop file exported from LECODE simulation'
    	print '-l : regular grid minimal horizontal value (left)'
    	print '-r : regular grid maximal horizontal value (right)'
    	print '-b : regular grid minimal vertical value (bottom)'
    	print '-t : regular grid maximal vertical value (top)'
    	print '-d : regular grid spacing'
    	print 'Note: you need to ensure that the vertical and horizontal are both multiple of the specify spacing'
        sys.exit()
    elif opt == "-s":
        ifile = arg
    elif opt == "-e":
        ofile = arg
    elif opt == "-l":
        hmin = float( arg )
    elif opt == "-r":
        hmax = float( arg )
    elif opt == "-b":
        vmin = float( arg )
    elif opt == "-t":
        vmax = float( arg )
    elif opt == "-d":
        dx = float( arg )
        
# Open and read file from std, and assign data (horizontal, vertical, mean grain size & porosity) 
# columns to an array.
x,y,mgz,poro = loadtxt( ifile, skiprows=0, delimiter=',', unpack = True)
            
# Extensions
hext = int( ( hmax - hmin ) / dx )
vext = int( ( vmax - vmin ) / dx )

# !!! Be careful the shape is defined as ( ny, nx )
shape = ( int( vext ) + 1, int( hext ) + 1 )

# Interpolation
grdx, grdy, gmgz = gridder.interp( x, y, mgz, shape, area=(hmin, hmax, vmin, vmax) , algorithm='cubic', extrapolate=True)
gx, gy, gporo = gridder.interp( x, y, poro, shape, area=(hmin, hmax, vmin, vmax) , algorithm='cubic', extrapolate=True)

# Open and read file from std, and assign first three (horizontal, low, high) columns to an array.
line,low,high = loadtxt( ofile, skiprows=0, delimiter=',', unpack = True)
xnew = linspace(hmin,hmax,int(rint((hmax-hmin)/dx+1)))

# Create the envelop bottom function
bfunc = interpolate.splrep(line,low,s=0)
envbot = interpolate.splev(xnew,bfunc,der=0)

# Create the envelop top function
tfunc = interpolate.splrep(line,high,s=0)
envtop = interpolate.splev(xnew,tfunc,der=0)

# Mask interpolated regular data using the envelope
k = 0
p = 0
vp = zeros(grdx.size)
vs = zeros(grdx.size)

# Here is a simple example which link porosity to velocity using Dominico formula 
for i in nditer(grdx):
    if k > hext:
        k = 0
    if grdy[ p ] < envbot[ k ] :		
        gmgz[ p ] = nan
        gporo[ p ] = nan
        vp[ p ] = nan
        vs[ p ] = nan
    elif grdy[ p ] > envtop[ k ] :		
        gmgz[ p ] = nan
        gporo[ p ] = nan
        vp[ p ] = nan
        vs[ p ] = nan
    else:
       	vp[ p ] = 1.0 / (0.163 + 0.365 * gporo[ p ] ) 
        vs[ p ] = 1.0 / (0.224 + 0.889 * gporo[ p ] )
    k = k + 1
    p = p + 1
       
savetxt('seismic_reg.csv',column_stack((grdx, grdy, gmgz, gporo)), fmt='%g', delimiter=',')
                
# Plot interpolated values
mpl.figure()

mpl.subplot(5, 1, 1)
mpl.title("Lecode line scattered points and cubic-spline envelop")
mpl.scatter(x, y, c=gmgz, label='points')
mpl.plot(xnew,envbot,c="cyan",linewidth=2.5)
mpl.plot(xnew,envtop,c="magenta",linewidth=2.5)
mpl.colorbar()

mpl.subplot(5, 1, 2)
mpl.title("Density map")
mpl.contourf(grdx, grdy, gmgz, shape, 200)
mpl.colorbar()

mpl.subplot(5, 1, 3)
mpl.title("Net 2 Gross map")
mpl.contourf(grdx, grdy, gporo, shape, 200)
mpl.colorbar()
mpl.axis('scaled')
mpl.plot(x, y, '.k', label='data points', alpha=0.4, ms=4 )
mpl.legend(loc='lower right', numpoints=1)

mpl.subplot(5, 1, 4)
mpl.title("S-wave velocity map")
mpl.contourf(grdx, grdy, vs, shape, 200)
mpl.colorbar()

mpl.subplot(5, 1, 5)
mpl.title("P-wave velocity map")
mpl.contourf(grdx, grdy, vp, shape, 200)
mpl.colorbar()

mpl.show()

'''
This Paraview Programmable Filter will produce a time series of 
cumulative elevation changes for your area (erosion/deposition area).
It could easily be adapted to look at the instantaneous elevation change.

Prior to applying the filter, you will need to use the - Merge Blocks -
filter from Paraview with the Points unchecked in the Properties list.

To use the script you will need to choose - vtkPolyData - in the
Output Data Set Type button.

After applying the filter you will have to use the - Delaunay2D -
filter from Paraview to visualise the cumulative elevation change.

'''
# Import Paraview Numpy Library
from paraview.vtk.dataset_adapter import numpyTovtkDataArray

data = inputs[0]
newPoints = vtk.vtkPoints()

# Offset should correspond to the initial deposit layer 
# It will work if you started with an uniform thickness
offset = -10000

# Definition of the current top layer
youngestIndex = max(data.PointData['Layer Index'])

# Find the Points within the stratigraphy corresponding to the 
# top layer
youngestLayer = data.PointData['Layer Index'] == youngestIndex

# Find the Points corresponding to the bottom layer
# Note in case you want to look at the instantaneous change
# replace it by youngestIndex - 1
oldestLayer = data.PointData['Layer Index'] == 1

# Define the X coordinates
youngestX = data.Points[:,0][youngestLayer].transpose()

# Define the Y coordinates
youngestY = data.Points[:,1][youngestLayer].transpose()
youngestZ = data.Points[:,2][youngestLayer].transpose()

# Define the Z coordinates
# The cumulative change of elevation will coincide with the
# current topography
oldestZ = data.Points[:,2][oldestLayer].transpose()


# Define the cumulative change of elevation
deltaZ = youngestZ - oldestZ + offset

# Create the VTK array grid
coords = hstack((youngestX, youngestY, youngestZ))
Points = numpyTovtkDataArray(coords,name="Points")
newPoints.SetData(Points)
output.SetPoints(newPoints)

# Add the cumulative elevation change property
output.PointData.append(deltaZ, "deltaZ")