If the stratigraphic structure is turned on in the XmL input file, Badlands produces sedimentary layers Hdf5 files. The stratigraphic layers are defined on a regularly spaced grid and a layer is recorded at each layer time interval given by the user.
Here we show how we can visualise quickly the structure of the stratigraphic layer in an IPython notebook.
In [ ]:
%matplotlib inline
import numpy as np
import colorlover as cl
# Import badlands grid generation toolbox
import pybadlands_companion.stratalAnalyse as strata
# display plots in SVG format
%config InlineBackend.figure_format = 'svg'
First we need to load one of the stratigraphic file. The files are located in the h5/ folder in the simulation main output folder and are named using the following convention:
sed.timeT.pX.hdf5with T the display time index and X the number of the partition (used in the parallel version). In cases where you ran your simulation in parallel you will also need to give the number of CPUs used (cpus).
To load a file you will need to give the folder name and the number of processors used in your simulation.
For more information regarding the function uncomment the following line.
In [ ]:
#help(strata.stratalSection.__init__)
In [ ]:
folder = '/workspace/volume/Examples/delta/output_0/h5/'
strat = strata.stratalSection(folder)
Then we need to load a particular output time interval (this is the T parameter in the hdf5 file name convention).
Note
This number is not always the number of sedimentary layers for this particular time step as you could have chosen in the input file to have more than 1 sedimentary layer recorded by output interval!
In [ ]:
#help(strat.loadStratigraphy)
In [ ]:
strat.loadStratigraphy(50)
We then slice the stratigraphic mesh to visualise the sedimentary architecture along a given cross-section.
To create the cross-section you will need to provide the position of the segment in the simulation space (xo,yo) and (xm,ym), a gaussian filter value for smoothing (gfilt a value of 0 can be used for non-smoothing) and the resolution of the cross-section (based on a number of points: nbpts).
In [ ]:
# Position of the cross-section
x1 = 0
y1 = 12000
x2 = 24000
y2 = 12000
# Interpolation parameters
gfilt = 5
nbpts = 1000
In [ ]:
#help(strat.buildSection)
In [ ]:
strat.buildSection(xo = x1, yo = y1, xm = x2, ym = y2, pts = nbpts, gfilter = gfilt)
In [ ]:
#help(strata.viewSection)
In [ ]:
strata.viewSection(width = 1000, height = 600, cs=strat,
rangeX=[2000,9000], rangeY=[-400,100],
linesize=0.5, title='stratigraphic stack coloured by time')
There are several models of systems tracts within depositional sequences, here we use the most simple one called the four systems tract model:
For each of these system tracts we define a given color. We use the 'colorlover' library link
In [ ]:
#colormap = cl.scales['9']['seq']['BuPu']
#colorrgb = cl.to_rgb( colormap )
#cLST = colormap[1]
#cTST = colormap[3]
#cHST = colormap[5]
#cFST = colormap[7]
cLST = 'rgb(220,213,166)' #colormap[1]
cTST = 'rgb(161,206,146)' #colormap[3]
cHST = 'rgb(72, 106,162)' #colormap[5]
cFST = 'rgb(119,11,116)' #colormap[7]
In [ ]:
seafile = '/workspace/volume/Examples/delta/data/sea.csv'
stime,sea = strata.viewSea(seafile)
Based on the sea-level curve, we define a series of system tracts periods.
In [ ]:
LST1 = np.array([0,90000],dtype=int)
TST1 = np.array([90000,160000],dtype=int)
HST1 = np.array([160000,220000],dtype=int)
FST1 = np.array([220000,340000],dtype=int)
LST2 = np.array([340000,420000],dtype=int)
TST2 = np.array([420000,500000],dtype=int)
The simulation ran for 500,000 years and each sedimentary layer is produced every 2500 years which give us 200 layers.
We can use the stratal mesh to retreive this information:
In [ ]:
LST1 = LST1/2500
TST1 = TST1/2500
HST1 = HST1/2500
FST1 = FST1/2500
LST2 = LST2/2500
TST2 = TST2/2500
print 'Number of layers recorded in the output :',strat.nz
We can now associate each layer to a given system tract by defining a specific color map:
In [ ]:
STcolors = []
for k in range(LST1[0],LST1[1]):
STcolors.append(cLST)
for k in range(TST1[0],TST1[1]):
STcolors.append(cTST)
for k in range(HST1[0],HST1[1]):
STcolors.append(cHST)
for k in range(FST1[0],FST1[1]):
STcolors.append(cFST)
for k in range(LST2[0],LST2[1]):
STcolors.append(cLST)
for k in range(TST2[0],TST2[1]):
STcolors.append(cTST)
Again we use plotly to visualise the cross-section and the associated system tracts:
In [ ]:
#help(strata.viewSectionST)
In [ ]:
strata.viewSectionST(width = 1000, height = 600, cs=strat, colors=STcolors,
rangeX=[2000,9000], rangeY=[-400,100],
linesize=0.5, title='system tract cross-section')
It is also possible to plot a wheeler-like diagram of the fluvial to marine depositional sequence using the previous cross-section (dip oriented stratigraphic succession). For that we need to define depositional environment types such as:
For each of these systems we will define their elevations range (this is relative to sea-level):
In [ ]:
envi = []
S = np.array([-100.,-60.])
envi.append(S)
LD = np.array([-60.,-30.])
envi.append(LD)
UD = np.array([-20.,-5.])
envi.append(UD)
MS = np.array([-5.,20.])
envi.append(MS)
BS = np.array([20.,100.])
envi.append(BS)
An associate a color value to each of them:
In [ ]:
WHcolors = []
cS = 'rgb(250,250,250)' # white
WHcolors.append(cS)
cLD = 'rgb(72, 106,162)' # blue
WHcolors.append(cLD)
cUD = 'rgb(158,202,225)' # light blue
WHcolors.append(cUD)
cMS = 'rgb(220,213,166)' # yellow
WHcolors.append(cMS)
cBS = 'rgb(133,193,143)' # green
WHcolors.append(cBS)
The wheeler diagram plot the extent of each of the depositional environments over time, which means that we need to define the time based on the simulation parameters. This is done by getting the following:
In [ ]:
times = np.linspace(0,500000,strat.nz)
In [ ]:
strata.viewWheeler(width = 800, height = 800, cs = strat, time = times, colors = WHcolors,
rangeE=envi, rangeX=[2000,9000], rangeY=[times[0],times[-1]], contourdx = 10,
title = 'wheeler-type diagram')
In [ ]: