Dear users,
I want find a way to insert into the CC header's files the pair stations distance for each CC.
In the file 'database_tools.py' I found in 'MISCS' this code:
############ MISCS ############
def azimuth(coordinates, x0, y0, x1, y1):
if coordinates == "DEG":
dist, azim, bazim = gps2DistAzimuth(y0, x0, y1, x1)
# print dist, azim, bazi
return azim
elif coordinates == 'UTM':
azim = 90. - np.arctan2((y1 - y0), (x1 - x0)) * 180. / np.pi
# print azim
return azim
else:
print "woooooow, please consider having a single coordinate system for all stations"
return 0
So my idea was insert this code inside the 'correlation' section to insert these information in the header's files like it happens with other parameters below:
if sac_format == "doublets":
tr.SetHvalue('A', 120)
else:
tr.SetHvalue('B', -maxlag)
tr.SetHvalue('DEPMIN', np.min(corr))
tr.SetHvalue('DEPMAX', np.max(corr))
tr.SetHvalue('DEPMEN', np.mean(corr))
tr.SetHvalue('SCALE', 1)
tr.SetHvalue('NPTS', len(corr))
tr.WriteSacBinary(filename)
del st, tr
return
I think also that insert this type of implementation in the code could be very helpful. I hope that someone can help me.
cheers,
Carmelo
-----------
Mr Carmelo Sammarco MSc BSc FGS
PhD student in Geology and Petroleum Geology
University of Aberdeen
The University of Aberdeen is a charity registered in Scotland, No SC013683.
Tha Oilthigh Obar Dheathain na charthannas clàraichte ann an Alba, Àir. SC013683.
MSNoise Workshop at AGU 2014
Dear All !
There are things coming every year, such as Thanksgiving, or Christmas…
And, before Christmas, there is the annual pilgrimage for some 20.000
geoscientists gathering in San Francisco for AGU. Just like last year,
we have booked a room at the Marriott Marquis to organise a MSNoise
workshop ! Last year an half-day was too short for going into detail, so
this year, we are pleased to announce a *1 day MSNoise workshop* on
*Sunday 14 December*!
The aim of this workshop is to demonstrate the powers (and limitations)
of the dv/v calculations using ambient seismic noise cross-correlation.
The MSNoise package (written in Python) will be introduced and then used
by participants on provided demo data. A general schedule will be:
08.30 – 08.45: Welcome
08.45 – 10.00: Introduction to ambient noise cross-correlation & dv/v –
Florent Brenguier (ISTerre)
10.00 – 12.00: MSNoise: general presentation - Thomas Lecocq (Royal
Obs. of Belgium)
12.00 – 12.45: Lunch Break (not included)
13.00 – 18.00: MSNoise practical session – Thomas & Corentin Caudron
(EOSingapore)
Because of high organisational costs, we can’t afford to organise this
workshop for free. We have fixed the fee to 250 USD (200 €) per person,
with a discount for students at 185 USD (150 €). Payments will be done
online via Credit Card (payment information will be sent by email).
To ensure the quality of the experience for everyone, we will not accept
more than 30 participants. Each participant is required to have a
sufficiently recent laptop with software requirements installed
beforehand (instructions will follow). Participants are not expected to
be Python experts, but to have at least a little experience in
programming. To register, use the form located here
<http://www.msnoise.org/msnoise-at-agu-2014-registration/>.
*Information on MSNoise:* MSNoise is the first complete software package
for computing and monitoring relative velocity variations using ambient
seismic noise. MSNoise is a fully-integrated solution that automatically
scans data archives and determines which jobs need to be done whenever
the scheduled task is executed (http://www.msnoise.org).
*Registration:* http://www.msnoise.org/msnoise-at-agu-2014-registration/
*Worth reading:* Lecocq, Thomas, Corentin Caudron, et Florent Brenguier
(2014), MSNoise, a Python Package for Monitoring Seismic Velocity
Changes Using Ambient Seismic Noise, /Seismological Research Letters/,
/85/(3), 715‑726, doi:10.1785/0220130073
<http://srl.geoscienceworld.org/content/85/3/715.full>.
http://srl.geoscienceworld.org/content/85/3/715.full*
*
*
Thomas Lecocq
*
--
Dr. Thomas Lecocq
Seismology - Gravimetry
Royal Observatory of Belgium
*
* * * * *
* * * *
---------
http://www.seismology.behttp://twitter.com/#!/Seismologie_be
Hello,
I was working with MSNoise without problems until I decided to change the
*mwcs_wlen* and *mwcs_step*. I ran *s05compute_mwcs.py* and
*s06compute_dtt.py* with the new parameters and then I ran
*s07plot_dtt.py* giving
me the next error message:
*oscar@terminus:~/Doctorado/TI_2/sweet_noise$ py s07plot_dtt.py *
*/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/setuptools-3.6-py2.7.egg/pkg_resources.py:1045:
UserWarning: /home/oscar/.python-eggs is writable by group/others and
vulnerable to attack when used with get_resource_filename. Consider a more
secure location (set with .set_extraction_path or the PYTHON_EGG_CACHE
environment variable).*
*Traceback (most recent call last):*
* File "s07plot_dtt.py", line 96, in <module>*
* data = detrend(py1_wmean)*
* File
"/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/statsmodels/tsa/tsatools.py",
line 230, in detrend*
* beta = np.linalg.lstsq(trends, x)[0]*
* File
"/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py",
line 1837, in lstsq*
* nlvl = max( 0, int( math.log( float(min(m, n))/2. ) ) + 1 )*
*ValueError: math domain error*
*oscar@terminus:~/Doctorado/TI_2/sweet_noise$ vi s06compute_dtt.py *
*oscar@terminus:~/Doctorado/TI_2/sweet_noise$ py s07plot_dtt.py *
*/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/setuptools-3.6-py2.7.egg/pkg_resources.py:1045:
UserWarning: /home/oscar/.python-eggs is writable by group/others and
vulnerable to attack when used with get_resource_filename. Consider a more
secure location (set with .set_extraction_path or the PYTHON_EGG_CACHE
environment variable).*
*Traceback (most recent call last):*
* File "s07plot_dtt.py", line 96, in <module>*
* data = detrend(py1_wmean)*
* File
"/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/statsmodels/tsa/tsatools.py",
line 230, in detrend*
* beta = np.linalg.lstsq(trends, x)[0]*
* File
"/home/oscar/bin/python/anaconda/lib/python2.7/site-packages/numpy/linalg/linalg.py",
line 1837, in lstsq*
* nlvl = max( 0, int( math.log( float(min(m, n))/2. ) ) + 1 )*
*ValueError: math domain error*
Does anyone can help me? Thank you very much in advance!
Oscar
Hi Thomas,
Thanks again for your information. I have already added 1bit
normalization by myself and am checking if it maked any differences or
not (or, whether I made any mistakes or not). So far, changes between on
& off 1 bit normalizations are subtle except for some combinations of
stations. I will let you know if I find anything on this.
Also, thank you very much for your application of additional signal
processing to avoid frequency leakage. That helps me to understand your
code. As for difference of amplitude, I also find slight difference
between your MSNoise and Florent's system on cirteria around RMS. I will
let you know if it make systematic difference of amplitude of not.
Sincerely,
Run
(11/11/14 1:24 PM), Thomas Lecocq wrote:
>
> Hello Run,
>
>
> Le 4/11/2014 11:01, Kohtaro R. Araragi a e'crit :
>> Hi Thomas,
>>
>> I have a quick question for you. Do you apply 1 bit normalization
>> (Bensen et al., 2007, GJI) at MSNoise?
> No, only 3*RMS clipping (windsorization) for each window.
>
>> I have been looking for your documentation and MSNoise code, but I
>> could not find or I may not understand your code well.
>>
>> At least, I can find that you apply similar processes around
>> "Pre-Whitening Traces"(below). If you actually do not apply 1 bit
>> normalization, should I apply by myself or do you intend to add it to
>> regular MSNoise package?
>>
>> If I just misunderstand your code, I'm sorry and appreciate if you
>> clarify this part.
>>
>> FYI, as long as I compare output from MSNoise and Florent's system,
>> they match very well if we analyze data that do not have many
>> earthquakes (attached picture below. Blue line from MSNoise, Red line
>> from Florent's system).
>
> Are you using MSNoise:master or MSNoise:releases:1.2.5 form github ?
> The master contains a bug correction in the maths of the whiten
> function, which changes the output slightly (but is correct !!). BTW,
> I also observed that my MSNoise shows more amplitude (but hopefully
> same phase) as Florent's matlab edition. This is because I care a
> little more about the signal processing (cos tapering windows, etc...,
> to avoid frequency leakage, ....etc)... :-) :-)
>
> Best
>
> Thomas
>
>>
>> Thanks in advance,
>> Run
>> --
>> logging.info("Pre-Whitening Traces")
>> whitened_slices = np.zeros((len(stations), len(get_filters(db, all=False)), tranches, int(Nfft)), dtype=np.complex)
>> for istation, station in enumerate(stations):
>> for itranche in range(tranches):
>> tmp = tramef_Z[istation, itranche * int(min30):(itranche + 1) * int(min30)]
>> rmsmat = np.std(np.abs(tmp))
>> indexes = np.where(
>> np.abs(tmp) > (windsorizing * rmsmat))[0]
>> tmp[indexes] = (tmp[indexes] / np.abs(
>> tmp[indexes])) * windsorizing * rmsmat
>> tmp *= cp
>> for ifilter, filter in enumerate(get_filters(db, all=False)):
>> whitened_slices[istation, ifilter, itranche,:] = whiten(tmp, Nfft, dt, float(filter.low), float(filter.high), plot=False)
>> del tmp
>> del tramef_Z
>>
>>
>> ____________________________________________________
>> K. Run Araragi
>> PhD studnet/M.Sc.
>> ISTerre, France &
>> Earthquake Research Institute
>> 1-1-1 Yayoi, Bunkyo-ku
>> Tokyo, Japan 113-0032
>>
>> Skype name: kalessinlord
>> Tel: 81-3-5841-5697
>> Cel: 81-90-1166-3091
>> 33-970-461-805
>> mail:kohtaro.araragi@ujf-grenoble.fr
>> ____________________________________________________
>
--
____________________________________________________
K. Run Araragi
PhD studnet/M.Sc.
ISTerre, France &
Earthquake Research Institute
1-1-1 Yayoi, Bunkyo-ku
Tokyo, Japan 113-0032
Skype name: kalessinlord
Tel: 81-3-5841-5697
Cel: 81-90-1166-3091
33-970-461-805
mail:kohtaro.araragi@ujf-grenoble.fr
____________________________________________________
Hi Thomas,
Thank you very much, Thomas. I should have downloaded from link at
https://github.com/ROBelgium/MSNoise/releases NOT from
https://github.com/ROBelgium/MSNoise. I will check the release version
and let you know if you find anything around "keep_all" options or other
things.
Sincerely,
Run
(11/11/14 1:20 PM), Thomas Lecocq wrote:
> Hello Run,
>
> Please send your emails to the mailing list, your remarks could be
> useful for many users !
>
>
>
> Le 22/10/2014 11:47, Kohtaro R. Araragi a e'crit :
>> Hi Thomas,
>>
>> Thank you very much for your help the other day. I had to change a
>> little bit after that, but your e-mail with other MSNoise users in
>> the past can help me to resolve all prolbem and now I can use all
>> function until s06compute_dtt.py.
>>
>> Today, I would like to you let you know about what I find in MSNoise
>> function.
>>
>> In "s001configurator.py", we can set "keep_all" as Y if we want to
>> keep each corr_duration CCF (or all 30 seconds cross-corr). However,
>> when I set Y in "keep_all", "s03compute_cc.py" did not work. When I
>> check the code and error message, it seems that current version of
>> s03compute_cc.py use "thisdate" for processing(see bold part in the
>> source below) of "keep_all" before value is assigned. When I set
>> "keep_all" as "N", everything is working and I can compute C.C. I may
>> make some mistakes, and I would be happy if you let me know where I
>> can change.
> You are using the "master" version of MSNoise, which contains this bug
> ! use "releases" on github to be sure to use the bug-free versions !
>
> Thomas
>
>>
>> You may have updated this part and if so, I appreciate if you let me
>> know.
>>
>> Sincerely,
>> Run
>>
>> ----------------------------------------------------------------
>> *if keep_all:**
>> **add_corr(db, station1.replace('.', '_'), station2.replace(**
>> **'.', '_'), filter.ref, thisdate, thistime, min30 / fe, "ZZ", corr,
>> fe)**
>> * if keep_days:
>> if not np.any(np.isnan(corr)) and not np.any(np.isinf(corr)):
>> daycorr += corr
>> ndaycorr += 1
>>
>> if keep_days:
>> thisdate = time.strftime(
>> "%Y-%m-%d", time.gmtime(basetime))
>> thistime = time.strftime(
>> "%H_%M", time.gmtime(basetime))
>> add_corr(
>> db, station1.replace(
>> '.', '_'), station2.replace('.', '_'), filter.ref,
>> thisdate, thistime, min30 / fe, 'ZZ', daycorr, fe, day=True,
>> ncorr=ndaycorr)
>> update_job(db, goal_day, orig_pair, 'CC', 'D')
>> logging.info("Job Finished. It took %.2f seconds" % (time.time() - jt))
>> ----------------------------------------------------------------
>
--
____________________________________________________
K. Run Araragi
PhD studnet/M.Sc.
ISTerre, France &
Earthquake Research Institute
1-1-1 Yayoi, Bunkyo-ku
Tokyo, Japan 113-0032
Skype name: kalessinlord
Tel: 81-3-5841-5697
Cel: 81-90-1166-3091
33-970-461-805
mail:kohtaro.araragi@ujf-grenoble.fr
____________________________________________________
**
Dear Thomas and MSnoise users,
I just completed my first installation and processing on the test case
from AGU 2013 and I had a few observation and comments to improve the
existing documentation for future users.
In the Configurator:
*
I made the mistake to set keep_days as “No”, probably because I
copied the screenshot from the doc. This was, of course, a problem
when I computed the CC. I think it would be worth either to change
the screenshot, or to mention in the doc that this setting could be
problematic.
*
The ref begin and end are set by default in relative, potentially
causing errors if the REF folder is not created by s04stack.
*
I had no filter values by default, I think it would be good to have
suggested values for those who are quickly trying MSnoise.
*
I think it would also be nice to have more advice on the use of
sqlite and mysql. Maybe recommend a sqlite manager, the way it is
done with EasyPHP for mysql.
*
It was already mentioned in the mailing list but it would also be
nice to have suggestions on the best way to reset jobs from “D” to
“T”. It probably won’t be a problem when the toolbox will be ready.
*
I didn’t install wx, so to plot the dt/t I had to disable the line 3
“mpl.use('WxAgg')”. Also the red line didn’t show on the plot, but I
think this is a known problem.
I hope this will be useful for future users.
Cheers,
Raphael De Plaen
PhD Student
University of Luxembourg
**
Using MSNoise 1.24 there are two python scripts that plot; they are
07plot_dtt.py and s07plot_dtt.py. I get errors for either.
Using s07plot_dtt.py yields attachment s07plot_dtt.png. There are two
problems:
1) I entered mov_stack values 1,3,5, and 10 days into the configurator but
the graphical output includes only 1, 3, and 5 days moving window graphs,
and the only graph with any data plotted onto it is the 5 days moving
window.
2) the graph x axes automatically include a much larger time period than
covered by the data plotted. This includes dates back to 2013 and forward
into 2016, while all the data I have is from 2014. I am unsure if this
problem is related to some of my settings in the configurator. My settings
that I can see might be related are: enddate=2100-01-01, ref_begin=-100,
ref_end=0, startdate=2014-05-04.
Using 07plot_dtt.py yields attachment image 07.plot_dtt.png the following:
>>> runfile('/home/devora/MSNoise-master/07.plot_dtt.py',
wdir=r'/home/devora/MSNoise-master')
UMD has deleted: database_tools, msnoise_table_def
loading 1 days
loading 3 days
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
File
"/usr/lib/python2.7/site-packages/spyderlib/widgets/externalshell/sitecustomize.py",
line 540, in runfile
execfile(filename, namespace)
File "/home/devora/MSNoise-master/07.plot_dtt.py", line 128, in <module>
plt.fill_between(ALL.index,ALL[dttname]-ALL[errname],ALL[dttname]+ALL[errname],lw=1,color='red',zorder=-1,alpha=0.3)
File "/usr/lib64/python2.7/site-packages/matplotlib/pyplot.py", line 2757,
in fill_between
interpolate=interpolate, **kwargs)
File "/usr/lib64/python2.7/site-packages/matplotlib/axes.py", line 6988, in
fill_between
x = ma.masked_invalid(self.convert_xunits(x))
File "/usr/lib64/python2.7/site-packages/numpy/ma/core.py", line 2239, in
masked_invalid
condition = ~(np.isfinite(a))
TypeError: ufunc 'isfinite' not supported for the input types, and the
inputs could not be safely coerced to any supported types according to the
casting rule ''safe''
>>>
Has anyone else encountered these problems before? Does anyone have any
suggestions for what I might do to overcome them?
Thankyou
Josiah
