This notebooks looks at the main striplog object. For the basic objects it depends on, see Basic objects.
First, import anything we might need.
In [1]:
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
import striplog
striplog.__version__
Out[1]:
In [2]:
from striplog import Legend, Lexicon, Interval, Component
In [3]:
legend = Legend.builtin('NSDOE')
lexicon = Lexicon.default()
In [4]:
from striplog import Striplog
imgfile = "M-MG-70_14.3_135.9.png"
In [5]:
strip = Striplog.from_img(imgfile, 14.3, 135.9, legend=legend)
strip
Out[5]:
In [6]:
strip.plot(legend, ladder=True, aspect=3)
This striplog doesn't have any gaps...
In [7]:
strip.find_gaps() or "No gaps!"
Out[7]:
But we can make some by deleting indices:
In [8]:
del strip[[2, 7, 12]]
strip.find_gaps()
Out[8]:
We can also get a list of the indices of intervals that are followed by gaps (i.e. are directly above gaps in 'depth' order, or directly below gaps in 'elevation' order).
In [9]:
strip.find_gaps(index=True)
Out[9]:
In [10]:
strip.thinnest()
Out[10]:
In [11]:
strip.thickest(n=5).plot(legend=legend, aspect=2)
In [12]:
strip
Out[12]:
In [13]:
strip = strip.prune(limit=1)
strip.plot(legend=legend, aspect=5)
In [14]:
strip = strip.anneal()
strip.plot(legend=legend, aspect=5)
In [15]:
strip.find_gaps() or "No gaps!"
Out[15]:
In [16]:
strip.plot(label='lithology', legend=legend, aspect=5)
In [17]:
fig, axs = plt.subplots(ncols=2, sharey=True)
axs[0] = strip.plot(ax=axs[0], legend=legend)
axs[1] = strip.plot_tops(axs[1], field='lithology', )
axs[1].axis('off')
plt.show()
In [18]:
strip[0]
Out[18]:
In [19]:
round(3.1415936, 2)
Out[19]:
In [20]:
import random
for iv in strip:
iv.data['por'] = round(random.random() / 4, 3)
strip[0]
Out[20]:
In [21]:
strip.plot(style='field', field='por', aspect=3)
In [22]:
strip[0].middle
Out[22]:
In [23]:
strip[0].primary['lithology']
Out[23]:
In [ ]:
In [24]:
s = strip[[1, 3]]
s.plot(legend=legend, aspect=1)
In [25]:
# Default behaviour: meet in middle.
fig, axs = plt.subplots(ncols=2, figsize=(4,2))
s.plot(legend=legend, aspect=1, ax=axs[0])
s.anneal(mode='middle').plot(legend=legend, ax=axs[1])
plt.show()
In [26]:
fig, axs = plt.subplots(ncols=2, figsize=(4,2))
s.plot(legend=legend, aspect=1, ax=axs[0])
s.anneal(mode='up').plot(legend=legend, ax=axs[1])
plt.show()
In [27]:
fig, axs = plt.subplots(ncols=2, figsize=(4,2))
s.plot(legend=legend, aspect=1, ax=axs[0])
s.anneal(mode='down').plot(legend=legend, ax=axs[1])
plt.show()
In [28]:
lappy = Striplog([Interval(**{'top': 50, 'base': 60, 'components':[Component({'lithology': 'dolomite'}),]}),
Interval(**{'top': 55, 'base': 75, 'components':[Component({'lithology': 'limestone'}),]}),
Interval(**{'top': 75, 'base': 80, 'components':[Component({'lithology': 'volcanic'}),]}),
Interval(**{'top': 78, 'base': 90, 'components':[Component({'lithology': 'anhydrite'}),]})
])
In [29]:
lappy.find_overlaps(index=True)
Out[29]:
In [30]:
overlaps = lappy.find_overlaps()
fig, (ax0, ax1) = plt.subplots(1, 2, sharey=True, figsize=(3,5))
# Use alpha to highlight overlaps.
ax0 = lappy.plot(legend, ax=ax0, alpha=0.75, lw=1)
ax1 = overlaps.plot(ax=ax1)
ax0.set_title('lappy')
ax1.set_title('overlaps')
ax0.set_ylim(100, 40)
plt.show()
The merge_overlaps() method operates in place and returns nothing.
In [31]:
lappy.merge_overlaps()
Now there are no overlaps!
In [32]:
lappy.find_overlaps()
In [33]:
lappy.plot(legend, aspect=3, alpha=0.75, ec='k', lw=0.5)
The merged intervals have mixed components:
In [34]:
lappy[1].components
Out[34]:
Note that the description is rather garbled.
In [35]:
lappy[1].description
Out[35]:
In [36]:
lappy = Striplog([Interval(**{'top': 50, 'base': 60, 'components':[Component({'lithology': 'dolomite'}),]}),
Interval(**{'top': 55, 'base': 75, 'components':[Component({'lithology': 'limestone'}),]}),
Interval(**{'top': 75, 'base': 80, 'components':[Component({'lithology': 'volcanic'}),]}),
Interval(**{'top': 78, 'base': 90, 'components':[Component({'lithology': 'anhydrite'}),]})
])
lappy.plot(legend, aspect=2, alpha=0.75)
You have to provide an Interval attribute (like 'top', 'base', or 'thickness') or a component attribute to merge on (if you use a component attribute, merge gets this attribute from the primary component to decide what to keep in each zone). The attribute must suport ordering — striplog keeps the thing with the larger value (or the thing which comes last in the ordering), e.g. the thickest, deepest, etc. If you want the reverse behaviour (keep the thinnest, shallowest, etc, then pass the reverse=True.
In [37]:
lappy.plot(legend, aspect=2, alpha=0.75)
In [38]:
lappy.merge('top').plot(legend, aspect=2, lw=1)
In [39]:
lappy.merge('top')[0].base.z
Out[39]:
In [40]:
lappy.merge('top', reverse=True)[0].base.z
Out[40]:
In [41]:
fig, axs = plt.subplots(ncols=2)
lappy.plot(legend, aspect=2, alpha=0.75, ax=axs[0])
lappy.merge('top').plot(legend, aspect=2, alpha=0.75, ax=axs[1])
plt.show()
In [42]:
fig, axs = plt.subplots(ncols=2)
lappy.plot(legend, aspect=2, alpha=0.75, ax=axs[0])
lappy.merge('top', reverse=True).plot(legend, aspect=2, alpha=0.75, ax=axs[1])
plt.show()
It is not currently possible to blend or combine units in the merge zones; you have to choose one.
This results in a new Striplog, contianing only the intervals requested.
In [43]:
strip.find('sandstone')
Out[43]:
In [44]:
strip.find('sandstone').unique
Out[44]:
In [45]:
strip.find('sandstone').cum
Out[45]:
In [46]:
strip.find('sandstone').plot(aspect=3)
Let's ask for the rock we just found by seaching.
In [47]:
rock = strip.find('sandstone')[1].components[0]
rock
Out[47]:
We can also search for a rock...
In [48]:
strip.find(rock).plot(legend, aspect=3)
In [49]:
rock in strip
Out[49]:
And we can ask what is at a particular depth.
In [50]:
strip.read_at(90).primary
Out[50]:
In [51]:
chrono = Striplog([Interval(**{'top': 0, 'base': 60, 'components':[Component({'age': 'Holocene'})]}),
Interval(**{'top': 60, 'base': 75, 'components':[Component({'age': 'Palaeogene'})]}),
Interval(**{'top': 75, 'base': 100, 'components':[Component({'age': 'Cretaceous'})]}),
])
In [52]:
time = Legend.default_timescale()
In [53]:
time[2]
Out[53]:
In [54]:
chrono.plot(time, aspect=3)
In [55]:
sands = strip.find('sandstone')
cretaceous = chrono.find('Palaeogene')
fig, (ax0, ax1, ax2) = plt.subplots(1, 3, sharey=True)
ax0 = sands.plot(legend, ax=ax0)
ax1 = chrono.plot(time, ax=ax1)
ax2 = sands.intersect(cretaceous).plot(legend, ax=ax2)
ax0.set_title('Sands')
ax1.set_title('Chrono')
ax2.set_title('Tert. sands')
ax0.set_ylim(100, 40)
plt.show()
In [56]:
for i in sands:
a = chrono.read_at(i.top.z)
i.data = {'age': a.primary['age']}
In [57]:
sands[1]
Out[57]:
Striplogs are just lists of Intervals, so you can use Python's functional programming patterns on them quite easily. For example, you can map functions and lambdas onto striplogs:
In [58]:
tops_in_feet = map(lambda i: i.top.z/0.3048, strip)
list(tops_in_feet)[:5] # First 5 only
Out[58]:
Don't forget the humble list comprehension though...
In [59]:
[i.thickness for i in strip][:5]
Out[59]:
Add all the thicknesses of intervals with a top depth > 100 m:
In [60]:
from functools import reduce
def sumr(a, b): return a + b
reduce(sumr, [i.thickness for i in strip if i.top.z > 100])
Out[60]:
To go even further, let's add a porosity array to each interval's primary component:
In [61]:
import random
for iv in strip:
iv.data['porosity'] = np.random.random(3)/4
strip[4]
Out[61]:
We can also write a function that returns True for some condition, and then filter intervals on that condition:
In [62]:
def porous(component):
return component.data['porosity'].mean() > 0.15
Striplog(list(filter(porous, strip)))
Out[62]:
It's a bit clunky now that filter returns an iterator. But it's also clunky because you can't pass arguments to the function you're giving filter — so you can't set the porosity to compare against when you call it, you have to edit the function itself.
To pass another argument to the filter function, you'll have to use a closure:
In [63]:
def min_porosity(x):
def compare(component):
return component.data['porosity'].mean() > x
return compare
Striplog(list(filter(min_porosity(0.15), strip)))
Out[63]:
In [64]:
liths = strip.to_log()
plt.plot(liths)
plt.show()
Pass a legend to get the ordering from the legend ('1' is given to the first component in the legend, '2' to the next, and so on).
In [65]:
liths = strip.to_log(legend=legend)
plt.plot(liths)
plt.show()
Recall that we added porosity to the components in this striplog:
In [66]:
strip[4]
Out[66]:
In [67]:
# I have broken this.
# por = strip.to_log(field='porosity', field_function=np.mean)
# plt.plot(por)
# plt.show()
We can also export any value corresponding to components from the legend, for example a 'width' log:
In [68]:
w, z, table = strip.to_log(legend=legend, legend_field='width', return_meta=True, step=0.1)
w
Out[68]:
...and we can make a composite plot in matplotlib:
In [69]:
width = 3
fig = plt.figure(figsize=(1,10))
ax = fig.add_axes([0, 0, 1, 1])
ax = strip.plot_axis(ax, legend, default_width=width+1)
plt.plot(w, z, color='white')
plt.fill_betweenx(z, w, width+1, edgecolor='white', facecolor='white', zorder=2)
ax.set_xlim([0, width+1])
ax.set_ylim([strip.stop.z, strip.start.z])
plt.show()
In [70]:
import lasio
l = lasio.read("P-129_out.LAS")
z, gr = l['DEPT'], l['GR']
In [71]:
strip.find('sandstone').plot(aspect=2) # There are actually 4 intervals here; 2 are touching
In [72]:
sand = strip.find('sandstone').to_flag(basis=z)
gr[sand].mean()
Out[72]:
Which is just a convenience; it does the same as:
In [73]:
sand = strip.find('sandstone').to_log(basis=z).astype(bool)
gr[sand].mean()
Out[73]:
Do we need to make this easier?
In [74]:
strip.plot(aspect=3)
In [75]:
strip[-1]
Out[75]:
In [76]:
s3 = strip.invert(copy=True)
In [77]:
s3.plot(aspect=3)
In [78]:
s3[0]
Out[78]:
In [79]:
l = """base,top,description
101,120,Till
100,101,Gypsum
50,100,Limestone Formation
28,50,Shale Formation
13,28,Granite Wash
0,13,Basement"""
In [80]:
log = Striplog.from_csv(text=l, lexicon=Lexicon.default())
In [81]:
log.read_at(30)
Out[81]:
In [82]:
log.order
Out[82]:
In [83]:
log.plot(aspect=3)
I recommend treating tops as intervals, not as point data.
In [84]:
tops_csv = """top,formation
100, Escanilla Fm.
200, Sobrarbe Fm.
350, San Vicente Fm.
500, Cretaceous
"""
In [85]:
tops = Striplog.from_csv(text=tops_csv)
In [86]:
print(tops)
In [87]:
tops.read_at(254.0)
Out[87]:
Some things really are point data. Sort of like a log, but irregular, more discrete. Here are some lab measurements...
In [88]:
data_csv = """depth, bodacity
1200, 6.4
1205, 7.3
1210, 8.2
1250, 9.2
1275, 4.3
1300, 2.2
"""
You must specify points=True otherwise Striplog will 'fill in' and create the bases for you, based on the next top.
In [89]:
points = Striplog.from_csv(text=data_csv, points=True)
In [90]:
print(points)
One day, when we have a use case, we can do something nice with this, like treat it as numerical data, and make a plot for it. We need an elegant way to get that number into a 'rock', like {'x': 6.4}, etc.
In [91]:
points.order
Out[91]:
We can read a log from an LAS file with lasio:
In [92]:
import lasio
Read a gamma-ray log.
In [93]:
l = lasio.read("P-129_out.LAS")
z, gr = l['DEPT'], l['GR']
In [94]:
z[-2000]
Out[94]:
Next we make a list of components to pass into the new striplog. The order must match the values you pass in the to_log() function:
In [95]:
comps = [Component({'lithology': 'sandstone'}),
Component({'lithology': 'greywacke'}),
Component({'lithology': 'shale'}), ]
Make a striplog from the GR curve, using the cutoffs given as cutoff = [10, 50]. These cutoffs define 3 lithologies, whichi is what we're passing in as comps. There must be enough components for the intervals you're defining.
If you don't provide components, you can provide legend instead; the components will be drawn from that. If you pass 'too many' components, they will be used in order and the 'extra' ones ignored.
You have to pass in the depth/elevation basis as well, because no assumptions are made about the log's extent.
In [96]:
s = Striplog.from_log(gr, cutoff=[10, 50], components=comps, basis=z)
s
Out[96]:
Now we can, say, remove the thin beds:
In [97]:
s.prune(limit=5)
s.anneal()
s
Out[97]:
And then read the log back into the intervals, 'reducing' with a function if we want:
In [98]:
s.extract(gr[2000:-2000], basis=z[2000:-2000], name='GR', function=np.mean)
In [99]:
s[20]
Out[99]:
Now close the loop by exporting these values as a new log and comparing to the original. Since we reduced with np.mean, we will get a blocked log...
In [100]:
g, gz, _ = s.to_log(field="GR", start=500, stop=1500, return_meta=True)
In [101]:
g2, gz2, _ = s.to_log(field="GR", return_meta=True)
In [102]:
plt.figure(figsize=(16,3))
plt.plot(z, gr, color='lightblue')
plt.plot(gz, g, lw=3, color='red')
plt.plot(gz2, g2, lw=1, color='black')
plt.show()
In [103]:
a = np.array([1,1,1,1,1,3,2,2,2,2,3,2,2,2,2,2,2,2,1,1,1,1,0,0,0,0,0,2,3,3,3,3,3,3,10,2,2,2,2,2,2,10,10,10,10,2,2,2,2,2])
In [104]:
z = np.linspace(100,200,50)
In [105]:
s = Striplog.from_log(a, legend=legend[:5], basis=z)
In [106]:
s[1]
Out[106]:
In [107]:
fig, (ax0, ax1) = plt.subplots(1, 2, sharey=True, figsize=(3,20))
# Use alpha to highlight overlaps.
ax0 = s.plot(ax=ax0)
ax1.plot(a, z, 'o-')
ax0.set_title('Striplog')
ax1.set_title('Log')
ax1.set_ylim(200, 100)
plt.show()
In [164]:
_ = hist(strip, legend=legend, rotation=-45, ha='left')
In [95]:
strip.crop((20, 100), copy=True)
Out[95]:
In [96]:
strip.crop((20, 100)) # in place
strip[0]
Out[96]:
Limitation — right now you cannot 'crop' to an extent larger than the current striplog. Maybe we should allow that, with padding...
In [97]:
strip.crop((0, 200))
# This should result in an error:
You should be able to natively save any format. If matplotlib complains, try replacing your usual import with
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
Or in Jupyter Notebook...
%matplotlib Agg
Then carry on as normal.
You need the figure object to save the striplog plot, so set return_fig to True:
In [98]:
fig = strip.plot(return_fig=True)
fig.savefig('test.png')
fig.savefig('test.pdf')
fig.savefig('test.svg')
To find out which backend you are using:
In [99]:
import matplotlib
matplotlib.get_backend()
Out[99]:
©2015 Agile Geoscience. Licensed CC-BY. striplog.py