In this tutorial, we will see how to use ResIPy API to time-lapse inversion of the problem in Figure 5.21 Note that more datasets are included here (as per Musgrave and Binley, 2011, doi:10.1016/j.jhydrol.2010.11.008
from resipy import R2
k = R2(typ='R2')
import matplotlib.pyplot as plt # used for more control of plots
Time-lapse datasets are stored in a sub-folder
k.createTimeLapseSurvey(r'C:\ResIPy\Notebook examples\Time-lapse inversion\data', ftype='ProtocolDC')
k.importElec(r'C:\ResIPy\Notebook examples\Time-lapse inversion\electrodes.csv') # read electrode geometry
k.createMesh(typ='trian', cl=0.25, fmd=5) # triangular mesh with element size near electrodes equal to 0.25m. fmd sets depth for fine mesh
k.zlim = [-5, 0]
k.showMesh(xlim=[0,63],zlim=[-5,0])
k.param['b_wgt'] = 0.02 # we are now informing R2 that the relative error is 0.02 for inversion
k.param['a_wgt'] = 0.001 # the offset error is 0.001 ohms
k.invert(parallel=True) # run time-lapse difference inversion in parallel once reference dataset is inverted
k.meshResults[0].df.columns # this will display the attributes that can be plotted
Plot log resistivity image for first (reference) dataset and show sensitivity pattern (see Figure 5.23b)
fig, axs = plt.subplots(2, 1, figsize=(18,6))
ax=axs[0]
k.showResults(ax=ax, index=0, attr='Resistivity(log10)', contour=False, vmin=1.5, vmax=2.4, color_map='jet', cropMaxDepth=True)
ax.set_title('09 Feb 2005 log10 resistivity')
# We can now plot the sensitivity map (as shown in Figure 5.23b)
ax=axs[1]
k.showResults(ax=ax, index=0, attr='Sensitivity(log10)', contour=True, vmin=-3, vmax=0, color_map='jet', cropMaxDepth=True)
ax.set_title('Sensitivity plot')
Now plot resistivity difference for the subsequent datasets relative to first dataset
headings = ['10 Mar 2005','07 Apr 2005','21 Apr 2005','03 May 2005','17 May 2005', '07 Jun 2005', '28 Jun 2005', '12 Jul 2005', '30 Aug 2005', '20 Sep 2005','30 Nov 2005', '18 Jan 2006', '09 Feb 2006']
fig, axs = plt.subplots(13, 1, figsize=(18,30))
for i in range(len(headings)):
ax=axs[i]
k.showResults(ax=axs[i], index=i+1, attr='difference(percent)', contour=False, vmin=-30, vmax=30, color_map='seismic', cropMaxDepth=True)
ax.set_title(headings[i])
if i < 12:
ax.set_xlabel('')