eqcorrscan / EQcorrscan

Earthquake detection and analysis in Python.
https://eqcorrscan.readthedocs.io/en/latest/
Other
166 stars 86 forks source link

[HELP] A tribe with P and S pick from Sfile #480

Closed tonino0013 closed 2 years ago

tonino0013 commented 2 years ago
Python 3.9.7 --> Eqcorrscan EQcorrscan version: 0.4 # What do you need help with? Creating a Tribe from Sfile with P ans S picks. # Provide an example so that we can reproduce your problem So far I have the following code: ``` stream = read(r'D:\Inves\ML\2021-05-23-2339-06S.NSN___029') Sfile = read_events(r'D:\Inves\ML23-2340-02L.S202105', format='NORDIC') print (Sfile) print(stream.__str__(extended=True)) ``` ``` # ----- Selecting the P pick from one station Sfile[0].picks = [p for p in Sfile[0].picks if ((p.waveform_id.station_code == "SOET") and (p.waveform_id.channel_code == "HZ"))] # to work with miniseed data Sfile[0].picks[0].waveform_id.channel_code='HHZ' ``` If I work with only one station and one channel it is ok, however I am stuck if I want to select all EHZ (short period), HHZ and HHN channels from the sfile (nordic format) So I tried: `SfileP[0].picks = [p for p in Sfile[0].picks if ((p.waveform_id.channel_code == "HZ"))]` result: `[Pick(resource_id=ResourceIdentifier(id="smi:local/67219537-509f-49bc-9895-b926c7db93b7"), time=UTCDateTime(2021, 5, 23, 23, 40, 19, 770000), waveform_id=WaveformStreamID(network_code='NA', station_code='SOEA', channel_code='HZ'), onset='emergent', phase_hint='Pg', evaluation_mode='manual'), Pick(resource_id=ResourceIdentifier(id="smi:local/d7b9ea0b-8810-4d1d-91e8-e7d7e35c098f"), time=UTCDateTime(2021, 5, 23, 23, 40, 19, 850000), waveform_id=WaveformStreamID(network_code='NA', station_code='AOEA', channel_code='HZ'), onset='emergent', phase_hint='Pg', evaluation_mode='manual'),....` `SfileS[0].picks = [p for p in Sfile[0].picks if ((p.waveform_id.channel_code == "HN"))]` result: `[Pick(resource_id=ResourceIdentifier(id="smi:local/8e6a5468-fc55-449a-87c3-a88771a9ad77"), time=UTCDateTime(2021, 5, 23, 23, 40, 22, 170000), waveform_id=WaveformStreamID(network_code='NA', station_code='SOET', channel_code='HN'), onset='emergent', phase_hint='Sg', evaluation_mode='manual'), Pick(resource_id=ResourceIdentifier(id="smi:local/19be6d3a-6810-4b9b-85c4-c9d6195ef820"), time=UTCDateTime(2021, 5, 23, 23, 40, 23, 250000), waveform_id=WaveformStreamID(network_code='NA', station_code='SOET', channel_code='HN'), phase_hint='IAML', evaluation_mode='manual'),` Then I try to create the template: ``` template_base = Tribe().construct(method='from_meta_file', st=stream, meta_file=Sfile, length=25, lowcut=0.8, highcut=3.5, filt_order=4, process_len=600, swin='all', prepick=0.35, samp_rate=50, process=True) ``` but I get the following: ``` 2021-12-06 22:08:37,704 eqcorrscan.core.template_gen INFO Pre-processing data 2021-12-06 22:08:39,344 eqcorrscan.core.match_filter.tribe ERROR Empty Template ``` # What help would you like? How can I create a tribe with P and S pick from Sfile (Nordic format), just in case I upload the waveform and Sfile Thanks in advance for your time. [Sfile_and_Waveform.zip](https://github.com/eqcorrscan/EQcorrscan/files/7664951/Sfile_and_Waveform.zip) # What is your setup? (please complete the following information):** - Operating System: MS Window10 - Python version: 3.9.7 - EQcorrscan version 0.4
flixha commented 2 years ago

This appears to be most likely a problem of no matching channel names between picks and waveforms. When using the old Nordic format used in Seisan, then as you saw your channel codes are just two characters while the channel codes in the miniseed waveform files are 3 characters. So they don't match up right now, but there's a parameter that helps you work around that:

template_base = Tribe().construct(method='from_meta_file', st=stream,
                                  meta_file=Sfile, length=25,
                                  lowcut=0.8, highcut=3.5, filt_order=4,
                                  process_len=600, 
                                  swin='all', prepick=0.35, samp_rate=50, 
                                  process=True,
                                  seisan_chan_names=True
)

Alternatively (if you find the different channel names and missing location and network codes annoying), you can start using the new Nordic format in Seisan's latest version, for which obspy I/O-support is available in this PR .

tonino0013 commented 2 years ago

@flixha thank you, is this new feature included on eqcorrscan v0.4.3?, the command to update will be: conda update eqcorrscan? The new version of Seisan will be v12.0 rigth? Thank you

flixha commented 2 years ago

@tonino0013 Yes, the newest Seisan version is 12.0, though you will still need to switch it manually to use the new format. You'd also need to convert all old files then with a tool. The New Nordic I/O is part of obspy, not EQcorrscan, and it will be part of it once the changes are merged - right now if you need it you'd have to install the branch https://github.com/flixha/obspy/tree/new_nordic for obspy.

So - as the quickest solution I hope that seisan_chan_names=True solves the empty-template problem for you?

tonino0013 commented 2 years ago

@flixha , thank you again, so if I understand well:

Once I complete the update procedure I pointed, I will try to create the tribe. Thank you

