metadata (or the lack thereof) in getwaveform from the iris client.

Users Mailing List Archive
1

Dear all,

I am using the iris client getwaveform to plot an event on all available stations from the IU network. However, the (mseed) trace object does not come with the station location information. (I would like to import this info to compute epicentral distance, so I can sort the waveforms from close to far.)

I see the arclink version of getwaveform comes with an optional argument metadata=true to get this info, but I have been unsuccessful to use the arclink client for the IU stations.

Any suggestions on how to either merge iris client waveform data with IU metadata, or use the arclink client for IU stations?

Cheers,
Kasper

2

Hello Kasper,

I have been using the dataless to find the station latitude and longitude by scanning through to find blockette 50 for a given station at a given time. I apologize for this not being optimal, looking fairly clunky, and my poor coding!

Note that the below code only gets the stations latitude and longitude and it doesn’t get the latitude and longitude at the sensor level (blockette 52).

If it helps this function sits in the following code: https://github.com/aringler-usgs/syncomppython/blob/master/syncomp.py

Best,
Adam


def getlatlon(sta,etime,xseedval):

#A function to get the lat and lon of a station at a given time

        for cursta in xseedval.stations:

#As we scan through blockettes we need to find blockettes 50

                for blkt in cursta:

                        if blkt.id == 50:

#Pull the station info for blockette 50
                                stacall = blkt.station_call_letters.strip()

                                if stacall == sta:

                                        lat = blkt.latitude

                                        lon = blkt.longitude        

                                        if type(blkt.end_effective_date) is str:

                                                curdoy = strftime("%j",gmtime())

                                                curyear = strftime("%Y",gmtime())

                                                curtime = UTCDateTime(curyear + "-" + curdoy + "T00:00:00.0") 

                                                if blkt.start_effective_date <= etime:

                                                        lat = blkt.latitude

                                                        lon = blkt.longitude

                                        elif blkt.start_effective_date <= etime and blkt.end_effective_date >= etime:

                                                lat = blkt.latitude

                                                lon = blkt.longitude        

        return lat,lon
3

Hi all,

I would like to retrieve mseed records for various streams from the SeedLink of a SeisComP3 server.
(Actually something like slinktool.)
Is there a way to use obspy for that purpose? I saw that there is a seedlink library in obspy, but there are no examples.
Could someone help?

Thank you in advance

4

The groundwork of obspy.seedlink came as an external contribution and
the documentation could indeed use a lot more detail/examples etc.
Anybody who is using it is welcome to improve on it.

Meanwhile, I can point you at the following example project that uses
obspy.seedlink:
https://github.com/bonaime/seedlink_plotter

The latest version has some complexity due to multiple threads for
different tasks, so you might want to look at earlier/less sophisticated
versions first:
https://github.com/bonaime/seedlink_plotter/blob/08f3460d2a5b0c55b0219cfbd85b8dc02e463924/seedlink_plotter.py

best,
Tobias

5

Hello,

The obspy.seedlink package is a Python port of libslink and JSeedLink (). obspy.seedlink.slclient.SLClient is a general purpose tool similar to slinktool. obspy.seedlink.tests.test_slclient.py shows very basic use of SLClient. Best regards, Anthony

6

Hi Adam,

the obspy.xseed Parser object also has a getCoordinates() method which does the same as your code but it looks for blockettes 52:

http://docs.obspy.org/packages/autogen/obspy.xseed.parser.Parser.getCoordinates.html

@Kasper: You need to use the obspy.iris.Client.station() method to get access to station metainformation. But the IRIS service will shut down soon anyways. The upcoming ObsPy version will support the FDSN webservices!

Cheers!

Lion

7

Adam and Lion,

Thanks for the feedback. I ended up doing things in a hybrid between your suggestions. See below. Note that in the plotting with the section command, I did some ugly things to get the xlim to be 180 degrees. It may be nice to set a max in the st.plot(‘section’) that offset cannot exceed 180 or $\sim$20,000 km.

I got most my info from the advanced plotting demo, but it may also be nice to have an example (like the one below) in the gallery?!

Anyway, I look forward to getting better at this, and one day making a contribution!
Kasper

#!/usr/bin/python
from obspy import UTCDateTime
from obspy.core.util import gps2DistAzimuth
from datetime import date
from obspy.iris import Client
client = Client()
from obspy import UTCDateTime
import matplotlib.pyplot as plt

read in the dataless information for all stations in the IU network

from here http://www.iris.edu/data/DatalessRequest.htm

from obspy.xseed import Parser
parser = Parser(“IU94.885435.dataless”)

bolivia EQ:

starttime = UTCDateTime(“1994-06-09T00:30:00.000”)
endtime = UTCDateTime(“1994-06-10T00:00:00.000”)

find all events between start and endtime above minmag:

cat = client.getEvents(starttime=starttime, endtime=endtime,minmag=7)
print(cat)

epicenter of the first event:

lat=cat[0].origins[0].latitude
lon=cat[0].origins[0].longitude

get MSEED data from network IU, component BHZ:

st = client.getWaveform(“IU”, “”, “”, “BHZ”, starttime, starttime+60606)

for tr in st:

get lat and lon from dataless:

seedid = tr.getId()
tr.stats.coordinates = parser.getCoordinates(seedid,starttime)

compute epicentral distance (not really used in the plot, though):

tr.stats.distance = gps2DistAzimuth(tr.stats.coordinates.latitude,tr.stats.coordinates.longitude,lat,lon)[0]

re-sample (otherwise too many points to plot):

tr.decimate(factor=2)

remove a trace with bad data from the stream:

for tr in st.select(station=“YAK”):
st.remove(tr)

filter the data to highlight the surface waves:

st.filter(‘lowpass’, freq=0.005)

#fig = plt.figure()
fig, ax = plt.subplots()
st.plot(fig=fig, type=‘section’,equal_scale=False, time_down=True, ev_coord=(lat,lon),linewidth=1,alpha=1,scale=1,dist_degree=True)

ugly (manual) way to set limit: 1 would be tight, but I want 180 to be max):

ax.set_xlim(xmax=1.1)
plt.xlabel(‘Epicentral distance [degrees]’)
fig.savefig(‘bolivia.pdf’, transparent=True, bbox_inches=‘tight’, pad_inches=0)
fig.show()