obspy / obspy

ObsPy: A Python Toolbox for seismology/seismological observatories.
https://www.obspy.org
Other
1.17k stars 536 forks source link

Problem of Vidale method of Polarization analysis #2441

Open xichaoqiang opened 5 years ago

xichaoqiang commented 5 years ago

a=polarization_analysis(st,20.0,1.0,10,12,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False) b=polarization_analysis(st,20.0,1.0,10,12,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False,"flinn") c=polarization_analysis(st,20.0,1.0,10,12,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False,"vidale") print(len(a['azimuth'])) print(len(b['azimuth'])) print(len(c['azimuth'])) Result is here.


55
55
569774

1.fcnt.zip

ThomasLecocq commented 5 years ago

Hi,

I just tested with your file and one thing could be the issue:

In [2]: st = read("1.fcnt/1.fcnt", format="rg16")

In [3]: st
Out[3]:
114 Trace(s) in Stream:

2001.5121.1.DP3 | 2018-05-29T09:42:46.550000Z - 2018-05-29T09:43:16.548000Z | 500.0 Hz, 15000 samples
...
(112 other traces)
...
2001.5121.1.DP2 | 2018-05-29T10:01:16.550000Z - 2018-05-29T10:01:46.548000Z | 500.0 Hz, 15000 samples

[Use "print(Stream.__str__(extended=True))" to print all Traces]

In [4]: st.merge()
Out[4]:
3 Trace(s) in Stream:
2001.5121.1.DP2 | 2018-05-29T09:42:46.550000Z - 2018-05-29T10:01:46.548000Z | 500.0 Hz, 570000 samples
2001.5121.1.DP3 | 2018-05-29T09:42:46.550000Z - 2018-05-29T10:01:46.548000Z | 500.0 Hz, 570000 samples
2001.5121.1.DP4 | 2018-05-29T09:42:46.550000Z - 2018-05-29T10:01:46.548000Z | 500.0 Hz, 570000 samples

Your data file doesn't contain Z, N or E components, I don't understand how the method could actually work (except if you forgot to copy something you do to the code above).

ThomasLecocq commented 5 years ago

but ok, even with the demo file (3000 samples @100 Hz), it gives different lengths:

st = read()
a=polarization_analysis(st,5.0,1.0,8,10,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False)
b=polarization_analysis(st,5.0,1.0,8,10,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False,"flinn")
c=polarization_analysis(st,5.0,1.0,8,10,st.traces[0].meta.starttime,st.traces[0].meta.endtime,False,"vidale")
print(len(a['azimuth']))
print(len(b['azimuth']))
print(len(c['azimuth']))

Out:

4
4
2944

not an expert, but @jwassermann can maybe explain why/if it's normal that the vidale method returns the azimuth for each sample (except some start/end buffer I guess)?

in the code:

if method.lower() == "vidale":
        res = vidale_adapt(stream, var_noise, fs, frqlow, frqhigh, spoint,
                           stime, etime)

doesn't use the wlen and win_frac indeed. If this is normal, it should be reflected in the Documentation.

obspy-bot commented 2 years ago

This issue has been mentioned on ObsPy Forum. There might be relevant details there:

https://discourse.obspy.org/t/polarization-examples/1388/8