flixha commented 2 years ago

@tonino0013 Sorry for mayeb confusing you. For now, you don't need to update anything. Just try to run:

template_base = Tribe().construct(method='from_meta_file', st=stream,
                                  meta_file=Sfile, length=25,
                                  lowcut=0.8, highcut=3.5, filt_order=4,
                                  process_len=600, 
                                  swin='all', prepick=0.35, samp_rate=50, 
                                  process=True,
                                  seisan_chan_names=True)

Note that I added the last line with seisan_chan_names=True. Does that work for you?

tonino0013 commented 2 years ago

Dear @flixha, Not worked, same error as I mentioned:

2021-12-07 09:31:58,154 eqcorrscan.core.template_gen    INFO    Pre-processing data
2021-12-07 09:31:59,662 eqcorrscan.core.match_filter.tribe  ERROR   Empty Template
Tribe of 0 templates
tonino0013 commented 2 years ago

Dear @flixha, I tried so far this option.


Sfile[0].picks = [p for p in Sfile[0].picks if ((p.waveform_id.channel_code == "HZ"))]
for ii in range(len(Sfile[0].picks)):
    Sfile[0].picks[ii].waveform_id.channel_code='HHZ'

and I obtained the tribe with Z channel (figure attach), I can do the same for S picks at HHN channel, but the how can I merge the two tribes to create one with P and S picks?

Tribe_test

Thanks in advance

calum-chamberlain commented 2 years ago

Hi all,

  1. I don't think seisan_chan_names is a supported kwarg in any current version of EQcorrscan.
  2. There isn't a great way to add templates together - you could just add the streams together and ensure that the template.event contains the full event
  3. The way I would do this is to just match your SEISAN picks to the stream and run that using a function like:

    def match_pick_to_channel(pick, st):
    """
    Try and match a SEISAN pick ID to a full seed ID in a stream.
    
    Parameters
    ----------
    pick:
        An obspy Pick object with SEISAN (old) format station and 2-char channel code
    st:
        Stream assciated with the pick.
    
    Returns
    -------
        Pick with full seed id.
    """
    possible_traces = [
        tr for tr in st if tr.stats.station == pick.waveform_id.station_code
        and f"{tr.stats.channel[0]}{tr.stats.channel[-1]}" == pick.waveform_id.channel_code]
    if len(possible_traces) == 0:
        print(f"No match found for {pick}")
        return None
    elif len(possible_traces) > 1:
        print(f"Found multiple possible matches for {pick}! I don't know what to do with this.")
        return None
    pick_back = pick.copy()
    tr = possible_traces[0]
    pick_back.waveform_id.channel_code = tr.stats.channel
    pick_back.waveform_id.network_code = tr.stats.network
    pick_back.waveform_id.location_code = tr.stats.location
    return pick_back

    called like:

corrected_picks = []
for pick in event.picks:
    corrected_pick = match_pick_to_channel(pick, st)
    if corrected_pick:
        corrected_picks.append(corrected_picks)
event.picks = corrected_picks
tonino0013 commented 2 years ago

Dear @calum-chamberlain, thanks for the code suggestions. I tried to implement the code you share, while I try to do the for-loop I got: AttributeError: 'Catalog' object has no attribute 'pick' But when I check the catalog object with var(Sfile) I have the following result:


{'events': [Event:  2021-05-23T23:39:57.700000Z | -17.387,  -66.105 | 3.5 ML

                resource_id: ResourceIdentifier(id="smi:local/2c557252-6c32-4ff4-ad46-6feee9400e47")
              creation_info: CreationInfo(agency_id='OSC')
        preferred_origin_id: ResourceIdentifier(id="smi:local/a2b9d5a0-3bbf-429d-a1c1-1a2dbceff8de")
     preferred_magnitude_id: ResourceIdentifier(id="smi:local/fe688ff2-7505-4570-9234-2c4166b9542e")
                       ---------
         event_descriptions: 1 Elements
                      picks: 20 Elements
                 amplitudes: 2 Elements
                    origins: 1 Elements
                 magnitudes: 1 Elements],
 'comments': [],
 'resource_id': smi:local/8f2bc85c-c49e-4f7e-ab55-e1e21ec77e23,
 'description': '',
 'creation_info': CreationInfo()}

I attach part of the code with your suggestions with sfile and waveform, what am I doing wrong?, or the implementation will be placed after the tribe generation?, thanks in advance.

Sfle_Wavs.zip

PandS.txt

calum-chamberlain commented 2 years ago

You have to loop through the picks of an individual event - a Catalog is a container of multiple events - you probably need a loop over your Catalog, or if you only have one Event in your Catalog then just index the zeroth element of the Catalog (e.g. for pick in catalog[0]:)

tonino0013 commented 2 years ago

Dear @calum-chamberlain, ok I will loop over the sfile in oder to get the picks (station, channel, time) and move them to obspy object catalog.

calum-chamberlain commented 2 years ago

One s-file usually only contains one event - a collect file will contain multiple events. ObsPy reads any event file into a Catalog object. For example:

from obspy import read_events

cat = read_events("insert-your-sfile-here")
print(len(cat))
# Will probably print "1"
event = cat[0]
for pick in event.picks:
# As before.

Your error before shows that you already have a Catalog object and you do not need to do any conversion (your Sfile variable is a Catalog because you read your nordic formatted file using obspy's read_events). You may want to read up on ObsPy Catalogs though and go through some of the obspy tutorials to get used to their objects and usage.

tonino0013 commented 2 years ago

Dear @calum-chamberlain, thank you for the suggestions, I will keep posted the solution, thanks for the time and patience, stay safe and best regards