Induced Polarisation Ratio of Two Separate Inversions/Datasets

[JMcKercher]

Problem description

Hello, I have been loving this software and starting to get a good understanding of it. My knowledge of coding is still quite limited and I am unsure if this has already been covered in an example or a tutorial, but I would like to know if it is possible to create an IP ratio (graph and/or inversion plot) from two different inversions or datasets.

Problem: I have collected two sets of data each at different frequencies from the same field transect. I would like to test or evaluate the IP difference between the two frequencies to ultimately determine which frequency is better at detecting a signal from the field. Please let me know if you require further information.

Your environment

            OS : Windows
        CPU(s) : 8
       Machine : AMD64
  Architecture : 64bit
           RAM : 7.7 GiB
   Environment : Jupyter

Python 3.9.18 | packaged by conda-forge | (main, Dec 23 2023, 16:29:04) [MSCv.1929 64 bit (AMD64)]
pygimli : 1.4.6
pgcore : 1.4.0
numpy : 1.21.6
matplotlib : 3.7.2
scipy : 1.11.2
IPython : 8.18.1
pyvista : 0.43.2

Steps to reproduce

Below is the basic inversion of two data sets showing the IP results. I am unsure how to go about what I would like to achieve, but I believe it to be like ratio.

import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
import pygimli.meshtools as mt
from pygimli.physics import ert
%matplotlib inline

lowfdat = ert.load('C:/PATHTODATA/Rest point 03-04 cross 0.52.tx0',
                  load=True, 
                  verbose=True);
modfdat = ert.load('C:/PATHTODATA/Rest point 03-04 cross 4.16.tx0',
                  load=True, 
                  verbose=True);
print("Low frequency data=", lowfdat)
print("Moderate frequency data =", modfdat)
ert.show(lowfdat, label="Low frequency")
ax, cb = ert.show(modfdat, label="Moderate frequency")

#Low Frequency Data
lowfdat.remove(lowfdat["rhoa"] < 0.01)
lowfdat.remove(lowfdat["rhoa"] > 2500)
ert.show(lowfdat, label="Low frequency")
#Moderate Frequency Data
modfdat.remove(modfdat["rhoa"] < 0.01)
modfdat.remove(modfdat["rhoa"] > 2500)
ax, cb = ert.show(modfdat, label="Moderate frequency");

ipmgr = ert.ERTIPManager(sr=False);

lowfdat['k'] = ert.createGeometricFactors(lowfdat, numerical=True)
k1 = ert.createGeometricFactors(lowfdat, numerical=True)
ert.show(lowfdat, 
         vals=k1/lowfdat['k'], 
         label='Topography effect - Low Hz',
         cMap="bwr", 
         logScale=True)
modfdat['k'] = ert.createGeometricFactors(modfdat, numerical=True)
k2 = ert.createGeometricFactors(modfdat, numerical=True)
ax, cb = ert.show(modfdat,
                  vals=k2/modfdat['k'],
                  label='Topography effect - Mod Hz',
                  cMap="bwr", 
                  logScale=True);

_ = lowfdat.show("ip", 
                 cMap="inferno", 
                 label="Low Hz")
ax, cb = modfdat.show("ip", 
                      cMap="BrBG", 
                      label="Mod Hz");

# Low Frequency Data
lowfdat.remove(lowfdat["ip"] < 0.01)
lowfdat.remove(lowfdat["ip"] > 200)
lowfdat.show("ip",
             cMap="inferno_r",
             label="Low Hz")
# Moderate frequency data
modfdat.remove(modfdat["ip"] < 0.01)
modfdat.remove(modfdat["ip"] > 200)
ax, cb = modfdat.show("ip",
                      cMap="BrBG",
                      label="Mod Hz");

# Low frequency data
lowfdat['err'] = ert.estimateError(lowfdat, absoluteUError=5e-5, relativeError=0.03); #Estimates the errors
ert.show(lowfdat, lowfdat['err']*100, label="Low Hz Error [%]")
# Moderate frequency data
modfdat['err'] = ert.estimateError(modfdat, absoluteUError=5e-5, relativeError=0.03); #Estimates the errors
ert.show(modfdat, modfdat['err']*100, label="Moderate Hz Error [%]");

# Low frequency
lowipmgr = ert.ERTIPManager(lowfdat)
lowip = lowipmgr.invert(lowfdat,
                     verbose=True,
                     paraMaxCellSize=10,
                     paraDepth=30,
                     quality=33.6)
# Moderate frequency
modipmgr = ert.ERTIPManager(modfdat)
modip = modipmgr.invert(modfdat,
                     verbose=True,
                     paraMaxCellSize=10,
                     paraDepth=30,
                     quality=33.6);

# Somehow model the ratio  of the two seperate datasets?
ipratio = pg.log(lowipmgr.showIPModel / modipmgr.showIPModel)

Expected behavior

Create a inversion graphical output displaying the ratio between the two datasets IP models

Actual behavior

None so far, is this possible?

Data to reproduce

Rest point 03-04 cross 0.52.txt
Rest point 03-04 cross 4.16.txt

[halbmy]

Sorry for the very late response. I've had a look into your scripts. There might be some further modifications necessary (filtering ranges, error estimation) but as to the technical point, you can simple access the IP model by mgr.modelIP.

However, I suggest NOT using the ratio of the IP values but the difference. If you consider resistivity as complex number of amplitude and phase A*exp(w P), the ratio of two values exhibits the phase difference, i.e.
A2 exp(w P2) / (A1 exp(w P1)) = A2/A1 exp(w(P2-P1))

Done by

lowipmgr.showResult(lowipmgr.modelIP - modipmgr.modelIP,
                     cMap="bwr", logScale=False)

I ended up in this image

[JMcKercher]

Thank you for getting back to me. I am able to return the same response using the provided code.

[halbmy]

I consider this issue solved. If not, reopen and describe the problem in more detail.