Problem with Seismic refraction Inversion example with different Data set

[oubaaid]

Inversion for seismic refraction profile cannot execute due to "getPhiD".

I am following the "plot04_koenigsee_jupyer_notebook" example. I have a seismic refraction data files. These files were not internaly organized as the example i followed. So i re-organized a test file to try and use it to troubleshoot all of the steps in the process.

This is the example of file that i have followed : https://github.com/gimli-org/example-data/blob/master/traveltime/koenigsee.sgt I followed the example using this data file :

aconda Environnement, Python3, Pygimli

Operating system: Windows Python version: 3.8.12 pyGIMLi version: 1.2.3 Way of installation: Conda package

Steps to reproduce

import numpy as np
import pygimli as pg
import pygimli.physics.traveltime as tt
!pip install pandas
import pandas as pd
from pygimli.physics import TravelTimeManager

data = tt.load("test.sgt", load=True, verbose=True)
print(data)

fig, ax = pg.plt.subplots()
pg.physics.traveltime.drawFirstPicks(ax, data)

mgr = TravelTimeManager()

# Alternatively, one can plot a matrix plot of apparent velocities which is the
# more general function also making sense for crosshole data.
ax, cbar = mgr.showData(data)

mgr.invert(data, secNodes=3, paraMaxCellSize=5.0,
           zWeight=0.2, vTop=500, vBottom=5000,
           verbose=1)

ax, cbar = mgr.showResult(logScale=True)
mgr.drawRayPaths(ax=ax, color="w", lw=0.3, alpha=0.5)

Expected behavior

after print(data) : after pg.physics.traveltime.drawFirstPicks(ax, data) : after vest = mgr.invert(...) :

Actual behavior

It shows me missing hodochrone of the reversed offset shot then the velocity graph:

And the last 2 lines of code "vest = mgr.invert(....) : It shows the following error code

