thomaslecocq / msnoise-tomo Goto Github PK

I am working on trying to get through the whole tomography process with this code, but having some slight confusion and not sure how to fix because documentation doesn't explain.

When I run msnoise p tomo ftan -s 0 it runs through all of my data and returns the following files.

GP.01101_GP.01102_MEAN.sac_FP.txt
GP.01101_GP.01102_MEAN.sac_TV.txt
GP.01101_GP.01102_MEAN.sac_amp.txt
GP.01101_GP.01102_MEAN.sac_disp.txt
GP.01101_GP.01102_MEAN.sac_ph.txt

These files are created for each station pair; the example above is for stations 01101 and 01102 in the GP network. I ask for the following periods

0.0333
0.0345
0.0357
0.0370
0.0385
0.0400
0.0417
0.0435
0.0455
0.0476
0.0500
0.0526
0.0556
0.0588
0.0625
0.0667
0.0714
0.0769
0.0833
0.0909
0.1000
0.1111
0.1250
0.1429
0.1667
0.2000

using the ftan_periods field in the tomo config. These are frequencies 5-30 Hz, spaced every 1 Hz. I have also set ftan_nfreq=41, which as I understand it means that the FTAN matrix will have 41 columns spanning ftan_fmin=5 to ftan_famx=30 in my particular case.

One station output then looks like

fmin=5, fmax=30, nfreqs=41
bmin=0.02, bmax=2
ampMin=0.05
vgMin=0.07 vgMax=0.4
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 5.22906
write amplitude
write phase
write FT
write TV
fmin=5, fmax=30, nfreqs=41
bmin=0.02, bmax=2
ampMin=0.05
vgMin=0.07 vgMax=0.4
disp: cont,  finit=5.22906 tginit=4.67366e-310
diag: Period - Vgroup
out: enable (matrix)
max: 5.22906
freq near 5.22906 = 5.22906
write amplitude
write phase
write FT
write TV
[ 0.033  0.034  0.036  0.037  0.038  0.04   0.042  0.043  0.045  0.048
  0.05   0.053  0.056  0.059  0.062  0.067  0.071  0.077  0.083  0.091
  0.1    0.111  0.125  0.143  0.167  0.2  ]
[        nan  0.08543644         nan         nan         nan         nan
         nan         nan         nan         nan         nan         nan
         nan         nan         nan         nan         nan         nan
         nan         nan         nan]

That makes sense. I understand that nans. That can be because of how close the stations are or the minimum amplitude parameter, ftan_ampmin=0.05.

What I do not understand is the resulting TOMO_DISP/ files. For example, for the stations above, I get the file TOMO_DISP/01/ZZ/GP_01101_GP_01102_MEAN.csv, which contains

,GP.01101_GP.01102_MEAN.sac
0.0,
0.1,0.192776595044
0.2,
0.30000000000000004,
0.4,
0.5,
0.6000000000000001,
0.7000000000000001,
0.8,
0.9,
1.0,
1.1,
1.2000000000000002,
1.3,
1.4000000000000001,
1.5,
1.6,
1.7000000000000002,
1.8,
1.9000000000000001,
2.0,

Now, I would have thought this contains two columns, Period and Velocity. But it does not appear to be the case. Can someone please explain this file to me? And I thought that runningmsnoise p tomo prepare_tomo combines the information in the TOMO_DISP folder into tomography input files on a per period basis in TOMO_FILES, but I don't think that can be because of the information contained in the DISP file above and the information we see in the TOMO_FILE below.

