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