20/02/22` - 15:50:13 - pyGIMLi - ERROR - <class pygimli.physics.traveltime.TravelTimeManager.TravelTimeManager'>.checkError(C:\Users\dell\Anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py:87) DataContainer has no "err" values. Fallback to 3% 20/02/22 - 15:50:13 - pyGIMLi - INFO - Starting inversion. 20/02/22 - 15:50:13 - pyGIMLi - INFO - Create gradient starting model. 500: 5000 20/02/22 - 15:50:13 - pyGIMLi - INFO - Created startmodel from forward operator: [0.0003432 0.00056897 0.00177862 ... 0.00049306 0.00141945 0.0007828 ] fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x000001207FB3B2C0> Data transformation: <pgcore.pygimli.RTrans object at 0x000001207F9A4B80> Model transformation: <pgcore.pygimli.RTransLog object at 0x000001207F9A4720> min/max (data): 3.10/37.20 min/max (error): 3%/3% min/max (start model): 2.0e-04/0.0020

---

RuntimeError Traceback (most recent call last) ~\AppData\Local\Temp/ipykernel_2460/2457910414.py in 4 ax, cbar = mgr.showResult(logScale=True) 5 mgr.drawRayPaths(ax=ax, color="w", lw=0.3, alpha=0.5)""" ----> 6 vest = mgr.invert(data, secNodes=2, paraMaxCellSize=15.0, 7 maxIter=10, verbose=True) 8 np.testing.assert_array_less(mgr.inv.inv.chi2(), 1.1)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py in invert(self, data, useGradient, vTop, vBottom, secNodes, kwargs) 239 self.fop._useGradient = None 240 --> 241 slowness = super().invert(data, mesh, kwargs) 242 velocity = 1.0 / slowness 243 self.fw.model = velocity

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py in invert(self, data, mesh, zWeight, startModel, kwargs) 747 748 self.preRun(kwargs) --> 749 self.fw.run(dataVals, errorVals, kwargs) 750 self.postRun(kwargs) 751 return self.paraModel(self.fw.model)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\inversion.py in run(self, dataVals, errorVals, **kwargs) 462 self._preStep(0, self) 463 --> 464 self.inv.start() 465 self.maxIter = maxIterTmp 466

RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/inversion.cpp:95 double GIMLI::RInversion::getPhiD(const Vec&) const getPhiD == inf

"question" , "Seismic refraction" , "inversion" , "Pygimly"

[halbmy]

I suspect the error in the inversion could be due to zero error. I recommend doing an error estimate with tt.estimateError using both an absolute and a relative error. Try removing zero traveltime (zero offset) data.

[oubaaid]

Thank you for your fast answer @halbmy . I tried what you have suggested : TravelTimeManager().estimateError(data, errLevel=0.01, absError=None)

And it gave me this error which i honestly couldn't relate to.

---

TypeError Traceback (most recent call last) ~\AppData\Local\Temp/ipykernel_2460/1682671606.py in 1 #TravelTimeManager.estimateError(data, absoluteError=0.001, relativeError=0.001) ----> 2 TravelTimeManager().estimateError(data, errLevel=0.01, absError=None) 3 #tt.(data, errLevel=0.01, absError=None)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py in estimateError(self, data, errLevel, absError) 286 return absError + data errLevel 287 --> 288 return np.ones(len(data)) errLevel 289 290 def simulate(self, model, **kwargs):

TypeError: object of type 'DataContainerTT' has no len()

[carsten-forty2]

data need to be of type iterable not of type DataContainerTT.. so please provide the method your traveltime values like data['t']

[oubaaid]

I changed to data["t"], i think it works fine. And this is what i got, i don't understand what to do now in order to eliminate what's stopping me from moving to the inversion.

[carsten-forty2]

The array are a data error of 1e-3 ms + 1%(data) . You can take these estimated error values and put them into your data container before you put them into the inversion.

data['err'] = mgr.estimateError(data['t'], errLevel=0.01, absError=1e-3)

[oubaaid]

Thank you for your fast answers, Now i have err column in my data

But i am still encountering the same error 22/02/22 - 22:28:50 - pyGIMLi - INFO - Starting inversion. 22/02/22 - 22:28:50 - pyGIMLi - INFO - Create gradient starting model. 160: 7740 22/02/22 - 22:28:50 - pyGIMLi - INFO - Created startmodel from forward operator: [0.0002128 0.00031391 0.00372424 ... 0.00013228 0.00013672 0.00012945] fop: <pygimli.physics.traveltime.modelling.TravelTimeDijkstraModelling object at 0x000001207F64AAE0> Data transformation: <pgcore.pygimli.RTrans object at 0x000001207F64A090> Model transformation (cumulative): 0 <pgcore.pygimli.RTransLogLU object at 0x0000012069F01D60> min/max (data): 3.10/37.20 min/max (error): 1%/1.03% min/max (start model): 1.3e-04/0.0062

---

RuntimeError Traceback (most recent call last) ~\AppData\Local\Temp/ipykernel_2460/343464553.py in ----> 1 mgr.invert(data, secNodes=3, paraMaxCellSize=5.0, 2 zWeight=0.2, vTop=160, vBottom=7740, 3 verbose=1) 4 5 ax, cbar = mgr.showResult(logScale=True)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\physics\traveltime\TravelTimeManager.py in invert(self, data, useGradient, vTop, vBottom, secNodes, kwargs) 239 self.fop._useGradient = None 240 --> 241 slowness = super().invert(data, mesh, kwargs) 242 velocity = 1.0 / slowness 243 self.fw.model = velocity

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\methodManager.py in invert(self, data, mesh, zWeight, startModel, kwargs) 747 748 self.preRun(kwargs) --> 749 self.fw.run(dataVals, errorVals, kwargs) 750 self.postRun(kwargs) 751 return self.paraModel(self.fw.model)

~\Anaconda3\envs\pg\lib\site-packages\pygimli\frameworks\inversion.py in run(self, dataVals, errorVals, **kwargs) 462 self._preStep(0, self) 463 --> 464 self.inv.start() 465 self.maxIter = maxIterTmp 466

RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/inversion.cpp:95 double GIMLI::RInversion::getPhiD(const Vec&) const getPhiD == inf

[carsten-forty2]

would you mind providing the file test.sgt?

[oubaaid]

i made it work, only after i introduced this line of code data["t"] = data["t"]>0 i had this result

and when i changed the parameters of the invert function to these, i had this other weird result

I think there is a lot of tuning that needs to be done. this is the file "test.sgt" link: https://we.tl/t-UjMG3uvc7A

[carsten-forty2]

data["t"] = data["t"]>0 will set all data values to the value True (interpreted as 1.0) if they are greater than 0 and to False in the other case. I guess its not you really want.

If you want to remove data values you can try:

data.remove(data['t'] < 1e-6) # don't compare with zero this may lead to tolerance issues

Its always a good idea to check data range after manipulations. e.g.:

pg.info(f"data.min = {min(data['t'])}, data.max = {max(data['t'])})"

[halbmy]

Please not that in your file

...
#s   g   t
1    3    22.2

the traveltime is considered 22.2 seconds, but I guess you mean milliseconds?

As a consequence, your starting model not matching your data.

[oubaaid]

@carsten-forty2 i passed a long time trying your suggestion but did not work for me, so i skipped it. @halbmy Time in milliseconds yes that's right

[florian-wagner]

@oubaaid Have you tried again with the traveltime given in seconds?

[halbmy]

In view of missing response I am closing the issue. If the problem persists, please repeat this with the newest pygimli version 1.3 and reopen the issue.