We had other issues (issue #4) with the small values of the periods (i.e. < 1.0). I wonder if there is another issues here. For example, let's look at the top three lines of the file for 0.2s period.

$ more TOMO_FILES/01/ZZ/TestGroupVel_0.2sGLISN.dat

GP.01101 GP.01102 0.2 0.192776595044 0 15.9273541759
GP.01101 GP.01103 0.2 0.26905235881 0 35.7949056343
GP.01101 GP.01104 0.2 0.108233877741 0 55.4774217922

We can see that stations GP.01101-GP.01102 at 0.2 s have a velocity of 0.192776595044 km/s, 0 is for the std. deviation and 15.9273541759 is the distance between stations.

Where did this information come from if not the TOMO_DISP/ files? And what is the point of the TOMO_DISP files? I confused as to the linearity of steps in the process. The DISP file above shows a velocity of this amount but for 0.1 s period (if I can interpret the first column as the period).

I am happy to dig through the source code to figure this out, but was hoping Tom or Aurelien can give a quicker answer. Thanks! (And we should add some documentation once this is figured out...which I can help with.)

I'm getting weird behavior with extremely small step sizes (this example has x and y steps of 0.0002 degrees). This is presumably not expected, and I suspect the culprit is the hard coded lat/long padding in prepare_tomo.py. I suggest padding with step size.

We think it would be good to add a t-max variable to the iftan GUI, similar to the normed FTAN box. This way the user can zoom in on the trace plot and see what the trace looks likes.

The setuptools.Feature class was removed from setuptools v46.0.0. For now, that case could be detected in setup.py at line 70 (adding the extra condition "or not hasattr(setuptools, 'Feature')" ) just like proposed by elygeo for obspy (issue #2578).

Hi Thomas,

I had created the tomo_sac files, populated the station table with correct coordinates and elevations. then I ran
msnoise p tomo
I was trying to plot the ftan curves and I got this error-

(base) C:\Users\akash\msnoise__testing>msnoise p tomo iftan

it opened a gui, but when I selected the TOMO_SAC files I got this error-

Running FTAN the first time: C:\Users\akash\msnoise__testing\TOMO_SAC\03\ZZ\AK_GREN_AK_GRES_MEAN.sac

fmin=0.1, fmax=1, nfreqs=40
bmin=0.0022, bmax=0.025
ampMin=0.05
vgMin=0.5 vgMax=5
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 0.1
write amplitude
write phase
write FT
write TV
C:\Users\akash\anaconda3\lib\site-packages\msnoise_tomo-0.1b0-py3.7-win-amd64.egg\msnoise_tomo\ftan_call.py:28: UserWarning: loadtxt: Empty input file: "write_amp.txt"
amp = np.loadtxt('write_amp.txt') # this the FTAN matrix
Exception in Tkinter callback
Traceback (most recent call last):
File "C:\Users\akash\anaconda3\lib\tkinter_init_.py", line 1705, in call
return self.func(*args)
File "C:\Users\akash\anaconda3\lib\site-packages\msnoise_tomo-0.1b0-py3.7-win-amd64.egg\msnoise_tomo\iftan.py", line 119, in process
ydata)
File "C:\Users\akash\anaconda3\lib\site-packages\msnoise_tomo-0.1b0-py3.7-win-amd64.egg\msnoise_tomo\ftan_call.py", line 32, in pickgroupdispcurv
index = np.unravel_index(np.argmax(amp), amp.shape)
File "<array_function internals>", line 6, in argmax
File "C:\Users\akash\anaconda3\lib\site-packages\numpy\core\fromnumeric.py", line 1186, in argmax
return _wrapfunc(a, 'argmax', axis=axis, out=out)
File "C:\Users\akash\anaconda3\lib\site-packages\numpy\core\fromnumeric.py", line 61, in _wrapfunc
return bound(*args, **kwds)
ValueError: attempt to get argmax of an empty sequence

Here is the output of conda list -
list.txt

I think there is an error in the t_obs calculation in ANSWT.py. This value is currently computed as

dist=(np.hypot(x1-x2, y1-y2))
t_obs = dist/Vg1

where dist is being computed from x1,x2,y1,y2, which are in units of degrees. I spoke with Aurelien today and he suggested using the rows in the G-matrix. This approach worked, but only after converting x1,x2,y1,y2 to a local cartesian grid. I tried many different things prior to this, and I could never get a good solution using units of degrees. Right now I am working on a grid that is 200m on either side.

I propose this solution.

# Compute the observed times in whatever units the grid is.
dist = GG.toarray().sum(axis=1) # distance in units of the grid
t_obs = dist / Vg1 # [??] arrival time data

But as I said, this has only worked for cartesian coordinates, and I think it has to do with src/mk_MatPaths.c. This is again just using the hypotenuse to compute distances in each grid cell and build the G-matrix.

Regardless of units, t_obs and t_calc should match in the way they are computed and right now they are not.

I also propose that we convert station coordinates to a local cartesian system inside of ANSWT.py and use a local grid. I can work on this. We can write the output back to actual coordinates, but I can't get the code to work right now without converting units...I even tried tracking down which units are causing the problem because Aurelien thought G is normalized and distance units would not matter. But I cannot find the problem yet.

Thoughts?

The descriptions of FTAN parameters on the tomoconfig page are not correct. This should be fixed.

Adding this issue because I am about to submit a commit to fix this. Now we know what this commit is for.

It seems like the minimum wavelength criteria is only actually function in iftan.py (e.g. line 166 idxx = np.where( float( _minWL.get() ) < dist/(disper*per) )[0]), but not in ftan.py. Therefore, the cutoff based on minimum wavelength is not being propagated into the data for ANSWT.py. I do not see any sort of wavelength filtering in prepare_tomo.py either. Should this be added to ftan.py?

Should add a warning or error when running prepare_ccf if the station coordinates are set to the default values 0. Otherwise, the distance between the stations is set to zero in the SAC files and the FTAN output empty files because the velocities are infinite.

There is a going to be problems with folks trying to use the tomo plugin for high frequency data. For example, when running msnoise p tomo prepare_tomo, the following file is written:
TestGroupVel_0.1sGLISN.dat.

(First, the word GLISN is in there for some reason. I guess this a hold over from Aurelien's original version. I don't see a need for that. The files are likely already in a project folder.)

Second, and more important, the files are named based on period with only one decimal. I am trying to do a tomography from 10 to 30 Hz, with the following periods:
0.033,0.034,0.036,0.037,0.038,0.040,0.042,0.043,0.045,0.048,0.050,0.053,0.056,0.059,0.062,0.067,0.071,0.077,0.083,0.091,0.100,0.111,0.125,0.143,0.167,0.200

When I run the prepare_tomo command, I get three files:

TestGroupVel_0.0sGLISN.dat
TestGroupVel_0.1sGLISN.dat  
TestGroupVel_0.2sGLISN.dat

Each one related to whichever was the last file written in this period band.

0.091
0.167
0.200

Because of the 1 decimal naming convention, the files just get overwritten. It would be very useful to have at least 3 decimals in the file names if we are going to use the period in the file name.

We are getting some strange behavior (or we think it is strange) regarding alpha_2. When alpha_2 is zero, we get a smoother model than when alpha_2 is 10^-6. See attached image. We are setting the _1 variables to 0 so that as many data as possible are kept. Then we vary alpha_2 to see how the final model is smoothed. Any ideas about what is going on?

when I do msnoise p tomo install this gave me the following error:

Usage: msnoise p [OPTIONS] COMMAND [ARGS]...
Try "msnoise p --help" for help.

Error: No such command "tomo".

while this command msnoise config set plugins=msnoise_tomo gave me this output: Successfully updated parameter plugins = msnoise_tomo

so I wanted to install https://github.com/ThomasLecocq/msnoise-tomo but it also doesn't install. my python version Python 3.8.8
it gives UnsatisfiableError: The following specifications were found to be incompatible with each other:
Package expat conflicts for ....

Please help me on this issue.
Regards.

I am posting about retrieving plots of the model resolution in msnoise-tomo. It appears this is done within the ANSWT.py script (Lines 464-532). However, these lines cannot be reached because they are under a if False: statement (Line 464) that does not appear to point to a variable being True or False. It also has a module/def that's not imported in (measure.find_contours, line 483) and a variable (pasgrille, line 523) that is not defined. Has this been corrected in any recent versions? or do you have any insight on how this should be corrected?

Also, it appears the lambda parameter input for the matrix H in the second iteration is lambda1, shouldn't it be lambda2 (line 306). This also leads me to ask, is there any built-in way in the MsNoise-tomo package to produce L-curves for the range of smoothing/damping parameter values?

Hi Thomas,

I am trying to install msnoise_tomo, to begin with -
whenever I type
python setup.py install in the main directory, I get following error -

Hi Thomas!

Thanks for your work. I'm testing tomo on ambient noise data. Can you share the workflow for this? The documentation seems to be missing it.

Thanks a lot, in advance!

Linus

I noticed that the tomo plugin fails with the newest numpy (1.19). This has already been noticed in Obspy (issue 2643). I plan to fix this after I submit a different pull request to address some FTAN.py issues. Will try to get to it later today.

Example of the current error thrown is below. Current work around is to use numpy 1.18.1.

python setup.py install

Traceback (most recent call last):
  File "setup.py", line 19, in <module>
    from numpy.distutils.core import DistutilsSetupError, setup
ImportError: cannot import name 'DistutilsSetupError' from 'numpy.distutils.core' (/Users/dmikesell/anaconda3/envs/tomo/lib/python3.7/site-packages/numpy/distutils/core.py)

In 'msnoise p tomo plot' the dispersion curves are concatenated then plotted. It would be nice to add NaNs in between each dispersion curve concatenation so that the endpoint of one curve is not linked to the first point of the next curve. That would make the plot much nicer.

Hi Thomas!
Recently I was ready to start learning MSNoise-TOMO, but there was very little in Documentation. Can you provide more detailed introductory documentation? I would be grateful if I could get your help.
Thanks a lot, in advance!

Hi @ThomasLecocq

I am planning to help my students and make some modifications to msnoise-tomo source code...and fix our own issues.

The previous changes (issue #4) were simple and I didn't need to test things. I am now trying to do an installation of msnoise-tomo from scratch in my forked version of the repo.
I am running python 3.8, obspy 1.2.1, and setuptools 47.3.1
I cannot get around this error, and I wonder if you have dealt with this already.
I am running

$ python setup.py install 

and getting the following output

/Users/dmikesell/GIT/msnoise-tomo
/Users/dmikesell/GIT/msnoise-tomo
Traceback (most recent call last):
  File "setup.py", line 229, in
    setupPackage()
  File "setup.py", line 200, in setupPackage
    features=add_features(),
  File "setup.py", line 73, in add_features
    class ExternalLibFeature(setuptools.Feature):
AttributeError: module 'setuptools' has no attribute 'Feature'

I tried running $ python setup.py clean --all and then rerunning the install. Same problem. Have you had this problem with setupPackage() before? I have tried a number of install things so far and none have worked.

Any thoughts? I found this on the setuptools package page: https://setuptools.readthedocs.io/en/latest/history.html#v46-0-0.

We cannot just remove features=add_features(), on line 200 in setup.py. I tried commenting this line out, but still got an error when trying to do the gcc compilation. Wondering if I need to go back to python 3.6 perhaps? add_features() is your location function that sets up the external libraries.

I am creating an issue for the diagram type plotting in both the iftan.py and ftan.py call. I will submit a commit to fix this issue, but wanted an issue to reference in the commit.

I have been working to modify iftan.py and ftan.py to work with any of the combination (PV, FV, PT, FT). I have everything working now, but the velocity seeding is broken. I am comparing my updated code to the version in the WHL file. Below is the iftan output for a SAC file.

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.5
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
write amplitude
write phase
write FT
write TV
fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.5
disp: cont,  finit=13.2006 vginit=0.185903
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
freq near 13.2006 = 13.2007
write amplitude
write phase
write FT
write TV

This is the output from running ftan() twice. ftan() is the external C routine. The first time the ftan() function is ran, you can see that the the "disp" is set to "no". The second time, the ftan() function is ran with the seed "finit" and vginit". The ftan() function is actually invoke in ftan_call.py. See below

ftan(filename, fmin, fmax, vgmin, vgmax, bmin, bmax,
         diagramtype, nfreq, ampmin, dist, disp="cont", tinit=pinit, vginit=vinit)

Now, in the modified version of my code, which is compiled on my mac, I cannot get the seeding to work properly. Here is the output for the exact same SAC file and db parameters. You can see I have added a couple of outputs, and a print of the vinit value before the second call to ftan().

Running FTAN the first time: /Users/dmikesell/Desktop/test_for_john/msnoise_geopark/TOMO_SAC/01/ZZ/GP_01101_GP_01205_MEAN.sac

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.35
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
write amplitude
write phase
write FT
write TV
Seed velocity: 0.185903 [km/s]
Seed period: 0.075754 [s]
vinit: 0.185903

Running FTAN a second time with the SEED

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.35
disp: cont,  finit=13.2006 tginit=2.24583e-314
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
freq near 13.2006 = 13.2007
write amplitude
write phase
write FT
write TV

I am using the same diagramtype (i.e. PV), but the C code always spits out tginit=0, rather than vinit, which is what I am passing to it. My call is identical (I think) to that in the WHL file.

ftan(filename, fmin, fmax, vgmin, vgmax, bmin, bmax,
     diagramtype, nfreq, ampmin, dist, disp="cont", tinit=pinit, vginit=vinit)

However, the output is not the same. How can I tell which version of the code was used to create the WHL file? Or can someone tell me why the C code keeps looking for tginit (the time seed) instead of the velocity seed (vginit) when I use the PV diagram type?

Export format has to be mseed in MSNoise otherwise the hardcoded values in export single sided fail